Speed Optimization for Container Ship Fleet Deployment Considering Fuel Consumption

Abstract

In recent years, low energy consumption has become the common choice of economic development in the world. In order to control energy consumption, shipping line speed optimization has become strategically important. to reduce fuel consumption, this study optimizes the container ship fleet deployment problem by adopting the strategy of adjusting each leg of each route’s sailing speed. To calculate fuel consumption more accurately, both sailing speed and the ship’s payload are considered. A multi-objective mixed integer nonlinear programming model is established to optimize the allocation of liner routes with multiple ship types on multiple routes. A linear outer-approximation algorithm and an improved piecewise linear approximation algorithm are used for linearization. If segments of an interval increase, the results will be more accurate but will take more time to compute. As fuel prices increase, to make trade-offs among economic and environmental considerations, the shipping company is adopting the “adding ship and slow down its speed” strategy, which verifies the validity and applicability of the established model.

1. Introduction

Although maritime transport is considered to be an environment-friendly mode of transport, due to over 80% of the world’s trade being carried by maritime transport [1], maritime transport has become a significant contributor to global greenhouse gas (GHG) emissions [2]. As GHG emissions of a ship are strongly related to its sailing speed [3], they can be reduced by optimizing sailing speed [4]. Hence, speed optimization has been proven to be one of the most effective operational measures to reduce GHG emissions.

In the context of GHG emissions reduction, shipping companies need to consider GHG emissions reduction when they make decisions on container ship fleet deployment (CSFD). For example, Maersk attaches great importance to reducing emissions [5]. GHG emissions from ships are directly proportional to fuel consumption [6]. The fuel consumption rate is a nonlinear function of sailing speed [7]. Thus, sailing speed optimization is one of the possible solutions to reduce ship fuel consumption. At the same time, the deployed number of ships on the route is dependent on sailing speed. The faster the speed, the fewer ships are configured in the route, and the slower the speed is, the more ships the route needs to be equipped. Therefore, shipping companies adopt the strategy of speed optimization to achieve GHG emission reduction and reduce ship operating costs, which is an effective means for the survival and profitability of ship owners and operators.

Liner shipping companies operate several container transportation routes and provide container transportation services for designated ports regularly. Container ships return to and from designated ports. When taking minimum cost as the objective function, the shipping companies are determined to reduce speed for each leg in the route so that the fuel consumption of their ships is minimized. However, this speed reduction causes the sailing time between ports to increase and requires more ships to be deployed on the route, which increases ship operating costs. Accordingly, the best way to balance operating costs and fuel consumption costs depends on sailing speed.

To bring flexibility to fleet deployment decisions, this paper considers sailing speed as a decision variable to make trade-offs among economic and environmental considerations. Regarding literature studies on fuel consumption, most studies only considered the effect of ship sailing speed, and few studies considered influence by payload (weight of containers on the ship) of a ship. In our work, to minimize the total costs, a mixed-integer nonlinear program (MINP) with container transshipment is proposed, which takes fuel consumption as the objective function and considers the influence of ship payload and ship speed. Meanwhile, for each route, we optimize each leg sailing speed, and we allocate appropriate ship types and deploy suitable number.

The remainder of this paper is organized as follows. Section 2 makes reviews the previous works. Section 3 develops a variation of an existing container ship deployment model, and the objective of reducing emissions. Section 4 proposes a series of techniques to a convert nonlinear model. Section 5 reports on computational experiments. Section 6 is the conclusion of the paper.

2. Literature Review

Container fleet deployment problem has become pressing for shipping lines. Christiansen et al. [8] claimed that shipping lines need to deploy suitable numbers and types of ships in order to maximize profits. Therefore, there are more and more researches focusing on container fleet deployment problems. Jaramillo and Perakis [9] addressed the fleet deployment problem by developing an integer linear programming model, but they did not consider container transshipment operations. As port-to-port cargo transportation may have many routes, it is necessary to consider container transshipment for the container fleet deployment problem. Gelareh and Pisinger [10] considered container transshipment but did not consider routing optimization because there is only one leg for each port pair. Wang and Meng [11] established a mixed-integer programming model, which incorporated container transshipment and routing optimization to realize container transshipment at any port and any number of times. Considering container transshipment operations will increase computational complexity, this is very realistic model.

In most of the existing container fleet deployment [12,13] models, it is assumed that the sailing speed is independent of the optimal fleet deployment decisions and set as constant, but in reality, the sailing speed is variable. Gelareh and Meng [14] studied the CSFD problem and proposed a MIP model with speed as a decision variable. Andersson et al. [15] studied CSFD problem with optimizing sailing speed simultaneously and established an integrated optimization model for RoRo shipping lines. Wang et al. [16] considered market uncertainty optimally regarding the CSFD problem with sailing speed optimization and developed a two-stage stochastic programming model. Zhen et al. [17] showed that to make ensure regular service frequency, shipping lines need to deploy more containerships when reducing sailing speed. In other words, the speed of the ship determines the number of ships deployed. Norstad et al. [18] stated that sailing speed is an important planning decision regarding the CSFD problem. Hence, it is essential to consider speed optimization for fleet deployment of a shipping network.

The speed of ships impacts fuel consumption costs, which in turn affects GHG emissions. Zacharioudakis et al. [19] proposed a practical solution to solve CSFD problem, which considered the effect of speed on fuel consumption. Du et al. [20] aimed to minimize fuel consumption by optimizing ship sailing speed and trim. Most studies assume that ship fuel consumption only depends on speed. Ronen [21] researched the daily fuel consumption function approach to the cubic of sailing speed. Wang and Meng [22] found the relationship between speed and fuel consumption by linear regression and calculated the fuel consumption function, which is entirely related to shipping speed (the index is between 2.7 and 3.3). Hence, it is essential to consider speed optimization for reducing fuel consumption.

However, in practice, various other factors affect a ship’s fuel consumption, such as ship hull conditions and payload. Ignoring the payload may lead to significant fuel consumption estimation errors. Research by Gkonis and Psaraftis [23] showed that at the same speed, ship fuel consumption was different between laden and ballast. The difference is approximately 25–30%. N.Psaraftis and A.Kontovas [24] provided a literature review of sailing speed optimization issues and showed that when speed was given, the fuel consumption of a ship varies with its payload. For example, the fuel consumption is different when the ship is full, in an intermediate loading condition, or empty. Therefore, it is necessary to comprehensively consider the influence of sailing speed and ship payload for reducing fuel consumption.

To accurately estimate fuel consumption, Xia et al. [25] used a nonlinear function to represent the relationship between ship payload and ship fuel consumption. However, there is a significant approximation error in the nonlinear approximation method, which may have seriously affected the quality of their solution. Wang et al. [26] studied shipping revenue management by adopting the fuel consumption rate function depending on sailing speed and ship payload, which developed a nonlinear fuel consumption rate function. However, they did not consider selecting the appropriate ship type according to the container flow. Zhen et al. [17] developed a multi-objective optimization model to solved CSFD problem considering the influence of container transshipment and fuel consumption. However, they assumed that fuel consumption cost only depended on sailing speed, and each port operation time was given.

To overcome the aforementioned limitations, our work is based on Wang and Meng’s [11] proposed model, which assumed that each leg of a ship’s sailing speed of route is deterministic. To relax this assumption and optimize each leg sailing speed of route, we put sailing speed optimization incorporate into CSFD problem models as a decision variable. As the number of deployed ships and fuel consumption depends on sailing speed, in order to balance economic and environmental considerations, we proposed a multi-objective mixed integer programming model to optimize sailing speed, which minimizes the total costs and total fuel consumption cost. The proposed model was nonlinear, and we use liner techniques to linearize it.

In this article, we primarily contribute to the relevant literature in the following ways. First, we investigate the theoretical and practical implications of fleet deployment that considers container transshipment and optimizes ships’ speeds. We propose a mixed-integer nonlinear program and linearize the nonlinear expressions in the developed models. Second, a linear outer-approximation algorithm and an improved piecewise linear approximation algorithm are applied to further approximate nonlinear functions. Thus, we develop a mixed-integer linear program whose problems can be solved by on-the-shelf solvers.

3. Mathematical Model Formulation

