Ship Emission Mitigation Strategies Choice Under Uncertainty

Abstract

Various mitigation strategies have been proposed to reduce the CO₂ emissions from ships, which have become a major contributor to global emissions. The fuel consumption under different mitigation strategies can be evaluated based on two data sources, real data from the real ship systems and simulated data from the simulation models. In practice, the uncertainties in the obtained data may have non-negligible impacts on the evaluation of mitigation strategies. In this paper, a Gaussian process metamodel-based approach is proposed to evaluate the ship fuel consumption under different mitigation strategies. The proposed method not only can incorporate different data sources but also consider the uncertainties in the data to obtain a more reliable evaluation. A cost-effectiveness analysis based on the fuel consumption prediction is then applied to rank the mitigation strategies under uncertainty. The accuracy and efficiency of the proposed method is illustrated in a chemical tanker case study, and the results indicate that it is critical to consider the uncertainty, as they can lead to suboptimal decisions when ignored. Here, trim optimisation is ranked more effective than draft optimisation when the uncertainty is ignored, but the reverse is the case when the uncertainty in the estimations are fully accounted for.

1. Introduction

The transportation sector contributes to nearly one-quarter of the total global CO₂ emission in 2018 [1]. Among them, road transport is the largest emission source and is responsible for 75% of the transportation sector emission. The aviation and the rail account for around 2.5% and 0.3% of the global CO₂ emission, respectively. International shipping transports 90% of the international traded goods and contributes around 693 million metric tons of CO₂ in 2018, accounting for about 2.1% of total global CO₂ emissions. Although the international shipping is considered as the most energy-efficient way for freight transport, the carbon emission from international shipping is likely to increase by 50–250% in 2050 under current business trends [2].

In an attempt to reduce CO₂ emissions, the International Maritime Organization (IMO) [3] identified more than 50 technical and operational mitigation strategies for improving the efficiency of international shipping and, hence, reducing CO₂ emissions. However, due to economic and other constraints, it is not practical to implement all relevant mitigation strategies at the same time. Very often, selected strategies are evaluated and ranked according to costs and CO₂ emission reduction potentials so as to minimise impacts on business competitiveness and the environment [4,5].

Fuel consumptions account for the bulk of the life cycle CO₂ emissions of a ship [6,7]. Reduction in fuel consumptions after implementing the mitigation strategies can be used to compute the CO₂ emission reductions [8]. The cost of a mitigation strategy typically can be calculated from the investment cost, operational cost, opportunity cost and savings of fuel cost resulting from reduction in fuel consumptions after implementing the mitigation strategy [9]. Since emissions reductions and costs usually represent two conflicting objectives, it is necessary to address both objectives when selecting the most cost-effective strategies for improving efficiency [10].

In the case of international shipping, there are three main approaches in estimating reductions in fuel consumptions by different mitigation strategies. The first approach is based on expert judgement. The IMO Marine Environment Protection Committee (MEPC) 62 report [3] examined reductions in fuel consumptions for 28 mitigation strategies, including 8 operational strategies and 20 technical strategies. This approach is usually considered to be subjective [2]. The second approach is based on data from experiments conducted on real ships. Wärtsilä [11] investigated reductions in fuel consumptions for several mitigation strategies, such as lightweight construction, optimisation of hull dimensions, propulsion upgrade and optimisation of the trim and ballast through experiments on real ships. This approach is usually reliable and accurate but can be expensive and time-consuming. The third approach is to use stochastic simulation models of an energy system of a ship. Baldi and Gabrielii [12] employed a stochastic simulation model to evaluate the feasibility of a waste heat recovery system for a ship. Simulation models are much more economical and efficient in deriving insights on mitigation strategies and are thus more appealing as compared to experiments on a real ship. However, simulation models can be expensive to develop and require sufficient expert knowledge.

Recently, there is an increasing trend to use metamodels, simpler and cheaper statistical approximation models to simulation models, to evaluate fuel consumptions in the literature [13,14]. Among the various approaches, Gaussian process (GP)-based metamodels [15,16] are widely applied and shown to have advantages in representing ship’s energy systems [17,18]. Such metamodels can be built using real and/or simulated data on fuel consumptions. Real data obtained through on-board measurements are usually more reliable than simulated data, but the quantity of real observed data is usually much lower than that of simulated data [12]. In the attempt to leverage on the advantages of both real and simulated data, Yuan et al. [19] proposed a GP-based model that uses both data sources simultaneously. Results from their study suggest that the simultaneous use of different data sources can lead to more reliable predictions on fuel consumptions as compared to the use of only one data source. However, their model did not consider the uncertainties associated with the data.

There are two types of inputs influencing the fuel consumptions of a ship, namely, control and environmental inputs. Control inputs literally refer to variables that can be controlled by human intervention, such as the ship speed and trim. In existing studies, control inputs are always assumed to be known and/or fixed when analysing the fuel consumptions of ships. However, the observed “control inputs” can be subject to the influence of external uncontrollable factors, leading to uncertainties in the observed control inputs [20,21]. For example, signal errors can cause uncertainties in the speed over ground (SOG) obtained from the Automatic Identification System (AIS), which is a system that uses transponders on ships for the automatic tracking of vessels. In the past, the SOG obtained from the AIS was always used directly as the true control input. Neglecting the uncertainties associated with observed SOG can lead to erroneous decisions in mitigation strategies, as discussed in [21].

The environmental inputs refer to inputs that cannot be controlled by humans but have an influence on the ship’s fuel consumptions [22]. For example, the distribution of the random electricity demand on a ship is considered as an environmental input. Very often, only a limited set of historical data can be obtained and used to estimate the distribution of the electricity demand, leading to uncertainty in the estimation [23,24]. Neglecting the uncertainties associated with the environmental inputs can lead to suboptimal or misinformed decisions [25]. As such, there is a need to consider the uncertainty associated with the environmental inputs [26,27] to hedge against those uncertainties in a metamodel [28] and an optimisation model [22]. Ignoring those uncertainties can lead to non-negligible impacts on the ranking results of the mitigation strategies.

In this work, the GP metamodel-based method is proposed to investigate the fuel consumptions of ships. The proposed model is capable of accounting for the uncertainties in the input data source while utilising both real and simulated data. The ability to account for uncertainties in the model formulation is considered an important extension from the previous model, as described in [19], to more realistically reflect the characteristics of the data collected. Findings from the case study suggest that the proposed GP metamodel is advantageous over models that ignore these input uncertainties. Based on the predicted fuel consumptions by the proposed model, the CO₂ emission reductions and costs of mitigation strategies are further assessed. These selected mitigation strategies are ranked according to their marginal cost-effectiveness (MCE) in reducing CO₂ emissions, as proposed in [9].

The main contributions of this study can be summarised as follows:

• A data-driven approach utilising both real data and simulated data is proposed to evaluate the fuel consumptions of a ship under different mitigation strategies.
• The proposed GP model is able to account for the uncertainty associated with control inputs and environmental inputs and propagates such input uncertainty when predicting the distribution of the fuel consumptions.
• A more reliable and robust ranking of mitigation strategies is developed by accounting for the uncertainty in CO₂ emission reductions and costs.

The rest of the paper is organised as follows: The proposed GP metamodel used to evaluate the fuel consumptions is presented in Section 2. The CO₂ emission reduction and cost assessments, including the MCE criterion for ranking the mitigation strategies, are presented in Section 3. The case study is presented in Section 4. Conclusions and recommendations for future research are presented in Section 5.

2. Methodology

2.1. Energy System of a Ship

Fuel consumed in a ship is converted to different energy forms, such as mechanical power, electricity power and thermal power, for different purposes. From an energy perspective, the entire ship can be treated as an energy system. An example of a ship energy system for a chemical tanker is shown in Figure 1.

In this ship energy system, fuel is converted to mechanical energy through the main engines. The mechanical energy is mainly used for propulsion. The propulsion is mostly influenced by resistance, which depends on various control input factors, such as speed, trim and draft. As mentioned earlier, control inputs refer to those inputs that can be controlled or designed in practice [23]. Electricity power can be generated from both auxiliary engines and main engines. Electricity is the most common form of energy used for powering on-board facilities such as the HVAC (acronym for heating ventilation and air-conditioning), cargo pumps and nitrogen compressor. Due to the random nature of demand for an HVAC and other electrical systems, the distribution of electricity demand is treated as an environmental input. Thermal power can be generated from the boilers, auxiliary engines and main engines. The thermal power is mainly used for fuel heating, accommodation and tank cleaning. Thermal energy is not modelled in this study due to its complexity and lack of data.

2.2. Data Sources

In general, fuel consumption data can be obtained through two sources, namely real fuel consumption data observed from on-board measurements and simulated fuel consumption data simulated from simulation models. As explained earlier, the real fuel consumption data are low in quantity due to limited reporting but tend to be more accurate. In comparison, the simulated fuel consumption data can be much higher in quantity but heavily rely on the assumptions of the simulation model.

The control inputs for the real fuel consumption data can be obtained directly from both the on-board measurements and AIS. More specifically, control inputs such as the fan and pump speed can be obtained from on-board measurements. SOG, trim and draft values and the navigating course can be obtained from the AIS. The environmental input considered in this paper is the distribution of the electricity demand. Electricity demand can be estimated based on historical data obtained from on-board measurements. These data are summarised in

