Optimal Scheduling of an Electric Propulsion Tugboat Considering Various Operating Conditions and Navigation Uncertainties

Abstract

The operating conditions of all-electric tugboats are flexible and changeable. They are more complicated than conventional vessels in terms of joint voyages and power generation scheduling. To guarantee the reliable operation of the ship, a new coordinated optimization scheme that combines economy and operational reliability is proposed. It is based on the various operating conditions of the tugboat during its voyage, taking into account the random outages of equipment and load fluctuations due to speed and wave uncertainties. Due to the difficulty of implementing a stochastic sampling method with space-time coupling constraints (e.g., the voyage is related to propulsion load), an analytical approach is needed to transform the model into a readily solvable mixed-integer linear program (MINP) which attributes risk scenarios to load fluctuations under various conditional probabilities. In addition, this paper proposes an improved piecewise linearization method based on a differential evolutionary algorithm to speed up the solution process and improve computational accuracy. Meanwhile, the energy storage loss cost due to battery degradation is added to the optimization target. The battery’s cycle life is extended by rational scheduling of charging and discharging. Simulations validate this paper’s joint scheduling optimization scheme in multiple comparison experiments. The results show that it can effectively balance the economic and reliability levels under various risk scenarios and improve the environmental energy efficiency indicators.

1. Introduction

All-electric ships (AESs) combine propulsion and service load through an energy management system [1]. The flexibility, safety, and energy utilization efficiency are improved compared with conventional mechanical propulsion ships [2]. Furthermore, for China to achieve the ambitious target of carbon peaking by 2030 and carbon neutrality by 2060, the research and promotion of AES are crucial.

In [3,4,5,6], the unit combination, generation scheduling, sizing of the energy storage system, and energy management of the ship power system have been intensively studied. As the above study, the influence of the voyage information on the optimization objective is ignored. Being a mobile microgrid, the propulsion load of the AES accounts for more than 70% of the total load demand [7]. Therefore, the optimization of ship speed by scheduling propulsion load and the generation scheduling by variable load demand need to be considered simultaneously to constitute joint generation and voyage scheduling [8,9].

To solve this problem, many improved optimization models have been proposed to jointly optimize AES generation and load dispatch in recent years [7,8,9,10,11,12,13]. The authors of [10] integrated the voyage scheduling and power system dispatching of AESs to develop the carbon price, and the method can better incentivize AESs to achieve emission control. In [7], a two-layer robust optimization model is used to minimize the fuel consumption and energy efficiency operational indicator by considering the ship speed loss due to wave and wind uncertainties for joint power generation and voyage schedules. In [11], a stochastic programming model is adopted to avoid the operational risk caused by renewable energy generation and load-side uncertainty. They introduced the conditional value-at-risk indicator in the objective function. Joint optimization of the voyage and multi-objective energy management is achieved by coordinating the AES with the hybrid energy storage system [12]. The study in [13] tries to extend the GHG emissions as an objective function to the AES joint scheduling model, which is solved using the non-dominated sorting genetic algorithm II algorithm. However, the studies mentioned above for AES voyage optimization only focus on economic and GHG emission targets. The dynamic reliability level of ship operation optimization has not been focused on.

The voyage conditions of tugboats are significantly different from those of conventional ships, where the power demand is high during towing, with a total load factor of more than 90% [14,15]. In contrast, the load level is lower during berthing or cruising. Therefore, it challenges the optimal economic dispatch of power under guaranteed reliability conditions when the operating conditions are variable. The studies show that resistance generated by wind and waves causes loss of ship speed [16,17]. The above study considered the effect of ship speed loss in AES operation optimization due to wind and wave variation. However, the propulsion power variation required to maintain the predefined voyage is not considered, which is significant for the reserve power in the ship’s power system at each moment. Additional power reserves are required to maintain the preset speed, considering the added resistance caused by waves. It is necessary to avoid power shortages under heavy load conditions. Therefore, the chance-constrained model [18,19,20] is used in this study. Stochastic outages of equipment components such as generators and batteries under load uncertainty are considered to avoid overly optimistic scheduling decisions.

Solving the voyage and power generation scheduling model with wind and wave uncertainty and random component outages has several issues to be addressed. First, the typical equation of ship speed and propulsion power used in the literature [11,13,21] does not apply to tugboats. When a tug is in dragging condition, its speed will change due to the additional resistance of the towed vessel. Therefore, it is necessary to model the propulsion system under the influence of the meteorological environment according to different working conditions of the tugboat [17]. Secondly, the piecewise linearization methods mentioned in the above studies convert nonlinear functions into MILP by equally spaced segmentation. However, introducing too many integers increases the burden of the solution. In addition, probabilistic constraints in the model need to be approximately transformed while considering balancing accuracy and computational efficiency [22,23]. The study in [14,19] used battery balancing for power fluctuations in the system, but the effect of cyclic and irregular charging/discharging strategies on the battery life was not considered. Consequently, the battery loss model needs to be established so that the control strategy can be optimized [24,25,26].