Shipping lines regularly serve a set of ports. The ports’ order was given and formed a closed loop, which provided customers with circular liner transportation services and ensured weekly flights. The shipping lines allocated the appropriate number and types of ships on the route. Figure 1 shows the shipping network of composition by three routes. Containers attached to ports on different routes require transshipment at public ports to meet the transportation demand of OD flows. For example, container transshipment operations of route 1 and route 3 in Figure 1 can be performed at ports 1 and 3. However, the transshipment operation at the transit port will increase the handling capacity of containers, thus increasing the transshipment cost.

3.1. Parameter and Variable Definition

The parameter and variable definitions in this study are introduced in Table 1.

3.2. Fuel Consumption Cost

N. Psaraftis and A. Kontovas [24] proposed the fuel consumption non-linear function, which depends on the ship sailing speed and ship payload, in which the daily bunker consumption rate is expressed as $f\left(s,w\right)$, $s$ is the sailing speed, and $w$ is the ship payload of each segment. A realistic approximation of fuel consumption function is $f\left(s,w\right)=c_1^r\left(q+s^{c_2^r}\right){\left(w+A\right)}^{c_3^r}$ with $c_1^r,c_2^r,c_3^r$ and $q$ constants such as $c_1^r>0,c_2^r\ge 3,c_3^r\ge 0$ and $q\ge 0$, and $A$ is the weight of the empty container ship. As the sailing speed and payload of different segments of each route are different, we assume that $q=0$ and $A=0$; the daily fuel consumption function of leg $i$ on route $r$ is equal to $f\left(s_{ri},w_{ri}\right)=c_1^r{\left(s_{ri}\right)}^{c_2^r}{\left(w_{ri}\right)}^{c_3^r}$. As far as we know, the leg $i$ length of route $r$ is $L_{ri}$, Then, the leg $i$ sailing time of route $r$ can be calculated by $L_{ri}/s_{ri}$. $\beta$ is the fuel price, 24 is the hour of the day, and the ship fuel consumption cost expressed as:

(1)$f^{FC}=\beta ·{{\displaystyle \sum }}_{\begin{array}{c}r\in R\\ i\in I_r\end{array}}{c}_1^r{\left(s_{ri}\right)}^{c_2^r-1}{\left(w_{ri}\right)}^{c_{\mathrm{a}}^r}·L_{ri}/24$

The developed model in this study optimizes the sailing speed of each leg. As the number of deployed ships and fuel consumption depend on sailing speed, it is beneficial for container shipping lines to optimize the sailing speed of each leg when they decide to minimize the total fuel consumption cost and total operating costs. This study adopts a realistic function of ship fuel consumption rate; optimizing the sailing speed of each leg allows the fuel consumption cost of ships in sailing to be more precisely estimated, which improves the management accuracy for liner shipping companies.

3.3. Container Ships’ Fleet Deployment Model

The objective function of the model is to minimize the cost of ship fuel consumption and the total operating cost of the fleet in the weekly plan, which is defined as a Constraint (2)–(8). Where $f^{CF}$ is expressed as a fixed cost of ship navigation, $f^{CB}$ is expressed as occupancy cost of berthing berth in port, $f^{CT}$ is expressed as total container transshipment cost, $f^{CO}$ is expressed as container handling cost, $f^{CI}$ is expressed as the cost of renting ships, and $f^{RO}$ is expressed as proceeds from charter-owned ships.

(2) $\mathrm{Minimize}f=f^{FC}+\left(f^{CF}+f^{CB}+f^{CT}+f^{CO}+f^{CI}-f^{RO}\right)$

(3) $f^{CF}={{\displaystyle \sum }}_{r\in R,v\in V_r}\left(m_{r,v}·C_v^{opr}+C_{r,v}^{fix}·x_{r,v}\right)$

(4) $f^{CB}={{\displaystyle \sum }}_{r,i\in I_r,v\in V_r}\left(C_{P_{r,i},v}^{ber}·T_{P_{r,i},v}^{cont}·z_{r,v,i}\right)$

(5) $f^{CT}={{\displaystyle \sum }}_{p\in P}\left(C_p^{tran}·z_p^{tran}\right)$

(6) $f^{CO}={{\displaystyle \sum }}_{o,d\in P}\left(\left(C_o^{load}+C_d^{disc}\right)·D_{o,d}\right)$

(7) $f^{CI}={{\displaystyle \sum }}_{v\in V}\left(C_v^{in}·n_v^{in}\right)$

(8) $f^{RO}={{\displaystyle \sum }}_{v\in V}\left(C_v^{out}·n_v^{out}\right)$

s.t.

(9) ${{\displaystyle \sum }}_{v\in V_r}{x}_{r,v}=1,\forall r$

(10) $m_{r,v}\le M^1·x_{r,v},\forall r,v$

(11) $168·m_{r,v}+M^2·\left(1-x_{r,v}\right)\ge {{\displaystyle \sum }}_{i\in I_r}{L}_{ri}/s_{ri}+{{\displaystyle \sum }}_{i\in I_r}\left(T_{P_{r,i},v}^{cont}·z_{r,v,i}\right),\forall r,v\in V_r$

(12) $z_{v,r,i}\le M^3·x_{r,v},\forall r,i\in I_r,v\in V_r$

(13) $z_{v,r,i}+M^4·\left(1-x_{r,v}\right)\ge {{\displaystyle \sum }}_o\left(z_{o,r,i}^{load}+z_{o,r,i}^{disc}\right),\forall r,i\in I_r,v\in V_r$

(14) ${{\displaystyle \sum }}_o{z}_{o,r,i}^{flow}-{{\displaystyle \sum }}_{v\in V_r}{Q}_v·x_{r,v}\le 0,\forall r,i\in I_r$

(15) $z_{o,r,i-1}^{flow}+z_{o,r,i}^{load}-z_{o,r,i}^{disc}=z_{o,r,i}^{flow},\forall r,i\in I_r,o$

(16) ${{\displaystyle \sum }}_{r\in R_i,i\in I_{r,i}}\left(z_{o,r,i}^{disc}-z_{o,r,i}^{load}\right)=D_{o,d},\forall o\ne d$

(17) $z_{o,r,i}^{flow}=0,\forall r,i\in I_r,o=P_{r,i+1}$

(18) $z_{o,r,i}^{disc}=0, ,\forall r,i\in I_r,o=P_{r,i}$

(19) ${{\displaystyle \sum }}_{r\in R_p,i\in R_p}{m}_{r,v}-n_v^{in}+n_v^{out}=N_v^{own},\forall v$

(20) $n_v^{in}\le N_v^{in},\forall v$

(21) $z_p^{tran}={{\displaystyle \sum }}_{r\in R_r,i\in I_{r,p},o\in P}\left(z_{o,r,i}^{load}+z_{o,r,i}^{disc}\right)-{{\displaystyle \sum }}_d{D}_{p,d}-{{\displaystyle \sum }}_o{D}_{o,p},\forall p$

(22) $w_{r,i}^{in}\le w_{r,i}\le w_{r,i}^{max},\forall r,i\in I_r$

(23) ${{\displaystyle \sum }}_o{z}_{o,r,i}^{flow}+e_{r,i}=w_{r,i},\forall r,i\in I_r$

(24) $s_{r,i}^{min}\le s_{r,i}\le s_{r,i}^{max},\forall r,i\in I_r$

Constraint (9) ensures that each route deployed only one type of ship; Constraints (10) and (11) indicate that the number of ships configured for each route can meet the service demand in the weekly plan, 168 represents the number of hours in a week, set ${\mathrm{M}}^1$ and ${\mathrm{M}}^2$ to $9999$; Constraints (12) and (13) define the relationship between variables, ${\mathrm{M}}^3$ and ${\mathrm{M}}^4$ set to $9999$; in Constraint (14) the container transportation volume of each leg cannot exceed the ship’s carrying capacity; and Constraint (15) represents the flow of conservation. The carrying capacity before berthing at port $i$, plus the amount of containers loading and subtract the amount of unloading, is equal to the shipping capacity leaving the port; Constraint (16) ensures that all container transportation needs between ports are served; Constraint (17) indicates containers from the port $o$ will not be shipped back to the port $o$; Constraint (18) indicates containers from the port $o$ will not be unloaded at the port $o$; Constraint (19) indicates conservation of the number of ships; Constraint (20) indicates that the number of chartered ships cannot exceed the upper limit; Constraint (21) indicates the weight of ballast water required to ensure the stability of the ship’s navigation; Constraint (22) is that for safety or other purposes, the ship’s payload should be within the range; Constraint (23) the weight of ballast water needed to ensure the stability of ship navigation; and Constraint (24) is the range of ship speed.