Table 1. Sources and variables from real data. AIS: Automatic Identification System and SOG: speed over ground.

	Data Source	Variable
Output	On-board measurements	Real fuel consumption
Control input	On-board measurements	Speed of pumps and fans
	AIS	SOG, trim, draft and course
Environmental input	On-board measurements	Distribution of electricity demand

.

There are uncertainties associated with control inputs obtained from the AIS mainly due to signal errors caused by transponders. On the other hand, the uncertainties associated with control inputs obtained from on-board measurements are caused by inaccuracies in the ship’s control system. For example, the speed of pumps and fans can be set using the monitor on-board, but the actual speed of pumps and fans may differ slightly due to engine, weather and sea conditions. The distribution of electricity demand as the environmental input is also subject to uncertainties due to the fact that only limited historical data of on-board measurements can be obtained to estimate the distribution of the electricity demand.

For the simulated data, the control and environmental inputs can be directly obtained from the simulation model. In addition, as the stochastic simulation model can be evaluated at every precise control and environmental input, there is no control and environmental input uncertainty associated with the inputs of the simulation data.

Both real and simulated data are used simultaneously to build a GP model in this study in two steps. The first step is a GP model without considering the uncertainties, and the next step is a GP metamodel incorporating uncertainties in the formulation.

2.3. Model Formulation without Input Uncertainties

In this model, $\mathit{x}\in R^d$ is defined as a vector of control inputs in the ship energy system, such as the speed, course and trim. $\mathit{\lambda }\in R^l$ is defined as a vector of the environmental input, which is the distribution of electricity demand in our setting. A known parametric form for the distribution of the environmental input is assumed, and the associated parameter $\mathit{\lambda }$ is used as the environmental input. This parametric assumption is commonly used in studies dealing with environmental inputs [26,27]. In this study, the electricity demand is assumed to follow an exponential distribution exp($\mathit{\lambda })$, where $\mathit{\lambda }$ is the rate of the exponential distribution.

$y_r\left(\mathit{x},\mathit{\lambda }\right)$ is defined as the real system output (real fuel consumptions observed from the on-board measurements), and $y_s\left(\mathit{x},\mathit{\lambda }\right)$ is defined as the simulation output (simulated fuel consumptions obtained from the simulation model) for a specific control input $\mathit{x}$ and environmental input $\mathit{\lambda }$. Since the simulation model is an approximation to the real ship, the relationship between the expected simulation output and the real system output can be represented as

$$y_r\left(\mathit{x},\mathit{\lambda }\right)=z\left(\mathit{x},\mathit{\lambda }\right)+e,$$

(1)

where $z\left(\mathit{x},\mathit{\lambda }\right)$ represents the true system output (true fuel consumption), and $e$ represents the observation error.

The model as described in Equation (1) has been commonly used to represent the relationship between different outputs [29,30,31]. For stochastic simulation models, there is also a stochastic error. The relationship between the observed simulation output of fuel consumption (with noise) and the expected simulation output of fuel consumption can be expressed as

$$y_s\left(\mathit{x},\mathit{\lambda }\right)=S\left(\mathit{x},\mathit{\lambda }\right)+ϵ=z\left(\mathit{x},\mathit{\lambda }\right)-\delta \left(\mathit{x},\mathit{\lambda }\right)+ϵ,$$

(2)

where $S\left(\mathit{x},\mathit{\lambda }\right)$ represents the expected simulation output of fuel consumption, $\delta \left(\mathit{x},\mathit{\lambda }\right)$ represents the discrepancy between the true fuel consumption and the expected simulation output of fuel consumption and $\epsilon$ represents the stochastic error.

As explained earlier, there are two types of data sources for a ship’s fuel consumption, namely observed fuel consumption data from on-board measurements and estimated data from a simulation model. The objective is to develop a model based on these two sources of data to predict the true fuel consumption $z\left(\mathit{x},\mathit{\lambda }\right)$ of a ship under input control $\mathit{x}$ and environmental conditions $\mathit{\lambda }$.

A GP is a stochastic process in which a finite collection of random variables has a multivariate normal distribution [22]. Due to its flexibility and efficiency, the GP has been widely used as a metamodel to study complex simulation models and real systems. We adopt this form and model both the true fuel consumption $z\left(\mathit{x},\mathit{\lambda }\right)$ and the discrepancy $\delta \left(\mathit{x},\mathit{\lambda }\right)$ in Equation (2) with independent GPs. The same assumption is also used in [19,23] and has been shown to be reasonable in practice [23]. More specifically, the true fuel consumption $z\left(\mathit{x},\mathit{\lambda }\right)$ is assumed to be a GP with a mean function ${\beta }_z\left(\mathit{x},\mathit{\lambda }\right)$ and covariance function ${\sigma }_z^2{R}_z\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)$. The discrepancy $\delta \left(\mathit{x},\mathit{\lambda }\right)$ is assumed to be a GP with a mean function ${\beta }_{\delta }\left(\mathit{x},\mathit{\lambda }\right)$ and a covariance function ${\sigma }_{\delta }^2{R}_{\delta }\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)$.

With reference to [23], several mean and covariance functions can be used to specify the GP. In this case, the mean functions for $z\left(\mathit{x},\mathit{\lambda }\right)$ and $\delta \left(\mathit{x},\mathit{\lambda }\right)$ are assumed to be unknown constants ${\beta }_z$ and ${\beta }_{\delta }$, respectively. The correlation function is assumed to be Gaussian. The covariance function for the true system output $z\left(\mathit{x},\mathit{\lambda }\right)$ can be expressed as

$${\sigma }_z^2{R}_z\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)={\sigma }_z^2\exp \left(-\frac{1}{2}{\left(\mathit{x}-{\mathit{x}}^{\prime }\right)}^{\mathrm{T}}{\mathit{W}}_{z,x}{}^{-1}\left(\mathit{x}-{\mathit{x}}^{\prime }\right)\right)\exp \left(-\frac{1}{2}{\left(\mathit{\lambda }-{\mathit{\lambda }}^{\prime }\right)}^T{\mathit{W}}_{z,\lambda }{}^{-1}\left(\mathit{\lambda }-{\mathit{\lambda }}^{\prime }\right)\right).$$

(3)

The covariance function for the discrepancy $\delta \left(\mathit{x},\mathit{\lambda }\right)$ can be expressed as

$${\sigma }_{\delta }^2{R}_{\delta }\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)={\mathsf{\sigma }}_{\mathsf{\delta }}^2\exp \left(-\frac{1}{2}{\left(\mathit{x}-{\mathit{x}}^{\prime }\right)}^{\mathrm{T}}{\mathit{W}}_{\delta ,x}{}^{-1}\left(\mathit{x}-{\mathit{x}}^{\prime }\right)\right)\exp \left(-\frac{1}{2}{\left(\mathit{\lambda }-{\mathit{\lambda }}^{\prime }\right)}^{\mathrm{T}}{\mathit{W}}_{\delta ,\lambda }{}^{-1}\left(\mathit{\lambda }-{\mathit{\lambda }}^{\prime }\right)\right),$$

(4)

where ${\mathit{W}}_{z,x},{\mathit{W}}_{z,\lambda },{\mathit{W}}_{\delta ,x},{\mathit{W}}_{\delta ,\lambda }$ represent the diagonal matrices. The covariance functions ${\sigma }_z^2{R}_z\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)$ and ${\sigma }_{\delta }^2{R}_{\delta }\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)$ are used to capture the spatial correlation of the true fuel consumptions and the discrepancy between the real fuel consumption and simulated consumption under different control and environmental input pairs. A smaller Euclidean norm between $\left(\mathit{x},\mathit{\lambda }\right)$ and $\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)$ signifies a higher correlation. A higher value along the diagonal of the diagonal matrices signifies a higher spatial correlation in that dimension of control input and environmental input.

Based on Equations (1) and (2), there is also an observation error term $e$ for the real system output and a stochastic error term $\epsilon$ for the stochastic simulation model output. It is assumed that these two terms follow a normal distribution with a mean of zero and variances of ${\sigma }_e^2$ and ${\sigma }_{ϵ}^2$, respectively. These assumptions are reasonable according to the central limit theorem [32]. Under these assumptions, the unknown parameters can be represented as $\left\{\beta ,\mathit{W},{\sigma }^2\right\},$ where $\beta =\left\{{\beta }_z,{\beta }_{\delta }\right\}$ represent the unknown constant mean values for the GP, $\mathit{W}=\left\{{\mathit{W}}_{z,x},{\mathit{W}}_{z,\lambda },{\mathit{W}}_{\delta ,x},{\mathit{W}}_{\delta ,\lambda }\right\}$ represent the unknown decaying parameters in the Gaussian correlation functions and ${\sigma }^2=\left\{{\sigma }_z^2,{\sigma }_{\delta }^2,{\sigma }_e^2,{\sigma }_{ϵ}^2\right\}$ represent the unknown variances in the model.

