Simultaneous Optimization of the Liner Shipping Route and Ship Schedule Designs with Time Windows

(is paper proposes to address the problem of the simultaneous optimization of the liner shipping route and ship schedule designs by incorporating port time windows. Amathematical programmingmodel was developed tominimize the carrier’s total operating cost by simultaneously optimizing the port call sequence, ship arrival time per port of call, and sailing speed per shipping leg under port time window constraints. In view of its structure, the nonlinear nonconvex optimization model is further transformed into a mixed-integer linear programming model that can be efficiently solved by extant solvers to provide a global optimal solution. (e results of the numerical experiments performed using a real-world case study indicated that the proposed model performs significantly better than the models that handle the design problems separately. (e results also showed that different time windows will affect the optimal port call sequence. Moreover, port time windows, bunker price, and port efficiency all affect the total operating cost of the designed shipping route.


Introduction
Container lines transport containers following fixed sequences of port calls with a regular service frequency [1,2]. A port, the core of the shipping logistics chain, needs to service a large number of ships on different routes according to predetermined schedules. In recent years, with the fast growth of container trade, ports are becoming more and more congested. In April 2017, Shanghai port, the world's largest container port, was not immune from port congestion, which surprised the industry. Meanwhile, slow port expansion projects in some nondeveloped countries and regions have exacerbated this phenomenon. is shows that the daily service capacity of a port is limited, and, therefore, it cannot guarantee that a ship can be serviced as soon as it arrives at the port. As a result, the availability of ports, i.e., port time windows, is the first factor to be considered when designing a ship schedule [3].
However, in practice, it is not feasible to design a ship schedule under port time windows constraints for a given route without considering the optimization of the port rotation. We take the S2 service route operated by SITC Container Lines as an example to illustrate. e route, consisting of four ports, provides services following the port call sequence of Taicang-Shanghai-Osaka-Kobe-Taicang. Without the constraints of time windows, only one container ship is needed to maintain weekly service frequency on the route. Now, we assume that the time windows of each port call are as follows: Taicang (Tuesday), Shanghai (Monday), Osaka (Friday), and Kobe ( ursday). As shown in Figure 1(a), owing to the constraints of time windows, a ship needs to wait in the anchorage for nearly 6 days before arriving at Shanghai as well as at Kobe following the original port call sequence on the route. (It takes only a few hours for a ship from Taicang to Shanghai as well as from Osaka to Kobe.) erefore, it would take 3 weeks for a ship to perform a round-trip journey, which means that three ships need to be deployed on the route. As a result, the weekly operating cost of the route will increase significantly. In addition, the longer transit time leads to delay in delivery, which also results in loss of sales, goodwill, and service satisfaction to the carrier. However, if we change the port call sequence to Shanghai-Taicang-Kobe-Osaka-Shanghai, as shown in Figure 1(b), it only takes 1 week for a ship to perform a round-trip journey. From the above comparison, it is easy to see that the adjusted port call sequence can greatly reduce the operating cost compared with the original one, subject to the constraints of time windows. e reason why we can make such an adjustment is the operational flexibility of the route. Taicang and Shanghai are both ports in China. For the carriers operating foreign trade, there is no cargo flow carried by the ships between the two adjacent O-D port pairs in the same country except for special reasons (empty container repositioning, transshipment operations due to skipping port calls, etc.). Moreover, the distance between the two ports is very short. In fact, such similar port pairs or groups are fairly common, such as the Kansai port group, which includes the Osaka and Kobe ports; the Kanto port group, which includes the Tokyo, Yokohama, and Kawasaki ports; the Manila north and south ports; and the Laem Chabang and Bangkok ports.
Hence, liner shipping route design and ship schedule design with time windows should be carried out at the same time rather than one after another. Otherwise, the designed shipping route and schedule will not be feasible in reality.
is study provides a decision support system for container lines to improve their services and reduce their operating costs.

Literature Review.
is study is related to two research topics: liner shipping network design and ship scheduling. e former has been a hot topic in the field of liner shipping research over recent decades. e work of Rana and Vickson [4] was the early contribution in this field. ey formulated a mathematical model to determine the port rotation, the number of trips each ship makes in a planning horizon, and the amount of transported cargo for multiple ships by employing a Lagrangian relaxation method with the objective of maximizing the total profit. Shintani et al. [5] proposed a two-stage programming model to design a single shipping route by explicitly considering empty container repositioning. A genetic algorithm-based heuristic was employed to solve the problem. e work by Agarwal and Ergun [6] was the first academic study to consider the weekly frequency constraint in the design of liner transport service networks composed of two or more routes. However, the sailing speed in their model was predetermined, not variable. Alvarez [7] jointly optimized the optimal routing and fleet deployment by using a combined tabu search and a column generation-based heuristic method. Wang and Meng [8] proposed a novel aspect in the liner shipping network design problem by reversing the port rotation direction. ey demonstrated that the operating cost of a container liner shipping network could be significantly reduced by reversing the port rotation direction of liner service routes. Song and Dong [9] proposed a shipping route design problem that considers ship deployment and empty container repositioning with the objective of minimizing the total cost. A three-stage solution method was proposed to deal with the problem. However, they failed to consider the speed variation in each leg. Wang and Meng [10] considered the delivery deadlines in the liner shipping network design problem. e problem was proved to be NP-hard and solved by a column generation-based heuristic method. Karsten et al. [11] examined the liner shipping network design problem that considers the coordination between vessels and cargo transit time restrictions. Wang et al. [12] designed a single intercontinental service by simultaneously optimizing the sailing speed, the port rotation of selected ports in a given set of ports, and the shipping demand for each type. A mixed-integer linear programming model was developed and solved by using exact algorithms. Wang et al. [13] proposed a port call adjustment problem to determine the new additional port calls to be inserted and the existing port calls to be removed for the routes in a shipping network. A two-phase solution method was proposed using two heuristic methods.
A considerable number of papers on vessel scheduling have been published in the past. Meng and Wang [14] jointly determined the sailing speed, service frequency, and fleet deployment for an operating strategy problem. An exact ε-optimal algorithm was proposed to solve the problem. Wang and Meng [15] investigated the optimal sailing speed on each shipping leg and the optimal fleet size in a liner shipping network. ey reformulated the bunker consumption function from historical operating data, which was applied to a real-world case study. Wang et al. [16] simultaneously optimized the schedules of shipping routes and cargo allocation with the objective of reducing the demurrage cost. Cheaitou and Cariou [17] tried to optimize the sailing speed under a semielastic demand. eir results showed that the slow steaming strategy may not be optimal. Hvattum et al. [18] developed an exact algorithm to optimize the sailing speed for a fixed sequence of port calls with hard time windows. ey proved that the optimal speeds can be obtained in quadratic time. Wang et al. [19] studied the liner ship schedule design problem with hard time windows. Wang et al. [3] studied the same problem with consideration of the availability of each berth with the objective of minimizing the bunker cost, operating cost, and inventory cost for a single shipping route. Alharbi et al. [20] examined the ship schedule design problem with time windows for a liner shipping network. In their model, the vessel could utilize a premium berth with high penalty when violating the berth time window. Aydin et al. [21] focused on the sailing speed optimization problem with time windows on a single voyage that considers the waiting cost of early arrivals. Dulebenets [22] investigated the vessel scheduling problem with a heterogeneous fleet for a single shipping route. Zhen et al. [23] proposed an integrated planning model to simultaneously determine the fleet size, ship schedule, sailing speed, and cargo allocation for a given shipping network. Kim et al. [24] designed a simple single route by considering the variable sailing speed on each leg and the fleet size to maximize a carrier's profit. Koza et al. [25] proposed a mixed-integer programming model for the integrated liner shipping network design and scheduling problem that incorporates the transshipment times, together with a solution method based on column generation. Jiang et al. [26] developed a mixed-integer linear programming model for the near-sea liner shipping schedule design problem considering big customers' time preferences.
From the abovementioned literature, we can see that the two key and interrelated problems, namely, liner shipping route and ship schedule, are usually treated separately by most existing studies. Although a few articles combined the two, they did not consider the port time window constraints and the variable sailing speed per shipping leg. is study tries to fill the gap. As mentioned in Section 1 of this study, it is highly desirable to simultaneously optimize port call sequence and ship schedule subject to the constraints of time windows.

Objectives and Contributions.
is study focuses on jointly optimizing the port rotation, sailing speed on each leg, and ship arrival time at each port of call for a given set of ports with time windows while minimizing the sum of the vessel operating cost, fuel cost, and waiting cost. is problem has practical significance because it can help container lines improve service levels and reduce total operating costs. e contribution of this study is threefold. First, it simultaneously determines the optimal port rotation, the optimal sailing speed per shipping leg, and the optimal arrival time per port of call under port time window constraints. Second, a tailored mixed-integer programming model was developed that considers the special structure of the problem. Finally, some useful management insights obtained from case studies can provide support for container lines.

Problem Description and Mathematical Formulation
We consider a group of ports, denoted by a set I. e set I has N number of physical ports indexed by i and j. e voyage from port i to port j is called the shipping leg (i, j).
We define the set of all shipping legs as A � (i, j)|i, j ∈ I . Since the shipping route is a closed circle, we can arbitrarily choose one port as the first port of call. In this study, we choose physical port 1 as the first port of call. e details on problem description and key constraints are provided in the following sections.

Port Rotation.
e port rotation of a shipping route forms a loop in practice. We use the binary variables x ij ∈ 0, 1 { } to indicate whether the shipping leg (i, j) is used. In this study, we assume that each physical port can be visited only once during a round-trip journey; then, the constraint j∈I,j≠i x ij � j∈I,j≠i x ji � 1 should hold for each port i.

Bunker Consumption.
A ship's unit bunker consumption largely depends on its sailing speed. We denote by the variable v ij the sailing speed of a ship on leg (i, j), and we define f ij (v ij ) as the bunker consumption in tons per nautical mile at the sailing speed v ij on leg (i, j). Based on the results of existing studies [15,27], the fuel consumption is a power function of the sailing speed: where a ij and b ij are two coefficients estimated from the practical data. Let β be the unit bunker price (in US dollars per ton), and let L ij be the length of shipping leg (i, j); then, the weekly fuel cost of the designed shipping route can be calculated by β i∈I j∈I x ij L ij a ij (v ij ) b ij .

Port Time Windows.
Since the service provided by the ports has time windows, we set Ω i to represent the time windows at port i. We define the time 00 : 00 of a certain Sunday as time 0. Let t a i be the ship arrival time at port i, and let Ω i ⊆[0, 168] be the time windows at port i in the first cycle. For example, Ω i � [24, 48]⋃ [96, 168] means that port i is available on Monday, ursday, Friday, and Saturday, and a ship can call the port for service in these periods. Otherwise, the ship must wait in the anchorage until the time windows open again. Owing to the characteristics of weekly service frequency of liner shipping, ships call port i at the same time in different weeks. Hence, an arrival time is feasible if and only if (t a i mod168) ∈ Ω i . For example, a ship can call port i for service at t a i � 298 because (298mod168 � 130) ∈ Ω i . It should be noted that, as a ship has a minimum sailing speed, it needs to wait in the anchorage when arriving at the port outside the time windows. On the one hand, it is well known that, when a ship is waiting in the anchorage, although its main engine is not running, its auxiliary engine still does. us, there is still fuel consumption on the ship. On the other hand, the inventory cost of containers also rises with the increase in waiting time. Hence, the weekly waiting cost of the designed route can be calculated by i∈I j∈I t w ij c w j , where the variable t w ij is the waiting time for the ship in the anchorage at port j after using shipping leg (i, j), and c w j is the fixed waiting cost per hour in the anchorage at port j.

Fixed Operating Cost.
A fleet of homogeneous ships is deployed on the designed shipping route to achieve a weekly service frequency for each port of call. Let the integer variable m represent the number of ships deployed on the route. us, the weekly operating cost of the designed route can be calculated by C op m, where C op is the fixed weekly operating cost in US dollars per week per ship.

Mathematical Model.
e problem of the simultaneous optimization of the liner shipping route and ship schedule designs with time windows aims to determine the optimal port call sequence, the optimal sailing of each shipping leg, and the optimal ship arrival time that satisfies the time window constraints at each port to minimize the total cost including the vessel operating cost, fuel cost, and waiting cost. Before presenting the model, we list the variables and parameters below.
Decision variables: x ij : binary variables indicating whether shipping leg (i, j) is used on the designed shipping route subject to e objective function (1) minimizes the total weekly operating cost, which amounts to the sum of the vessel operating cost, fuel cost, and waiting cost. Constraint (2) ensures that each port can be visited only once. Constraint (3) is the Miller-Tucker-Zemlin (MTZ) subtour elimination constraints. Constraints (4) and (5) jointly guarantee that if leg (i, j) is used, T ij � (L ij /v ij ) or T ij � 0 otherwise. Constraint (6) defines the lower and upper bounds of the sailing speed on each shipping leg. Constraints (7) and (8) jointly impose the port time window restrictions at each port. Constraints (9) and (10) establish the relationship between the port call arrival times: if leg (i, j) is used, the ship arrival time at port j is equal to the sum of the ship arrival time at port i, the fixed dwell time at port i, the sailing time on shipping leg (i, j), and the waiting time in the anchorage at port j. Constraints (11) and (12) define the time when the ship returns to the first port of call after a round-trip journey. Constraints (13) and (14) indicate that if leg (i, j) is not used, the waiting time for the ship in the anchorage at port j should be 0. Constraint (15) enforces a weekly service frequency on the designed shipping route. Constraint (16) indicates that the number of ships deployed on the shipping route is a positive integer. Constraints (17) and (18) indicate that all container shipping demands should be satisfied while ensuring the container flow balance. Constraint (19) indicates that the total number of containers loaded on each shipping leg should not exceed the capacity of one container ship.

Approximate MILP Model
To take advantage of state-of-the-art MILP solvers such as CPLEX to solve the developed model [M1], we need to deal with three difficulties: (i) the bunker consumption, which contains two variables and power functions in the objective function, is nonlinear; (ii) constraint (5) contains the reciprocal of the sailing speed; and (iii) both the "mod" operation and the set Ω i lead to a nonconvex domain. ese three difficulties make the proposed model challenging to solve using the MILP solvers. Next, we deal with these difficulties and transform the model [M1] into an approximate mixed-integer linear programming model.

Linearization of the Bunker Consumptions.
First, we define a new variable u ij to be the reciprocal of the sailing speed v ij en, constraints (5) and (6) can be transformed into linear constraints associated with the variable u ij : We introduce an auxiliary variable B ij to replace the nonlinear combinational items x ij and v ij in the objective function. We define By introducing a very large number M 4 , we can replace the above if-then complementary conditions equivalently with the big-M reformulations of the following constraints: By introducing the auxiliary variables g ij , we define g ij � a ij (u ij ) − b ij . Objective function (1) can be transformed by constructing its epigraph form as follows: subject to and to constraints (23) and (24).
where M is a given integer number. As long as M is large enough, the error of the approximate model can be controlled within an acceptable range. We use u m ij to represent any value of the M equal segments. From the property of the first-order condition of the convex function, we can obtain where h ij ′ (u m ij ) represents the tangent line of h ij (u ij ) at the point u m ij . us, the function h ij (u ij ) � a ij (u ij ) − b ij can be converted into the following piecewise function: By combining constraints (26) and (29) e penalty cost is high enough to limit a ship to call the port only within the time windows. For example, as shown in Figure 2, ships can only call port i on Tuesdays and ursdays. Next, we linearize the penalty function p i (t a i ). Since p i (t a i ) has n i turning points, let b ik be the turning points and let K i : � 1, 2, . . . , n i be the set of all turning points at port i. For example, as shown in Figure 2, there are 10 turning points: 9 � 120, and b i,10 � 168. By introducing the continuous auxiliary variables w ik and the binary auxiliary variables z ik , we can transform the piecewise function into the following linear constraints: where w ik can realize the function linearization and z ik can limit the value range of w ik to limit the interval of t a i . us, the total penalty costs of time windows at all ports are as follows: By introducing the auxiliary variables θ i ∈ Ζ + , we have It is noticeable that constraints (42)-(44) are equivalent to constraint (7).
[M2]min C op m + β i∈I j∈I

Case Study
is section presents the results of the numerical experiments we conducted on the PBT1 service operated by SITC Container Lines to assess the applicability of the proposed models and algorithms. We first tested the validity of the solution approach to determine how it can significantly reduce the total cost. en, we conducted a sensitivity analysis to test some key input parameters in the solution approach.
e operating cost C op � US, the maximum speed V max ij � 25knots, the minimum speed V min ij � 5knots, the bunker price β � US, the maximum number of ships deployed m max � 5, and the fixed waiting cost in the anchorage at each port c w j �US$100/h. For ease of calculation, we take the tightest upper bounds; then, the very large numbers (M 1 , M 2 , M 3 , and M 4 ) can be set as max(L ij )/V min ij , 168m max , 168, and max(L ij ) · a ij (V max ij ) b ij , respectively. To avoid violating the time window constraints, the very large number M 5 was set to 300,000. e number of segments for the interval [U min ij and U max ij ] was M � 200. e bunker consumption function of each leg was set to 0.001(v) 2 . e port-to-port distance, the dwell time at each port, and the time windows are shown in Tables 1-3, respectively. CPLEX 12.8.0, programmed by the MATLAB toolbox YALMIP running on a 2.5 GHz Intel Core Duo laptop with 16 GB of RAM, is called to solve the MILP model. e MILP model inputting the above parameters has 3,473 continuous variables, 207 integers (198 binary) variables, and 21,466 constraints, and the optimal solution with 0 gap among objective bounds can be found in about 1 hour.

Comparison with the Original Port Rotation.
In this section, we present the results of the validation of the proposed model and the solution approach. e results obtained using the parameters in Section 4.1 are shown in Tables 4 and 5. It is noticeable that the arrival times do not violate port time window constraints. e solution approach is feasible and optimal. We further compared our model [M2] with the model of [3], which optimized sailing speeds without optimizing port rotation, as shown in Table 6. We can see that it takes nearly 400 h for a container ship to wait in the anchorage in a round-trip journey using the model in [3]. Hence, five container ships need to be deployed on the route to meet the weekly service frequency, while only two container ships need to be deployed following the proposed model [M2]. e total cost is also reduced by 45% compared with that in [3]. erefore, the proposed model and the solution approach can greatly reduce the total cost.

Effects of Port Time Windows.
In this section, we discuss the impact of port time windows on the solution. First, we examine the effect of relaxing the time windows on the results. We set the time windows to three scenarios, where scenario 1 is the setting in section 4.1. As shown in Table 7, in scenario 2, we give an extra day for the time window of each port while keeping the time windows of scenario 1; we did the same in scenario 3 on the basis of scenario 2. e results are shown in Table 8. It is noticeable that the total cost decreases and that the port rotation is different.
is is because, with the relaxation of the time windows, there will be a better port rotation to reduce the total cost. en, we examine the effect of different time windows on the results. Scenarios 4 and 5 are set to different time windows from that of scenario 1, and each of them has only 1 day available in a week on each port, as shown in Table 9. From the results in Table 10, we can see that, under different time window restrictions, both the optimal port rotation and the total cost are different. e results also confirm the applicability and flexibility of our proposed models and algorithms, which always seek the optimal port rotation subject to the constraints of the time windows.

Effects of Bunker Prices.
e bunker price fluctuates, and the fuel cost constitutes the bulk of the total operating expenses. Hence, in this section, we examine the impact of bunker price on the solution. We set the bunker price from US$200 to US$600, and the other parameters are the same as those in Section 4.1. e results are shown in Table 11. We can see that, when the bunker price rises, the port rotation remains the same, while the total cost increases significantly. Moreover, with the increase in the bunker price, one more container ship is deployed to maintain the weekly service frequency. e reason is that when more ships are deployed, the round-trip time is longer, and hence ships can sail at a lower speed to reduce bunker consumption.

Effects of Port Efficiency.
In this section, we discuss the effect of port efficiency on the solution. We varied the dwell time at each port by −50%, −25%, +25%, +50%, and +100% of the benchmark value, but the other parameters remained the same. e results are shown in Table 12. We can see that, with the decrease in port efficiency, the total cost increased significantly. e port rotation changed only when the dwell time was +100% of the benchmark value, and in other cases, it remained the same. is indicates that the port rotation under the benchmark value cannot meet the limit of the time windows. Hence, the optimal solution approach gives a new port call order.

Conclusions and Future Work
In this paper, we investigated a new and practical tacticallevel decision support system that simultaneously optimizes the liner shipping route and ship schedule designs with time windows for the liner shipping industry. Owing to the constraint of the time windows, the liner shipping route design and the schedule design cannot be dealt with separately. Otherwise, the designed schedule may result in a delay in delivery and a huge increase in operating costs. To deal with the problem, we developed a mathematical programming model to simultaneously determine the optimal port rotation, the optimal arrival time at each port of call, and the optimal sailing speed on each leg, with the objective of minimizing the total cost of the shipping route.
e nonconvex and nonlinear model is then transformed into a mixed-integer linear programming model by reformulating the piecewise functions from the reconstruction of the hard time windows and by approximating the convex nonlinear function through an Mathematical Problems in Engineering 7 outer linear approximation technique. e new model can be solved efficiently by using the current mainstream optimization solver to obtain the optimal solution, which proves the applicability of the model to the problem. e model was applied to a real-world case study involving SITC Container Lines. A number of numerical experiments were conducted, and a sensitivity analysis was performed on the bunker price, port efficiency, and port time windows. First, a higher bunker price and a lower port efficiency result in a higher total cost, and the changes barely affect the port call sequence. Second, with the relaxation of the time windows, a better port rotation with a lower total cost is obtained.
Our model can also accommodate the possibility of violating time windows with only minor modifications. Actually, M 5 can be deemed as a high penalty cost using a premium berth. If M 5 is set to be an appropriate penalty cost, a ship can arrive at the port outside time windows by      paying a penalty cost that is equal to M 5 . en, the objective is to minimize the sum of the vessel operating cost, fuel cost, waiting cost, and penalty cost.
In the future, first, we will extend the proposed model by incorporating the inventory cost. Second, as this study focuses on a single shipping route, future research could extend the investigation of the problem to a liner shipping network composed of multiple shipping routes, on each of which a port can be visited twice a week.
ird, future research could incorporate the proposed problem into a fleet   deployment problem while considering decisions making for coordinating the contractual and spot voyages assignment [28,29]. Fourth, the convergent speed of the solvers is generally slow with large-scale cases, so we will attempt to use appropriate heuristic algorithms to accelerate solution efficiency.

Data Availability
All data included in this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this article.