4. Linearization of the Model

The above model for fleet deployment with sailing speed and routing optimization has a nonlinear objective (1) with nonlinear terms ${\left(s_{ri}\right)}^{c_2^r-1}{\left(w_{ri}\right)}^{c_3^r}$ and nonlinear Constraint (11) with nonlinear term $L_{ri}/s_{ri}$. Due to the following reasons, we linearize these formula expressions. First, although some mixed-integer program solvers (e.g., Gurobi) can handle nonlinear objective functions, it is challenging to solve models with complicated nonlinear objectives, especially nonlinear constraints. Second, customized linearization techniques are applicable, as studied by many pioneering studies in the literature. Third, after preliminary study and tests, the devised model in this study is not solvable directly when the nonlinear objective and constraints are considered. Therefore, we adopted the linearization technique to convert the model to a linear one. The logical structure is expressed as: first, the reciprocal of the speed to linearize nonlinear Constraint (11). Then, linearize the ship’s fuel consumption objective Function (1). Finally, perform bilinear optimization on the Constraint (30).

4.1. Linearization of the Reciprocal of Sailing Speeds

Due to $1/s_{ri}$ that is nonlinear, take its reciprocal as the decision variable $u_{ri}$. In general, the sailing speed is different on each leg of the route, and within the speed interval $\left. [ s_{ri}^{min},s_{ri}^{max}\right]$, this speed interval depends on mechanical power engine size and economic factors. Therefore, $u_{ri}$ lower bound is $u_{ri}^{min}=1/s_{ri}^{max}$ and upper bound is $u_{ri}^{max}=1/s_{ri}^{min}$. Then, the Constraint (11) becomes a linear constraint and can be rewritten by Constraint (27)

(25)$u_{ri}=1/s_{ri}$

(26)$u_{ri}^{min}\le u_{ri}\le u_{ri}^{max},r,i\in I_r$

(27)$168·m_{r,v}+M^2·\left(1-x_{r,v}\right)\ge {{\displaystyle \sum }}_{i\in I_r}{L}_{ri}·u_{ri}+{{\displaystyle \sum }}_{i\in I_r}\left(T_{P_{r,i},v}^{cont}·z_{r,v,i}\right),\forall r,v\in V_r$

4.2. Linearization of the Objective Function of Fuel Consumption Cost

The objective function of fuel consumption cost can be expressed as follows:

(28)$\beta ·{{\displaystyle \sum }}_{\begin{array}{c}r\in R\\ i\in I_r\end{array}}{c}_1^r{\left(u_{ri}\right)}^{1-c_2^r}·{\left(w_{ri}\right)}^{c_3^r}·L_{ri}/24$

The fuel consumption cost objective Function (28) can be further linearization by introducing the auxiliary variables: $B_{ri}$, ${\nabla }_{ri, }$ and ${\mathsf{\Lambda }}_{ri}$, transformed to the following objective Function (29) and Constraints (30)–(32).

(29)$\beta ·{{\displaystyle \sum }}_{r\in R,i\in I_r}{B}_{ri}$

(30)$24·B_{ri}/\left(c_1^r·L_{ri}\right)\ge {\nabla }_{ri}·{\mathsf{\Lambda }}_{ri}$

(31)${\nabla }_{ri}\ge {\left(u_{ri}\right)}^{1-c_2^r},\forall r,i\in I_r$

(32)${\mathsf{\Lambda }}_{ri}\ge {\left(w_{ri}\right)}^{c_3^r},\forall r,i\in I_r$

4.3. Underestimating Bilinear Terms

Al-Khayyal and Falk [27] showed that the bilinear term $x·y$ tightest convex lower bound is obtained via the following relationship $max\left\{x^{L}y+y^{L}x-x^L{y}^L, x^{U}y+y^{U}x-x^U{y}^U\right\}$ in the interval $\left. [x^L, x^U\right]\times \left. [y^L,y^U\right]$. Due to the fact that ${\nabla }_{ri}$ and ${\mathsf{\Lambda }}_{ri}$ have the lower limit and upper limit, Wang et al. [26] showed that ${\nabla }_{ri}·{\mathsf{\Lambda }}_{ri}$ term is bilinear. The lower bound for ${\nabla }_{ri}·{\mathsf{\Lambda }}_{ri}$ in the interval $\left. [{\left(u_{ri}^{max}\right)}^{1-c_2^r},{\left(u_{ri}^{min}\right)}^{1-c_2^r}\right]$ and $\left. [{\left(w_{ri}^{min}\right)}^{c_3^r},{\left(w_{ri}^{max}\right)}^{c_3^r}\right]$ is obtained by following Constraint (33)

(33)${\nabla }_{ri}·{\mathsf{\Lambda }}_{ri}=max\left. {\begin{array}{c}{\left(u_{ri}^{max}\right)}^{1-c_2^r}·{\mathsf{\Lambda }}_{ri}+{\left(w_{ri}^{min}\right)}^{c_3^r}·{\nabla }_{ri}-{\left(u_{ri}^{max}\right)}^{1-c_2^r}·{\left(w_{ri}^{min}\right)}^{c_3^r}\\ {\left(u_{ri}^{min}\right)}^{1-c_2^r}·{\mathsf{\Lambda }}_{ri}+{\left(w_{ri}^{max}\right)}^{c_3^r}·{\nabla }_{ri}-{\left(u_{ri}^{min}\right)}^{1-c_2^r}·{\left(w_{ri}^{max}\right)}^{c_3^r}\end{array}\right\}$

Therefore, Constraint (30) can be relaxed by adding linear Constraint (34) and Constraint (35):

(34)$24·B_{ri}/\left(c_1^r·L_{ri}\right)\ge {\left(u_{ri}^{max}\right)}^{1-c_2^r}·{\mathsf{\Lambda }}_{ri}+{\left(w_{ri}^{min}\right)}^{c_3^r}·{\nabla }_{ri}-{\left(u_{ri}^{max}\right)}^{1-c_2^r}·{\left(w_{ri}^{min}\right)}^{c_3^r}$

(35)$24·B_{ri}/\left(c_1^r·L_{ri}\right)\ge {\left(u_{ri}^{min}\right)}^{1-c_2^r}·{\mathsf{\Lambda }}_{ri}+{\left(w_{ri}^{max}\right)}^{c_3^r}·{\nabla }_{ri}-{\left(u_{ri}^{min}\right)}^{1-c_2^r}·{\left(w_{ri}^{max}\right)}^{c_3^r}$

5. Approximation Algorithm

The logical structure of the approximation algorithm is expressed as: first, a linear outer-approximation algorithm developed to handle nonlinear Constraint (31). Then, an improved piecewise linear approximation algorithm was applied to linearization nonlinear Constraint (32). Last, it is converted into a linear mixed-integer programming model.

5.1. Linear Outer-Approximation Algorithm

The principle of the algorithm is to divide the nonlinear characteristic curve into several linear segments and approximate it by a series of tangent lines. Therefore, the interval $\left. [u_{ri}^{min},u_{ri}^{max}\right]$ is subdivided into several segments. For each segment, the function ${\left(u_{ri}\right)}^{1-c_2^r}$ is approximated by a tangent line.

Proposition 1.

${\left(u_{ri}\right)}^{1-c_2^r}$ is convex in $u_{ri}$

Proof. 

The first-order derivatives of ${\left(u_{ri}\right)}^{1-c_2^r}$ are denoted by ${\left({\left(u_{ri}\right)}^{1-c_2^r}\right)}^{\prime }$ and can be calculated as:

${\left({\left(u_{ri}\right)}^{1-c_2^r}\right)}^{\prime }= \left(1-c_2^r\right) · {\left(u_{ri}\right)}^{-c_2^r}$

The second-order derivative of ${\left(u_{ri}\right)}^{1-c_2^r}$ with respect to $u_{ri}$ is denoted by ${\left({\left(u_{ri}\right)}^{1-c_2^r}\right)}^{″}$ and can be calculated as:

${\left({\left(u_{ri}\right)}^{1-c_2^r}\right)}^{″}=-\left(1-c_2^r\right) · c_2^r · {\left(u_{ri}\right)}^{-c_2^r-1}$

According to Wang and Meng [22], the coefficient $c_2^r$ is greater than 2, hence the second-order derivative ${\left({\left(u_{ri}\right)}^{1-c_2^r}\right)}^{″}>0$. When on the interval $\left. [u_{ri}^{min},u_{ri}^{max}\right]$, ${\left(u_{ri}\right)}^{1-c_2^r}$ is convex. □

In order to approximate ${\left(u_{ri}\right)}^{1-c_2^r}$, Figure 2 shows the tangent lines. To control the approximation error ${\epsilon }_1$, some tangent lines are generated a priori. Then, whether the calculation gap is larger than ${\epsilon }_1$ is determined. If so, the interval is further subdivided to obtain more approximate tangent lines. Otherwise, the approximation gap is acceptable, and the tangent point is recorded in the set ${\mathsf{\Omega }}_{ri}$.

At tangent point $u_{ri}^k$, the tangent line of the function ${\left(u_{ri}\right)}^{1-c_2^r}$ can be derived as:

(36)${\left(u_{ri}^k\right)}^{1-c_2^r}+{\left({\left(u_{ri}^k\right)}^{1-c_2^r}\right)}^{\prime }\times \left(u_{ri}-u_{ri}^k\right),\forall r\in R,i\in I_r,u_{ri}^k$

where $\mathsf{\Omega }\left({\epsilon }_1\right)$ is the tangent points set of function ${\left(u_{ri}\right)}^{1-c_2^r}$, when given ${\epsilon }_1$, the procedure to obtain set $\mathsf{\Omega }\left({\epsilon }_1\right)$ is shown in Algorithm 1. Then, the nonlinear Constraint (31) can be expressed as:

(37)${\nabla }_{ri}\ge \left(1-c_2^r\right){\left(u_{ri}^k\right)}^{-c_2^r}{u}_{ri}+c_2^r{\left(u_{ri}^k\right)}^{1-c_2^r},\forall r\in R, i\in I_r,u_{ri}^k\in {\mathsf{\Omega }}_{ri}\left({\epsilon }_1\right)$

Algorithm 1. The procedure of linear outer-approximation algorithm to obtain tangent point sets.
Input:	Convex function ${\left(u_{ri}\right)}^{1-c_2^r}$, the tangent point set ${\mathsf{\Omega }}_{ri}$. The lower limit and upper limit of $u_{ri}$ is $\left. [ u_{ri}^{min},u_{ri}^{max}\right]$, the approximation relative error ${\epsilon }_1$
Output:	the tangent point set ${\mathsf{\Omega }}_{ri}\left({\epsilon }_1\right)$
Step 1	$\mathrm{Calculate}\mathrm{the}\mathrm{midpoint}{u}_{ri}^{mid}$ of the interval $\left. [u_{ri}^{min},u_{ri}^{max}\right]$; $u_{ri}^{mid}$ can be carried out by the bisection search method
Step 2	At tangent point $u_{ri}^{mid}$, the tangent line is defined as ${\left(u_{ri}^{mid}\right)}^{1-c_2^r}+{\left({\left(u_{ri}^{mid}\right)}^{1-c_2^r}\right)}^{\prime }\times \left(u_{ri}-u_{ri}^{mid}\right)$
Step 3	Calculate the relative approximation error of points $u_{ri}^{min}$ and $u_{ri}^{max}$ according to ${\epsilon }_l={\left(u_{ri}^{min}\right)}^{1-c_2^r}-{\left(u_{ri}^{mid}\right)}^{1-c_2^r}-{\left({\left(u_{ri}^{mid}\right)}^{1-c_2^r}\right)}^{\prime }\times \left(u_{ri}^{min}-u_{ri}^{mid}\right)$ ${\epsilon }_r={\left(u_{ri}^{max}\right)}^{1-c_2^r}-{\left(u_{ri}^{mid}\right)}^{1-c_2^r}-{\left({\left(u_{ri}^{mid}\right)}^{1-c_2^r}\right)}^{\prime }\times \left(u_{ri}^{max}-u_{ri}^{mid}\right)$
Step 4	If ${\epsilon }_l>\epsilon$ or ${\epsilon }_r>\epsilon$, then branch the feasible range of $u_{ri}$ is divided into two ranges: $\left. [ u_{ri}^{min},u_{ri}^{mid}\right]$ and $\left. [ u_{ri}^{mid},u_{ri}^{max}\right]$, and the tangent point set ${\mathsf{\Omega }}_{ri}\left({\epsilon }_1\right)={\mathsf{\Omega }}_{ri}\left({\epsilon }_1\right) \cup \left. {u_{ri}^{mid}\right\}$
	In one branch $u_{ri}\in \left. [ u_{ri}^{min},u_{ri}^{mid}\right]$, repeat the above step 1 to e step 4 until the stop criterion is reached
	In the other branch $u_{ri}\in \left. [ u_{ri}^{mid},u_{ri}^{max}\right]$, repeat the above step 1 to e step 4 until the stop criterion is reached
	Stop criterion check: if ${\epsilon }_{\mathrm{l}}<\epsilon$ and ${\epsilon }_r<\epsilon$, stop and output the current solution. Otherwise, go to Step4.
Step 5	Return output ${\mathsf{\Omega }}_{ri}\left({\epsilon }_1\right)$

5.2. Improved Piecewise Linear Approximation Algorithm

The principle of the algorithm is to divide the nonlinear characteristic curve into several linear segments and approximate the characteristic curve in each segment with a straight-line section. Therefore, the interval $\left. [w_{ri}^{min},w_{ri}^{max}\right]$ is subdivided into several segments; for each segment, the function ${\left(w_{ri}\right)}^{c_3^r}$ is approximated by the linear function $l_K\left(w_{ri}\right)$.

Proposition 2.

$f\left(w_{ri}\right)={\left(w_{ri}\right)}^{c_3^r}$is a continuous functionon the interval$\left. [w_{ri}^{min},w_{ri}^{max}\right]$, then the piecewise linear function$l_K\left(w_{ri}\right)$is uniformly approximated to$f\left(w_{ri}\right)$.

Proof. 

Since $f\left(w_{ri}\right)$ is continuous on $\left. [w_{ri}^{min},w_{ri}^{max}\right]$, then $f\left(w_{ri}\right)$ is uniformly continuous on $\left. [w_{ri}^{min},w_{ri}^{max}\right]$. That is, for any given $\epsilon >0$, there exists $\delta \left(\epsilon \right)>0$; for any $w_{ri}^{\prime },w_{ri}^{″}\in \left. [w_{ri}^{min},w_{ri}^{max}\right]$, as long as $\left|w_{ri}^{\prime }-w_{ri}^{″}\right|<\delta \left(\epsilon \right)$, there is $\left. |f\left(w_{ri}^{\prime }\right)-f\left(w_{ri}^{″}\right)\right|<\epsilon$.

The interval $\left. [w_{ri}^{min},w_{ri}^{max}\right]$ is divided into $K$ segments; $a_j,\left(j=1,\cdots ,K,K+1\right)$ is expressed as the breakpoints of the piecewise linear function. Command $w_{ri}^{min}=a_1<a_2<a_3<\cdots <a_K<a_{K+1}=w_{ri}^{max}$. When $K>1/\delta \left(\epsilon \right)$, for $\forall w_{ri}^1,w_{ri}^2\in \left. [a_j,a_{j+1}\right]$, there is $\left. |f\left(w_{ri}^1\right)-f\left(w_{ri}^2\right)\right|<\epsilon$.