The purpose is to derive the predictive distribution of the true fuel consumption $\mathrm{z}\left(\mathit{x},\mathit{\lambda }\right)$ given the available data sources. ${\mathit{Y}}_r$ is defined as the real fuel consumption observed from the ship on-board measurements at the input settings $({\mathit{x}}_r,{\mathit{\lambda }}_r)$ in the ship energy system. ${\mathit{Y}}_s$ is defined as the simulated fuel consumption at the input settings $({\mathit{x}}_s,{\mathit{\lambda }}_s)$ specified in the simulation model. Based on the characteristics of the multivariate normal distribution, the predictive distribution of the true fuel consumption $\mathrm{z}\left(\mathit{x},\mathit{\lambda }\right)$ for any unknown input settings $\left(\mathit{x},\mathit{\lambda }\right)$ can be derived as a conditional normal distribution with the following form:

$$\mathrm{z}\left(\mathit{x},\mathit{\lambda }\right)|{\mathit{Y}}_r,{\mathit{x}}_r,{\mathit{\lambda }}_r,{\mathit{Y}}_s,{\mathit{x}}_s,{\mathit{\lambda }}_s~\mathrm{N}\left(\mathrm{m}\left(\mathit{x},\mathit{\lambda }\right),{\mathsf{\sigma }}^2\left(\mathit{x},\mathit{\lambda }\right)\right),$$

(5)

in which

$$\mathrm{m}\left(\mathit{x},\mathit{\lambda }\right)={\beta }_z+\mathsf{\Sigma }{\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right)\right)}^{\mathrm{T}}\mathsf{\Sigma }{\left(\left({\mathit{x}}_D,{\mathit{\lambda }}_D),({\mathit{x}}_D,{\mathit{\lambda }}_D\right)\right)}^{-1}\left({\mathit{Y}}_D-[\begin{array}{c}{\beta }_z{\mathbf{1}}_m\\ ({\beta }_z-{\beta }_{\delta }){\mathbf{1}}_n\end{array}]\right),$$

(6)

$${\mathsf{\sigma }}^2\left(\mathit{x},\mathit{\lambda }\right)=\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left(\mathit{x},\mathit{\lambda }\right)\right)-\mathsf{\Sigma }{\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right)\right)}^{\mathrm{T}}\mathsf{\Sigma }{\left(\left({\mathit{x}}_D,{\mathit{\lambda }}_D),({\mathit{x}}_D,{\mathit{\lambda }}_D\right)\right)}^{-1}\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right)\right),$$

(7)

where $\mathrm{m}\left(\mathit{x},\mathit{\lambda }\right)$ represents the predictive mean of fuel consumptions, and ${\mathsf{\sigma }}^2\left(\mathit{x},\mathit{\lambda }\right)$ represents the predictive variance. $({\mathit{x}}_D,{\mathit{\lambda }}_D)=\left[\begin{array}{c}({\mathit{x}}_r,{\mathit{\lambda }}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_s)\end{array}\right]$ represents the concatenated input settings, including both real system inputs $({\mathit{x}}_r,{\mathit{\lambda }}_r)$ and simulation inputs $({\mathit{x}}_s,{\mathit{\lambda }}_s)$, where $({\mathit{x}}_r,{\mathit{\lambda }}_{\mathrm{r}})$ is a subset of $({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})$. ${\mathit{Y}}_{\mathrm{D}}=\left({\mathit{Y}}_r,{\mathit{Y}}_s\right)$ represents the outputs, including both real-system observations and simulated outputs. ${\mathbf{1}}_m$ represents an m-dimensional vector of 1, where $m$ is the number of real observations, and n is the number of simulation outputs. $\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right)\right)={\sigma }_z^2{R}_z\left(\left(\mathit{x},\mathit{\lambda }\right)\right),\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right))$ represents the covariance between $\left(\mathit{x},\mathit{\lambda }\right)$ and $({\mathit{x}}_D,{\mathit{\lambda }}_D)$ with a Gaussian correlation. $\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left(\mathit{x},\mathit{\lambda }\right)\right)={\sigma }_z^2{R}_z\left(\left(\mathit{x},\mathit{\lambda }\right)\right),\left(\mathit{x},\mathit{\lambda }\right))$ represents the variance at $\left(\mathit{x},\mathit{\lambda }\right)$. $\mathsf{\Sigma }\left(\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right),\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right)\right)$ represents the covariance of $\left({\mathit{x}}_D,{\mathit{\lambda }}_D\right)$, which can be represented as

$$\mathsf{\Sigma }\left(({\mathit{x}}_D,{\mathit{\lambda }}_D),({\mathit{x}}_D,{\mathit{\lambda }}_D)\right)=\left(\begin{array}{cc}{\sigma }_z^2{R}_z\left(({\mathit{x}}_r,{\mathit{\lambda }}_r),\left({\mathit{x}}_{r,}{\mathit{\lambda }}_r\right)\right)+{\sigma }_e^2{I}_m& {\sigma }_z^2{R}_z\left(({\mathit{x}}_{r,}{\mathit{\lambda }}_r),\left({\mathit{x}}_{s,}{\mathit{\lambda }}_s\right)\right)\\ {\sigma }_z^2{R}_z\left(({\mathit{x}}_s,{\mathit{\lambda }}_s),\left({\mathit{x}}_{r,}{\mathit{\lambda }}_r\right)\right)& {\sigma }_z^2{R}_z\left(\left({\mathit{x}}_{s,}{\mathit{\lambda }}_s\right),\left({\mathit{x}}_{s,}{\mathit{\lambda }}_s\right)\right)+{\sigma }_{\delta }^2{R}_{\delta }\left(\left({\mathit{x}}_{s,}{\mathit{\lambda }}_s\right),\left({\mathit{x}}_{s,}{\mathit{\lambda }}_s\right)\right)+{\sigma }_{\epsilon }^2{I}_n\end{array}\right),$$

(8)

where ${\mathit{I}}_m$ represents an $m\times m$ identity matrix, and ${\mathit{I}}_n$ represents an $n\times n$ identity matrix.

The MLE estimation method [32] can be used to estimate the parameters. The loglikelihood function of these unknown parameters, $l\left(\beta ,{\sigma }^2,\mathit{W}\right)$, can be expressed as:

$$l\left(\beta ,\mathit{W},{\sigma }^2\right)=-\ln \left[{\left(2\pi \right)}^{\frac{d}{2}}\right]-\frac{1}{2}\ln \left[\left|\mathsf{\Sigma }\left(({\mathit{x}}_D,{\mathit{\lambda }}_D),({\mathit{x}}_D,{\mathit{\lambda }}_D)\right)\right|\right]-\frac{1}{2}{\left({\mathit{Y}}_D-[\begin{array}{c}{\beta }_z{\mathbf{1}}_m\\ ({\beta }_z-{\beta }_{\delta }){\mathbf{1}}_n\end{array}]\right)}^T{\left[\mathsf{\Sigma }\left(({\mathit{x}}_D,{\mathit{\lambda }}_D),({\mathit{x}}_D,{\mathit{\lambda }}_D)\right)\right]}^{-1}\left({\mathit{Y}}_D-[\begin{array}{c}{\beta }_z{\mathbf{1}}_m\\ ({\beta }_z-{\beta }_{\delta }){\mathbf{1}}_n\end{array}]\right),$$

(9)

The parameters $\beta$, $\mathit{W}$ and ${\sigma }^2$ can be estimated by maximising $l\left(\beta ,\mathit{W},{\sigma }^2\right)$. Various standard optimisation packages can be used for this estimation, such as “mle.tools” in the R package used in this study.

2.4. Model Formulation with Input Uncertainties

2.4.1. Uncertainties in Observed Inputs

As mentioned earlier, in reality, the observed control inputs such as SOG in a real system suffer from control input uncertainty due to AIS noise. The uncertainty in the observed control input location, denoted as ${\mathit{u}}_i$, can be assumed to be represented by a distribution ${\mathit{u}}_i~P_{{\mathit{u}}_i}\triangleq N\left({\mathit{x}}_i,{\mathsf{\Sigma }}_i\right)$, $i=1\cdots n$, where ${\mathsf{\Sigma }}_i$ is assumed to be known. Besides, the uncertainties of the control inputs at different locations are independent.

Input uncertainty is also present in the environmental inputs. Typically, a parametric form is assumed for the environmental inputs, and for our ship electricity demand case, we assume a known parametric form with parameter $\mathit{\lambda }$ as the environmental input. The uncertainty associated with $\mathit{\lambda }$ happens due to the finiteness of the historical data used to estimate $\mathit{\lambda }$. In order to account for the environmental input uncertainty, a Bayesian perspective is taken by treating $\mathit{\lambda }$ as a random variable. More specifically, the Bayesian approach is to compute the posterior distribution of $\mathit{\lambda }$, defined as $p(\mathit{\lambda }|D_h)$, where $D_h$ is the historical electricity demand obtained from on-board measurements and $h$ is the data volume. Then, the uncertainty in the environmental input locations ${\mathit{\lambda }}_i,i=1\cdots n$ can be characterised as ${\mathit{\lambda }}_i~P_{{\mathit{\lambda }}_i}\triangleq p(\mathit{\lambda }|D_h)$.

2.4.2. Modelling Real Fuel Consumption $y_r$ with Input Uncertainty

With reference to Equation (1), when uncertainties are not considered, the real fuel consumption $z\left(\mathit{x},\mathit{\lambda }\right)$ can be modelled using GP with a mean ${\beta }_z$ and a covariance ${\sigma }_z^2{R}_z\left(\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right),\left({\mathit{x}}_j,{\mathit{\lambda }}_j\right)\right)$. Thus,

