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

Authors: Xi Jiang; Haijun Mao; Hao Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper proposes to address the problem of the simultaneous optimization of the liner shipping route and ship schedule designs by incorporating port time windows. A mathematical programming model was developed to minimize 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. The 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. The 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.


---

## Body

## 1. 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. This 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. The 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 (Thursday). As shown in Figure1(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.) Therefore, 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. The 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.Figure 1
S2 service provided by SITC. (a) Original sequence. (b) Adjusted sequence.
(a)(b)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. This study provides a decision support system for container lines to improve their services and reduce their operating costs.


### 1.1. Literature Review

This study is related to two research topics: liner shipping network design and ship scheduling. The former has been a hot topic in the field of liner shipping research over recent decades. The work of Rana and Vickson [4] was the early contribution in this field. They 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. The 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. Álvarez [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. They 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. The 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. They 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. Their 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. They 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. This study tries to fill the gap. As mentioned in Section1 of this study, it is highly desirable to simultaneously optimize port call sequence and ship schedule subject to the constraints of time windows.


### 1.2. Objectives and Contributions

This 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. This problem has practical significance because it can help container lines improve service levels and reduce total operating costs.The 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.


## 1.1. Literature Review

This study is related to two research topics: liner shipping network design and ship scheduling. The former has been a hot topic in the field of liner shipping research over recent decades. The work of Rana and Vickson [4] was the early contribution in this field. They 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. The 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. Álvarez [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. They 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. The 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. They 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. Their 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. They 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. This study tries to fill the gap. As mentioned in Section1 of this study, it is highly desirable to simultaneously optimize port call sequence and ship schedule subject to the constraints of time windows.


## 1.2. Objectives and Contributions

This 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. This problem has practical significance because it can help container lines improve service levels and reduce total operating costs.The 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.


## 2. Problem Description and Mathematical Formulation

We consider a group of ports, denoted by a setI. The set I has N number of physical ports indexed by i and j. The voyage from port i to port j is called the shipping leg i,j. We define the set of all shipping legs as Α=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. The details on problem description and key constraints are provided in the following sections.


### 2.1. Port Rotation

The port rotation of a shipping route forms a loop in practice. We use the binary variablesxij∈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≠ixij=∑j∈I,j≠ixji=1 should hold for each port i.


### 2.2. Bunker Consumption

A ship’s unit bunker consumption largely depends on its sailing speed. We denote by the variablevij the sailing speed of a ship on leg i,j, and we define fijvij as the bunker consumption in tons per nautical mile at the sailing speed vij on leg i,j. Based on the results of existing studies [15, 27], the fuel consumption is a power function of the sailing speed: fijvij=aijvijbij, where aij and bij are two coefficients estimated from the practical data. Let β be the unit bunker price (in US dollars per ton), and let Lij 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∈IxijLijaijvijbij.


### 2.3. 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 tia 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, Thursday, 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 tiamod168∈Ωi. For example, a ship can call port i for service at tia=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. Thus, 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∈Itijwcjw, where the variable tijw is the waiting time for the ship in the anchorage at port j after using shipping leg i,j, and cjw is the fixed waiting cost per hour in the anchorage at port j.


### 2.4. 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 variablem represent the number of ships deployed on the route. Thus, the weekly operating cost of the designed route can be calculated by Copm, where Cop is the fixed weekly operating cost in US dollars per week per ship.


### 2.5. Mathematical Model

The 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:xij: binary variables indicating whether shipping leg i,j is used on the designed shipping routeηi: auxiliary variables for subtour elimination constraintsvij: sailing speed of a container ship on leg i,jTij: sailing time of a container ship on leg i,jtijw: waiting time for the ship to be in the anchorage at port j after using leg i,jm: integer variable representing the number of ships deployed on the designed shipping routetia: ship arrival time at port it˜1a: the time when the ship returns to the first port of call after visiting all ports of call of a round-trip journeyt^ia: ship arrival time at port i in the first cyclefijh: numbers of containers of container shipping demand h∈H on leg i,jSets:I: set of all portsA: set of all shipping legsH: set of all container shipping demandsΖ+: set of nonnegative integersParameters:N: number of ports on the designed shipping routeLij: length of shipping leg i,jVijmin: minimum sailing speed of a container ship on leg i,jVijmax: maximum sailing speed of a container ship on leg i,jΩi: port time windows at port icjw: fixed waiting cost per hour in the anchorage at port jtip: fixed dwell time at port iCop: fixed weekly operating cost of one container shipmmax: maximum number of ships deployed on the designed shipping routeDh: volume of container shipping demand hohanddh: port of loading and port of discharging for container shipping demand hM1,M2, andM3: very large numbersβ: unit bunker priceai: a coefficient estimated from the practical databi: a coefficient estimated from the practical dataCap: capacity of one container shipThus, the problem of the simultaneous optimization of the liner shipping route and ship schedule designs with time windows can be formulated as follows:(1)M1minCopm+β∑i∈I∑j∈IxijLijaijvijbij+∑i∈I∑j∈Itijwcjw,subject to(2)∑j∈I,j≠ixij=∑j∈I,j≠ixji=1,∀i∈I,(3)ηi−ηj+Nxij≤N−1,∀i,j∈I,i,j≠1,i≠j,(4)0≤Tij≤M1xij,∀i,j∈I,(5)Lijvij−M11−xij≤Tij≤Lijvij,∀i,j∈I,(6)Vijmin≤vij≤Vijmax,∀i,j∈I,(7)t^ia=tiamod168,i∈I,(8)t^ia∈Ωi,i∈I,(9)tia+tip+Tij+tijw−tja≤M21−xij,∀i∈I,j∈I\1,(10)tia+tip+Tij+tijw−tja≥−M21−xij,∀i∈I,j∈I\1,(11)tia+tip+Ti1+ti1w−t˜1a≤M21−xi1,∀i∈I,(12)tia+tip+Ti1+ti1w−t˜1a≥−M21−xi1,∀i∈I,(13)0≤tijw≤M3xij,∀i,j∈I,(14)0≤ti1w≤M3xi1,∀i∈I,(15)t˜1a≤t1a+168m,(16)m∈1,2,3,……,mmax,(17)∑j∈Ifjih−∑j∈Ifijh=0,∀i≠oh,dh,Dh,∀i=dh,∀h∈H,(18)∑j∈Ifijh=Dh,∀i∈oh,∀h∈H,(19)∑h∈Hfijh≤Cap⋅xij,∀i,j∈A.The 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, Tij=Lij/vij or Tij=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.


## 2.1. Port Rotation

The port rotation of a shipping route forms a loop in practice. We use the binary variablesxij∈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≠ixij=∑j∈I,j≠ixji=1 should hold for each port i.


## 2.2. Bunker Consumption

A ship’s unit bunker consumption largely depends on its sailing speed. We denote by the variablevij the sailing speed of a ship on leg i,j, and we define fijvij as the bunker consumption in tons per nautical mile at the sailing speed vij on leg i,j. Based on the results of existing studies [15, 27], the fuel consumption is a power function of the sailing speed: fijvij=aijvijbij, where aij and bij are two coefficients estimated from the practical data. Let β be the unit bunker price (in US dollars per ton), and let Lij 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∈IxijLijaijvijbij.


## 2.3. 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 tia 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, Thursday, 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 tiamod168∈Ωi. For example, a ship can call port i for service at tia=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. Thus, 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∈Itijwcjw, where the variable tijw is the waiting time for the ship in the anchorage at port j after using shipping leg i,j, and cjw is the fixed waiting cost per hour in the anchorage at port j.


## 2.4. 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 variablem represent the number of ships deployed on the route. Thus, the weekly operating cost of the designed route can be calculated by Copm, where Cop is the fixed weekly operating cost in US dollars per week per ship.


## 2.5. Mathematical Model

The 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:xij: binary variables indicating whether shipping leg i,j is used on the designed shipping routeηi: auxiliary variables for subtour elimination constraintsvij: sailing speed of a container ship on leg i,jTij: sailing time of a container ship on leg i,jtijw: waiting time for the ship to be in the anchorage at port j after using leg i,jm: integer variable representing the number of ships deployed on the designed shipping routetia: ship arrival time at port it˜1a: the time when the ship returns to the first port of call after visiting all ports of call of a round-trip journeyt^ia: ship arrival time at port i in the first cyclefijh: numbers of containers of container shipping demand h∈H on leg i,jSets:I: set of all portsA: set of all shipping legsH: set of all container shipping demandsΖ+: set of nonnegative integersParameters:N: number of ports on the designed shipping routeLij: length of shipping leg i,jVijmin: minimum sailing speed of a container ship on leg i,jVijmax: maximum sailing speed of a container ship on leg i,jΩi: port time windows at port icjw: fixed waiting cost per hour in the anchorage at port jtip: fixed dwell time at port iCop: fixed weekly operating cost of one container shipmmax: maximum number of ships deployed on the designed shipping routeDh: volume of container shipping demand hohanddh: port of loading and port of discharging for container shipping demand hM1,M2, andM3: very large numbersβ: unit bunker priceai: a coefficient estimated from the practical databi: a coefficient estimated from the practical dataCap: capacity of one container shipThus, the problem of the simultaneous optimization of the liner shipping route and ship schedule designs with time windows can be formulated as follows:(1)M1minCopm+β∑i∈I∑j∈IxijLijaijvijbij+∑i∈I∑j∈Itijwcjw,subject to(2)∑j∈I,j≠ixij=∑j∈I,j≠ixji=1,∀i∈I,(3)ηi−ηj+Nxij≤N−1,∀i,j∈I,i,j≠1,i≠j,(4)0≤Tij≤M1xij,∀i,j∈I,(5)Lijvij−M11−xij≤Tij≤Lijvij,∀i,j∈I,(6)Vijmin≤vij≤Vijmax,∀i,j∈I,(7)t^ia=tiamod168,i∈I,(8)t^ia∈Ωi,i∈I,(9)tia+tip+Tij+tijw−tja≤M21−xij,∀i∈I,j∈I\1,(10)tia+tip+Tij+tijw−tja≥−M21−xij,∀i∈I,j∈I\1,(11)tia+tip+Ti1+ti1w−t˜1a≤M21−xi1,∀i∈I,(12)tia+tip+Ti1+ti1w−t˜1a≥−M21−xi1,∀i∈I,(13)0≤tijw≤M3xij,∀i,j∈I,(14)0≤ti1w≤M3xi1,∀i∈I,(15)t˜1a≤t1a+168m,(16)m∈1,2,3,……,mmax,(17)∑j∈Ifjih−∑j∈Ifijh=0,∀i≠oh,dh,Dh,∀i=dh,∀h∈H,(18)∑j∈Ifijh=Dh,∀i∈oh,∀h∈H,(19)∑h∈Hfijh≤Cap⋅xij,∀i,j∈A.The 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, Tij=Lij/vij or Tij=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.


## 3. 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. These 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.


### 3.1. Linearization of the Bunker Consumptions

First, we define a new variableuij to be the reciprocal of the sailing speed vij(20)uij=1vij,∀i,j∈I.Then, constraints (5) and (6) can be transformed into linear constraints associated with the variable uij:(21)Lijuij−M11−xij≤Tij≤Lijuij,∀i,j∈I,(22)Uijmin=1Vijmax≤uij≤1Vijmin=Uijmax,∀i,j∈I.We introduce an auxiliary variableBij to replace the nonlinear combinational items xij and vij in the objective function. We define Bij=xijLijaijuij−bij. Hence, we have Bij=Lijaijuij−bij if xij=1, i.e., leg i,j is used; otherwise, Bij=0. By introducing a very large number M4, we can replace the above if-then complementary conditions equivalently with the big-M reformulations of the following constraints:(23)0≤Bij≤M4xij,∀i,j∈I,(24)Lijaijuij−bij−M41−xij≤Bij≤Lijaijuij−bij,∀i,j∈I.By introducing the auxiliary variablesgij, we define gij=aijuij−bij. Objective function (1) can be transformed by constructing its epigraph form as follows:(25)minCopm+β∑i∈I∑j∈IBij+∑i∈I∑j∈Itijwcjw,subject to(26)gij≥aijuij−bij,∀i,j∈I,and to constraints (23) and (24).We definehijuij=aijuij−bij. Because −bij≤−2<0 and aij>0, hijuij is convex in uij∈UijminandUijmax. We divide the interval UijminandUijmax into M equal segments:(27)Uijmin=uij1<uij2<⋯<uijM=Uijmax,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 uijm to represent any value of the M equal segments. From the property of the first-order condition of the convex function, we can obtain(28)hijuij≥hijuijm+hij′uijmuij−uijm,m=1,2,…,M,∀i,j∈I,where hij′uijm represents the tangent line of hijuij at the point uijm. Thus, the function hijuij=aijuij−bij can be converted into the following piecewise function:(29)hijuij=maxm=1,…,Mhijuijm+hij′uijmuij−uijm,∀i,j∈I.By combining constraints (26) and (29), we can obtain M separate inequalities(30)gij≥hijuijm+hij′uijmuij−uijm,m=1,…,M,∀i,j∈I.


### 3.2. Transformation of the Time Windows

Since the setΩi may be a discrete domain, by transforming the hard time windows into soft time windows, we can obtain piecewise functions that can eliminate the discrete domain and can be linearized. For example, we assume that the time windows at port i are Tuesday and Thursday. By introducing a very large number M5, we can express the penalty function pit^ia at port i as follows:(31)pit^ia=M5,0≤t^ia<48,0,48≤t^ia<72,M5,72≤t^ia<96,0,96≤t^ia<120,M5,120≤t^ia<168.The penalty cost is high enough to limit a ship to call the port only within the time windows. For example, as shown in Figure2, ships can only call port i on Tuesdays and Thursdays.Figure 2
Penalty functional image of soft time windows on Tuesdays and Thursdays.Next, we linearize the penalty functionpit^ia. Since pit^ia has ni turning points, let bik be the turning points and let Ki:=1,2,…,ni be the set of all turning points at port i. For example, as shown in Figure 2, there are 10 turning points: bi,1=0,bi,2=48,bi,3=48,bi,4=72,bi,5=72,bi,6=96,bi,7=96,bi,8=120,bi,9=120,andbi,10=168. By introducing the continuous auxiliary variables wik and the binary auxiliary variables zik, we can transform the piecewise function into the following linear constraints:(32)Pi=∑k∈Kiwikpibik,∀i∈I,(33)t^ia=∑k∈Kiwikbik,∀i∈I,(34)wi1≤zi1,∀i∈I,(35)wik≤zi,k−1+zik,k=2,…,ni−1,∀i∈I,(36)wi,ni≤zi,ni−1,∀i∈I,(37)∑k∈Kiwik=1,∀i∈I,(38)∑k=1ni−1zik=1,∀i∈I,(39)wik≥0,k=1,…,ni,∀i∈I,(40)zik∈0,1,k=1,…,ni−1,∀i∈I,where wik can realize the function linearization and zik can limit the value range of wik to limit the interval of t^ia.Thus, the total penalty costs of time windows at all ports are as follows:(41)∑i∈IPi=∑i∈I∑k∈Kiwikpibik.By introducing the auxiliary variablesθi∈Ζ+, we have(42)t^ia=tia−168θi,∀i∈I,(43)0≤t^ia≤168,∀i∈I,(44)θi∈Ζ+,∀i∈I.It is noticeable that constraints (42)–(44) are equivalent to constraint (7).Hence, an approximated mixed-integer linear optimization model [M2] can be obtained, and it can be solved efficiently by using state-of-the-art solvers, subject to constraints (2)–(4), (9)–(19), (21)–(24), (30), and (32)–(44).(45)M2minCopm+β∑i∈I∑j∈IBij+∑i∈I∑j∈Itijwcjw+∑i∈I∑k∈Kiwikpibik.


## 3.1. Linearization of the Bunker Consumptions

First, we define a new variableuij to be the reciprocal of the sailing speed vij(20)uij=1vij,∀i,j∈I.Then, constraints (5) and (6) can be transformed into linear constraints associated with the variable uij:(21)Lijuij−M11−xij≤Tij≤Lijuij,∀i,j∈I,(22)Uijmin=1Vijmax≤uij≤1Vijmin=Uijmax,∀i,j∈I.We introduce an auxiliary variableBij to replace the nonlinear combinational items xij and vij in the objective function. We define Bij=xijLijaijuij−bij. Hence, we have Bij=Lijaijuij−bij if xij=1, i.e., leg i,j is used; otherwise, Bij=0. By introducing a very large number M4, we can replace the above if-then complementary conditions equivalently with the big-M reformulations of the following constraints:(23)0≤Bij≤M4xij,∀i,j∈I,(24)Lijaijuij−bij−M41−xij≤Bij≤Lijaijuij−bij,∀i,j∈I.By introducing the auxiliary variablesgij, we define gij=aijuij−bij. Objective function (1) can be transformed by constructing its epigraph form as follows:(25)minCopm+β∑i∈I∑j∈IBij+∑i∈I∑j∈Itijwcjw,subject to(26)gij≥aijuij−bij,∀i,j∈I,and to constraints (23) and (24).We definehijuij=aijuij−bij. Because −bij≤−2<0 and aij>0, hijuij is convex in uij∈UijminandUijmax. We divide the interval UijminandUijmax into M equal segments:(27)Uijmin=uij1<uij2<⋯<uijM=Uijmax,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 uijm to represent any value of the M equal segments. From the property of the first-order condition of the convex function, we can obtain(28)hijuij≥hijuijm+hij′uijmuij−uijm,m=1,2,…,M,∀i,j∈I,where hij′uijm represents the tangent line of hijuij at the point uijm. Thus, the function hijuij=aijuij−bij can be converted into the following piecewise function:(29)hijuij=maxm=1,…,Mhijuijm+hij′uijmuij−uijm,∀i,j∈I.By combining constraints (26) and (29), we can obtain M separate inequalities(30)gij≥hijuijm+hij′uijmuij−uijm,m=1,…,M,∀i,j∈I.


## 3.2. Transformation of the Time Windows

Since the setΩi may be a discrete domain, by transforming the hard time windows into soft time windows, we can obtain piecewise functions that can eliminate the discrete domain and can be linearized. For example, we assume that the time windows at port i are Tuesday and Thursday. By introducing a very large number M5, we can express the penalty function pit^ia at port i as follows:(31)pit^ia=M5,0≤t^ia<48,0,48≤t^ia<72,M5,72≤t^ia<96,0,96≤t^ia<120,M5,120≤t^ia<168.The penalty cost is high enough to limit a ship to call the port only within the time windows. For example, as shown in Figure2, ships can only call port i on Tuesdays and Thursdays.Figure 2
Penalty functional image of soft time windows on Tuesdays and Thursdays.Next, we linearize the penalty functionpit^ia. Since pit^ia has ni turning points, let bik be the turning points and let Ki:=1,2,…,ni be the set of all turning points at port i. For example, as shown in Figure 2, there are 10 turning points: bi,1=0,bi,2=48,bi,3=48,bi,4=72,bi,5=72,bi,6=96,bi,7=96,bi,8=120,bi,9=120,andbi,10=168. By introducing the continuous auxiliary variables wik and the binary auxiliary variables zik, we can transform the piecewise function into the following linear constraints:(32)Pi=∑k∈Kiwikpibik,∀i∈I,(33)t^ia=∑k∈Kiwikbik,∀i∈I,(34)wi1≤zi1,∀i∈I,(35)wik≤zi,k−1+zik,k=2,…,ni−1,∀i∈I,(36)wi,ni≤zi,ni−1,∀i∈I,(37)∑k∈Kiwik=1,∀i∈I,(38)∑k=1ni−1zik=1,∀i∈I,(39)wik≥0,k=1,…,ni,∀i∈I,(40)zik∈0,1,k=1,…,ni−1,∀i∈I,where wik can realize the function linearization and zik can limit the value range of wik to limit the interval of t^ia.Thus, the total penalty costs of time windows at all ports are as follows:(41)∑i∈IPi=∑i∈I∑k∈Kiwikpibik.By introducing the auxiliary variablesθi∈Ζ+, we have(42)t^ia=tia−168θi,∀i∈I,(43)0≤t^ia≤168,∀i∈I,(44)θi∈Ζ+,∀i∈I.It is noticeable that constraints (42)–(44) are equivalent to constraint (7).Hence, an approximated mixed-integer linear optimization model [M2] can be obtained, and it can be solved efficiently by using state-of-the-art solvers, subject to constraints (2)–(4), (9)–(19), (21)–(24), (30), and (32)–(44).(45)M2minCopm+β∑i∈I∑j∈IBij+∑i∈I∑j∈Itijwcjw+∑i∈I∑k∈Kiwikpibik.


## 4. Case Study

This 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. Then, we conducted a sensitivity analysis to test some key input parameters in the solution approach.


### 4.1. Parameter Settings

The PBT1 service includes nine ports of call, and its port rotation is Qinhuangdao (QHD) ⟶ Xingang (TXG) ⟶ Dalian (DAL) ⟶ Nagoya (NAG) ⟶ Yokkaichi (YCI) ⟶ Shimizu (SZU) ⟶ Tokyo (TOK) ⟶ Yokohama (YOK) ⟶ Kawasaki (KAW) ⟶ Qinhuangdao (QHD). We assumed that 2000-TEU ships were deployed on the designed route, i.e.,cap=2,000. The operating cost Cop=US$200,000/week, the maximum speed Vijmax=25knots, the minimum speed Vijmin=5knots, the bunker price β=US$300/t, the maximum number of ships deployed mmax=5, and the fixed waiting cost in the anchorage at each port cjw=US$100/h. For ease of calculation, we take the tightest upper bounds; then, the very large numbers (M1,M2,M3, and M4) can be set as maxLij/Vijmin, 168mmax, 168, and maxLij⋅aijVijmaxbij, respectively. To avoid violating the time window constraints, the very large number M5 was set to 300,000. The number of segments for the interval UijminandUijmax was M=200. The bunker consumption function of each leg was set to 0.001v2. The 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. The 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.Table 1
Port-to-port distance matrix (in nautical miles).
PortQHDTXGDALNAGYCITOKYOKKAWSZUQHD0131.4167.21,169.51,165.21,297.81,281.61,288.21,212.7TXG131.40218.71,221.11,216.71,349.41,333.21,339.71,264.3DAL167.2218.701,058.11,053.71,186.41,170.21,176.71,101.3NAG1,169.51,221.11,058.1014.1225.6209.4215.5139.1YCI1,165.21,216.71,053.714.10221.5205.2211.3135TOK1,297.81,349.41,186.4225.6221.5021.218135.9YOK1,281.61,333.21,170.2209.4205.221.2011.6119.7KAW1,288.21,339.71,176.7215.5211.31811.60125.8SZU1,212.71,264.31,101.3139.1135135.9119.7125.80Source: result of the Port Distance Desktop software.Table 2
Dwell time at each port (in hours).
PortQHDTXGDALNAGYCITOKYOKKAWSZUDwell time20132412420886Table 3
Available time at each port.
PortSunMonTueWedThuFriSatQHDFreeTXGFreeDALFreeNAGFreeYCIFreeTOKFreeYOKFreeKAWFreeSZUFreeThe blanks indicate that the port is busy on this day.


### 4.2. Comparison with the Original Port Rotation

In this section, we present the results of the validation of the proposed model and the solution approach. The results obtained using the parameters in Section4.1 are shown in Tables 4 and 5. It is noticeable that the arrival times do not violate port time window constraints. The 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]. The total cost is also reduced by 45% compared with that in [3]. Therefore, the proposed model and the solution approach can greatly reduce the total cost.Table 4
Optimal sailing speed on each leg.
LegQHD⟶TXG⟶DAL⟶YCI⟶NAG⟶YOK⟶KAW⟶TOK⟶SZU⟶Speed (knots)10.7319.8821.9510.7310.7310.7310.7310.7310.68Table 5
Optimal arrival time at each port.
Port callQHDTXGDALYCINAGTOKKAWTOKSZUQHDArrival time87.75120.00144.00216.00221.31252.83261.91271.59304.25423.75Table 6
Comparison between the model [M2] and the model of Wang et al. [3].
ModelTotal cost (×105 US$)Fleet sizeTotal waiting time in the anchorage (h)The model [M2]6.3820Route: QHD ⟶ TXG⟶DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHDThe model of [3]11.65400.76Route: QHD ⟶ TXG ⟶ DAL ⟶ NAG ⟶ YCI ⟶ SZU ⟶ TOK ⟶ YOK ⟶ KAW ⟶ QHD


### 4.3. 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 section4.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. The results are shown in Table 8. It is noticeable that the total cost decreases and that the port rotation is different. This is because, with the relaxation of the time windows, there will be a better port rotation to reduce the total cost.Table 7
Relaxations of the time windows at each port.
Port callScenarioSunMonTueWedThuFriSatQHD2FreeFree3FreeFreeFreeTXG2FreeFree3FreeFreeFreeDAL2FreeFree3FreeFreeFreeNAG2FreeFree3FreeFreeFreeYCI2FreeFree3FreeFreeFreeTOK2FreeFree3FreeFreeFreeYOK2FreeFree3FreeFreeFreeKAW2FreeFree3FreeFreeFreeSZU2FreeFree3FreeFreeFreeTable 8
Optimal solutions with different relaxations of the time windows.
ScenarioFleet sizeTotal cost (×105 US$)126.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD225.67Route: QHD ⟶ DAL ⟶ YCI ⟶ NAG ⟶ SZU ⟶ YOK ⟶ KAW ⟶ TOK ⟶ TXG ⟶ QHD325.66Route: QHD ⟶ DAL ⟶ YCI ⟶ NAG ⟶ SZU ⟶ KAW ⟶ TOK ⟶ YOK ⟶ TXG ⟶ QHDThen, 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 Table9. 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. The 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.Table 9
Different time windows at each port.
Port callScenarioSunMonTueWedThuFriSatQHD4Free5FreeTXG4Free5FreeDAL4Free5FreeNAG4Free5FreeYCI4Free5FreeTOK4Free5FreeYOK4Free5FreeKAW4Free5FreeSZU4Free5FreeTable 10
Optimal solutions with different time windows.
ScenarioFleet sizeTotal cost (×105 US$)126.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD426.13Route: QHD ⟶ NAG ⟶ TOK ⟶ YOK ⟶ KAW ⟶ SZU ⟶ YCI ⟶ DAL ⟶ TXG ⟶ QHD526.29Route: QHD ⟶ DAL ⟶ TXG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ NAG ⟶ YCI ⟶ QHD


### 4.4. Effects of Bunker Prices

The 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 Section4.1. The 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. The 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.Table 11
Optimal solutions with different bunker prices.
Bunker price (US$/t)Fleet sizeTotal cost (×105 US$)20025.59Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD30026.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD40036.86Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD50037.08Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD60037.29Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD


### 4.5. 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. The results are shown in Table12. We can see that, with the decrease in port efficiency, the total cost increased significantly. The port rotation changed only when the dwell time was +100% of the benchmark value, and in other cases, it remained the same. This 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.Table 12
Optimal solutions with different port efficiencies.
Dwell time (%)Fleet sizeTotal cost (×105 US$)−5025.47Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD−2525.83Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD026.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD2536.79Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD5037.04Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD10038.22Route: QHD⟶TXG⟶KAW⟶TOK⟶SZU⟶YCI⟶NAG⟶YOK⟶DAL⟶QHD


## 4.1. Parameter Settings

The PBT1 service includes nine ports of call, and its port rotation is Qinhuangdao (QHD) ⟶ Xingang (TXG) ⟶ Dalian (DAL) ⟶ Nagoya (NAG) ⟶ Yokkaichi (YCI) ⟶ Shimizu (SZU) ⟶ Tokyo (TOK) ⟶ Yokohama (YOK) ⟶ Kawasaki (KAW) ⟶ Qinhuangdao (QHD). We assumed that 2000-TEU ships were deployed on the designed route, i.e.,cap=2,000. The operating cost Cop=US$200,000/week, the maximum speed Vijmax=25knots, the minimum speed Vijmin=5knots, the bunker price β=US$300/t, the maximum number of ships deployed mmax=5, and the fixed waiting cost in the anchorage at each port cjw=US$100/h. For ease of calculation, we take the tightest upper bounds; then, the very large numbers (M1,M2,M3, and M4) can be set as maxLij/Vijmin, 168mmax, 168, and maxLij⋅aijVijmaxbij, respectively. To avoid violating the time window constraints, the very large number M5 was set to 300,000. The number of segments for the interval UijminandUijmax was M=200. The bunker consumption function of each leg was set to 0.001v2. The 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. The 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.Table 1
Port-to-port distance matrix (in nautical miles).
PortQHDTXGDALNAGYCITOKYOKKAWSZUQHD0131.4167.21,169.51,165.21,297.81,281.61,288.21,212.7TXG131.40218.71,221.11,216.71,349.41,333.21,339.71,264.3DAL167.2218.701,058.11,053.71,186.41,170.21,176.71,101.3NAG1,169.51,221.11,058.1014.1225.6209.4215.5139.1YCI1,165.21,216.71,053.714.10221.5205.2211.3135TOK1,297.81,349.41,186.4225.6221.5021.218135.9YOK1,281.61,333.21,170.2209.4205.221.2011.6119.7KAW1,288.21,339.71,176.7215.5211.31811.60125.8SZU1,212.71,264.31,101.3139.1135135.9119.7125.80Source: result of the Port Distance Desktop software.Table 2
Dwell time at each port (in hours).
PortQHDTXGDALNAGYCITOKYOKKAWSZUDwell time20132412420886Table 3
Available time at each port.
PortSunMonTueWedThuFriSatQHDFreeTXGFreeDALFreeNAGFreeYCIFreeTOKFreeYOKFreeKAWFreeSZUFreeThe blanks indicate that the port is busy on this day.


## 4.2. Comparison with the Original Port Rotation

In this section, we present the results of the validation of the proposed model and the solution approach. The results obtained using the parameters in Section4.1 are shown in Tables 4 and 5. It is noticeable that the arrival times do not violate port time window constraints. The 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]. The total cost is also reduced by 45% compared with that in [3]. Therefore, the proposed model and the solution approach can greatly reduce the total cost.Table 4
Optimal sailing speed on each leg.
LegQHD⟶TXG⟶DAL⟶YCI⟶NAG⟶YOK⟶KAW⟶TOK⟶SZU⟶Speed (knots)10.7319.8821.9510.7310.7310.7310.7310.7310.68Table 5
Optimal arrival time at each port.
Port callQHDTXGDALYCINAGTOKKAWTOKSZUQHDArrival time87.75120.00144.00216.00221.31252.83261.91271.59304.25423.75Table 6
Comparison between the model [M2] and the model of Wang et al. [3].
ModelTotal cost (×105 US$)Fleet sizeTotal waiting time in the anchorage (h)The model [M2]6.3820Route: QHD ⟶ TXG⟶DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHDThe model of [3]11.65400.76Route: QHD ⟶ TXG ⟶ DAL ⟶ NAG ⟶ YCI ⟶ SZU ⟶ TOK ⟶ YOK ⟶ KAW ⟶ QHD


