Distributionally Robust Optimal Scheduling of Hybrid Ship Microgrids Considering Uncertain Wind and Wave Conditions

Abstract

A hybrid ship uses integrated generators, an energy storage system (ESS), and photovoltaics (PV) to match its propulsion and service loads, and together with optimal power and voyage scheduling, this can lead to a substantial improvement in ship operation cost, ensuring compliance with the environmental constraints and enhancing ship sustainability. During the operation, significant uncertainties such as waves, wind, and PV result in considerable speed loss, which may lead to voyage delays and operation cost increases. To address this issue, a distributionally robust optimization (DRO) model is proposed to schedule power generation and voyage. The problem is decoupled into a bi-level optimization model, the slave level can be solved directly by commercial solvers, the master level is further formulated as a two-stage DRO model, and linear decision rules and column and constraint generation algorithms are adopted to solve the model. The algorithm aims at minimizing the operation cost, limiting greenhouse gas (GHG) emissions, and satisfying the technical and operational constraints considering the uncertainty. Extensive simulations demonstrate that the expected total cost under the worst-case distribution is minimized, and compared with the conventional robust optimization methods, some distribution information can be incorporated into the ambiguity sets to generate fewer conservative results. This method can fully ensure the on-time arrival of hybrid ships in various uncertain scenarios while achieving expected operation cost minimization and limiting greenhouse gas (GHG) emissions.

1. Introduction

About 3–5% of global greenhouse gas (GHG) emissions are produced by marine transport [1,2]. If no mitigation measures are taken, GHG emissions from marine transport will rise to 8% of global GHG emissions by 2050 [3,4], so marine transport needs to reduce GHG emissions significantly. To meet the demand for emission reduction, hybrid ships have gained industry attention. Hybrid ships refer to ships driven by different types of engines or electric motors and contain multiple energy sources [5]. In the past, hybrid ships were usually referred to as diesel–electric hybrid ships. With the development of new energy technology, new energy technology represented by photovoltaics (PV) has begun to be applied on ships, and PV has been used when carrying out certain research studies on hybrid ships composed of diesel generators (DGs), energy storage systems (ESSs), and PV [6,7]. In [8], it is verified by the results of actual navigation experiments that the participation of photovoltaic power generation in reasonable scheduling can help the shipping industry meet the “low-carbon navigation” policy.

Energy management is the key to achieving economic objectives and emission reduction on the premise of ensuring power performance. In [9], a deterministic energy management approach is proposed for a cruising ship, to minimize the fuel consumption cost by effectively scheduling the power generation units and voyage. In [10], a coordinated optimal power management method is studied according to the different objectives. The management model is solved by particle swarm optimization. The system is initialized with a population of random solutions and searches for optima by updating generations. At each time step, the velocity of each particle is changed as it heads toward its p-best and l-best locations. Ref. [11] proposes a scheduling model to optimize the economic and environmental objectives for the ship integrated with hybrid BS and energy storage.

However, due to the influence of weather, photovoltaic power generation has fluctuations and uncertainties. Compared to land, ships may encounter more uncertainties when navigating at sea, such as ship pitching and rolling, which greatly increases the prediction error of photovoltaic power generation and affects the accuracy of overall power generation and navigation scheduling. Based on this, in recent years, stochastic optimization and robust optimization have been proposed to deal with uncertain power grid optimization problems. In [12], a stochastic optimization method is adopted to solve the coordinated operation optimization problem of a multiple-energy microgrid under various uncertainties; however, many scenarios need to be generated according to the probability density function in advance, which increases the calculation burden. Compared with stochastic optimization, robust optimization only needs to construct an uncertain set and does not need the probability distribution of uncertain variables which are difficult to obtain, and the solution of robust optimization is robust to all possible scenarios. In [13], a robust optimization method is adopted to address the uncertainty of photovoltaic power generation on a ship to minimize the operating cost under the worst uncertainty scenario while satisfying all the operational constraints.

In addition, most existing studies have not considered the speed loss caused by water waves or wind during ship navigation. Although the speed loss models for different ships may be diverse, various tests have shown that the speed loss induced by water and waves for most ship prototypes is significant and may usually exceed 10% of the nominal speed even in moderate weather conditions. Since the propulsion load has a cubic relationship with the cruising speed, the propulsion power loss led by speed loss is much more significant. In this view, speed loss is always an important issue, and in conventional maritime applications, this problem is referred to as “weather routing” [14]. Therefore, in [14], a robust optimization method is proposed for the joint scheduling of power generation and voyage of a ship considering the uncertain wind and wave, and the results show that the proposed robust optimization model can ensure the punctuality of ships in various uncertain situations. In [15], hydrogen fuel cells are integrated into ship microgrids, and a two-stage robust optimization method is proposed. A two-level column constraint generation algorithm was used to solve the optimization model, aiming to minimize operating costs for the entire voyage in the worst-case scenario.

The above robust optimization methods only consider the worst-case scenario where uncertainty occurs. However, the probability of uncertainty occurring at the extremes is small, so the results of robust optimization are too conservative and often require higher economic costs. Combining the characteristics of stochastic and robust optimization, distributionally robust optimization (DRO) has been proposed in recent years to solve the uncertainty problem [16,17]. DRO is an effective way to handle diverse uncertainties, which incorporates the available probability distribution information into an ambiguity set to characterize the true probability distribution of uncertainties and reduce the solution conservatism. At present, this method has not been used to deal with the uncertainty of ship power systems, which is worthy of further study.