$$\mathrm{E}({\mathrm{y}}_{\mathrm{r}}\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right))={\beta }_z,$$

(10)

$$\mathrm{Cov}[{\mathrm{y}}_{\mathrm{r}}\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right),{\mathrm{y}}_{\mathrm{r}}\left({\mathit{x}}_j,{\mathit{\lambda }}_j\right)]={\sigma }_z^2{R}_z\left(\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right),\left({\mathit{x}}_j,{\mathit{\lambda }}_j\right)\right)+{\gamma }_{ij}{\sigma }_e^2,$$

(11)

where ${\gamma }_{ij}$ is the Kronecker delta.

When the input uncertainties are considered in the modelling, the input locations are now ${\left\{\left(P_{{\mathit{u}}_i},P_{{\mathit{\lambda }}_i}\right)\right\}}_{\mathrm{i}=1}^n$ instead of the fixed input locations ${\left\{\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right)\right\}}_{\mathrm{i}=1}^n$. The resulting real system output (fuel consumption) can be defined as follows:

$${\mathrm{y}}_{\mathrm{r}}\left(P_{{\mathit{u}}_i},P_{{\mathit{\lambda }}_i}\right)=E_{{\mathit{u}}_i~P_{{\mathit{u}}_i}}{E}_{{\mathit{\lambda }}_i~P_{{\mathit{\lambda }}_i}}\left[{\mathrm{y}}_{\mathrm{r}}\left({\mathit{u}}_i,{\mathit{\lambda }}_i\right)\right]=\int \int {\mathrm{y}}_{\mathrm{r}}\left({\mathit{u}}_i,{\mathit{\lambda }}_i\right)p\left({\mathit{u}}_i\right)p\left({\mathit{\lambda }}_i\right)d{\mathit{u}}_{i}d{\mathit{\lambda }}_i.$$

(12)

${\mathrm{y}}_{\mathrm{r}}\left(P_{{\mathit{u}}_i},P_{\mathit{\lambda }}{}_i\right)$ can be modelled as a GP, as it can be rewritten as a Riemann sum of Gaussian random variables [33]. The mean and the covariance function of ${\mathrm{y}}_{\mathrm{r}}\left(P_{{\mathit{u}}_i},P_{\mathit{\lambda }}{}_i\right)$ can be derived to obtain its distribution:

$$\mathrm{E}\left({\mathrm{y}}_{\mathrm{r}}\left(P_{{\mathit{u}}_i},P_{\mathit{\lambda }}{}_i\right)\right)={\mathrm{E}}_{{\mathit{u}}_i~P_{{\mathit{u}}_i}}{\mathrm{E}}_{{\mathit{\lambda }}_i~P_{{\mathit{\lambda }}_i}}\left[\mathrm{E}\left[{\mathrm{y}}_{\mathrm{r}}\left({\mathit{u}}_i,{\mathit{\lambda }}_i\right)\right]\right]={\mathrm{E}}_{{\mathit{u}}_i}{\mathrm{E}}_{{\mathit{\lambda }}_i}[{\mathsf{\beta }}_{\mathrm{z}}]={\mathsf{\beta }}_{\mathrm{z}},$$

(13)

$$\mathrm{Cov}\left(y_r\left(P_{{\mathit{u}}_i},P_{\mathit{\lambda }}{}_i\right),y_r\left(P_{{\mathit{u}}_j},P_{\mathit{\lambda }}{}_j\right)\right)={\sigma }_z^2{\mathrm{E}}_{{\mathit{u}}_i~P_{{\mathit{u}}_i}}{\mathrm{E}}_{{\mathit{u}}_j~P_{{\mathit{u}}_j}}\left[R_{z,\mathit{x}}\left({\mathit{u}}_i,{\mathit{u}}_j\right)\right]{\mathrm{E}}_{{\mathit{\lambda }}_i~P_{{\mathit{\lambda }}_i}}{\mathrm{E}}_{{\mathit{\lambda }}_j~P_{{\mathit{\lambda }}_j}}\left[R_{z,\mathit{\lambda }}\left({\mathit{\lambda }}_i,{\mathit{\lambda }}_j\right)\right]+{\gamma }_{ij}{\sigma }_e^2,$$

(14)

where $R_{z,\mathit{x}}\left({\mathit{u}}_i,{\mathit{u}}_j\right)=\exp \left(-\frac{1}{2}{\left({\mathit{u}}_i-{\mathit{u}}_j\right)}^{\mathrm{T}}{\mathrm{W}}_{\mathrm{z},\mathrm{x}}{}^{-1}\left({\mathit{u}}_i-{\mathit{u}}_j\right)\right)$ and $R_{z,\mathit{\lambda }}\left({\mathit{\lambda }}_i,{\mathit{\lambda }}_j\right)=\exp \left(-\frac{1}{2}{\left({\mathit{\lambda }}_i-{\mathit{\lambda }}_j\right)}^T{\mathit{W}}_{z,\mathit{\lambda }}{}^{-1}\left({\mathit{\lambda }}_i-{\mathit{\lambda }}_j\right)\right)$.

Suppose that the control input uncertainty follows a normal distribution ${\mathit{u}}_i~\mathrm{N}\left({\mathit{x}}_i,{\mathsf{\Sigma }}_i\right)$, ${\mathit{u}}_j\sim \mathrm{N}\left({\mathit{x}}_j,{\mathsf{\Sigma }}_j\right)$ and ${\mathsf{\Sigma }}_i,{\mathsf{\Sigma }}_j$ are assumed to be known, then ${\mathrm{E}}_{{\mathit{u}}_i}{\mathrm{E}}_{{\mathit{u}}_j}[R_{z,\mathit{x}}\left({\mathit{u}}_i,{\mathit{u}}_j\right)]$ can be computed analytically, as shown in [20]:

$${\mathrm{E}}_{{\mathit{u}}_i}{\mathrm{E}}_{{\mathit{u}}_j}[R_{z,x}\left({\mathit{u}}_i,{\mathit{u}}_j\right)]=\frac{\exp \left(-\frac{1}{2}{\left({\mathit{x}}_i-{\mathit{x}}_j\right)}^T{\left({\mathit{W}}_{z,x}+{\mathsf{\Sigma }}_i+{\mathsf{\Sigma }}_j\right)}^{-1}\left({\mathit{x}}_i-{\mathit{x}}_j\right)\right)}{{\left|I+{\mathit{W}}_{z,x}{}^{-1}\left({\mathsf{\Sigma }}_i+{\mathsf{\Sigma }}_j\right)\right|}^{\frac{1}{2}}}.$$

(15)

In addition, as $P_{{\lambda }_i}$= $P_{{\lambda }_j}$ = $p(\mathit{\lambda }|D_h)$, we have

$$E_{{\mathit{\lambda }}_i~P_{{\mathit{\lambda }}_i}}{E}_{{\mathit{\lambda }}_j~P_{{\mathit{\lambda }}_j}}\left[R_{z,\lambda }\left({\mathit{\lambda }}_i,{\mathit{\lambda }}_j\right)\right]=R_{z,\lambda }\left(P_{{\lambda }_i},P_{{\lambda }_j}\right)=R_{z,\lambda }\left(P_{{\lambda }_i},P_{{\lambda }_i}\right)=1.$$

(16)

Then, with the above results, Equation (14) can be calculated as follows:

$$\mathrm{Cov}[y_r\left(P_{{\mathit{u}}_j},P_{{\mathit{\lambda }}_i}\right),y_r\left(P_{{\mathit{u}}_j},P_{{\mathit{\lambda }}_j}\right)]={{\mathsf{\sigma }}^{\prime }}_{\mathrm{z}}^2\exp \left(-\frac{1}{2}{\left({\mathit{u}}_i-{\mathit{u}}_j\right)}^T{\left({\mathit{W}}_{z,x}^{\prime }\right)}^{-1}\left({\mathit{u}}_i-{\mathit{u}}_j\right)\right).$$

(17)

Equation (17) is still a Gaussian covariance function with ${\mathsf{\sigma }\prime }_{\mathrm{z}}^2={\mathsf{\sigma }}_{\mathrm{z}}^2{\left|\mathit{I}+{\mathit{W}}_{z,x}{}^{-1}\left({\mathsf{\Sigma }}_i+{\mathsf{\Sigma }}_j\right)\right|}^{-\frac{1}{2}}$, $W_{z,x}^{\prime }$ = ${\mathit{W}}_{z,x}+{\mathsf{\Sigma }}_i+{\mathsf{\Sigma }}_j$. Furthermore, the covariance between $z\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right),y_r\left(P_{{\mathit{u}}_j},P_{{\mathit{\lambda }}_j}\right)$ and the covariance between $y_s\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right),y_r\left(P_{{\mathit{u}}_j},P_{{\mathit{\lambda }}_j}\right)$ can be derived as

$$\begin{gathered} \mathrm{Cov}\left[z\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right),y_r\left(P_{{\mathit{u}}_j},P_{\mathit{\lambda }}{}_j\right)\right]=\mathrm{Cov}\left[y_s\left({\mathit{x}}_i,{\mathit{\lambda }}_i\right),y_r\left(P_{{\mathit{u}}_j},P_{\mathit{\lambda }}{}_j\right)\right] \\ ={\sigma }_z^2{E}_{{\mathit{u}}_j~P_{{\mathit{u}}_j}}\left[R_{z,\mathit{x}}\left({\mathit{x}}_i,{\mathit{u}}_j\right)\right]E_{{\mathit{\lambda }}_j~P_{\mathit{\lambda }}{}_j}\left[R_{z,\mathit{\lambda }}\left({\mathit{\lambda }}_i,{\mathit{\lambda }}_j\right)\right], \end{gathered}$$