$min\left. {f\left(a_j\right),f\left(a_{j+1}\right)\right\}\le l_{K,j+1}\left(w_{ri}\right)\le max\left. {f\left(a_j\right),f\left(a_{j+1}\right)\right\}$,$l_{K,j+1}\left(w_{ri}\right)$ indicates the function section for the interval $\left. [a_j,a_{j+1}\right]$, that is, the $j+1$ segment linear function. So $min\left\{f_{j+1}\left(w_{ri}\right)-f\left(a_j\right),f_{j+1}\left(w_{ri}\right)-f\left(a_{j+1}\right)\right\}\le f_{j+1}\left(w_{ri}\right)-l_{K,j+1}\left(w_{ri}\right)\le max\left. {f_{j+1}\left(w_{ri}\right)-f\left(a_j\right),f_{j+1}\left(w_{ri}\right)-f\left(a_{j+1}\right)\right\}$, $f_{j+1}\left(w_{ri}\right)$ represents the corresponding function of $f\left(w_{ri}\right)$ on interval $\left. [a_j,a_{j+1}\right]$, which is $\left|f_{j+1}\left(w_{ri}\right)-l_{K,j+1}\left(w_{ri}\right)\right|<max\left. {\left. |f_{j+1}\left(w_{ri}\right)-f\left(a_j\right)\right|,\left. |f_{j+1}\left(w_{ri}\right)-f\left(a_{j+1}\right)\right|\right\}<\epsilon$, then for $\forall w_{ri}\in \underset{j=0}{\overset{n-1}{{\displaystyle \cup }}}\left. [a_j,a_{j+1}\right]$, there is $\left|f\left(w_{ri}\right)-l_K\left(w_{ri}\right)\right|<\epsilon$. Therefore, on interval $\left. [w_{ri}^{min},w_{ri}^{max}\right]$, there is a piecewise linear function $l_K\left(w_{ri}\right)$ that is uniformly approximated to the continuous function $f\left(w_{ri}\right)$. □

Therefore, this paper proposes an improved function piecewise linear approximation algorithm, which can complete the optimal selection of the number of pieces under the premise of known accuracy, so that the calculation function error is within a specific range. The main idea of the algorithm is: further subdivide the function interval with more significant errors, avoid the waste of segmentation of the smooth part of the curve, and reduce the computational complexity. The algorithm flow implemented in this paper is shown in Figure 3. The calculation of the segment interval is completed by an iterative method.

On each segment interval, the concavity of the approximated function and whether the function has an inflection point in the interval are judged. For interval $\left. [w_{ri}^{min},w_{ri}^{max}\right]$, the approximation lines can be generated as follows:

Function Generate approximate straight lines$l_K\left(w_{ri}\right)$ two endpoints of the curve $f\left(w_{ri}\right)={\left(w_{ri}\right)}^{c_3^r}$ in interval $\left. [w_{ri}^{min},w_{ri}^{max}\right]$ are $A=\left(w_{ri}^{min},{\left(w_{ri}^{min}\right)}^{c_3^r}\right)$ and $B=\left(w_{ri}^{max},{\left(w_{ri}^{max}\right)}^{c_3^r}\right)$, the equation of a straight line is found through these two points. Assuming the point $C=\left(w_{ri},l_K\left(w_{ri}\right)\right)$ is on a straight line, hence the slope of the line AC and the line AB are the same, expressed as:

$\left(l_K\left(w_{ri}\right)-{\left(w_{ri}^{min}\right)}^{c_3^r}\right)/\left(w_{ri}-w_{ri}^{min}\right) = \left({\left(w_{ri}^{max}\right)}^{c_3^r}-{\left(w_{ri}^{min}\right)}^{c_3^r}\right)/\left(w_{ri}^{max}-w_{ri}^{min}\right)$

By transforming the formula, we can get an approximately straight line:

$l_K\left(w_{ri}\right) = \frac{{\left(w_{ri}^{max}\right)}^{c_3^r}-{\left(w_{ri}^{min}\right)}^{c_3^r}}{w_{ri}^{max}-w_{ri}^{min}}·w_{ri}+\frac{w_{ri}^{max}·{\left(w_{ri}^{min}\right)}^{c_3^r}-w_{ri}^{min}·{\left(w_{ri}^{max}\right)}^{c_3^r}}{w_{ri}^{max}-w_{ri}^{min}}$

The approximation error between straight line and function is expressed as $d= \left. |f\left(w_{ri}\right) - l_K\left(w_{ri}\right)\right|$. Then, the approximation error is calculated between the approximate straight line $l_K\left(w_{ri}\right)$ and function ${\left(w_{ri}\right)}^{c_3^r}$ at the midpoint $w_{ri}^{mid}$. The concavity of the function ${\left(w_{ri}\right)}^{c_3^r}$ is judged; the following lemma proves the function property is nonconvex:

Proposition 3.

${\left(w_{ri}\right)}^{c_3^r}$ is concave in $w_{ri}$

Proof. 

The second-order derivative of ${\left(w_{ri}\right)}^{c_3^r}$ with respect to $w_{ri}$ denoted by ${\left({\left(w_{ri}\right)}^{c_3^r}\right)}^{″}$ can be calculated by constraint

${\left({\left(w_{ri}\right)}^{c_3^r}\right)}^{″}=c_3^r\left(c_3^r-1\right){\left(w_{ri}\right)}^{c_3^r-2}$

According to the N.Psaraftis and A.Kontovas [24], the coefficient $c_3^r$ is between 0 and 1; hence, the second-order derivative ${\left({\left(w_{ri}\right)}^{c_3^r}\right)}^{″}<0$. So, the function ${\left(w_{ri}\right)}^{c_3^r}$ is concave when $w_{ri}$ takes value in interval $\left. [w_{ri}^{min},w_{ri}^{max}\right]$. □

Therefore, the function ${\left(w_{ri}\right)}^{c_3^r}$ is concave, then the approximation line is moved to $l_K\left(w_{ri}\right) + d/2$, shown in Figure 4. The maximum error is found between the function and the new approximation line on the interval $\left. [ w_{ri}^{min},w_{ri}^{max}\right]$, expressed as $d_{max} = \left. |f\left(w_{ri}\right)-\left(l_K\left(w_{ri}\right) + d/2\right)\right|$. Then, the derivative of the function $d_{max}$ is taken and the derivative is made equal to 0, to obtain the value of $w_{ri}^{*}$. At last, $w_{ri}^{*}$ is substituted into the function $d_{max}$. The maximum error between the approximation line and function is obtained and then compared with the given error ${\epsilon }_2$: if the maximum error $d_{max}>{\epsilon }_2$, the interval of the approximation is further subdivided and the breakpoints of the piecewise linear function set ${\mathrm{W}}_{ri}={\mathrm{W}}_{ri}\cup \left. {w_{ri}^{mid}\right\}$, then we branch the feasible range of $w_{ri}$ into two ranges: $\left. [ w_{ri}^{min},w_{ri}^{mid}\right]$ and $\left. [ w_{ri}^{mid},w_{ri}^{max}\right]$. If the maximum error $d_{max}<{\epsilon }_2$, the continuous subdivision of the segment is stopped.

$w_{ri}$ is divided into $K$ segments in the interval $\left. [w_{ri}^{min},w_{ri}^{max}\right]$, which is determined by improved piecewise linear approximation algorithm, the breakpoints of the piecewise linear function set expressed as $W_{ri}\left(K\right) = \left. {w_{r,i,k},k\in M\right\}$, where $M= \left. {1,2,\cdots \left. |K\right|,\left. |K+1\right|\right\}$; by introducing the auxiliary variable ${\eta }_{rik}$ and ${\pi }_{rik}$, then, the nonlinearity Constraints (32) can be rewritten as:

(38)${\mathsf{\Lambda }}_{ri}\ge {{\displaystyle \sum }}_{k\in M}{\eta }_{r,i,k}·{\left(w_{r,i,k}\right)}^{c_3^r},\forall w_{r,i,k}\in W_{ri}\left(K\right),r\in R, i\in I_r$

(39)${\mathrm{w}}_{ri}={{\displaystyle \sum }}_{k\in M}{\eta }_{r,i,k}·w_{r,i,k},\forall w_{r,i,k}\in W_{ri}\left(K\right),r\in R, i\in I_r$

(40)${{\displaystyle \sum }}_{k\in M}{\eta }_{r,i,k}=1,\forall r\in R, i\in I_r$

(41)${{\displaystyle \sum }}_{k\in M-\left. |K+1\right|}{\pi }_{r,i,k}=1,\forall r\in R, i\in I_r$

(42)${\eta }_{r,i,1}\le {\pi }_{r,i,1},\forall r\in R, i\in I_r$

(43)${\eta }_{r,i,1}\le {\pi }_{r,i,1},\forall r\in R, i\in I_r$