This paper comprehensively considers the uncertainties of photovoltaic output power and wind and wave events and proposes a DRO method to jointly optimize the scheduling of power and voyage speed to minimize the ship operation cost.

This paper contributes to this field of research in the following aspects:

(1)

To address diverse uncertainties from photovoltaic systems and wind and wave events, a distributionally robust optimization (DRO) model is proposed to schedule power generation and voyage. With the proposed method, ships can arrive at each port on time, while ensuring a lower operation cost compared to existing methods.

(2)

The original model is decoupled into a bi-level optimization model, the slave level can be solved directly by commercial solvers, the master level is further formulated as a two-stage DRO framework, and linear decision rules are adopted to solve the model, which is suitable for practical applications.

2. Deterministic Optimization Model

As shown in Figure 1, a hybrid ship microgrid consists of DGs, ESSs, and PV for propulsion and service loads.

2.1. Objective

In the optimal scheduling model constructed in this paper, it is assumed that the operation cost of PV and ESSs is zero, and the objective function includes the start–stop cost and the operation cost of the DG units, and the expression is shown as follows:

$$\min (F1+F2)$$

(1)

where $F1$ is the startup and shutdown cost of the DG units; and $F2$ is the operation cost of the DG units.

$F1$ and $F2$ are calculated in (2) and (3), respectively:

$$F1={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i=1}^N(C_t^{ion}}+}{C}_t^{ioff})$$

(2)

$$F2={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i=1}^N{C}_t^i}}$$

(3)

$$C_t^i=a_{0i}{P}_t^{i2}+a_{1i}{P}_t^i+a_{2i}$$

(4)

where $T$ is the time period during navigation; $t$ is the time interval; $N$ is the number of onboard DGs; $C_t^{ion}$ denotes the startup cost of the $i$th generator; $C_t^{ioff}$ denotes the cost of the shutdown of the $i$th DG; $C_t^i$ denotes the fuel consumption cost of the $i$th generator; $a_{0i}$, $a_{1i}$, and $a_{2i}$ are cost coefficients for the $i$th generator, respectively; and $P_t^i$ is the output power of the $i$th DG.

For the convenience of subsequent solving, segmental linearization is used to transform the fuel consumption cost function into a linear function in this paper. By dividing the $i$th DG to the output $P_t^i$ at the $t$th time interval into $K^i$ segments, the fuel consumption cost of the DG can be approximated as a linear function in each of the $k$ segments. The expression is shown in Equation (5):