The following are the main contributions of this paper. This paper proposes a joint optimal scheduling model for navigation and power generation considering tugboat operation’s reliability to solve the above-proposed problem. The main contributions are as follows.

(1) System operational reliability has not been a concern in previous AES operational optimization problems. In this study, a combined navigation and power generation scheduling model is constructed for different operating conditions of tugboats. A probabilistic constraint on operational reliability is adopted instead of the traditional deterministic constraint to avoid over-optimistic optimization schemes.

(2) The typical ship speed and propulsion power formula do not apply to tugboats. Nonlinear expressions for the tugboat’s speed in calm water and the propulsion power for the variable operating conditions of the tugboat departure, docking, dragging, cruising, and berthing are developed. Due to the increased resistance caused by wind speed and waves, the preset speed of the ship is effectively maintained by dispatching the corresponding reserve power. This way, ships’ operational reliability and arrival rate can be enhanced.

(3) For nonlinear and chance constraints in the model, the optimal adaptive piecewise linearization and the probability distribution discretization method are employed to transform them into linear deterministic constraints. Previous studies have not been concerned with the piecewise linearization method’s effect on the calculation’s accuracy. The improved model used in this study significantly reduces the number of decision variables and decreases the computational burden.

(4) The energy storage system consisting of battery units combined with generator units can improve energy utilization and flatten power fluctuations. Conventional battery charging and discharging strategies do not consider the economic loss due to degradation. In this study, the loss of energy storage lifetime is quantified and translated into an economic loss. It effectively mitigates the loss of battery life.

The rest of this article is organized as follows: Section 2 presents the mathematical model of the problem. Section 3 illustrates the optimization problem formulations. The solution method for the optimal operation of tugboats is introduced in Section 4. Then, Section 5 provides case studies and a discussion of the results. Section 6 concludes the work.

2. Mathematical Modeling

2.1. Problem Description

Tugboats are vessels that carry out towing, transportation operations, and rescue missions and play an essential role in economic production. The primary operating area of the port tugboat studied in this paper is about 12 sea miles from the shore. In this study, an electric propulsion tugboat with an energy storage system is used as the research object, considering that the ship, as a mobile microgrid, needs to regard the reasonable distribution of power and is constrained by the voyage. Therefore, a mathematical model for the joint dispatch of power generation and navigation schedule under multiple working conditions of a tugboat is established.

2.2. Tugboat Propulsion System Modeling

The system topology of an AES is shown in Figure 1. Three diesel generators provide the necessary power for the propulsion and service load of the entire ship, and a storage system consisting of two lithium batteries. The batteries can be used to achieve peak shaving, energy saving, and emission reduction depending on the load profile of the generator under various operating conditions. Since the ship propulsion load accounts for more than 90% of the total load while towing, the power fluctuations due to wind and wave resistance cannot be ignored. The following equation usually expresses the relationship between propulsion power and ship speed, and $h_1$ and $h_2$ are fixed parameters after fitting [11,13].

$$P_{pl}^t=h_1·{(V^t)}^{h_2}$$

(1)

The formula does not apply to tugboats. The sailing resistance increases significantly under towing conditions, and the power required for propulsion to overcome the resistance increases. The power increment due to wind and wave resistance must also be solved for large vessels. Therefore, the model of propulsion power of an all-electric tugboat under multiple working conditions and different meteorological conditions is established for the above problems, as in (1)–(6). In the towing condition, the required power of the tugboat and the dragged vessel need to be summed. In the following equations, $\mathcal{L}$, $\mathcal{B}$, $\mathcal{D}$, and $C_B$ denote the ship’s length, breadth, draft, and block coefficient, respectively. $H_w$ represents the wave height. $c_1-c_5$ denote the different resistance coefficients, respectively, and they are derived from [15,16,27,28]. Equations (3) and (4) are used to calculate the frictional resistance and residual resistance of a ship in calm water, respectively [29]. Equation (2) denotes the desired propulsion power of an AES in calm water conditions. Wave height $H_w$ is expressed as a quadratic function of wind speed as in Equation (5). The additional power required to maintain the preset speed is solved by (6) [30], which describes the effect of varying wind and wave conditions on the propulsion power.

$$P_{pl}^t=(R_F^t+R_R^t)·V^t/{\eta }_{pl}$$

(2)

$$R_F^t=c_1·\mathcal{L}·(\mathcal{B}+2\mathcal{D})·{(V^t)}^{c_2}$$

(3)

$$R_R^t=c_3·C_B·\mathcal{B}·\mathcal{D}·{(V^t)}^{(c_4+c_5{V}^t)}$$

(4)

$$H_w=0.02·{(V_s^t)}^2$$

(5)