(44)${\eta }_{r,i,k}\le {\pi }_{r,i,k-1}+{\pi }_{r,i,k},\forall r\in R, i\in I_{r}k=2,\cdots ,K$

(45)${\pi }_{r,i,k}=\left. {0,1\right\},k\in M-\left. |K+1\right|$

(46)${\eta }_{r,i,k}\ge 0,k\in M$

5.3. Mixed Integer Linear Programming Model

In summary, through the above series linearization operations, the nonlinear MINP model can become an equivalent mixed-integer linear programming MIIP.

(47) $\mathrm{Minimize}f=\beta ·{{\displaystyle \sum }}_{r\in R,i\in I_r}{B}_{ri}+\left(f^{CF}+f^{CB}+f^{CT}+f^{CO}+f^{CI}-f^{RO}\right)$

s.t.  Constraints (3)–(10), Constraints (12)–(24), Constraints (25)–(27), Constraints (34)–(35), and Constraints (37)–(46).

6. Numerical Experiments

6.1. Parameter Setting

The numerical experiment data are derived from reference Wang and Meng [11], which contains 46 ports and 12 routes. Table 2 shows the values of parameters for different ships. Table 3 shows parameters related to the 12 ship routes. Table 4 shows the ship fixed operating costs of 12 ship routes. The value range of container flow between ports is $U\left. [0,30\right]$. The loading cost $C_p^{load}$ and discharge cost $C_p^{disc}$ are 150; the container transshipment cost is 200. For the value of $c_1^r$,$c_2^r$, and $c_3^r$ from literature Wang et al. [26], the coefficient $c_1^r$ is 0.0006, $c_2^r$ is 2.5, and $c_3^r$ is 0.56. The model was implemented in Python version 3.8 and solved by using GUROBI version 8.1.1. We carry out calculation experiments run by an Intel 2.7 GHz core i3 CPU and 16GB RAM machine with Windows 10.

Wang et al. ’s [26] research shows a small approximation error means a large number of segments and more computation time, so they set minimum linear outer-approximation error at $1\times {10}^{-3}$. Yang and Xing [28] set the number of line fragments for linear outer-approximation to be 100. We set the number of line fragments for linear outer-approximation to be 68 and 135. For piecewise linear function approximation error, Wang et al. [26] set the number of line segments 8. To obtain higher accuracy, we set the number of line segments to 11 and 22, providing a more accurate optimal solution. We composed four scenarios according to the combination of the different approximation error ${\epsilon }_1$ and piecewise linear function approximating error ${\epsilon }_2$, as shown in Table 5.

6.2. Sensitivity Analysis of Various Fleet Costs

According to the four scenarios composed of different error accuracy, the total fleet operation cost is calculated, along with the ship consumption cost, the total cost of berthing operation, and the total cost of container transshipment under four scenarios. Results of Scenario 2 compared to Scenario 1 and Scenario 4 compared to Scenario 3 are shown in Table 6:

The negative values of ${∆}_1$ and ${∆}_2$ show that the number of segments within the range of ship’s payload increased, as shown in Table 6. The same as ${∆}_3$ and ${∆}_4$. If the ship’s fuel consumption estimation is accurate, then the total cost of the fleet and fuel consumption cost in the weekly plan is reduced, which shows the effectiveness of the model. More segments will incur more accurate results, but the calculating time will be longer. For example, the calculation time of scenario 2 is greater than scenario 1. The fluctuation range of container transshipment cost is small, which indicates that when the OD flow of container is determined, the container transshipment network is formed and the container transit port is determined, so the reduction of calculation accuracy has little effect on it.

6.3. Analyze the Relationship between Ship Deployment and Sailing Speed

As the unit fuel consumption cost of ships changes, so does the type and quantity of ships deployed on routes 2 and 6, as well as the optimal sailing speed for different legs of the two routes.

Table 7 shows that with the increase of unit fuel consumption cost of ships, the types and number of ships deployed have changed, which means shipping companies need to adjust the deployment of appropriate ship types and determine the appropriate number to reduce the total operation cost of fleet deployment. For the ships deployed in route 2, when the unit fuel consumption cost changes from 100 to 200, one medium-sized ship is replaced by two small-sized ships. It can be seen from Table 7 that the average speed of the ship’s navigation decreases, which is called “adding ship and slow down its speed”. When the unit fuel consumption cost changes from 200 to 300, two small-sized ships are replaced by one medium-size ship. It can be seen from Table 7 that the average speed of the ship’s navigation has increased, which is called “reducing ship and accelerating its speed”.

With the rise of fuel price, the legs of route 2 sailing speed changed alternately because the type and number of ships changed alternately, which can be seen from Table 7. Simultaneously, the last leg had the same sailing speed as leg 4 of route 2 and leg 8 of route 6. The reason is that in the last leg, with the reduction of ship’s payload, the weight of the cargo was negligible relative to the weight of the ship. Therefore, the ship sailed at the optimal economic speed.

6.4. Analysis of the Relationship between Loading Rate and Sailing Speed

The optimal speed and loading rate of route 10 are selected as the research objects in order to study the relationship between loading rate and sailing speed and then get the following trend (Figure 5).

It can be seen from Figure 5 that there is a negative correlation between the optimal speed and the loading rate. With the increase of ship loading rate, the ship’s sailing speed will decrease; on the contrary, the ship’s loading rate will decrease, and the sailing speed will increase. It can be seen from the fuel consumption formula that the fuel consumption will increase with the increase of the ship’s payload or sailing speed. Hence, the increase of ship’s payload and the decrease of sailing speed can improve the energy conservation of ship operation.

7. Conclusions

In this paper, a speed optimization model is proposed for container ship fleet deployment to minimize the total operating costs and fuel consumption costs over the shipping network. Due to container transshipment, containers can be transported from origin port to destination port through many different routes. Although this provides flexibility for the operation of shipping companies, it brings challenges to the deployment of ship type and quantity.

From industry and practical perspectives, this study provides a multi-objective mixed-integer program that potentially can be applied to ship fleet deployment by decision-makers (e.g., the shipowner or the charterer) to reduce fuel consumption costs and operating costs. This paper serves as a valuable guide to determine the fleet deployment measures for liner shipping companies and make trade-offs among economic and environmental considerations.

The critical point of this paper is to balance operating costs and fuel consumption costs, which means that when studying the container ship fleet deployment, it is necessary to calculate fuel consumption according to on-ship payload and ship speed. On this basis, it can determine the ship’s fuel consumption more accurately by optimizing the speed of each leg to confirm the allocation of appropriate ship type and number for the route.

Numerical experiments were performed on 46 ports and 12 routes. Three important conclusions are summarized as follows. First, the number of segments within the range of ship’s payload is increased from 11 to 22, the average total cost is reduced 0.84%, and the average fuel consumption cost is reduced 4.41%. Thus, to obtain the more accurate solution, more piecewise numbers will be needed, but the calculating time will longer. Second, when the unit fuel consumption cost of the ship changes, the shipping company adopts “adding ship and slow down its speed” and “reducing ship and accelerate its speed” strategies to optimize the sailing speed. Third, from the analysis of the relationship between loading rate and sailing speed, it can be seen that there is a negative correlation between sailing speed and loading rate.

The research overview of this article is as follows: first, a mixed-integer nonlinear programming model is transformed into a mixed linear programming model through a series of linearization, taking the speed as the decision variable. Second, the calculation of bunker consumption is more scientific by considering the influence of ship’s payload and speed through nonlinear function. Third, the piecewise linear approximation algorithm is used to optimize the selection of segments on the premise of known accuracy so that the error of the calculation function can find a reasonable balance between the accuracy and the number of segments within a specific range. Under the condition of uncertain demand and speed change of each segment on the route, fleet deployment is the next research direction.

Author Contributions

Conceptualization and methodology, C.-F.G. and Z.-H.H.; writing—original draft preparation, C.-F.G.; writing—review and editing, Z.-H.H.; software, C.-F.G.; supervision, Z.-H.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 71871136.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors gratefully acknowledge the editors and anonymous reviewers for their valuable comments that helped improve this paper. This study is partially supported by the National Nature Science of China (71871136).

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Shipping network of composition by three routes.

Figure 2. The linear outer-approximation technique.

Figure 3. Improved piecewise linear approximation algorithm flow.

Figure 4. The piecewise linear function to approximate.