$$\left\{\begin{array}{c}{C}_t^i=C_0^i+{\displaystyle \sum _{k=1}^{K^i}{c}_{ik}{P}_t^{ik}}\\ P_t^i=P_{\min }^i+{\displaystyle \sum _{k=1}^{K^i}{P}_t^{ik}}\\ 0\le P_t^{ik}\le (P_t^i-P_{\min }^i)\end{array}\right.$$

(5)

where $c_{ik}$ is the generation cost coefficient of the $i$th DG after linearization in section $k$; $P_t^{ik}$ is the power increment of the $i$th DG in section $k$; and $C_0^i$ is the $i$th DG cost at minimum output power $P_{\min }^i$, which is calculated by the SFC curve and fuel price.

Since the start–stop state of the DG is a discrete variable and the on/off s cost of the unit at each time interval depends on the on/off s state, the start–stop cost is also a discrete function that needs to be linearized as follows:

$$\left\{\begin{array}{l}{C}_t^{ion}\ge C_{on}^i(o_{t+1}^i-o_t^i)\\ C_t^{ion}\ge 0\\ C_t^{ioff}\ge C_{off}^i(o_t^i-o_{t+1}^i)\\ C_t^{ioff}\ge 0\end{array}\right.$$

(6)

where $o_t^i$ is the start–stop state of the $i$th DG at time $t$, with 0 indicating that the generator is in the shutdown state and 1 indicating that the generator is in the startup state.

2.2. Constraints

2.2.1. Generation Constraints

• Power Balance

In a power system, the total power output of the generating unit must be equal to the load power, i.e., the power balance is maintained as follows:

$${\displaystyle \sum _{i=1}^N{P}_t^i+}{P}_t^{Ess}+P_t^{Pv}=P_t^{lp}+P_t^{ls}$$

(7)

where $P_t^{Ess}$ is the charging/discharging power of the storage battery $t$; $P_t^{Pv}$ is the photovoltaic power; $P_t^{lp}$ is the load power of the electric propulsion device; and $P_t^{ls}$ is the load power of the electric service device.

The relationship between ship speed and propulsive load power is as follows:

$$P_t^{lp}=c_1{v}^{c_2}$$

(8)

where $c_1$ and $c_2$ are the proportional and exponential coefficients, respectively; generation scheduling and voyage scheduling are coupled by Equation (8).

2.

Generator

The power of DGs is limited by its lower and upper limits, as shown in (9):

$$o_t^i{P}_{\min }^i\le P_t^i\le o_t^i{P}_{\max }^i$$

(9)

where $P_{\min }^i$ is the minimum power; and $P_{\max }^i$ is the rated output power.

The maximal/minimal ramp rates of DGs are limited, as shown in (10):

$$\left\{\begin{array}{c}{P}_{t+1}^i-P_t^i\le R_{\max }^i\\ P_t^i-P_{t+1}^i\le D_{\max }^i\end{array}\right.$$

(10)

3.

ESS

State of charge (SOC) refers to the proportion of the remaining capacity to the total capacity, and SOC at any moment is determined by the remaining capacity of the battery at the previous moment and the charging/discharging power, so the dynamic model of the SOC of ESSs at each time step can be described by (12) and (13):

$$SO{C}_t^{Ess}={\displaystyle \frac{E_t^{Ess}}{E^{Ess}}}$$

(11)

$$E_t^{Ess}=\left\{\begin{array}{c}{E}_t^{Ess}-P_t^{Ess}/{\eta }_{ch}\mathsf{\Delta }tt\ne 1,P_t^{Ess}<0\\ E_t^{Ess}-P_t^{Ess}/{\eta }_{dis}\mathsf{\Delta }tt\ne 1,P_t^{Ess}\ge 0\end{array}\right.$$

(12)

$E_t^{Ess}$ and $E^{Ess}$ are the remaining capacity and rated capacity of the ESS at the t-th moment, respectively; $P_t^{Ess}$ is the ESS output power at the t − 1 moment; and η_ch and η_dis are charging/discharging efficiencies.

The SOC limits during the voyage are expressed as follows:

$$SO{C}_{\min }\le SO{C}_t\le SO{C}_{\max }$$

(13)

where $SO{C}_{\min }$ and $SO{C}_{\max }$ represent the minimal and maximal $SOC$s of an $ESS$, respectively.

The power of an $ESS$ is limited by its lower and upper limits, as shown in (14):

$$\left\{\begin{array}{c}{P}_{ch.\max }<P_{t-1}^{Ess}<0P_{t-1}^{Ess}<0\\ 0\le P_{t-1}^{Ess}<P_{dis.\max }{P}_{t-1}^{Ess}\ge 0\end{array}\right.$$

(14)

where $P_{ch.\max }$ indicates the maximum charging power of the battery and $P_{dis.\max }$ indicates the maximum discharging power of the battery.

4.

Minimum On/Off Time

Once a DG is turned off (or on), it must remain in the same state for a time interval until the next time interval when it is turned on (or off); this is described by (15):

$$t_2-t_1\ge T_{on/off}$$

(15)

where t₂ and t₁ are any two adjacent on/off state switching (on to off or off to on) times.

5.

EEOI

To improve the energy efficiency of ship operation and reduce emissions, EEOI is limited during ship operation as shown in (16) and (17):

$$EEOI={\displaystyle {\sum }_{i\in N}{\displaystyle {\sum }_{t\in T}\left(\frac{o_i^t·G_{DG}^{i,t}·\mathsf{\Delta }t}{m_{AES}Dis{t}_t}\right)}}$$

(16)

$$EEOI\le EEO{I}_{\max }$$

(17)

where $G_{DG}^{i,t}$ is the carbon dioxide emissions during the operation of the DG, for which the formula is shown in Equation (18); $m_{AES}$ is the weight of the ship cargo, in knots; $Dis{t}_t$ is nautical miles of voyage; and $EEO{I}_{\max }$ is the maximum allowed.

$$G_{DG}^{i,t}(t)=g_{0,i}{(r_{DG}^{i,t}(t))}^2+g_{1,i}(r_{DG}^{i,t}(t))+g_{2,i}$$

(18)

Here, $g_{1,i}$, $g_{2,i}$, and $g_{3,i}$ are the calculation factors for CO₂ emissions from DGs.

6.

Action Spinning Reserves Constraints

To ensure the stable operation of the ship’s power system, it is necessary to consider the active spinning reserve; the active spinning reserve mathematical model for DGs and $ESS$s is shown in (19) and (20), and the total active spinning reserve is shown in Equation (21):

$$R_{DG}(t)={\displaystyle {\sum }_{n\in N}(P_{DG,\max }-}{P}_{DG,i}(t)),\forall t\in T$$

(19)

$$R_{ESS}(t)=\left\{\begin{array}{c}{P}_{ESS}^{\max }-P_{ESS}(t),P_{ESS}(t)\ge 0\\ P_{ESS}^{\max },P_{ESS}(t)<0\end{array},\forall t\in T\right.$$

(20)

$$\left\{\begin{array}{c}{R}_R(t)=R_{DG}(t)+R_{ESS}(t)\\ R_R(t)\ge {\delta }_R·N_{DG}·P_{DG,\max }\end{array},\forall t\in T\right.$$

(21)

where $R_{DG}$, $R_{ESS}$, and $R_R$ are the active spinning reserves of DGs and ESSs and the total active spinning reserve, respectively; ${\delta }_R$ is the spinning reserve coefficient; and $N_{DG}$ is the number of diesel generators placed.

2.2.2. Voyage Constraints

The ships sail between multiple ports, which contain a variety of operating conditions such as berthing, cruising, and docking, with different speed requirements, and the speed constraints are shown in Figure 2 [14]. The black arrows indicate the speed limits; in each interval, the cruising speeds should be within the upper and lower bounds. The blue arrows indicate the voyage distance limit, meaning that the ship must arrive at the next port within the specified time from the previous port. Equations (22)–(24) are the cruising speed constraints during the cruising intervals and the partial speed intervals.

$$(1-{\delta }_v)v_n\le v_t\le (1+{\delta }_v)v_n\begin{array}{c}\forall t\in T_c\end{array}$$

(22)

$${\eta }_v(1-{\delta }_v)v_n\le v_t\le {\eta }_v(1+{\delta }_v)v_n\begin{array}{c}\forall t\in T_d\end{array}$$

(23)

$$v_t=0\begin{array}{c}\forall t\in T_b\end{array}$$

(24)

Here, $T_c$ is the cruising intervals, $T_d$ is the partial speed interval, $T_b$ is the berthing interval; $v_n$ is the rated speed; ${\delta }_v$ is the rate of change of $v_n$, and when the ship is cruising, the ship speed is shown in (22); and ${\eta }_v$ is the ratio between $v_t$ and $v_n$, and when the ship is sailing at a partial speed, the constraint is shown in (23).

Since the ship needs to arrive at the port on time, the voyage distance to the intermediate and terminal ports is required, respectively, as shown in (25) and (26):

$$(1-{\delta }_{mid})Dis{t}_n^t\le Dis{t}^t\le (1+{\delta }_{mid})Dis{t}_n^t$$

(25)

$$Dis{t}_n^{nT}\le Dis{t}^{nT}\le (1+{\delta }_{term})Dis{t}_n^{nT}$$

(26)

where $Dist$ is the actual sailing distance of the ship; and ${\delta }_{mid}$ and ${\delta }_{term}$ are the maximum permissible distance error for intermediate and terminal ports, respectively.

2.3. Wave and Wind Resistance

Wind resistance and wave resistance, combined with the calm water resistance of the ship, can be used to obtain the speed loss of the ship in different sea conditions.

Calm water resistance is mainly composed of frictional resistance $R_f$, residual resistance $R_r$, and attached resistance $R_{ap}$ [18,19], as shown in (27):

$$R_T=R_f+R_r+R_{ap}$$

(27)

The International Organization for Standardization (ISO) has constructed a mathematical model for calculating wind resistance and wave resistance based on the data obtained from wind tunnel model tests [20,21], as follows:

$$R_{wind}={\displaystyle \frac{1}{2}}{\rho }_{air}{S}_t{C}_{air}[{(v_t+v_w\cos {\theta }_t)}^2-v_t^2]$$

(28)

where $R_{wind}$ is the wind resistance; ${\rho }_{air}$ is the air density, kg/m³; $S_t$ is the lateral projected area above the waterline, m²; $C_{air}$ is the wind resistance coefficient; $v_w$ is the wind speed, m/s; and ${\theta }_t$ is the wind direction.

The wave resistance is calculated using (29):

$$R_{wave}={\displaystyle \frac{1}{L_s}}{\rho }_{water}g{h}_t^2{B}_{\mathrm{int}}^2{C}_{D,wat}({\tau }_t,{\theta }_t)$$

(29)

where $g$ is the acceleration of gravity, m/s²; $h_t$ is the wave height, m; $B_{\mathrm{int}}$ is the ship’s width, m; and $C_{D,wat}$ is the wave resistance coefficient, with the wavelength ${\tau }_t$ and the wind direction ${\theta }_t$. When the ship is in a certain sea area, the wind direction and wave direction can be approximated to be the same, as both are ${\theta }_t$.

It should be noted that the variables $v_w$ and $h_t$ also depend on the season and geographic position. For example, in the South China Sea region, $v_w$ may be large and $h_t$ may be tall in winter, thus $R_{wind}$ and $R_{wave}$ are large.

The total resistance ($R_c$) due to wind and waves can be estimated using Equation (30):

$$R_c=R_T+R_{wind}+R_{wave}$$

(30)

The ship speed loss is estimated based on the assumption that the required power $P_t$ at the ship speed in calm water is the same as the required power $P_t^{\prime }$ in the specific sea condition, as follows:

$$P_t={\displaystyle \frac{R_T·v_t}{{\eta }_d{\eta }_s}},P_t^{\prime }={\displaystyle \frac{R_c·v_t^{\prime }}{{\eta }_d{\eta }_s}}{P}_t=P_t^{\prime }$$

(31)

where ${\eta }_d$ and ${\eta }_s$ are the propulsion and transmission efficiencies.

Then, the cruising speed with speed loss can be calculated using (32):

$$v_t^{\prime }={\displaystyle \frac{R_T}{R_c}}{v}_t={\displaystyle \frac{R_T}{R_c}}\sqrt[c_2]{{\displaystyle \frac{P_t}{c_1}}}$$

(32)

3. Uncertainty Model

3.1. Uncertainty Variables

In this paper, the uncertain PV power generation $P_t^{Pv}$ and propulsion load power $P_t^{lp}$ associated with speed are expressed as the sum of the expected value and a random term for forecast error ${\tilde{P}}_t$, given as

$$P_t^{Pv}={\overline{P}}_t^{Pv}+{\tilde{P}}_t^{Pv}$$

(33)

$$P_t^{lp}={\overline{P}}_t^{lp}+{\tilde{P}}_t^{lp}$$

(34)

where ${\overline{P}}_t^{Pv}$ and ${\overline{P}}_t^{lp}$ denote the expected value of PV power generation and propulsion load power at $t$ time interval, respectively, and ${\tilde{P}}_t^{Pv}$ and ${\tilde{P}}_t^{lp}$ denote the random terms of the forecast error.

The next section will describe the process of a fuzzy set for ${\tilde{P}}_t^{Pv}$ and ${\tilde{P}}_t^{lp}$. For the convenience of the study, the uncertainty parameter matrix $\theta$ is used, where $\theta =[{\tilde{P}}_t^{Pv},{\tilde{P}}_t^{lp}]$.

3.2. Fuzzy Set Construction

This section focuses on modeling the random term of PV power and the propulsive load power. According to the historical data of the random variable, as expected, the upper bounds and lower bounds of a random variable are known. The mathematical expression of the fuzzy set $\mathsf{\Omega }$ constructed in this paper based on the first-order deviation moment function is as follows:

$$\mathsf{\Omega }=\left\{P\in \Im |\begin{array}{c}{\mathrm{E}}_P\left\{\mathsf{\theta }\right\}=\mu \\ P\left\{\theta \in D\right\}=1\\ {\mathrm{E}}_P\left\{g_j(\mathsf{\theta })\right\}\le {\sigma }_j,\forall j\in [j]\end{array}\right\}$$

(35)

P is the distribution of a random variable, and $\Im$ denotes the set of all probability distributions.

The first line suggests that the expectation of random variable $\theta$ is $\mu$, which can be obtained from the historical data; the second line suggests that all the values of random variable $\theta$ are within the uncertainty set $D$, defined the same as the conventional robust optimization. In this paper, the box-type uncertainty set is adopted to limit the range of the fluctuation of the random variable; for $D=\{\underset{¯}{\theta }\le \theta \le \overline{\theta }\}$, $\overline{\theta }$ and $\underset{¯}{\theta }$ are the upper and lower bounds, respectively. $P(g)$ is the probability function; the third line represents the uncertainty of the first-order moment information of $\theta$ by adding the first-order deviation moment function $g_j(\theta )$ to the fuzzy set. $g_j(\theta )$ can be approximated to a segmented linear form, $g_j(\theta )=\max \{\theta -{\phi }_j,0\}$, which means the first-order statistical information deviates from a certain given value ${\phi }_j$, ${\sigma }_j$ is a cut-off constant, which can be estimated from historical statistics, $j$ is the index of segments describing the distribution of each random variable, and $[j]$ is a set of segments describing the distribution of each random variable.

Since it is difficult to directly evaluate the expectations for each $g_j(\theta )$ under an uncertain distribution, to facilitate the equivalent transformation of the model, the auxiliary variable $\mathrm{u}$ is introduced into the model in this paper [17]. With the extension of the original fuzzy set, this is expressed as follows:

$$\mathrm{G}=\left\{\mathrm{Q}\in \mathrm{B}|\begin{array}{c}{\mathrm{E}}_Q\left\{\theta \right\}=\mu \\ P\{(\theta ,\mathrm{u})\in H\}=1\\ {\mathrm{E}}_Q\left\{{\mathrm{u}}_j\right\}\le {\sigma }_j,\forall j\in [j]\end{array}\right\}$$

(36)

$ℍ$ is the new extended uncertainty set, as shown in (37):

$$\left.ℍ=\left\{\begin{array}{c}\overline{\theta }\le \theta \le \underset{\_}{\theta }\\ u_j\ge \theta -{\phi }_j,\forall j\in [j]\\ u_j\ge 0,\forall j\in [j]\end{array}\right.\right\}$$

(37)

4. Solution Method

4.1. Matrix Form of the Model

To clearly demonstrate the solution of the proposed distributionally robust nonlinear programming (DRONLP), the compact forms are shown in Equations (38)–(45):

$$Q(x_t,\theta )=\begin{array}{c}\max \\ \theta \end{array}\begin{array}{c}\min \\ x_t\end{array}\left[{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{n=1}^N(C_t^i+C_t^{ion}+C_t^{ioff})}}\right]$$

(38)

$$s.t.{\mathcal{M}}_{IEQ}·x_t=m_{IEQ},\forall t\in \mathrm{T}$$

(39)

$${\mathcal{M}}_{IEQ}·x_t\le m_{IEQ},\forall t\in \mathrm{T}$$

(40)

$$C\theta +Du\le d$$

(41)

$$1_a^T·x_t=c_1·{(y_t)}^{c_2},\forall t\in \mathrm{T}$$

(42)

$${\mathcal{Y}}_{min}\le \mathcal{Y}={[y_1,\cdots ,y_{\mathrm{T}}]}^T\le {\mathcal{Y}}_{max}$$

(43)

$${\mathrm{D}}_{\min }\le {\mathrm{V}}_T·\left[\begin{array}{ccc}{k}_v(u_1)& \mathrm{L}& 0\\ \mathrm{M}& \mathrm{O}& \mathrm{M}\\ 0& \mathrm{L}& k_v(u_T)\end{array}\right]·\mathrm{Y}\le {\mathrm{D}}_{\max }$$

(44)

$$(\theta ,u)\in ℍ$$

(45)

where $x_t$ = [${(o_{1,t}^{DG},\cdots ,o_{N,t}^{DG})}^T$, ${(P_{1,t}^{DG},\cdots ,P_{N,t}^{DG})}^T$, $P_t^{ESS}$, $P_t^{Pv}$, $P_t^{PL}$] is the generation variable; $y_t={\nu }_t^c$ is the speed variable; and $Q(x_t,\theta )$ is the generator fuel consumption cost and the on/off cost. Equations (39) and (40) represent linear equation and inequality constraints; Equation (41) is the matrix form of the extended fuzzy set; Equation (42) is the relationship between propulsion load and speed; Equation (43) is the sailing speed constraint; Equation (44) is the sailing distance constraint; and Equation (45) is the uncertain parameter set constraint. All of the above coefficient matrices are derived from the original model.

4.2. Bi-Level Formulation of Proposed Model

Equations (38)–(43) are two scheduling problems, i.e., power generation scheduling and voyage scheduling. The two scheduling problems are coupled through Equation (42), i.e., the relationship between propulsion load and cruising speed. If this coupling constraint is removed, the original problem can be decomposed into two levels as shown in Equations (46)–(48) [14]:

$$master:Q(x_t,\theta )$$

(46)

$$\begin{gathered} s.t.{\mathcal{M}}_{IEQ}·x_t\le m_{IEQ},\forall t\in \mathrm{T} \\ \hspace{1em}\hspace{1em}C\theta +Du\le d \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{1}_a^T·x_t=c_1·{(y_t)}^{c_2},\forall t\in \mathrm{T} \end{gathered}$$

(47)

$$slave:y_t^{\ast }=\arg \min \left(\left|\begin{array}{c}{D}_{\max }-{\mathrm{V}}_T·K_v·\mathrm{Y}\\ {\mathrm{V}}_T·K_v·\mathrm{Y}-D_{\min }\end{array}\right|:(44),(45)\right)$$

(48)

In the above formulation, $x_t$ and $y_t$ are determined in different levels, and $Q(x_t,\theta )$ becomes related to $x_t$ only, while $y_t$ is considered as a constant in the master problem, which is updated by the slave level. For the slave level, $y_t$ is obtained by minimizing the voyage deviation. With this decomposition, Equations (46) and (47) are quadratic programming, and Equation (48) is linear programming. The slave level can be solved directly by commercial solvers; the master level is described in the next section.

4.3. Linearization for the Master Level

4.3.1. Two-Stage Optimization Model

The master level model can be transformed into a two-stage optimization problem, and the variables in the master level are decomposed into first-stage variables and second-stage variables. The first-stage variables are decision variables that do not change with uncertainties, denoted by $x_1$, and before the power of PV and propulsion loads are determined, the on/off plan of generating units will be formulated first; the second-stage variables will be affected by uncertainties, denoted by $x_1$, and the power of each generating unit will be optimized according to the actual power of PV and propulsion loads under the constraints.

The first-stage optimization problem can be formulated as shown in (49) and (50):

$$\min c^T{x}_1+Q(x_1,\theta )$$

(49)

$$s.t.A·x_1\le b$$

(50)

The second-stage objective function $Q(x_1,\theta )$ can be represented by the second-stage variable $x_2$, as follows:

$$Q\left(x_1,\theta \right)={\min \mathrm{q}}^T{x}_2$$

(51)

$$s.t.N{x}_1+M{x}_2\le h(\theta )$$

(52)

The function $h(\theta )$ contains the random variable $\theta$, which will be affected by uncertainty, as follows:

$$h(\theta )=h^0+\sum _{s=1}^s{h}_s^{\theta }{\theta }_s$$

(53)

where $h^0$ denotes the constant coefficient column vector and $h_s^{\theta }$ denotes the column vector associated with ${\theta }_s$.

4.3.2. Two-Stage Distributional Robust Optimization Model

Based on the fuzzy sets, the previous two-stage optimization model can be written as a two-stage distributionally robust optimization model, as follows:

$$\begin{gathered} \min c^T{x}_1+\underset{\mathbb{Q}\in \mathbb{G}}{\max }{\mathbb{E}}_{\mathbb{Q}}\left\{Q(x_1,\theta )\right\} \\ s.t.A·x_1\le b \end{gathered}$$

(54)

The objective function minimizes the expectation of function $Q(x_1,\theta )$ under a distribution $P$. The objective function minimizes the expectation of the function $Q$ under the probability distribution $\mathbb{Q}$, which is the worst distribution on the fuzzy set $\mathbb{G}$. The function $Q(x_1,\theta )$ represents the economic dispatch cost associated with $x_1$ under a random variable $\theta$, which can be solved by the second-stage problem.

The second-stage objective function $\underset{\mathbb{Q}\in \mathbb{G}}{\max }{\mathbb{E}}_{\mathbb{Q}}\left\{Q(x_1,\theta )\right\}$ can be written in a general matrix form, as follows:

$$\begin{gathered} \underset{\mathbb{Q}\in \mathbb{G}}{\max }{\mathbb{E}}_{\mathbb{Q}}\left\{Q(x_1,\theta )\right\}=\min \underset{\mathbb{Q}\in \mathbb{G}}{\max }{\mathbb{E}}_{\mathbb{Q}}\left\{q^T{x}_2\right\} \\ s.t.N{x}_1+M{x}_2\le h(\theta ) \end{gathered}$$

(55)

However, solving the second-stage decision variable $x_2$ based on all realizations $\theta$ within fuzzy set G is an NP-hard problem. A linear decision rule is an efficient way to deal with this y intractable problem [22,23], which restricts decision variable $x_2$ affine dependence on the uncertain parameter $\theta$.

In this paper, a linear decision rule model is adopted which approximates $x_2(\theta ,u)$ by a linear affine function of random variables $\theta$ and auxiliary variables u.

The linear affine function $x_2(\theta ,u)$ can be expressed as (56):

$$x_2(\theta ,u)={x_2}^0+\sum _{s=1}^S{x}_{2s}^{\theta }{\theta }_s+\sum _{j=1}^J\sum _{l=1}^L{x_2}_{jl}^u{u}_{jl}$$

(56)

By converting the second-stage decision variable $x_2$ to a linear affine function $x_2(\theta ,u)$, the original second-stage optimization problem can also be re-represented as follows:

$$\min \underset{\mathbb{Q}\in \mathbb{G}}{\max }{\mathbb{E}}_{\mathbb{Q}}\left\{q^T{x}_2(\theta ,u)\right\}$$

(57)

$$s.t.N{x}_1+M{x}_2(\theta ,u)\le h(\theta )$$

(58)

Instead of solving the recourse problem under all outcomes of uncertain parameters, the reformulated problem attempts to search for optimal coefficients of the decision rule $x_2(\theta ,u)$ such that the worst-case expectation is minimized under all second-stage constraints. A detailed proof is given in Appendix A in the literature [17].

In the second stage, the inner maximization problem $\underset{\mathbb{Q}\in \mathbb{G}}{\max }{\mathbb{E}}_{\mathbb{Q}}\left\{q^T{x}_2(\theta ,u)\right\}$ can be expressed as a semi-infinite optimization problem, the problem is further transformed into a finite-dimensional optimization problem by Lagrange duality, and the original two-layer min-max optimization problem is equivalently reduced to a minimization problem, as shown in Equations (59)–(62):

$$\min \eta +{\gamma }^T\lambda$$

(59)

$$\mathrm{s}.\mathrm{t}.\lambda \ge 0$$

(60)

$$\eta +{\theta }^T\rho +u^T\lambda \ge q^T{x}_2(\theta ,u)$$

(61)

$$N{x}_1+M{x}_2(\theta ,u)\le h(\theta )$$

(62)

After the equivalent transformation of the optimization model, constraints (61) and (62) are still infinite-dimensional linear constraints, they can be transformed based on dual theory, and the model of the second-stage problem is as shown in Equations (63)–(72):

$$\min \eta +{\gamma }^T\lambda$$

(63)

$$s.t.\lambda \ge 0$$

(64)

$${\pi }_0\ge 0$$

(65)

$${\pi }_m\ge 0$$

(66)

$$\eta -q^T{x_2}^0+{\pi }_0^{T}d\ge 0$$

(67)

$${\pi }_0^T{C}_s={\displaystyle \sum _{n\in N_s^v}{q}_n{x}_{2ns}^{\theta }-{\rho }_s=q^T{x}_{2s}^{\theta }}-{\rho }_s$$

(68)

$${\pi }_0^T{D}_i={\displaystyle \sum _{n\in N_i^u}{q}_n{x}_{2ni}^u-{\lambda }_i=q^T{x}_{2i}^u}-{\lambda }_i$$

(69)

$$N_k^T{x}_1+M_k^T{x_2}^0-h_k^0+d^T{\pi }_m\le 0$$

(70)

$${\pi }_m^T{C}_s={\displaystyle \sum _{n\in N_s^v}{W}_{mn}{x}_{2ns}^{\theta }-h_{ms}^{\theta }=W_m^T}{x}_{2s}^{\theta }-h_{ms}^{\theta }$$

(71)

$${\pi }_m^T{D}_i={\displaystyle \sum _{n\in N_i^u}{W}_{mn}{x}_{2ni}^u=W_m}{x}_{2i}^u$$

(72)

After the above series of transformations, the original two-stage distributional robust optimization model is transformed into a mixed-integer linear programming model, objective function is as (73), The constraints include (50) and (64)–(72).

$$\min c^T{x}_1+\eta +{\gamma }^T\lambda$$

(73)

As mentioned above, the original problem is decomposed into master and slave levels, where the master level is transformed into a mixed-integer linear programming problem by linear decision rule. The overall process of the proposed solution method is shown in Figure 3.

5. Case Study

5.1. Simulation Parameters

In this study, four DGs, PV, and ESSs were considered for simulations, and the cases had two intermediate ports and one terminal port, as shown in Figure 4. The total voyage is 24 h; ship cruising occurred during t = 1–7, 9–15, and 17–23, and staying in a port occurred during t = 8, 16, and 24, and furthermore when t = 1, 7, 9, 15, 17, and 23,to avoid outage events, a DG remained closed until it needed to be turned on.

The service load is shown in

Table 1. Electricity service load power.

Time interval	1	2	3	4	5	6	7	8	9	10	11	12
Service load	8.96	8.92	9.43	10.97	8.45	8.52	7.96	8.54	8.01	7.97	8.95	10.82
Time interval	13	14	15	16	17	18	19	20	21	22	23	24
Service load	8.91	7.93	8.99	8.35	8.9	7.46	9.67	10.34	8.51	7.62	7.38	9

and the basic parameters of DGs and $ESS$s are shown in

Table 2. Diesel generator specifications.

DG	${\mathit{P}}_{\mathit{D}\mathit{G}\mathit{n}}\mathbf{(}\mathbf{MW}\mathbf{)}$	${\mathit{P}}_{\mathit{D}\mathit{G}\mathbf{min}}\mathbf{(}\mathbf{MW}\mathbf{)}$	${\mathit{r}}_{\mathit{u}},{\mathit{r}}_{\mathit{d}}\mathbf{(}\mathit{\%}/\mathbf{\Delta }\mathit{t}\mathbf{)}$	${\mathit{a}}_{0,1,2}\mathbf{(}\mathit{m}.\mathit{u}./\mathit{p}.\mathit{u}.\mathbf{)}$	${\mathit{g}}_{0,1,2}\mathbf{(}\mathit{u}\mathit{C}{\mathit{O}}_2/\mathit{p}.\mathit{u}.\mathbf{)}$
DG1	15	4	+/−50%	13.5, 10, 450	386, −2000, 8383
DG2	15	4	+/−50%	13, 12, 430	386, −2000, 8383
DG3	15	4	+/−50%	13.5, 12, 460	363, −650, 950
DG4	15	4	+/−50%	5.6, 58, 390	125, 450, 850
ESS	${\eta }_{ch}=95\%$			${\eta }_{dis}=97\%$

[14,24]. The parameters of the voyage scheduling are shown in

Table 3. Parameters of voyage scheduling.

$l=175\mathrm{m}$	$B_{\mathrm{int}}=25.4\mathrm{m}$	$S_t=274\mathrm{m}$	${\eta }_d=98\%$	${\eta }_s=0.7$	$P_t^{lp}=0.003{v_t}^3$
$h_t\in [0,4.5] \mathrm{m}$	${\tau }_t\in [0,1.5] \mathrm{m}$	$v_w\in [0,15.5] \mathrm{m}/\mathrm{s}$	${\theta }_t\in [30,180] \mathrm{deg}$	$C_{air}=0.025$	${\rho }_{air}=1.29\times {10}^{-3}{\mathrm{g}/\mathrm{cm}}^3$
${\rho }_{water}=1.04{\mathrm{g}/\mathrm{cm}}^3$	$C_{D,wat}=0.041$ [26]	$g=9.8{\mathrm{m}/\mathrm{s}}^2$

. The predicted mean value and upper and lower bounds of PV power and propulsion power are shown in Figure 5. The capacity of the ESS is 30 MWh, and the rated power is 15 MW. The rated cruising speed is 19 knots. Environmental data along the cruise route were downloaded from the ERA5 global reanalysis data set [25].

5.2. Analysis of Results

5.2.1. On-Time Rates Under Different Methods

In this study, the speeds of the three methods are analyzed as shown in Figure 6, where RO is robust optimization, DRO is distributionally robust optimization, NR is without considering wind and wave resistance, and Nom denotes the nominal speed. In Figure 6, the speeds of both RO and DRO are significantly higher than those in the case of NR, so both can cover the speed loss caused by wind and waves and ensure that the ship arrives at the port on time. On the other hand, the speed in NR is obviously lower than the rated speed, which makes the ship fail to arrive at the harbor on time.

When t = 1~7 and 17~23, due to DRO not being as conservative as the RO for the uncertainty, the cruise speed is relatively low compared to the RO, and the power of the propulsion load is also lower at these intervals, so the generator output power is reduced accordingly, and so DRO improves the ship’s energy-saving and emission reduction capabilities.

5.2.2. Analysis of Economic Costs

Figure 7 shows the fuel consumption for each time step under four different methods, Due to PV being used as supplementary energy during t = 10–15 periods, fuel consumption during t = 10–15 is lower than at other times, so the integration of the PV system can reduce ship operating costs. In NR, the fuel consumption is lowest because the impact of wind and waves is not considered, but it also causes the ship to fail to arrive at the port on time. The RO and DRO consider the uncertainty; therefore, fuel consumption is higher than the Nom Speed scenario. It can also be seen from Figure 8 that the fuel consumption of DRO is always lower than that of RO.

Table 4. Total fuel consumption for different methods.

	Total Fuel Consumption (m.u.)
DRO	70,222
RO	79,950
NR	54,357
Nom Speed	64,890

shows the total fuel consumption for four different methods, and it can be clearly seen that DRO generally outperforms RO in yielding lower expected costs.

Figure 8 shows that two units are turned on during most of the navigation time for the DRO method since the two units are enough to follow wave, wind, and PV fluctuations at a minimum cost. But for the RO method, three units are turned on for most of the navigation time, which requires more economic costs. Therefore, compared to RO methods, DRO methods greatly reduce the conservatism of decision-making and lower scheduling costs.

5.2.3. Analysis of Sensitivity

When the uncertain variable fluctuates in the range of ±15~±35% of the expected value, it can be seen from Figure 9 that the total cost of DRO is generally lower than that of RO, and with an increase in interval fluctuation, the advantage of DRO is more obvious. This is because the decision made by the traditional RO depends on the worst case, while the DRO can simulate the distribution of uncertain variables according to the probability information contained in the fuzzy set, which reduces the conservative decision.

6. Conclusions

To address the effects of uncertain wave and wind conditions in the navigation route, this paper proposes a DRO formulation for solving the joint scheduling model of power and voyage, aiming at operation cost minimization and GHG emission limitation. The relationship between propulsion load and voyage as influenced by navigation speed and waves is first established, and then a DRO model is established. Due to the complexity of the examined model, the problem is decoupled into a bi-level optimization model, the slave level can be solved directly by commercial solvers, the master level is further formulated as a two-stage DRO model, and linear decision rules and column and constraint generation algorithms are adopted to solve the model. With the exploitation of the proposed method, joint scheduling is achieved with excellent results. The obtained simulation results confirm the effectiveness of the method, and compared with the conventional robust optimization methods, less conservative results are obtained. Expected operation costs are minimized and emissions are reduced while ensuring the on-time arrival of hybrid ships in various uncertain scenarios. Future efforts might be spent on the consideration of incorporating risk into the model to better guide operators.