(18)

where $E_{{\mathit{u}}_j~P_{{\mathit{u}}_j}}\left[R_{z,\mathit{x}}\left({\mathit{x}}_i,{\mathit{u}}_j\right)\right]$ = $\frac{\exp \left(-\frac{1}{2}{\left({\mathit{x}}_i-{\mathit{x}}_j\right)}^T{\left({\mathit{W}}_{z,x}+{\mathsf{\Sigma }}_j\right)}^{-1}\left({\mathit{x}}_i-{\mathit{x}}_j\right)\right)}{{\left|\mathit{I}+{\mathit{W}}_{z,x}{}^{-1}{\mathsf{\Sigma }}_j\right|}^{\frac{1}{2}}}$. $E_{{\mathit{\lambda }}_j~P_{{\mathit{\lambda }}_j}}\left[R_{z,\mathit{\lambda }}\left({\mathit{\lambda }}_i,{\mathit{\lambda }}_j\right)\right]$ can be hard to compute, even when the posterior distribution has an analytical form. In order to generalise our approach, we use a Laplace approximation for the posterior distribution; i.e., we use a Gaussian distribution to approximate the posterior distribution of $\mathit{\lambda }$. Assuming the approximated posterior distribution is $\mathrm{N}\left(\mathit{v},{\mathsf{\Sigma }}_{\lambda }\right)$, then $E_{{\mathit{\lambda }}_j~P_{\mathit{\lambda }}{}_j}\left[R_{z,\mathit{\lambda }}\left({\mathit{\lambda }}_i,{\mathit{\lambda }}_j\right)\right]$ = $\frac{\exp \left(-\frac{1}{2}{\left({\mathit{\lambda }}_i-\mathit{v}\right)}^T{\left({\mathit{W}}_{z,\lambda }+{\mathsf{\Sigma }}_{{\mathit{\lambda }}_j}\right)}^{-1}\left({\mathit{\lambda }}_i-\mathit{v}\right)\right)}{{\left|I+{\mathit{W}}_{z,\lambda }{}^{-1}{\mathsf{\Sigma }}_{\mathit{\lambda }}\right|}^{\frac{1}{2}}}$, and Equation (18) can be computed accordingly.

2.4.3. Predictive Distribution for the True Fuel Consumption $\mathrm{z}\left(\mathrm{x},\mathsf{\lambda }\right)$ with Input Uncertainty

Based on the characteristics of the multivariate normal distribution, the predictive distribution for $\mathrm{z}\left(\mathit{x},\mathit{\lambda }\right)$, the true fuel consumption at exact control input $\mathit{x}$ and exact environmental input $\mathit{\lambda }$ can be derived as

$$\mathrm{z}\left(\mathit{x},\mathit{\lambda }\right)|{\mathit{Y}}_r,P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r,{\mathit{Y}}_s,{\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}}~N(\mathrm{m}\left(\mathit{x},\mathit{\lambda }\right),{\mathsf{\sigma }}^2\left(\mathit{x},\mathit{\lambda }\right))$$

(19)

with

$$\mathrm{m}\left(\mathit{x},\mathit{\lambda }\right)={\beta }_z+\mathsf{\Sigma }{\left(\left(\mathit{x},\mathit{\lambda }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)}^{\mathrm{T}}\mathsf{\Sigma }{\left(\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right],\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)}^{-1}\left(\left[\begin{array}{c}{\mathit{Y}}_r\\ {\mathit{Y}}_s\end{array}\right]-[\begin{array}{c}{\beta }_z{\mathbf{1}}_m\\ ({\beta }_z-{\beta }_{\delta }){\mathbf{1}}_n\end{array}]\right),$$

(20)

$$\begin{array}{l}{\mathsf{\sigma }}^2\left(\mathit{x},\mathit{\lambda }\right)\\ =\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left(\mathit{x},\mathit{\lambda }\right)\right)\\ -\mathsf{\Sigma }{\left(\left(\mathit{x},\mathit{\lambda }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)}^{\mathrm{T}}\mathsf{\Sigma }{\left(\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right],\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)}^{-1}\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right),\end{array}$$

(21)

$$\begin{array}{l}\mathrm{k}\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)\\ =\left[\begin{array}{cc}\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left(\mathit{x},\mathit{\lambda }\right)\right)& \mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)\\ \mathsf{\Sigma }\left(\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right),\left(\mathit{x},\mathit{\lambda }\right)\right)& \mathsf{\Sigma }\left(\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)\end{array}\right]\\ -\left[\begin{array}{c}\mathsf{\Sigma }{\left(\left(\mathit{x},\mathit{\lambda }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)}^T\\ \mathsf{\Sigma }{\left(\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)}^T\end{array}\right]\mathsf{\Sigma }{\left(\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right],\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)}^{-1}\left[\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)\mathsf{\Sigma }\left(\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right),\right]\end{array}$$

(22)

where $\mathrm{m}\left(\mathit{x},\mathit{\lambda }\right)$ denotes the predictive mean for the true fuel consumption $z\left(\mathit{x},\mathit{\lambda }\right)$ and ${\mathsf{\sigma }}^2\left(\mathit{x},\mathit{\lambda }\right)$ denotes the predictive variance for $z\left(\mathit{x},\mathit{\lambda }\right)$. $\mathrm{k}\left(\left(\mathit{x},\mathit{\lambda }\right),\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)\right)$ denotes the predictive covariance between $z\left(\mathit{x},\mathit{\lambda }\right)$ and z$\left({\mathit{x}}^{\prime },{\mathit{\lambda }}^{\prime }\right)$. $\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]$ represents the input settings, including both real system inputs $(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)$ and simulation inputs $({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})$. $\left({\mathit{Y}}_r,{\mathit{Y}}_s\right)$ are the overall outputs, including both the real fuel consumption and simulated fuel consumption, and ${\mathbf{1}}_m$ is an m-dimensional vector of 1, where $m$ is the number of real observations and $n$ is the number of simulation outputs. $\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left(\mathit{x},\mathit{\lambda }\right)\right)={\sigma }_z^2{R}_z\left(\left(\mathit{x},\mathit{\lambda }\right)\right),\left(\mathit{x},\mathit{\lambda }\right))$ is the variance at $z\left(\mathit{x},\mathit{\lambda }\right)$. $\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)$ is the covariance between z$\left(\mathit{x},\mathit{\lambda }\right)$ and $\left[\begin{array}{c}{y}_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ y_s({\mathit{x}}_s,{\mathit{\lambda }}_s)\end{array}\right],$ which can be calculated as

$$\mathsf{\Sigma }\left(\left(\mathit{x},\mathit{\lambda }\right),\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_s)\end{array}\right]\right)=\left[\begin{array}{c}cov\left[z\left(\mathit{x},\mathit{\lambda }\right),y_r\left(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r\right)\right]\\ cov\left[z\left(\mathit{x},\mathit{\lambda }\right),y_s({\mathit{x}}_s,{\mathit{\lambda }}_s)\right]\end{array}\right]=\left[\begin{array}{c}cov\left[z\left(\mathit{x},\mathit{\lambda }\right),y_r\left(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r\right)\right]\\ cov\left[z\left(\mathit{x},\mathit{\lambda }\right),z\left({\mathit{x}}_s,{\mathit{\lambda }}_s\right)\right]\end{array}\right],$$

(23)