## 4.3. 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 section4.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. The results are shown in Table 8. It is noticeable that the total cost decreases and that the port rotation is different. This is because, with the relaxation of the time windows, there will be a better port rotation to reduce the total cost.Table 7
Relaxations of the time windows at each port.
Port callScenarioSunMonTueWedThuFriSatQHD2FreeFree3FreeFreeFreeTXG2FreeFree3FreeFreeFreeDAL2FreeFree3FreeFreeFreeNAG2FreeFree3FreeFreeFreeYCI2FreeFree3FreeFreeFreeTOK2FreeFree3FreeFreeFreeYOK2FreeFree3FreeFreeFreeKAW2FreeFree3FreeFreeFreeSZU2FreeFree3FreeFreeFreeTable 8
Optimal solutions with different relaxations of the time windows.
ScenarioFleet sizeTotal cost (×105 US$)126.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD225.67Route: QHD ⟶ DAL ⟶ YCI ⟶ NAG ⟶ SZU ⟶ YOK ⟶ KAW ⟶ TOK ⟶ TXG ⟶ QHD325.66Route: QHD ⟶ DAL ⟶ YCI ⟶ NAG ⟶ SZU ⟶ KAW ⟶ TOK ⟶ YOK ⟶ TXG ⟶ QHDThen, 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 Table9. 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. The 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.Table 9
Different time windows at each port.
Port callScenarioSunMonTueWedThuFriSatQHD4Free5FreeTXG4Free5FreeDAL4Free5FreeNAG4Free5FreeYCI4Free5FreeTOK4Free5FreeYOK4Free5FreeKAW4Free5FreeSZU4Free5FreeTable 10
Optimal solutions with different time windows.
ScenarioFleet sizeTotal cost (×105 US$)126.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD426.13Route: QHD ⟶ NAG ⟶ TOK ⟶ YOK ⟶ KAW ⟶ SZU ⟶ YCI ⟶ DAL ⟶ TXG ⟶ QHD526.29Route: QHD ⟶ DAL ⟶ TXG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ NAG ⟶ YCI ⟶ QHD