Figure 5. Relationship between sailing speed and loading rate of each leg of route 10.

Table 1

Description of symbols.

Sets	Description
$R$	$\mathrm{Ship}\mathrm{routes}\mathrm{set}\left(r\in R\right)$
$P$	Ports set $\left(p,o,d\in P\right)$
$V$	$\mathrm{Ship}\mathrm{types}\left(v\in V\right)$
Parameters
$V_r$	$\mathrm{Candidate}\mathrm{ship}\mathrm{types}\mathrm{set}\mathrm{of}\mathrm{route} r$
$I_r$	$\mathrm{Port}\mathrm{indices}\mathrm{for}\mathrm{route}r$
$P_{r,i}$	The ith port of call
$I_{r,p}$	$\mathrm{Port}p$ $\mathrm{indices}\mathrm{of}\mathrm{route}\mathrm{r}\left(p\in P,I_{r,p}\subseteq I_r\right)$
$R_p$	$\mathrm{Routes}\mathrm{set}\mathrm{of}\mathrm{berthing}\mathrm{port}p$
$C_p^{tran}$	Unit container transshipment cost
$C_p^{load}$	$\mathrm{Unit}\mathrm{loading}\cos \mathrm{t}\mathrm{at}\mathrm{port}p$
$C_p^{disc}$	$\mathrm{Unit}\mathrm{discharge}\cos \mathrm{t}\mathrm{at}\mathrm{port}p$
$D_{o,d}$	$\mathrm{Transport}\mathrm{demand}\mathrm{between}\mathrm{port}o$ $\mathrm{and}\mathrm{port} d$
$C_v^{opr}$	$\mathrm{Fixed}\cos \mathrm{ts}\mathrm{of}\mathrm{operating}\mathrm{a}\mathrm{ship}\mathrm{of}\mathrm{type}v$
$C_{p,v}^{ber}$	$\mathrm{Unit}\cos \mathrm{t}\mathrm{of}\mathrm{berthing}\mathrm{operation}\mathrm{at}\mathrm{port}p$ $\mathrm{for}\mathrm{ship}\mathrm{of}\mathrm{type}v$
$Q_v$	$\mathrm{Ship}\mathrm{of}\mathrm{type}v$ capacity
$C_v^{out}$	$\mathrm{Chartering}-\mathrm{out}\mathrm{profit}\mathrm{of}\mathrm{type}v$ ship
$C_{r,v}^{fix}$	Costs related to the voyage
$T_{p,v}^{cont}$	$\mathrm{Per}-\mathrm{container}\mathrm{operating}\mathrm{time}\mathrm{at}\mathrm{port}p$ $\mathrm{for}\mathrm{ship}\mathrm{type}v$
Decision Variables
$n_v^{in}$	$\mathrm{Number}\mathrm{of}\mathrm{ships}\mathrm{of}\mathrm{type}v$$\mathrm{charted}\mathrm{for}\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{type}v$ ships
$n_v^{out}$	$\mathrm{Chartering}-\mathrm{out}\mathrm{number}\mathrm{of}\mathrm{type}v$ ships
$m_{r,v}$	$\mathrm{Number}\mathrm{of}\mathrm{type}v$ $\mathrm{ships}\mathrm{deployed}\mathrm{on}\mathrm{ship}\mathrm{route}r$
$z_{r,v,i}$	$\mathrm{Number}\mathrm{of}\mathrm{containers}\mathrm{handled}\mathrm{for}\mathrm{a}\mathrm{ship}\mathrm{of}\mathrm{type}v$ $\mathrm{at}\mathrm{the}\mathrm{ith}\mathrm{port}\mathrm{of}\mathrm{call}\mathrm{on}\mathrm{ship}\mathrm{route}r$
$x_{r,v}$	$1,\mathrm{if}\mathrm{deployed}\mathrm{ships}\mathrm{of}\mathrm{type}v$$\mathrm{on}\mathrm{route}r$, and 0 otherwise
$z_{o,r,i}^{load}$	$\mathrm{Number}\mathrm{of}\mathrm{containers}\mathrm{loaded}\mathrm{at}\mathrm{the}\mathrm{i}-\mathrm{th}\mathrm{port}\mathrm{of}\mathrm{call}\mathrm{on}\mathrm{route}r$ $\mathrm{originating}\mathrm{from}\mathrm{port}o$
$z_{o,r,i}^{disc}$	$\mathrm{Number}\mathrm{of}\mathrm{containers}\mathrm{discharged}\mathrm{at}\mathrm{the}\mathrm{i}-\mathrm{th}\mathrm{port}\mathrm{of}\mathrm{call}\mathrm{on}\mathrm{route}r$ originating from port o
$z_{o,r,i}^{flow}$	$\mathrm{Number}\mathrm{of}\mathrm{containers}\mathrm{originating}\mathrm{from}\mathrm{port}o$ $\mathrm{sailing}\mathrm{on}\mathrm{the}\mathrm{ith}\mathrm{leg}\mathrm{of}\mathrm{route}r$
$z_p^{tran}$	$\mathrm{Number}\mathrm{of}\mathrm{containers}\mathrm{transshipped}\mathrm{at}\mathrm{a}\mathrm{particular}\mathrm{port}p$
$e_{r,i}$	Ballast water weight required for ship stability sailing on the i-th leg of route r
$s_{ri}$	$\mathrm{Sailing}\mathrm{speed}\mathrm{during}\mathrm{leg}i$ $\mathrm{of}\mathrm{route}r$
$w_{ri}$	$\mathrm{Ship}\mathrm{payload}\mathrm{during}\mathrm{leg}i$ $\mathrm{of}\mathrm{route}r$

Table 2

Basic ship parameters.

	Ship Type
	Small	Medium	Large	Giant
Capacity of different types of ships $Q_v$ (TEU)	1500	3000	5000	10,000
Fixed operating costs of different types of ships $C_v^{opr}$(week)	51,923	76,923	115,384	173,076
Unit cost of berthing operation of different types of ships $C_{p,v}^{ber }$(h)	500	1000	1666	3333
Per-container operating time of different types of ships $T_{p,v}^{cont}$(h)	0.025	0.012	0.011	0.008
Number of ships owned by the shipping lines $N_v^{own}$	20	20	20	20
Chartering-out profit of different types of ships $C_v^{out}$	52,500	77,000	98,000	140,000
Chartering-in cost of different types of ships $C_v^{in}$	66,500	94,500	122,500	175,000

Table 3

Parameters related to the 12 ship routes.

No.	Ports Corresponding Number, Ports Order (Distance between Ports)
1	1Yokohama (15)$\to$2Tokyo (177)$\to$3Nagoya (201)$\to$4Kobe (734)$\to$5Shanghai (745)$\to$6Hong Kong (1568) $\to$ 1Yokohama
2	7Ho Chi Minh (589)$\to$8Laem Chabang (755)$\to$9Singapore (187)$\to$10Port Klang (830) $\to$7Ho Chi Minh
3	11Brisbane (419)$\to$12Sydney (512)$\to$13Melbourne (470)$\to$14Adelaide (1325)$\to$15Fremantle (1733)$\to$16Jakarta (483)$\to$9Singapore (3649)$\to$11Brisbane
4	17Manila (527)$\to$ 18Kaohsiung (164)$\to$19Xiamen (260) $\to$6Hong Kong (15)$\to$20Yantian (19)$\to$21Chiwan (17)$\to$6Hong Kong (620)$\to$17Manila
5	22Dalian (187)$\to$23Xingang (379)$\to$24Qingdao (303)$\to$5Shanghai (93)$\to$25Ningbo (93)$\to$5Shanghai (383)$\to$26Kwangyang (72)$\to$27Busan (487)$\to$22Dalian
6	28Chittagong (872)$\to$29Chennai (573)$\to$30Colombo (306)$\to$31Cochin (723)$\to$32Nhava Sheva (723)$\to$31Cochin (306) $\to$30Colombo (573)$\to$29Chennai (872)$\to$28Chittagong
7	33Sokhna (265)$\to$34Aqabah (554)$\to$35Jeddah (1268)$\to$36Salalah (885)$\to$37Karachi (688)$\to$ 38Jebel Ali (862)$\to$36Salalah (1878)$\to$33Sokhna
8	39Southampton (165)$\to$40Thamesport (386)$\to$41Hamburg (82) $\to$42Bremerhaven (196)$\to$43Rotterdam (42)$\to$44Antwerp (51)$\to$45Zeebrugge (168)$\to$46Le Havre (103)$\to$39Southampton
9	10Port Klang (187)$\to$9Singapore (483)$\to$16Jakarta (1917)$\to$18Kaohsiung (904)$\to$27Busan (904) $\to$18Kaohsiung (342)$\to$6Hong Kong (17)$\to$21Chiwan (1597)$\to$10Port Klang
10	39Southampton (3162)$\to$33Sokhna (1878)$\to$36Salalah (1643)$\to$30Colombo (1560)$\to$9Singapore (1415)$\to$6Hong Kong (260) $\to$19Xiamen (486)$\to$5Shanghai (448)$\to$27Busan (487)$\to$22Dalian (187)$\to$23Xingang (379)$\to$24Qingdao (303) $\to$5Shanghai (745)$\to$6Hong Kong (1415)$\to$9Singapore (1560)$\to$30Colombo (1643)$\to$36Salalah (5029)$\to$39Southampton
11	11Brisbane (419)$\to$12Sydney (512)$\to$13Melbourne (470)$\to$14Adelaide (1325)$\to$15Fremantle (3148)$\to$30Colombo (1643)$\to$36Salalah (5244)$\to$43Rotterdam (5244)$\to$36Salalah (1643)$\to$30Colombo (5191)$\to$11Brisbane
12	20Yantian (9956)$\to$41Hamburg (3621)$\to$33Sokhna (620)$\to$35Jeddah (4156)$\to$10Port Klang (187)$\to$9Singapore (1309)$\to$17Manila (629)$\to$20Yantian