where $\mathrm{cov}\left[z\left(\mathit{x},\mathit{\lambda }\right),y_r\left(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r\right)\right]$ can be computed using Equation (18) and $\mathrm{cov}\left[z\left(\mathit{x},\mathit{\lambda }\right),z\left({\mathit{x}}_s,{\mathit{\lambda }}_s\right)\right]={\sigma }_z^2{R}_z\left(\left(\mathit{x},\mathit{\lambda }\right)\right),\left({\mathit{x}}_s,{\mathit{\lambda }}_s\right))$. Similarly, $\mathsf{\Sigma }\left(\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right],\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)$ is the covariance matrix of $\left[\begin{array}{c}{y}_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ y_s({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]$, which can be calculated as

$$\mathsf{\Sigma }\left(\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right],\left[\begin{array}{c}(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\\ ({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\end{array}\right]\right)=\left(\begin{array}{cc}\mathrm{Cov}\left[y_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\right]& \mathrm{Cov}{\left[y_s({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}}),y_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\right]}^T\\ \mathrm{Cov}\left[y_s({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}}),y_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\right]& \mathrm{Cov}\left[y_s({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\right]\end{array}\right),$$

(24)

where $\mathrm{Cov}\left[y_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\right]$ is the $m\times m$ covariance matrix for $y_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r),$ whose elements can be computed using Equation (17). $\mathrm{Cov}\left[y_s({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}}),y_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)\right]$ is the n × m covariance matrix between $y_s({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})$ and $y_r(P_{\mathit{u}}{}_r,P_{\mathit{\lambda }}{}_r)$, whose elements can be computed using Equation (18). $\mathrm{Cov}\left[y_s({\mathit{x}}_s,{\mathit{\lambda }}_{\mathrm{s}})\right]={\sigma }_z^2{R}_z\left(\left({\mathit{x}}_{s,}{\mathit{\lambda }}_s\right),({\mathit{x}}_s,{\mathit{\lambda }}_s)\right)+{\sigma }_{\delta }^2{R}_{\delta }\left(\left({\mathit{x}}_{s,}{\mathit{\lambda }}_s\right),({\mathit{x}}_s,{\mathit{\lambda }}_s)\right)+{\sigma }_{\epsilon }^2{I}_n$ is the $n\times n$ covariance matrix for $y_s({\mathit{x}}_s,{\mathit{\lambda }}_s).{\mathit{I}}_n$ is an $n\times n$ identity matrix.

The MLE estimation method can be used to estimate the parameters. With the mean function in Equation (20) as the covariance function in Equation (22), the resulting predicted distribution for the true fuel consumption $z\left(\mathit{x},\mathit{\lambda }\right)$ now takes both the control input uncertainty and the environmental input uncertainty into account and can be used for the prediction of fuel consumption.

2.4.4. Integrating out the Environmental Input Uncertainty in $\mathrm{z}\left(\mathrm{x},\mathsf{\lambda }\right)$

Using the above predictive distribution for $z\left(\mathit{x},\mathit{\lambda }\right)$ given by Equations (20)–(22), we can evaluate the fuel consumption at any control input setting $\mathit{x}$ given the environmental input setting $\mathit{\lambda }$. Clearly, for different values of $\mathit{\lambda }$, the fuel consumption evaluation may be quite different. Ideally, we want to evaluate the fuel consumption given the true environmental input, but unfortunately, the true value of the environmental input is unknown. As stated above, we use the Bayesian approach to quantify the environmental input uncertainty. The posterior distribution of the environmental input parameter is $p(\mathit{\lambda }|D_h)$ based on the historical electricity demand $D_h$ obtained from on-board measurements. Then, we can use

$$G\left(\mathit{x}\right)=\int z\left(\mathit{x},\mathit{\lambda }\right)p\left(\mathit{\lambda }|D_h\right)d\mathit{\lambda }$$

(25)

as the GP model for the true fuel consumption to hedge against the environmental input uncertainty. According to [22], $G\left(\mathit{x}\right)$ is also a GP, but there is usually no general analytical form for $G\left(\mathit{x}\right)$ due to the integration with respect to the posterior distribution. A Monte Carlo approach can be applied to approximate the mean and covariance of $G\left(\mathit{x}\right)$. Specifically, $N_{MC}$ samples$:\left\{{\mathit{\lambda }}_1,{\mathit{\lambda }}_2,\cdots ,{\mathit{\lambda }}_{N_{MC}}\right\}$ can be generated from $p(\mathit{\lambda }|D_h)$; then,

$$\widehat{G}\left(\mathit{x}\right)=\frac{1}{N_{MC}}{\displaystyle \sum }_{i=1}^{N_{MC}}z\left(\mathit{x},{\mathit{\lambda }}_i\right)$$

(26)

is used as the approximated GP model for $\widehat{G}\left(\mathit{x}\right)$. This $\widehat{G}\left(\mathit{x}\right)$ is also a GP, as it is a finite sum of a GP. Then, the mean and variance of $\widehat{G}\left(\mathit{x}\right)$ can be calculated as

$$\mathrm{E}\left[\widehat{G}\left(\mathit{x}\right)\right]=\frac{1}{{\mathrm{N}}_{\mathrm{MC}}}{\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{MC}}}\mathrm{m}\left(\mathit{x},{\mathit{\lambda }}_i\right),$$

(27)

$$\mathrm{VAR}\left[\widehat{G}\left(\mathit{x}\right)\right]=\frac{1}{{\mathrm{N}}_{\mathrm{MC}}^2}{\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{MC}}}{\displaystyle \sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{MC}}}\mathrm{k}((\mathit{x},{\mathit{\lambda }}_i),\left(\mathit{x},{\mathit{\lambda }}_j\right)).$$

(28)

To hedge against the environmental input uncertainty, the true fuel consumption at any control input location $\mathit{x}$ can be estimated with $\mathrm{E}\left[\widehat{G}\left(\mathit{x}\right)\right]$, and $\mathrm{VAR}\left[\widehat{G}\left(\mathit{x}\right)\right]$ can be used to measure the uncertainty of this estimation.

3. Assessment and Ranking of Mitigation Strategies

In this section, we will illustrate how the proposed GP model can be used to evaluate the fuel consumption and further used to assess the emission reduction and cost for various mitigation strategies.

3.1. Assessment of Emission Reduction

The emission reduction before and after applying the i-th mitigation strategy, denoted as $E{R}_i$, can be computed using the reduction of fuel consumption and the corresponding emission factor as follows:

$$E{R}_i=F{R}_i\times EF=\left(F{C}_i^{before}-F{C}_i^{after}\right)\times EF,$$

(29)

where $EF$ denotes the emission factor, and $F{R}_i$ denotes the reduction of fuel consumption for the i-th strategy. $F{C}_i^{before}$ and $F{C}_i^{after}$ denote the fuel consumption before and after implementing strategy i, respectively. Both $F{C}_i^{before}$ and $F{C}_i^{after}$ can be modelled by the predictive distribution $\widehat{G}\left(\mathit{x}\right)$ based on Equations (27) and (28). However, it is more accurate to use the real ship consumption, if available, to obtain $F{C}_i^{before}$. For the case study in Section 4, $F{C}_i^{before}$ are directly obtained from real ship on-board measurements.

By using the GP metamodel $\widehat{G}\left(\mathit{x}\right)$ to predict the fuel consumption, the control and environmental input uncertainties from the observed real data can be incorporated, and thus, the prediction results are more reliable. The reduction of fuel consumption is then obtained from the difference between the fuel consumption before (directly obtained from on-board measurements) and after (computed from E[$\widehat{G}\left(\mathit{x}\right)$]) a certain strategy is applied.

It is worthwhile to mention that, while the uncertainty in the observed control input due to the signal error from the AIS is incorporated when building the GP model, the prediction $\widehat{G}\left(\mathit{x}\right)$ is the fuel consumption evaluated at the exact control input $\mathit{x}$, which has no uncertainty with respect to $\mathit{x}$. This is because, in the ranking of mitigation strategies, we assume that the control on-board the vessel is completely deterministic; i.e., we can set the control inputs to exact values on-board so as to fairly compare different mitigation strategies. In reality, when there do exist control variability on-board the vessel, the practitioners can adjust the control inputs according to their experiences to hedge against the control variability, but this is out of the scope of our research.

3.2. Assessment of Cost

3.2.1. Cost Components

The cost of implementing a mitigation strategy includes the capital cost/investment cost, operational cost and opportunity cost. In addition, there is a fuel consumption saving cost from implementing a certain mitigation strategy $i$. The total cost for a mitigation strategy $i$ can be represented as follows:

$$T{C}_i=I{C}_i+R{C}_i+O{C}_i-S{C}_i,$$

(30)

where $I{C}_i$ is the annual investment cost, which can be obtained by annualising the initial investment cost over the remaining years. The initial investment cost is the cost invested at time zero, which can be calculated as the net present value of the total capital cost. The remaining years may be the lifetime of the strategy or the remaining lifetime of the ship, whichever is shorter. Let $I_{0,i}$ denote the initial investment for strategy $i$, $T_i$ the remaining years and $r$ the discount rate; then, the annualised investment cost is $I{C}_i=\frac{I_{0,i}·r}{1-1/{\left(1+r\right)}^{T_i}}$. $R{C}_i$ is the annual operational cost or recurring cost when a strategy $i$ is implemented. $O{C}_i$ is the annual opportunity cost due to the loss of service when implementing strategy $i$. $S{C}_i$ is the cost savings due to the fuel consumption saved when implementing strategy $i$. The fuel consumption savings can be obtained as $S{C}_i=\left(F{C}_i^{before}—F{C}_i^{after}\right)·FP=F{R}_i·FP$, where $FP$ is the fuel price. Similarly, $F{C}_i^{before}$ can be directly obtained from real ship on-board measurements if available (otherwise can be evaluated by $\widehat{G}\left(\mathit{x}\right)$), and $F{C}_i^{after}$ is evaluated using the predictive distribution of $\widehat{G}\left(\mathit{x}\right)$ given by Equations (27) and (28) in Section 2.4.4.

3.2.2. Cost Assessment with Uncertainty

There are various uncertainties in assessing the total cost $T{C}_i$ of implementing strategy $i$, including the uncertainties arising from the investment cost $I{C}_i$, the operational cost $R{C}_i$, the opportunity cost $O{C}_i$ and the fuel consumption saving cost $S{C}_i$. In addition to these, there are also uncertainties in some input factors when computing the total cost, such as the discount rate $r$ and the fuel price $FP$. It is important to consider all of these uncertainties in the total cost assessment [34].

