Speed Optimization for Container Ship Fleet Deployment Considering Fuel Consumption

: 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 ﬂeet 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 veriﬁes the validity and applicability of the established model.


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.

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 mixedinteger 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.

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.

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

Sets Description
Candidate ship types set of route r I r Port indices for route r P r,i The ith port of call I r,p Port p indices of route r p ∈ P, I r,p ⊆ I r R p Routes set of berthing port p C tran p Unit container transshipment cost C load p Unit loading cos t at port p C disc p Unit discharge cos t at port p D o,d Transport demand between port o and port d C Ballast water weight required for ship stability sailing on the i-th leg of route r s ri Sailing speed during leg i of route r w ri Ship payload during leg i of route r

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 (s, w), s is the sailing speed, and w is the ship payload of each segment. A realistic approximation of fuel consumption function is f (s, w) = c r 1 q + s c r 2 (w + A) c r 3 with c r 1 , c r 2 , c r 3 and q constants such as c r 1 > 0, c r 2 ≥ 3, c r 3 ≥ 0 and q ≥ 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 (s ri , w ri ) = c r 1 (s ri ) c r 2 (w ri ) c r 3 . 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 . β is the fuel price, 24 is the hour of the day, and the ship fuel consumption cost expressed as: 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.

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.
f CB = ∑ r,i∈I r ,v∈V r C ber 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 M 1 and M 2 to 9999; Constraints (12) and (13) define the relationship between variables, M 3 and 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.

Linearization of the Model
The above model for fleet deployment with sailing speed and routing optimization has a nonlinear objective (1) with nonlinear terms (s ri ) c r 2 −1 (w ri ) c r 3 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).

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 s min ri , s max ri , this speed interval depends on mechanical power engine size and economic factors. Therefore, u ri lower bound is u min ri = 1/s max ri and upper bound is u max ri = 1/s min ri . Then, the Constraint (11) becomes a linear constraint and can be rewritten by Constraint (27)

Linearization of the Objective Function of Fuel Consumption Cost
The objective function of fuel consumption cost can be expressed as follows: The fuel consumption cost objective Function (28) can be further linearization by introducing the auxiliary variables: B ri , ∇ ri, and Λ ri , transformed to the following objective Function (29) and Constraints (30)-(32). β· ∑ r∈R,i∈I r B ri

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 x L y + y L x − x L y L , x U y + y U x − x U y U in the interval x L , x U × y L , y U . Due to the fact that ∇ ri and Λ ri have the lower limit and upper limit, Wang et al. [26] showed that ∇ ri ·Λ ri term is bilinear. The lower bound for Therefore, Constraint (30) can be relaxed by adding linear Constraint (34) and Constraint (35):

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.

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 u min ri , u max ri is subdivided into several segments. For each segment, the function (u ri ) 1−c r 2 is approximated by a tangent line.
Proof. The first-order derivatives of (u ri ) 1−c r 2 are denoted by (u ri ) 1−c r 2 and can be calculated as: The second-order derivative of (u ri ) 1−c r 2 with respect to u ri is denoted by (u ri ) 1−c r 2 and can be calculated as: According to Wang and Meng [22], the coefficient c r 2 is greater than 2, hence the second-order derivative (u ri ) 1−c r 2 > 0. When on the interval u min ri , u max ri , (u ri ) 1−c r 2 is convex.
In order to approximate (u ri ) 1−c r 2 , Figure 2 shows the tangent lines. To control the approximation error ε 1 , some tangent lines are generated a priori. Then, whether the calculation gap is larger than ε 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 Ω ri . At tangent point u k ri , the tangent line of the function (u ri ) 1−c r 2 can be derived as: where Ω(ε 1 ) is the tangent points set of function (u ri ) 1−c r 2 , when given ε 1 , the procedure to obtain set Ω(ε 1 ) is shown in Algorithm 1. Then, the nonlinear Constraint (31) can be expressed as: Algorithm 1. The procedure of linear outer-approximation algorithm to obtain tangent point sets.