Table 4

The ship-fixed operating costs of 12 ship routes.

${\mathit{c}}_{\mathit{r},\mathit{v}}^{\mathit{f}}$	Ship Type					Ship Type
No.	Small	Medium	Large	Giant	No.	Small	Medium	Large	Giant
1	226,198	280,542	−	404,000	7	404,711	499,139	−	710,000
2	154,791	191,900	−	276,100	8	−	129,712	149,622	199,300
3	533,980	656,891	−	929,100	9	−	501,088	551,946	715,100
4	−	155,123	176,013	232,200	10	−	−	1,883,007	2,430,000
5	148,807	187,600	−	279,700	11	1,504,231	1,843,178	−	2,583,900
6	322,916	400,072	−	574,800	12	1,235,313	1,512,755	−	2,117,800

Table 5

Four scenarios composed of the different approximation error ${\epsilon }_1$ and ${\epsilon }_2$ combinations.

${\mathit{\epsilon }}_1$	Segments	${\mathit{\epsilon }}_2$	Segments	Scenario 1	Scenario 2	Scenario 3	Scenario 4
1.6 × 10−3	68	7.6 × 10−3	11	68	68	135	135
4.2 × 10−4	135	9.6 × 10−4	22	11	22	11	22

Table 6

Various cost comparisons of fleet deployment under different scenarios.

$\mathit{\beta }$	Scenario 1					Scenario 2					${\mathbf{\Delta }}_1\mathbf{\times }{\mathbf{10}}^{\mathbf{-}\mathbf{2}}$	${\mathbf{\Delta }}_{\mathbf{2}}\mathbf{\times }{\mathbf{10}}^{\mathbf{-}\mathbf{2}}$
	$\mathit{f}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{7}}\right)$	${\mathit{f}}^{\mathit{F}\mathit{C}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	${\mathit{f}}^{\mathit{C}\mathit{T}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	Gap	Time	$\mathit{f}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{7}}\right)$	${\mathit{f}}^{\mathit{F}\mathit{C}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	${\mathit{f}}^{\mathit{C}\mathit{T}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	Gap	Time
100	2.26345	1.456898	5.4227	4.5785%	157s	2.25179	1.305053	5.4227	4.1103%	1320s	−0.5	−1.0
200	2.38397	2.291148	5.4082	4.7686%	131s	2.37363	2.203396	5.4082	4.7209%	1693s	−0.4	−3.8
300	2.50588	2.978536	5.397	4.3182%	116s	2.4736	2.757697	5.397	4.7709%	1455s	−1.2	−7.4
400	2.61124	3.153422	5.3957	4.9317%	132s	2.55396	2.924008	5.3957	4.3785%	1791s	−2.1	−7.2
500	2.65745	3.470209	5.3562	4.9863%	170s	2.67637	3.780855	5.3562	4.7166%	1422s	−0.7	8.9
600	2.71938	3.725877	5.4509	4.8223%	176s	2.6914	3.686778	5.4509	4.7034%	1534s	−1.0	−1.0
$\mathit{\beta }$	Scenario 3					Scenario 4					${\mathbf{\Delta }}_{\mathbf{3}}\mathbf{\times }{\mathbf{10}}^{\mathbf{-}\mathbf{2}}$	${\mathbf{\Delta }}_{\mathbf{4}}\mathbf{\times }{\mathbf{10}}^{\mathbf{-}\mathbf{2}}$
	$\mathit{f}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{7}}\right)$	${\mathit{f}}^{\mathit{F}\mathit{C}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	${\mathit{f}}^{\mathit{C}\mathit{T}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	Gap	Time	$\mathit{f}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{7}}\right)$	${\mathit{f}}^{\mathit{F}\mathit{C}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	${\mathit{f}}^{\mathit{C}\mathit{T}}\left(\mathbf{\times }{\mathbf{10}}^{\mathbf{6}}\right)$	Gap	Time
100	2.25775	1.53476	5.4208	4.2802%	219s	2.2547	1.285647	5.3908	4.4730%	1663s	−0.1	−16
200	2.39083	2.230419	5.3724	4.8466%	182s	2.37366	2.208508	5.3239	4.6335%	1807s	−0.7	−0.9
300	2.50871	3.146133	5.408999	4.7001%	160s	2.47303	2.589513	5.352	4.9580%	2006s	−1.4	−18
400	2.59337	3.204796	5.3449	4.7656%	146s	2.56488	3.067922	5.3726	4.8849%	1514s	−1.1	−4.2
500	2.66613	3.596815	5.3561	4.9071%	159s	2.64155	3.536123	5.3642	4.9157%	1514s	−0.9	−1.6
600	2.6916	3.704831	5.346299	4.8691%	167s	2.69072	3.677827	5.3619	4.8956%	1606s	−0.0	−0.7

Notes: ${\Delta }_1=\left(f_{Scenario2}-f_{Scenario1}\right)/f_{Scenario1}$, ${\Delta }_2=\left(f_{Scenario2}^{FC}-f_{Scenario1}^{FC}\right)/f_{Scenario1}^{FC},{\Delta }_3=\left(f_{Scenario4}-f_{Scenario3}\right)/f_{Scenario3},{\Delta }_4=\left(f_{Scenario4}^{FC}-f_{Scenario3}^{FC}\right)/f_{Scenario3}^{FC}$, $Gap=\left|BP-BF\right|/\left|BP\right|$.

Table 7

Ship deployment and sailing speed changes of Route 2 and Route 6.

$\mathit{\beta }$	Ship Deployment		Each Legs Speed of Route 2				Ship Deployment		Each Legs Speed of Route 6
	Type	Num	1	2	3	4	Type	Num	1	2	3	4	5	6	7	8
100	3000	1	25.98	25.62	25.98	23.83	1500	3	16.26	16.00	16.00	16.26	16.21	16.26	16.26	24.45
200	1500	2	10.36	10.42	10.36	23.83	3000	2	21.06	20.62	20.62	20.89	20.62	20.62	20.62	24.45
300	3000	1	25.51	25.98	25.98	23.83	3000	2	20.62	21.06	20.62	20.89	20.62	21.06	20.62	24.45
400	1500	2	10.41	10.36	10.47	23.83	3000	2	21.06	20.62	20.62	20.62	20.71	21.06	20.62	24.45
500	3000	1	25.51	25.98	25.98	23.83	1500	3	16.26	16.07	16.00	16.26	16.26	16.00	16.26	24.45
600	1500	2	10.58	10.61	10.58	23.83	1500	3	11.24	11.17	11.17	11.17	11.17	11.17	11.17	24.45