$$△P_e^i=4.594·{(2H_w+0.152)}^2·C_B·{\mathcal{B}}^2/\mathcal{L}/{\eta }_{pl}$$

(6)

2.3. Wind Speed Probability Distribution Model

The above illustrates that the wind speed values are significant to solve for the power increment due to maintaining the pre-defined speed of the ship. Therefore, a probabilistic distribution model of wind speed is needed to further describe the propulsion power fluctuations during navigation. Many studies have shown that the wind speed follows the Weibull distribution [19,31]. Its probability density function and cumulative distribution function are (7) and (8), where $v_w$ denotes the wind speed, k is the shape factor, and $\eta$ is the scale factor [32,33]. Figure 2 shows the fitted probability distribution based on the annual wind speed data collected from the meteorological station in Yingkou.

$$f(v_s)={\displaystyle \frac{k}{\eta }}·{({\displaystyle \frac{v_s}{\eta }})}^{k-1}·exp[-{({\displaystyle \frac{v_s}{\eta }})}^k]$$

(7)

$$F(v_s)=1-exp[-{({\displaystyle \frac{v_s}{\eta }})}^k]$$

(8)

2.4. Battery Loss Assessment Model

The practical life of a battery is determined by a combination of factors such as temperature, the number of cycles, and depth of discharge [24]. In a hybrid energy power system, the scheduling strategy of energy storage directly affects the battery life loss [26]. Therefore, the battery losses throughout its lifetime need to be quantified and converted into economic costs to facilitate the joint scheduling optimization of ship energy [25]. The lifetime loss of battery energy storage is modeled by (9)–(11).

$$N_{cycle}^{n,t}=s_0+s_1·e^{s_2(1-SO{C}^{n,t})}+s_3·e^{s_4(1-SO{C}^{n,t})}$$

(9)

$$L_n^t=\frac{1}{N_{cycle}^{n,t}}$$

(10)

$$Los{s}^{n,w,t}=S_B^{n,w,t}·|L_n^t-L_n^{t-1}|$$

(11)

3. Problem Formulation

A mathematical model of the joint voyage and generation scheduling that includes chance constraints is developed to balance the economy and reliability of the voyage. The model simulates the navigation scenario and contains output constraints on the generators and batteries, voyage constraints, and reliability probability constraints. As in Figure 3, the tugboat receives the operation task and departs from port A, sailing in cruising mode to location B of the vessel to be towed. Finally, the tugboat drags the target vessel from point B to port C in dragging mode. When the ship is berthing in port C, the shore power charges the lithium battery to satisfy the power demand during berthing.

3.1. Objective Function

In this study, the expectation of the economic cost of one cruise is taken as the objective function of the optimization problem. As shown in (12), the total cost contains the generator’s fuel consumption cost $C_f^{w,t}$, start-up cost $C_s^{w,t}$, power reserve cost $R_G^{w,t}$, battery operation cost $C_B^{w,t}$, and power reserve cost $R_B^{w,t}$. In order to reduce the lifetime loss of battery energy storage due to random, irregular charging and discharging, the loss cost $L_B^{w,t}$ is added to the optimization objective. The specific modeling is shown in (13)–(20).

$$min\sum _{w=1}^W\sum _{t=1}^T{\rho }^w·(C_G^{w,t}+C_B^{w,t}+R_G^{w,t}+R_B^{w,t}+L_B^{w,t})$$

(12)

$$C_G^{w,t}=C_f^{w,t}+C_s^{w,t}$$

(13)

$$C_f^{w,t}=\sum _{k=1}^K{S}_G^{k,w,t}·u^{k,t}·[a_0^k+a_1^k·r_G^{k,t}+a_2^k·{(r_G^{k,t})}^2]·\Delta t$$

(14)

$$P_G^{k,t}=P_{G,max}^k·r_G^{k,t}$$

(15)

$$C_s^{w,t}=\sum _{k=1}^{K}max(0,u^{k,t}-u^{k,t-1})·F_{up}^k$$

(16)

$$C_B^{w,t}=\sum _{n=1}^N{S}_B^{n,w,t}·F_B^n·(P_c^{n,t}+P_{dc}^{n,t})·\Delta t$$

(17)

$$R_G^{w,t}=\sum _{k=1}^K{Q}_G^{k,w,t}·F_G^k·\Delta t$$

(18)

$$R_B^{w,t}=\sum _{n=1}^N{Q}_B^{n,w,t}·F_{Br}^n·\Delta t$$

(19)

$$L_B^{w,t}=\sum _{n=1}^{N}Los{s}^{n,w,t}·F_{Bl}^n$$

(20)

3.2. Problem Constrains

As a mobile microgrid, an AES is not only subject to the constraints of the system power balance and the physical constraints of equipment operation but also to the limitations of ship speed and voyage. Equation (21) represents the real-time power balance constraint in the system. Equations (22)–(25) shows the operating constraints of the generating units. Equations (23) and (24) are the minimum on and off time constraints, which represent the restrictions on the continuous on and off time of the generators.