Input:
Convex function (u ri ) 1−c r 2 , the tangent point set Ω ri . The lower limit and upper limit of u ri is u min ri , u max ri , the approximation relative error ε 1 Output: the tangent point set Ω ri (ε 1 ) Step 1 Calculate the midpoint u mid ri of the interval u min ri , u max ri ; u mid ri can be carried out by the bisection search method Step 2 At tangent point u mid ri , the tangent line is defined as u mid Step 3 Calculate the relative approximation error of points u min ri and u max ri according to Step 4 If ε l > ε or ε r > ε, then branch the feasible range of u ri is divided into two ranges: u min ri , u mid ri and u mid ri , u max ri , and the tangent point set Ω ri (ε 1 ) = Ω ri (ε 1 ) ∪ u mid ri In one branch u ri ∈ u min ri , u mid ri , repeat the above step 1 to e step 4 until the stop criterion is reached In the other branch u ri ∈ u mid ri , u max ri , repeat the above step 1 to e step 4 until the stop criterion is reached Stop criterion check: if ε l < ε and ε r < ε, stop and output the current solution. Otherwise, go to Step4.

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 w min ri , w max ri is subdivided into several segments; for each segment, the function (w ri ) c r 3 is approximated by the linear function l K (w ri ).
Proposition 2. f (w ri ) = (w ri ) c r 3 is a continuous function on the interval w min ri , w max ri , then the piecewise linear function l K (w ri ) is uniformly approximated to f (w ri ).
Proof. Since f (w ri ) is continuous on w min ri , w max ri , then f (w ri ) is uniformly continuous on w min ri , w max ri . That is, for any given ε > 0, there exists δ(ε) > 0; for any w ri , w ri ∈ w min ri , w max ri , as long as w ri − w ri < δ(ε), there is f w ri − f w ri < ε. The interval w min ri , w max ri is divided into K segments; a j , (j = 1, · · · , K, K + 1) is expressed as the breakpoints of the piecewise linear function. Command w min ri = a 1 < a 2 < a 3 < · · · < a K < a K+1 = w max ri . When K > 1/δ(ε), for ∀w 1 ri , w 2 ri ∈ a j , a j+1 , there is f w 1 ri − f w 2 ri < ε. min f a j , f a j+1 ≤ l K,j+1 (w ri ) ≤ max f a j , f a j+1 ,l K,j+1 (w ri ) indicates the function section for the interval a j , a j+1 , that is, the j + 1 segment linear function. So , there is a piecewise linear function l K (w ri ) that is uniformly approximated to the continuous function f (w ri ). 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 w min ri , w max ri , the approximation lines can be generated as follows: Function Generate approximate straight lines l K (w ri ) two endpoints of the curve f (w ri ) = (w ri ) c r 3 in interval w min ri , w max ri are A = w min ri , w min ri c r 3 and B = w max ri , w max ri c r 3 , the equation of a straight line is found through these two points.
Assuming the point C = (w ri , l K (w ri )) is on a straight line, hence the slope of the line AC and the line AB are the same, expressed as: By transforming the formula, we can get an approximately straight line: The approximation error between straight line and function is expressed as d = | f (w ri ) − l K (w ri )|. Then, the approximation error is calculated between the approximate straight line l K (w ri ) and function (w ri ) c r 3 at the midpoint w mid ri . The concavity of the function (w ri ) c r 3 is judged; the following lemma proves the function property is nonconvex: Proof. The second-order derivative of (w ri ) c r 3 with respect to w ri denoted by (w ri ) c r 3 can be calculated by constraint According to the N.Psaraftis and A.Kontovas [24], the coefficient c r 3 is between 0 and 1; hence, the second-order derivative (w ri ) c r 3 < 0. So, the function (w ri ) c r 3 is concave when w ri takes value in interval w min ri , w max ri . Therefore, the function (w ri ) c r 3 is concave, then the approximation line is moved to l K (w ri ) + d/2, shown in Figure 4. The maximum error is found between the function and the new approximation line on the interval w min ri , w max ri , expressed as d max = | f (w ri ) − (l K (w ri ) + d/2)|. 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 ε 2 : if the maximum error d max > ε 2 , the interval of the approximation is further subdivided and the breakpoints of the piecewise linear function set W ri = W ri ∪ w mid ri , then we branch the feasible range of w ri into two ranges: w min ri , w mid ri and w mid ri , w max ri . If the maximum error d max < ε 2 , the continuous subdivision of the segment is stopped. w ri is divided into K segments in the interval w min ri , w max ri , which is determined by improved piecewise linear approximation algorithm, the breakpoints of the piecewise linear function set expressed as W ri (K) = w r,i,k , k ∈ M , where M = {1, 2, · · · |K|, |K + 1|}; by introducing the auxiliary variable η rik and π rik , then, the nonlinearity Constraints (32) can be rewritten as: ∑ k∈M−|K+1| π r,i,k = 1, ∀r ∈ R, i ∈ I r (41) η r,i,1 ≤ π r,i,1 , ∀r ∈ R, i ∈ I r (42) η r,i,1 ≤ π r,i,1 , ∀r ∈ R, i ∈ I r (43) η r,i,k ≤ π r,i,k−1 + π r,i,k , ∀r ∈ R, i ∈ I r k = 2, · · · , K (44)

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.

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[0, 30]. The loading cost C load p and discharge cost C disc p are 150; the container transshipment cost is 200. For the value of c r 1 , c r 2 , and c r 3 from literature Wang et al. [26], the coefficient c r 1 is 0.0006, c r 2 is 2.5, and c r 3 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.    [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 × 10 −3 . Yang and Xing [28] set the number of line fragments for linear outerapproximation 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 ε 1 and piecewise linear function approximating error ε 2 , as shown in Table 5.

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.

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.

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.

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 mixedinteger 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.