## 4.4. Effects of Bunker Prices

The 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 Section4.1. The 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. The 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.Table 11
Optimal solutions with different bunker prices.
Bunker price (US$/t)Fleet sizeTotal cost (×105 US$)20025.59Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD30026.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD40036.86Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD50037.08Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD60037.29Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD


## 4.5. 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. The results are shown in Table12. We can see that, with the decrease in port efficiency, the total cost increased significantly. The port rotation changed only when the dwell time was +100% of the benchmark value, and in other cases, it remained the same. This 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.Table 12
Optimal solutions with different port efficiencies.
Dwell time (%)Fleet sizeTotal cost (×105 US$)−5025.47Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD−2525.83Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD026.38Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD2536.79Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD5037.04Route: QHD ⟶ TXG ⟶ DAL ⟶ YCI ⟶ NAG ⟶ YOK ⟶ KAW ⟶ TOK ⟶ SZU ⟶ QHD10038.22Route: QHD⟶TXG⟶KAW⟶TOK⟶SZU⟶YCI⟶NAG⟶YOK⟶DAL⟶QHD


## 5. Conclusions and Future Work

In this paper, we investigated a new and practical tactical-level 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. The 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 outer linear approximation technique. The 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.The 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,M5 can be deemed as a high penalty cost using a premium berth. If M5 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 M5. Then, 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. Third, 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.



---

*Source: 3287973-2020-12-22.xml*