To assess the total cost uncertainty, we use a general approach based on the Monte Carlo method [35]. Given the uncertain factors, one straightforward way to quantify these uncertainties is by using their density functions. Based on the distribution of the fuel consumption, discount rate and fuel price, the Monte Carlo method can be applied. This requires random samples to be generated from each of the distributions, and the sampled total cost is then computed. These random samples of the uncertain factors (the discount rate $r$, fuel consumption $F{C}_i^{after}$ and fuel price $FP$) are generated from their corresponding densities to obtain samples of the total cost. The uncertainty (distribution), including the sample mean and sample variance of the total cost $T{C}_i$, can then be assessed from the obtained total cost samples.

3.3. Ranking of Mitigation Strategies with Uncertainty

Given the uncertainties in the total cost and emission reduction for different strategies, the ranking under uncertainty method proposed in [9] can be used to rank the mitigation strategies. Let $T{C}_i$ and $E{R}_i$ denote the cost and emission reduction for strategy $i$ and $T{C}_j$ and $E{R}_j$ denote the cost and emission reduction for strategy $j$. Here, $T{C}_i$ is computed from Equation (30), and $E{R}_i$ is computed from Equation (29). To compare strategies $i$ and $j$, three situations must be considered:

• The first is that strategy $i$ dominates strategy $j$, for which the probability can be obtained as

  $$p\left(i>j\right)=p\left(T{C}_i<T{C}_j\right)\cap p\left(E{R}_i>E{R}_j\right),$$

  (31)

  which is the probability of strategy $i$ with less cost and larger emission reductions compared to strategy $j$.
• The second situation is that strategy $j$ dominates strategy $i$, for which the probability can be obtained as

  $$p\left(i<j\right)=p\left(T{C}_i>T{C}_j\right)\cap p\left(E{R}_i<E{R}_j\right),$$

  (32)

  which is the probability of strategy $j$ with less cost and larger emission reductions compared to strategy $i$.
• The third situation is nondominance between strategy $i$ and strategy $j$, for which the probability can be calculated as

  $$p\left(i\sim j\right)=1-p\left(i>j\right)-p\left(i<j\right).$$

  (33)

A further comparison of strategies $i$ and $j$ needs to be carried out in the nondominated case. In this situation, the MCE criterion proposed in [9,36] can be used. The MCE can be obtained as

$$MC{E}_{ij}=\frac{T{C}_i-T{C}_j}{E{R}_i-E{R}_j}$$

(34)

Given the uncertainties, the total probability of strategy $i$ being better than strategy $j$, also called preference probability (pp) of strategy $i$ over strategy $j$, can be obtained by

$$p\left(i\succ j\right)=p\left(i>j\right)+p\left(i\sim j\right)p\left(MC{E}_{ij}\le C{E}_0\right),$$

(35)

where $p\left(MC{E}_{ij}\le C{E}_0\right)$ is the probability of strategy $i$ being preferable to strategy $j$ in the nondominance situation. $C{E}_0$ denotes the accepted additional cost for one additional unit of emission reduction. Similarly, the probability of strategy $j$ being better than strategy $i$ can be obtained as

$$p\left(i\prec j\right)=p\left(i<j\right)+p\left(i\sim j\right)p(MC{E}_{ij}>C{E}_0).$$

(36)

All of these probability computations can be approximated by the Monte Carlo method described in Section 3.2.2.

4. Case Study

A chemical tanker is taken as an example here in order to evaluate several different mitigation strategies. The emission reductions for several mitigation strategies are evaluated using the developed GP metamodel based on the fuel consumption reductions. The predictive performances of the model are compared, with and without taking input uncertainty into consideration. An assessment of the cost-effectiveness of these different mitigation strategies is also carried out.

4.1. Mitigation Strategies

The chemical tanker studied here was built in 2014. It has a double hull, and the length and width of the ship are 181 and 31.3 m, respectively. The maximum capacity of the ship is 51,000 m³, and the maximum draft of the ship is 12.4 m. It has two main engines and two auxiliary engines. Based on the ship’s energy system and the available data, five mitigation strategies are evaluated: speed reduction (10%), trim optimisation, draft optimisation, autopilot adjustment and speed control of pumps and fans. The distribution of the electricity demand of this ship is also considered in this study and is treated as an environmental input. The corresponding control inputs for these mitigation strategies are shown in

Table 2. Control inputs for different mitigation strategies.

Mitigation Strategy	Control Input
Speed reduction (10%)	Vessel speed
Trim optimisation	Trim value
Draft optimisation	Draft value
Autopilot adjustment	Course
Speed control of pumps and fans	Speed of pumps and fans

.

From the real fuel consumption data collected from real ship on-board measurements and the simulated fuel consumption data obtained from a Simulink simulator, as in [19], a GP model can be developed (as described in Section 2.4) for all the strategies. As we consider five strategies, the GP model has control input $\mathit{x}\in {\mathrm{R}}^5$ and environmental input $\mathit{\lambda }\in {\mathrm{R}}^1$. Speed reduction has been comprehensively studied in recent years and is an important mitigation strategy, as it can provide large emission reductions. The recorded speed in the AIS data is the SOG, and this is taken as the input of the energy system. Here, the relationship between the SOG and the engine speed (the control input of the simulation model) is translated using an equation adopted from [37]. Here, a speed reduction of 10% from the design speed of the ship is considered. The reduction of 10% is acceptable by most ships in practice. However, the potential for speed reduction is limited. This is because the engines can be damaged if the ship is sailing in off-design conditions. The minimum load of the engine depends on the technical specification of each individual engine. Based on the developed GP model, the fuel consumption can be predicted for different vessel speeds. Trim and draft are two further control inputs that can have impacts on the fuel consumption, as both of them can influence the ship’s resistance. Again, with the developed GP model, the fuel consumptions for different values of trim and draft can be predicted, and the emission reductions can be evaluated for the optimal values. Autopilot adjustment can keep the ship on course and can save on fuel consumption by preventing the unnecessary use of the rudder. The course over ground (COG) is taken as the control input for this mitigation strategy, and the emission reduction for the optimal COG is evaluated. Speed control of the pumps and fans is another mitigation strategy considered in this paper, in which the speed of the pumps and fans is the corresponding control input, and the fuel consumption for different speeds is evaluated.

4.2. Data and Assumptions

Data are collected for both the emission reduction evaluation and the cost assessment. As explained in Section 2.2, there are two types of data: real data and simulation data. For real data, real fuel consumptions (outputs) are directly observed from on-board measurements. The control inputs such as SOG, trim, draft and COG for real data are obtained from the AIS, while speeds of the pumps and fans are observed from on-board measurements. The historical electricity demand is also observed from on-board measurements. The available on-board measurements and AIS data are from January 2017 to March 2018. For the simulation data, the simulated fuel consumptions (outputs), control inputs and environmental inputs are all obtained from the simulation model through experimental design. The simulation model is built using Simulink to represent the ship energy system, as in [18]. The relationship between the inputs (e.g., engine speed, trim and draft) and the output (i.e., fuel consumption) in the simulation model is either based on the physical process or the regression model. Based on the real data and simulation data, the GP metamodel can be developed, and the total emission reduction for different mitigation strategies can be predicted using this metamodel.

In the cost assessment, the total costs of different mitigation strategies are assessed. The total costs include nonrecurring costs, recurring costs and opportunity costs. Here, the estimated values of low and high nonrecurring and recurring costs are adopted from the IMO MEPC 62 report [3]. The nonrecurring costs of speed reduction also includes the additional cost arising from the fleet increase needed to transport in time to meet the demand, where the cost of building a new ship is estimated based on the Review of Maritime Transport series [38]. The lifetime of a new ship is assumed to be 30 years. The reference discount rate is estimated using the weighted average cost of capital over 12 major shipping companies, which is 7%. Opportunity costs occur due to the loss of service when the mitigation strategy is implemented. Most operational strategies do not have this opportunity cost, as they usually do not require extra dry docking days for updating. However, some technical strategies may have this cost, such as the speed control of pumps and fans [3]. Here, the opportunity cost is obtained by multiplying the term-chartered rate and the expected dry docking days for updating. The term-chartered rate is computed using the new build cost divided by the vessel’s life [38]. The extra dry docking days are adopted from the IMO MEPC 62 report [3].

When mitigation strategies are implemented, there are also fuel consumption saving costs due to the reduction in fuel consumption. The total fuel consumption reduction and the fuel price are required in order to compute this cost. The total fuel consumption reduction can be predicted using the developed GP model, and the fuel price is adopted from the IMO MEPC 62 report [3], where the low and high estimates for fuel price are US$500 and US$900, respectively. These prices are close to the estimates provided in the Annual Energy Outlook [39].

To estimate the distribution of the electricity demand, we assume the distribution of the electricity is exponential with rate parameter $\lambda$ and use the noninformative Jeffreys prior $p\left(\lambda \right)\propto 1/\mathit{\lambda }$ for the rate parameter. The posterior is calculated based on the recorded electricity demand in 2017 and the first quarter of 2018, adopted from on-board measurements.

4.3. Results and Discussion

The developed GP metamodel is first validated before it is applied for further analysis. The data obtained in 2017 is used for model development, and the data obtained in the first quarter of 2018 is used for model validation. The predictive performance of the two models is compared, i.e., with and without input uncertainty. The mean and 95% confidence interval (CI) of the predicted versus observed fuel consumptions for both training and validation data are shown in Figure 2 and Figure 3. The absolute error between the predicted and observed fuel consumptions versus the predictive error (half of the 95% CI) for both the training and validation data are also shown in Figure 2 and Figure 3. The figures indicate that the 95% CI of the predicted fuel consumptions can cover the observed fuel consumptions most of the time in the model with uncertainty, and that it is better than the model without uncertainty. The widths of the CIs for the model with input uncertainty are larger than those of the model without input uncertainty. The coverage rate indicates the number of real observed values that are within the 95% CI of the predicted values. The coverage rate of the 95% CI of the real observations for both models is also given in