Equation (25) indicates the constraint on the generator reserve power in case of an emergency w. Equations (26)–(30) represent the operational constraints of energy storage. Equation (31) describes the ship’s speed constraint. Equation (32) represents the distance constraint from port A to position B. Equation (33) represents the entire distance limitation of the voyage from port A to position B and then to port C. Equation (34) expresses the operational reliability constraint for the entire voyage, which is the form of the probability constraint, and the deterministic transformation of this formula is described in Section 4.3. To translate into a MILP problem, contingency events where generator or battery outages may occur are equated with loading surges.

Power Balance Constrains:

$$\sum _{n=1}^N(P_{dc}^{n,t}-P_c^{n,t})+\sum _{k=1}^K{P}_G^{k,t}·u^{k,t}=P_{pl}^t+P_{ser}^t$$

(21)

Generator Constraints:

$$P_{G,min}·u^{k,t}\le P_G^{k,t}\le P_{G,max}^k·u^{k,t}$$

(22)

$$(u^{k,t}-u^{k,t-1})(T{N}^{k,t-1}-T{N}_{min}^k)\le 0$$

(23)

$$(u^{k,t-1}-u^{k,t})(T{F}^{k,t-1}-T{F}_{min}^k)\le 0$$

(24)

$$Q_G^{k,w,t}\le S_G^{k,w,t}·u^{k,t}·(P_{G,max}^k-P_G^{k,t})$$

(25)

Battery Energy Storage Constraints:

$$P_{c,min}^n\le P_c^{n,t}\le P_{c,max}^n$$

(26)

$$P_{dc,min}^n\le P_{dc}^{n,t}\le P_{dc,max}^n$$

(27)

$$SO{C}_{min}^n\le SO{C}^{n,t}=\frac{E_B^{n,t}}{E_B^n}\le SO{C}_{max}^n$$

(28)

$$SO{C}^{n,t}=SO{C}^{n,t-1}+\frac{{\eta }_c^n·P_c^{n,t}·\Delta t}{E_B^n}-\frac{P_{dc}^{n,t}·\Delta t}{{\eta }_{dc}^n·E_B^n}$$

(29)

$$\begin{array}{c}\hfill Q_B^{n,w,t}\le S_B^{n,w,t}·min\left\{{\displaystyle \frac{(SO{C}^{n,t}-SO{C}_{min}^n)·E_B^n·{\eta }_{dc}^n}{\Delta t}},P_{dc,max}^n-P_{dc}^{n,t}\right\}\end{array}$$

(30)

AES Voyage Constraints:

$$\left\{\begin{array}{c}{V}_{doc}(1-{\mu }_1)\le V^t\le V_{doc}(1+{\mu }_1)t\in T_{doc}\hfill \\ V_{dep}(1-{\mu }_2)\le V^t\le V_{dep}(1+{\mu }_2)t\in T_{dep}\hfill \\ V_{cru}(1-\lambda )\le V^t\le V_{cru}(1+\lambda )t\in T_{cru}\hfill \\ V_{dra}(1-\epsilon )\le V^t\le V_{dra}(1+\epsilon )t\in T_{dra}\hfill \end{array}\right.$$

(31)

$$D_{AB}(1-{\delta }_{AB})\le \sum _{t=1}^{T_{AB}}{V}^t·\Delta t\le D_{AB}(1+{\delta }_{AB})$$

(32)

$$D_{AC}\le \sum _{t=1}^{T_{AC}}{V}^t·\Delta t\le D_{AC}(1+{\delta }_{AC})$$

(33)

Operational Reliability Constraints

$$\begin{array}{c}\hfill \sum _{w=1}^W{\rho }^w·Pr\left\{\sum _{k=1}^K{Q}_G^{k,w,t}+\sum _{n=1}^N{Q}_B^{n,w,t}\ge (1-S_G^{k,w,t})·P_G^{k,t}·u^{k,t}+(1-S_B^{n,w,t})·P_{dc}^{n,t}+△P_e^t\right\}\ge 1-Lolp\end{array}$$

(34)

4. Scheduling Optimization Model Conversion

The optimization problem described in Section 3 includes multiple nonlinear constraints such as (2), (9), (23), (24), and (34), which would result in an excessive computational burden. This section uses the corresponding linearization methods to reduce computational complexity and improve model accuracy. The converted MILP model is solved based on the CPLEX solver.

4.1. Improved Piecewise Linearization Model

In the above model, (2) and (9) are nonlinear functions, which are usually transformed into mixed integer forms using piecewise linearization. Equations (35)–(41) are specific linearization processes. The variable $x^t$ in the t period is split into the form of a summation of the product of multiple continuous variables $w_m^t$ and segmented points $x_m$. $\widehat{g}(·)$ is a linear function after segmentation, where $l_m^t$ is a binary variable.

$$\begin{array}{c}{x}^t=\sum _{m=1}^{M+1}{w}_m^t·x_m\end{array}$$

(35)

$$\begin{array}{c}{x}_m=(m-1)·\Delta x+x_{min}\end{array}$$

(36)

$$\begin{array}{c}\Delta x=(x_{max}-x_{min})/M\end{array}$$

(37)

$$\begin{array}{c}\widehat{g}(x^t)=\sum _{m=1}^{M+1}{w}_m^t·g(x_m)\end{array}$$

(38)

$$\begin{array}{c}\sum _{m=1}^{M+1}{w}_m^t=1w_m^t\ge 0\end{array}$$

(39)

$$\begin{array}{c}\sum _{m=1}^M{l}_m^t=1l_m^t=0or1\end{array}$$

(40)

$$\left\{\begin{array}{c}{w}_m^t\le l_m^{t}m=1\hfill \\ w_{m+1}^t\le l_m^{t}m=M\hfill \\ w_m^t\le l_{m-1}^t+l_m^{t}1<m<M\hfill \end{array}\right.$$

(41)

$$min\sum _{j\in Num}{(g(x^j)-\widehat{g}(x^j))}^2$$

(42)

In (36) and (37), the use of uniform segmentation affects the fitting accuracy of the curve. Therefore, a modified piecewise linearization model is proposed in this study to speed up the fitting speed and improve the estimation accuracy. A differential evolution method is used to optimize the location of the breakpoint that minimizes the sum of squared errors. The corresponding objective function is shown in (42). $Num$ is the number of samples. The fitted error values at x for the equally spaced segmentation function and the optimized segmentation function are shown in Figure 4. The red line represents the fitted function for each segment of equal length, denoted as the equally spaced segmentation function. The blue line indicates the fitted function after the optimized segmentation point by an improved piecewise linearization model, denoted as the optimized segmentation function. In addition, the yellow line represents a nonlinear function curve, which is denoted as a true nonlinear function. $Error1$, produced by the equally spaced linearization method, is significantly larger than $Error2$, produced by the proposed improved segmented linearization method.

4.2. Minimum on and off Time Linearization

The minimum on and off time constraints are expressed in (23) and (24). Taking the minimum operating time as an example, its mixed integer linearized expression is shown in (43)–(45). These equations are based only on the binary variables of the generator operating state.

$$\begin{array}{c}\sum _{t=1}^{G_k}(1-u^{k,t})=0\\ G_k=min(T,(T{N}_{min}^k-T{N}^{k,0})·u^{k,0})\end{array}$$

(43)

$$\begin{array}{c}\sum _{s=t}^{t+T{N}_{min}^k-1}{u}^{k,s}\ge T{N}_{min}^k·(u^{k,t}-u^{k,t-1})\\ \forall k\in K,\forall t=G_k+1\cdots T-T{N}_{min}^k+1\end{array}$$

(44)

$$\begin{array}{c}\sum _{s=t}^T{u}^{k,s}-(u^{k,t}-u^{k,t-1})\ge 0\\ \forall k\in K,\forall t=T-T{N}_{min}^k+2\cdots T\end{array}$$

(45)

where $G_k$ denotes the initial time period in which generator k must be online. $T{N}^{k,0}$ represents the number of periods when the initial moment is already online. The initial state of k is expressed by $u^{k,0}$.

4.3. Deterministic Transformation of Chance Constraints

Equation (34) is a chance constraint and needs to be transformed into a deterministic expression. The formula contains two kinds of random variables, and one is the discrete variable $S_G^{k,w,t}$ that expresses the occurrence of contingencies in the system with probability ${\rho }^w$. The occurrence of such random events is equated to sudden load increases. In this study, the critical event is set to be a single outage of the generator and battery. The other is a continuous variable with a sudden increase in power due to wind and waves. The cumulative distribution function (8) of the wind speed is discretized by (46)–(48).

$$\begin{array}{c}\sum _{z=0}^{Z}p(z)=1,p(z)\ge 0,z=0,1,2,...,Z\end{array}$$

(46)

$$\begin{array}{c}Z=V_{s,max}/d\end{array}$$

(47)

$$p(z)=\left\{\begin{array}{cc}F(d/2)-F(0)\hfill & z=0\hfill \\ F(z·d+d/2)-F(z·d-d/2)\hfill & z\ne Z\hfill \\ F(z·d)-F(z·d-d/2)\hfill & z=Z\hfill \end{array}\right.$$

(48)

where $p(z)$ is the discretized probability sequence and $V_{s,max}$ represents the maximum value in the historical wind speed statistics. The discrete step length of the sequence is denoted by d. Therefore, $△P_e^i$ can be transformed into the corresponding discrete values $△P_e^i(z·d)$.

$\mathcal{F}(t,s,w)$ of (49) represents the difference between the reserve power and the load increment at the moment t when contingency w occurs. The 0–1 variable ${\mathcal{Z}}_z^w$ is introduced and is recorded as 1 when $\mathcal{F}(t,s,w)\ge 0$ and 0 otherwise. Therefore, the operational reliability constraint of the system can be translated into (51), where $Lolp$ denotes the probability of load loss.

$$\begin{array}{c}\hfill \mathcal{F}(t,z,w)=\sum _{k=1}^K{Q}_G^{k,w,t}+\sum _{n=1}^N{Q}_B^{n,w,t}-(1-S_G^{k,w,t})·P_G^{k,t}·u^{k,t}-(1-S_B^{n,w,t})·P_{dc}^{n,t}-△P_e^t(z·d)\end{array}$$

(49)

$${\mathcal{Z}}_z^w=\left\{\begin{array}{cc}1\hfill & \mathcal{F}(t,z,w)\ge 0\hfill \\ 0\hfill & \mathcal{F}(t,z,w)<0\hfill \end{array}\right.$$

(50)

$$\sum _{w=1}^W\sum _{z=0}^Z{\rho }^w·{\mathcal{Z}}_z^w\ge 1-Lolp$$

(51)

The current model does not conform to the set rules of MILP and still needs further transformation. Equation (52) introduces a sizeable positive number $\mathcal{M}$ to screen for power reserve shortages in the system. In summary, the chance constraint (34) can be transformed into a deterministic constraint by (51) and (52).

$$\left\{\begin{array}{c}{\mathcal{Z}}_z^w\ge \mathcal{F}(t,z,w)/\mathcal{M}\hfill \\ {\mathcal{Z}}_z^w\le 1+\mathcal{F}(t,z,w)/\mathcal{M}\hfill \end{array}\right.$$

(52)

5. Case Study

The algorithms in this paper are programmed by Python 3.7. CPLEX solves the proposed mixed-integer quadratic programming model.

5.1. Test System and Parameters Introduction

This paper investigates the voyage and power generation scheduling of an all-electric tugboat equipped with battery energy storage. The target vessel supplies three generators and two battery units. The parameters of the generators and energy storage are set based on references [13,19], as shown in

Table 1. Economic and operating parameters of generators.

	${\mathit{P}}_{\mathbf{max}}$	${\mathit{P}}_{\mathbf{min}}$	${\mathit{a}}_0$	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{F}}_{\mathbf{up}}$	${\mathit{F}}_{\mathit{G}}$
	(MW)	(MW)	(USD/p.u.)	(USD/p.u.)	(USD/p.u.)	(USD)	(USD/MWh)
G1	1.4	0.3	111	533	14	25	140
G2	1.4	0.3	120	623	10	25	161
G3	1.4	0.4	116	530	14	25	140

. In addition, the generator reserve price is usually set to 25% of the most expensive block of the marginal generator function [34]. The parameters of the two battery storage with 980 kWh capacity are the same: ${\eta }_c$ = 95% and ${\eta }_{dc}$ = 97%. The maximum power of charging and discharging is 0.5 kW. The SOC at the initial moment is set to 100%. $F_B$ = 0.1 USD/kWh, $F_{B_r}$ = 0.035 USD/kWh [35]. The optimization parameters for this algorithm are set as: population size PN = 50, the maximum number of iterations MaxI = 1000, differential weight F = 0.5, crossover rate CR = 0.5, and absolute tolerance for convergence Atol = ${10}^{-4}$. The forced outage rate parameters for generators and batteries are referred to [36,37,38], and only the risk scenario of a single failure is considered in this paper for analysis.

The load side consists mainly of the propulsion load and service load. Based on the voyage of the towing vessel and some guiding documents, the Discharge Standard for Water Pollutants from Ships (GB 3552-2018) [39] and the Guidelines for Towage at Sea (GD 02-2012) [29], the voyage planning of the towing vessel is predetermined. For instance, the ship’s speed in and out of the port is more than 4 kn, and the tugboat’s speed in calm water is generally not less than 6 kn. Based on the above requirements, the details of an assumption voyage are shown in

Table 2. AES voyage model parameters.

AES route parameters
$△t$ = 15 min	T = 20	$T_{doc}$ = {1, 2}	$T_{dep}$ = {14, 15}
$T_{cru}$ = {3, 4,..., 7}	$T_{dra}$ = {8, 9,..., 13}	$T_{ber}$ = {16, 17,..., 20}	$V_{doc}$ = 6 kn
$V_{dep}$ = 6 kn	$V_{cru}$ = 10 kn	$V_{dra}$ = 8.3 kn	$V_{ber}$ = 0 kn
${\mu }_1$ = 0.1	${\mu }_2$ = 0.1	$\lambda$ = 0.2	$\epsilon$ = 0.2
${\delta }_{AB}$ = 0.01	${\delta }_{AC}$ = 0.005	$D_{AB}$ = 15.5 nm	$D_{AC}$ = 31 nm
$T_{AB}\in T_{doc}\cup T_{cru}$	$T_{AC}\in T_{dra}\cup T_{doc}$	$c_1=1.67\times {10}^{-3}$	$c_2$ = 1.83
$c_3$ = 0.147	$c_4$ = 1.74	$c_5$ = 0.15	${\eta }_{pl}$ = 0.97
AES main characteristics
${\mathcal{L}}_{tug}$ = 38.5 m	${\mathcal{B}}_{tug}$ = 10.6 m	${\mathcal{D}}_{tug}$ = 3.7 m	$C_{B,tug}$ = 0.97
${\mathcal{L}}_{ves}$ = 91.5 m	${\mathcal{B}}_{ves}$ = 24.5 m	${\mathcal{D}}_{ves}$ = 2.5 m	$C_{B,ves}$ = 0.95

. In the following, this study verifies the model’s validity through the four perspectives in Section 5.2, Section 5.3, Section 5.4 and Section 5.5 and analyzes the effect of different parameters on the degree of calculated results.

5.2. Voyage and Generation Scheduling Results Analysis

To test the validity of the above model, the navigation of the tugboat offshore of Yingkou is simulated. The economically optimal dispatch results under the given voyage constraint are shown in Figure 5. The optimized voyage speed is higher than the nominal speed in the docking and departure and slightly lower than in the towing condition with a heavy load. By constraining the generator on/off status, only one generator operates at the time of departure. The second generator is put into service as the tugboat enters the cruising stage. At the same time, the battery units are properly discharged to improve economic efficiency and reduce pollution emissions. When t ≥ 8, the tugboat performs towing duties and the load level rises significantly, requiring three generators to be put in simultaneously to meet the load demand. The ship’s docking speed decreases in t = 14–15, and the third generator is withdrawn from operation to maintain the economy. After that, the ship is at berth and the batteries are charged by shore power.

Three experimental cases are presented to compare the optimization results for different optimization objectives. As shown in

Table 3. Comparison of results for different optimization objectives.

	Operation Costs	Life Loss	$\mathit{Lolp}$
Case 1	5013	-	0.18
Case 2	4463	3.702 × ${10}^5$	0.02
Case 3	4473	3.427 × ${10}^5$	0.02

, Case 1 represents the scheduling result without energy storage. Case 2 represents the result, including energy storage but not considering its loss degradation. Case 3 represents the result with energy storage and adds the loss factor of energy storage degradation to the objective function. The results show that the $Lolp$ of the system is only 0.02 after adding the energy storage device, while the minimum load loss probability of the Case 1 system is 0.18. It follows that integrating energy storage can reduce operating costs and greenhouse gas emissions. Its flexible power supply can mitigate load fluctuations caused by wave resistance during the voyage and improve operational reliability. In addition, the operation results show that Case 3 has 7.43% lower life loss than Case 2.

5.3. Effect of Wind Speed Distribution

The effect of different wind speed distributions on the optimization results is discussed in this experiment. In the reliability constraint, the power increment $△P_e^t(V_s^t)$ is affected by the wind speed. Therefore, to make the results more consistent with practical applications, wind speed data from meteorological stations in three coastal cities, Weihai, Dalian, and Yingkou, fit cumulative distribution functions separately. The wind data for weather stations along the coast used herein are obtained from the SolarGIS dataset obtained by MTSAT and HIMAWARI satellites. The fitting results are shown in Figure 6 by setting $Lolp$ = 0.1. The results after discretizing the continuous distribution by (46)–(48) are shown in

Table 4. Probability sequences of wind speed at various locations.

	Weihai		Yingkou		Qinhuangdao
s	Speed	Probability	Speed	Probability	Speed	Probability
1	0.0	0.000	0.0	0.000	0.0	0.000
2	1.6	0.049	1.1	0.045	1.0	0.029
3	3.2	0.154	2.2	0.126	2.0	0.114
4	4.8	0.218	3.3	0.176	3.0	0.190
5	6.4	0.218	4.4	0.187	4.0	0.219
6	8.0	0.169	5.5	0.165	5.0	0.193
7	9.6	0.106	6.5	0.125	6.0	0.134
8	11.2	0.054	7.6	0.083	7.0	0.074
9	12.8	0.022	8.7	0.049	8.0	0.032
10	14.4	0.008	9.8	0.025	9.0	0.011
11	16.0	0.003	10.9	0.019	10.0	0.004

.

It can be seen that the probability of wind speed in the high-speed zone is higher in Weihai compared to Yingkou. Secondly, the maximum wind speed in Qinhuangdao is only 10.0 m/s, which is significantly lower than the other two locations. Equation (6) shows that the power increment required to maintain the sailing speed is exponentially related to the wind speed, i.e., the power demand is more significant at more substantial wind speeds. Therefore, the reserve power needs to be increased to respond quickly to fluctuations in propulsion load to ensure the reliability of navigation, and the operating costs will increase as a result. The results in

Table 5. Comparison of operating costs by locations.

	Weihai	Dalian	Yingkou
Total operation cost (USD)	4136	3782	3511

also confirm that the operating costs in Weihai are more expensive than those in Qinhuangdao and Yingkou.

5.4. Effect of $Lolp$ on Operational Costs

To study the impact of load loss probability $Lolp$ on operating costs and dispatch results, 11 indicators in the interval [0.02, 1] are selected for simulation. The calculation results are presented in Figure 7, and the system’s total operating cost decreases from USD 4281 to USD 3516 as the reliability level gradually decreases. The higher level of reliability imposes more stringent requirements on the spinning reserve of the system at different moments. It is equivalent to transforming the chance constraint into a deterministic constraint when $Lolp$ = 1, i.e., it does not consider the impact of load randomness and unexpected outage events on the system power balance.

For observation purposes, the variation of the system reserve power for six different loss-of-load probability indicators is analyzed. At $Lolp$ = 1, the reserve power at each moment is zero, and the results are displayed in Figure 8. At t > 15, the shore power supplies energy to the ship’s power system. Thus, it is considered that there is no possibility of load loss. As reliability levels increase, power reserves consequently rise to meet operational risk constraints. In addition, the towing condition requires a higher spinning reserve than cruising. Due to the heavy load factor under towing conditions, the risk of the system being affected by uncertainties increases, which requires focused attention.

5.5. Comparison of Different Piecewise Linearization Methods

The two linearization methods are tested for comparison. As speed variations will cause power fluctuations and thus affect the optimization results, pre-set nominal ship speeds are adopted for the experiments to exclude their effects. The actual values of the nonlinear function at nominal speed are used as the propulsion load, and the optimization result is recorded as the model baseline value. We compare the relative error between the optimization results of the two methods and the benchmark value under different numbers of segments from 3–10.

The results of the experimental comparison are shown in

Table 6. Comparison of results for different numbers of segments.

Segments	Optimized Spaced	Equally Spaced	Time	Variables
3	14.062%	38.244%	2.03	8145
4	9.337%	10.372%	2.13	8337
5	3.886%	9.813%	2.67	8529
6	3.383%	8.750%	5.02	8721
7	1.566%	6.206%	6.11	8913
8	1.510%	2.404%	9.52	9105
9	1.398%	2.153%	11.63	9297
10	0.587%	2.069%	15.83	9489

. The optimized spaced method in the table represents the relative error calculated using the improved piecewise linearization method. Equally spaced indicates the result of a piecewise linearization method using equally spaced segments. The accuracy of the fitting of the two methods increases with the number of segments. For convenience in measuring the experimental results, this paper uses a relative error of 2% as the standard. The table shows that the relative error of the improved piecewise linearization method is 1.566% at a segment number of 7. In contrast, the equally spaced method is still above 2% at a segment number of 10. With a segment number of 7, there are 576 fewer variables to optimize in the model than when the segment number is 10. Therefore, if the optimized spaced method with the number of segments set to 7 is used, the decision variables are reduced by 6.07% compared to the equally spaced method with the number of segments set to 10. Moreover, the former method is more precise in its calculation. The effectiveness of the presented approach for balancing computational resources and model accuracy is verified.

6. Conclusions

This paper combines the operating conditions of electric tugboats and proposes a new joint optimal scheduling model of navigation and power generation, considering the reliability of tugboat operation. The relationship between propulsion load and voyage as influenced by navigation speed and waves is first established. The risk scenarios (equipment outages and load fluctuations) that occur in the system are transformed into readily solvable mixed integer problems solved by probability distribution discretization and an improved piecewise linearization method. Finally, it is calculated by the CPLEX solver. The effect of various influencing factors on the scheduling results is verified and compared through four simulation experiments so that operators can adjust the navigation speed and generation scheme according to different demands.

This paper innovatively incorporates operational reliability opportunity constraints into the joint optimization to balance economy and reliability by adequately scheduling system reserve power. The simulation results show that the system can satisfy the power increment demand caused by wind and waves by enhancing the reserve power after adding the operational risk constraint and guaranteeing the operational safety of the ship throughout the voyage. Furthermore, the effect of segmentation points on precision is not studied in the traditional piecewise linearization method. In this paper, the improved linearization method with optimized segment points reduces the decision variables by 6.07% while ensuring the accuracy of the operation and improving the computational efficiency. Meanwhile, the economic indicator of battery degradation added to the optimization objective effectively avoids frequent charging and discharging, and the results show that the battery life loss is decreased.

The proposed scheduling model can be improved to consider other practical factors. The configuration and sizing of the energy storage system, combined with the tugboat’s working conditions, is essential. These will be further investigated.