Table 3. Coverage rate that the 95% confidence interval (CI) covers the observed value and the root mean square error (RMSE) for different models.

Model	Coverage Rate		RMSE
	Training Data	Validation Data	Training Data	Validation Data
Model without uncertainty	0.702	0.625	0.1496	0.2241
Model with uncertainty	0.992	0.937	0.1445	0.2173

, and it can be seen that the performance of the model with uncertainty is better, as its coverage rate is closer to the nominal 95% level as compared with the model without uncertainty. The root mean square errors (RMSEs) for both models are also computed in Table 3. The results show that the RMSE for the model with uncertainty is smaller than for the model without uncertainty, indicating that the former has better accuracy. In addition, a goodness-of-fit test is used to test the difference between the predicted and the observed values, and the results indicate that the difference is not significant at an α level of 0.05. The RMSEs are accepted for accuracy. Hence, the developed GP metamodel is valid and can accurately account for uncertainties in this case study. Furthermore, the GP metamodel with input uncertainty can provide more reliable and accurate results, especially when the effects of input uncertainty cannot be ignored.

The validated GP metamodel-incorporated input uncertainty is then used to evaluate the fuel consumption reduction for different mitigation strategies. The estimated annual fuel consumption for the chemical tanker in 2017 was 2846 metric tons (MT). The predicted annual fuel consumption reduction, its predictive variance and the corresponding mean percentage of the fuel reduction with its variance for different mitigation strategies are given in

Table 4. Annual fuel consumption reduction, emission reduction and percentage fuel reduction for different mitigation strategies.

Mitigation strategy	Speed Reduction (10%)	Trim Optimisation	Draft Optimisation	Autopilot Adjustment	Speed Control of Pumps and Fans
Annual fuel consumption reduction (MT), mean (variance)	523.47 (50.32)	48.56 (6.81)	46.76 (6.18)	33.14 (5.25)	18.08 (2.27)
Annual emission reduction (MT), mean (variance)	1669.90 (487.95)	151.22 (66.04)	145.61 (59.93)	103.20 (50.91)	56.30 (22.01)
Percentage of fuel reduction, mean (variance)	18.39% (1.77%)	1.71% (0.24%)	1.64% (0.22%)	1.16% (0.18%)	0.64% (0.08%)

. It can be seen that a speed reduction of 10% has the largest fuel consumption reduction (523.47 MT), which accounts for 18.39% of the annual fuel consumption. Compared to speed reduction, the other four mitigation strategies have lower fuel consumption reductions of less than 2%. This result is reasonable and is as expected, since the speed reduction has a significant influence on the fuel consumption. The emission reduction for different mitigation strategies can also be computed by multiplying the fuel consumption reduction and the emission factor. Here, the emission factor adopted from the third IMO GHG study for heavy fuel oil is 3.114 g/g fuel. In terms of both fuel consumption reduction and emission reduction, speed reduction is therefore the best choice. Compared to the results given in [18], where the input uncertainty is not considered, it can be found that the variances of the fuel consumption reduction and the emission reduction are larger using the GP metamodel-incorporated input uncertainty. This is expected as more uncertainties are taken into account. The results with input uncertainty can also provide more robust comparisons among different mitigation strategies.

In addition to the emission reduction evaluation, the cost of the different mitigation strategies is also assessed, where the cost data are described in Section 4.2. The 95% CIs of the assessed costs for the different mitigation strategies are given in

Table 5. Costs for different mitigation strategies.

Mitigation Strategy	Annual cost (US$)
	Investment Cost (min, max)	Operational Cost (min, max)	Opportunity Cost (min, max)	Cost Saved from Fuel Reduction (−95%, +95%)	Total Cost (min, max)
Speed reduction (10%)	(184,171, 225,097)	(80,289, 98,131)	(0, 0)	(367,109, 461,383)	(−43,881, −196,923)
Trim optimisation	(2392, 2924)	(806, 986)	(0, 0)	(22,486, 30,422)	(−18,577, −27,224)
Draft optimisation	(2146, 2687)	(753, 942)	(0, 0)	(24,973, 29,347)	(−20,169, −26,846)
Autopilot adjustment	(3339, 4081)	(0, 0)	(0, 0)	(16,899, 22,863)	(−12,818, −19,524)
Speed control of pumps and fans	(1396, 1706)	(0, 0)	(695, 849)	(8177, 11,063)	(−5622, −8972)

.

Based on the emission reduction and cost assessment, the mitigation strategies are then ranked according to different criteria, including the total emission reduction, cost, MCE without uncertainty and MCE with uncertainty. The ranking results are given in

Table 6. Ranking of the mitigation strategies. MCE: marginal cost-effectiveness.

Mitigation Strategy	Mean Emission Reduction (MT)	Rank	Mean Cost (US$)	Rank	MCE without Uncertainty	Rank	MCE with Uncertainty (pp)	Rank
Speed reduction (10%)	1669.90	1	−101,983	1	−52.20	1	1	1
Trim optimisation	151.22	2	−22,715	3	2.50	2	0.4816	3
Draft optimisation	145.61	3	−22,729	2	−148.36	3	0.9987	2
Autopilot adjustment	103.20	4	−16,437	4	−202.11	4	0.9938	4
Speed control of pumps and fans	56.30	5	−6,958	5	-	5	-	5

. It can be seen that the results are different for different criteria; for example, the trim optimisation is better than the draft optimisation in terms of emission reduction, although it is worse in terms of cost. It can also be seen that the ranking results may be different when the input uncertainty is taken into account. When the mitigation strategies are ranked using the MCE without considering the input uncertainty (see Section 3), the trim optimisation is better than the draft optimisation; when the input uncertainty is taken into account, the draft optimisation becomes better. It is therefore important to consider the input uncertainty in order to provide more reliable ranking results. The MCE with uncertainty ranks the speed reduction (10%) as the most effective mitigation strategy, followed by the draft optimisation and trim optimisation. Autopilot adjustment and the speed control of pumps and fans are ranked last. In summary, mitigation strategies may have different rankings using different criteria, such as emission reduction, cost and MCE. Therefore, which mitigation strategy is better depends on the priority of the target. In addition, the input uncertainty may have an impact on prioritising the mitigation strategies. This uncertainty is usually ignored in practice, which may result in suboptimal choices in mitigation strategies. To provide more reliable choices, it is better to consider the input uncertainty in mitigation strategies’ ranking and selection.

5. Conclusions

In this paper, a Gaussian process metamodel is developed to integrate real-world data and simulation data for evaluation of the emission reduction. As an important improvement to the existing methods in the literature, the uncertainties in the various sources of data are explicitly accounted for in the model. Furthermore, it is shown that neglecting such uncertainties by existing methods can lead to errors in evaluating mitigation strategies for emission reduction. The proposed improvement can significantly improve the reliability and robustness of GP-based metamodels in evaluating mitigation strategies for emission reductions. In turn, the method can improve the reliability in the ranking of mitigation strategies for better-informed decision-making. Specifically, in this paper, we show how this model can be applied into the marginal cost-effectiveness (MCE) approach for the ranking of mitigation strategies.

A case study of a chemical tanker is developed to demonstrate the advantages of the proposed improvements. The case study is developed based on multiple data sources, including both real and simulated data. The results show that the proposed GP model with input uncertainties has a better coverage rate and a smaller root mean square error than the model that does not account for these uncertainties. This highlights the importance of considering uncertainties in the model formulation. Similar to the GP modelling results on emission reductions, accounting for uncertainties in the model formulation also leads to different ranking results for the MCE criteria as compared to not accounting for uncertainties. The observed differences in results further demonstrate the risks of ignoring the uncertainties in the ranking of mitigation strategies by GP metamodels.

The proposed method can be particularly useful when attempting to evaluate mitigation strategies under different types and values of control inputs so as to choose the optimal one. Findings of this study can benefit both the policymakers and business owners in identifying the most cost-effective strategies in attaining decarbonisation goals and pathways for the shipping industry. The proposed method can provide an express alternative to expensive experimental investigations and time-consuming simulation exercises.

In this paper, all the considered control inputs are continuous variables. A possible direction for future work could be to consider strategies with discrete control inputs in a ship’s energy system. This direction would require further development of the proposed GP model’s mathematical formulations. In addition, only one specific type of ship’s energy system is considered in this paper. The other possible direction for future work could be to examine different types of ships so as to derive generalised conclusions and insights on mitigation strategies for the shipping industry.

In addition, the focus of the paper is to evaluate strategies with the use of fuel oil. In order to more efficiently reduce the emissions from international shipping, more advanced technology other than the traditional mitigation strategies will need to be implemented. For example, the designation of alternative propulsion systems, the use of renewable energy sources and the improvement of the cargo management system would have great potential in reducing emissions. The proposed method may also be applied to compare the performance of these alternative strategies in reducing emissions, and this is worth investigating in the future.