The SAS module describes how the continuous states of each ship evolve in different modes: sailing, berthing, and reset. Each ship is in a specific mode at any instant. In sailing mode, the ship’s position changes while the quantities of handled containers remain 0; however, the situation is reversed in the berthing mode. The ship enters reset mode when it accomplishes a round-circle voyage; then, its position is reset to 0 for a new voyage.

We consider a liner shipping system with a fleet of n v ships. The number of calling ports is n p and the total number of sailing legs is n l . The distance between adjacent ports in the rotation is called a leg.

(a) Position dynamics can be described as follows:

where T s is the sampling time; x i ( t ) (n mile) represents the position of ship i at instant t ; and v i ( t ) (knots) is the sailing speed of ship i at instant t . M s , M b , and M r represent three modes of sailing, berthing, and resetting, respectively.

(b) The dynamics of container handling can be expressed as follows:

where n i , j ( t ) denotes the quantities of handled containers of ship i at port j , which is the total number of containers being loaded or unloaded, or any other terminal operations; furthermore, n i , j ( t ) is set to 0 when the ship is sailing. h i , j ( t ) expresses the container handling efficiency of port j at instant t . The port handling efficiency is defined as an average container handling rate at port, which is calculated through dividing the total quantities of handled containers by the total observable time at port.

Binary events are triggered when the EG module detects that the continuous state variables reach certain ranges. In the liner shipping system, the events are triggered by the values of ships’ positions and the quantities of their handled containers.

(a) Event triggered by position dynamics is as follows:

where e i , j x ( t ) describes whether ship i has reached port j or not. At instant t , if ship i has reached port j , then e i , j x ( t ) is true; otherwise, e i , j x ( t ) is false.

(b) Event generated when ship i finishes a round-circle voyage is as follows:

where e i c i r ( t ) describes whether ship i has finished a circle voyage or not.

(c) Event triggered by container handling dynamics is as follows:

where e i , j h ( t ) describes whether ship i has finished container handling at port j or not. H i , j is the container quantities to be handled for ship i at port j , which is acquired through observation. e i , j h ( t ) is triggered when the already-handled container quantities exceed the due quantities.

In the MS module, different modes are activated by certain binary states and events, changing patterns as the continuous states evolve. The three modes of liner shipping operations are defined as follows:

For j = 1 ,   …   ,   n p , m = 1 ,   …   ,   n l

where binary state variable s i , m ( t ) represents at instant t that ship i is sailing on leg m ; furthermore, binary state variable b i , j ( t ) represents at instant t that ship i is berthing at port j . If ship i is sailing on a certain leg, then the sailing mode of this ship is activated; if ship i is berthing at any port, then the berthing mode is activated; and if ship i has finished a circle voyage, then the reset mode is activated. When ship i finishes sailing on the last leg and reaches the first port, its position will be reset to 0 at the same time it enters the berthing state and a new circular voyage begins. Therefore, sailing and berthing modes are mutually exclusive, while berthing and reset modes are compatible.

The FSM module describes the evolution of binary states. At every instant, the current binary states are calculated on former states and events generated in the EG module.

(a) The logical rules judging whether ship i is sailing on leg m are defined as follows:

where B i , m s ( t ) and C i , m s ( t ) are auxiliary Boolean variables and function as events in the FSM module. B i , m s ( t ) denotes that ship i begins sailing on leg m ; furthermore, C i , m s ( t ) denotes that ship i is still sailing on leg m . Ship i is deemed to be sailing on leg m if it finishes port handling at port j ( j = m ) or it is sailing on leg m and has not reached the next port.

(b) The logical rules judging whether ship i is berthing at port j are defined as follows:

where B i , j b ( t ) and C i , j b ( t ) are auxiliary Boolean variables and function as events in the FSM module. B i , j b ( t ) denotes that ship i reaches port j ( j = m + 1 ) and begins to berth at port j ; furthermore, C i , j b ( t ) denotes that ship i is still berthing at port j . Ship i is deemed to be berthing at port j if it reaches port j or it is still berthing at port j and port handling has not finished yet.

Figure 1 shows the state flow in the FSM module.

A liner shipping schedule is a plan that combines time and space factors. Ships serving a certain shipping route are required to reach the corresponding port within the specified period. The importance of schedule reliability is reflected in two ways: (a) Container shipping is a key link in international cargo transportation. A liner shipping schedule provides an important basis for planning by other participants in the global supply chain, such as shippers, inland carriers, and marine container terminal (MCT) operators; therefore, any schedule delays would affect the subsequent logistics chain. (b) Schedule reliability reflects the service level a shipping line can achieve for customers; therefore, maintaining high levels of schedule reliability is beneficial for improving customer satisfaction and gaining competitive advantages.

Schedules published by shipping lines contain estimated times of arrival (ETA) and departure (ETD) for all ports they serve. In the liner shipping schedule, ETA and ETD are sequences of days, meaning that ships have to arrive or leave a port on a certain date, with the exact time (hours) not announced. For example, if the ETA for a certain port is 10, then the ship is supposed to arrive at the port within the period from 00:00 to 24:00 on the 10th day after the voyage starts—there is no requirement regarding the exact time (hours) of arrival. Such a practice allows some schedule flexibility.

However, at an operational level, more detailed instructions for arrivals and departures should be made; therefore, the ETA needs to be accurate to hours. An accurate ETA is key to the coordination between ships and MCT operators (Tao et al., 2023). As such, we transform the schedule into a series of position–time coordinates to obtain a tractable reference trajectory—the time series is in hours and the ships have reference positions at any instant.

In the w t h voyage, the estimated time of arrival and departure of ship i at port j (denoted separately by E T A i , j w and E T D i , j w ) can be calculated by predetermined service speed v i s e r (knots), estimated average port handling efficiency h ^ j (TEU/hour), and estimated due quantities of containers to be handled H ^ i , j (TEU). Then, E T A i , j w and E T D i , j w are set as relay points; therefore, we can derive all the position–time coordinates. Algorithm 1 shows the detailed process for generating the reference trajectories. We minimized the difference between the reference trajectories and the ships’ actual trajectories, which is called a tracking error, using a receding horizon scheme in Section 2.3 so that the ships follow the reference trajectories.

Minimizing fuel consumption and carbon emissions are additional objectives when operating the liner shipping system in an energy-efficient way. Therefore, we calculate a formula for fuel consumption and carbon emissions as follows:

where N i e ( t ) is the indicated power (kW) of ship i ; D i is the ship’s deadweight tonnage (DWT); C is the admiralty coefficient; G i ( t ) is the fuel consumption per unit time (t/h); O i ( t ) and E i ( t ) are the fuel consumption (t) and carbon emissions (t) in the time period [ t ,   t + 1 ] ; and C F is the conversion factor between fuel consumption and C O 2 emissions.

For calculation convenience, we transformed the fuel consumption and carbon emissions into linear expressions of speed. We used the cubic root of fuel consumption and carbon emissions as performance indicators of ship energy efficiency. The formulas are as follows:

The EEXI standard is introduced in the 2021 Marine Environmental Protection Committee meeting MEPC76, which acts as an expansion of the Energy Efficiency Design Index (EEDI) and covers existing ships. The attained EEXI for ships in liner shipping lines should be calculated and meet the required EEXI (IMO, 2022). Our calculation formulas for attained EEXI are as follows:

For ship i ,

where R X is the required EEXI ( g / t · nm ), R D is the reference line of EEDI ( g / t · nm ), and A X is the attained EEXI ( g / t · nm ). Y is the reduction factor and Z is the safety margin, which are decided by the shipping lines themselves. V r e f is the service speed under the EEXI draft. f R is the correction factor and its value is 1. E M E and E A E are the carbon emissions of the main and auxiliary engines ( g ); P T I is the emissions of the shaft generator ( g ); and E f f A E and E f f M E are the carbon emissions reduced by the innovative energy efficient devices in the auxiliary and main engines ( g ). Our calculation formulas for the carbon emissions of various devices are as follows:

P M E ( i ) is the effective power of the i t h main engine (kW), which is 75% of the maximum continuous rating (MCR); P A E is the auxiliary engine power; P P T I ( i ) , P A E e f f ( i ) , and P M E e f f ( i ) are the effective power of the shaft motor, and auxiliary and main engine power reductions due to energy efficient technology. f i and f e f f ( i ) are correction factors and f i = f e f f ( i ) = 1 . C F M E ( i ) and C F A E are the carbon converting factors; furthermore, here, we assume C F M E ( i ) = C F A E = C F .

If the ship’s attained EEXI cannot meet the EEXI standard then adjusting measures should be carried out for EEXI compliance. Among various EEXI improvement solutions, EPL has been proven the most common choice by shipowners due to its effectiveness and economic feasibility. Therefore, our solution is to apply engine power limitation (EPL) to improve the ship’s EEXI.

Using the DHA model in Section 2.1 as the system plant, we created a DHA-DMPC framework, which is depicted in Figure 2. There is a local controller on each ship in the liner shipping fleet. The prediction model in the local MPC controller is a hybrid dynamical model built analogously to the DHA model from Section 2.1. We present the prediction model as follows:

where (29) is analogous to the switched affine system in (1) and (2); A M , B M , and C M are the corresponding coefficients in mode M; the continuous state vector x c ( k ) consists of the ship’s position x ( k ) at instant k and handled container quantities n j ( k ) at port j ; the control vector u ( k ) contains the ship’s sailing speed v ( k ) at instant k; and the disturbance vector d ( k ) contains the port handling efficiency h j ( k ) and due quantities of containers to be handled H j at port j . Formula (30) is analogous to the event generator in (3)-(5); event vector δ e ( k ) consists of e j x ( k ) (the ship reaches port j and container handling at port j finishes) and e c i r ( k ) (the ship accomplishes a round-circle voyage); and W e denotes the condition that triggers various events. Formula (31) is analogous to the mode selector in (6)-(8); the mode vector M ( k ) contains the sailing, berthing, and reset modes; the binary state vector x l ( k ) consists of s m ( k ) (the ship is sailing on leg m ) and b j ( k ) (the ship is berthing at port j ); and f M represents the Boolean operation between the binary variables. Formula (32) is analogous to the finite state machine in (9)-(14); furthermore, in addition to those events in (30), event vector δ e ( k ) in (32) also consists of B m s ( k ) (the ship begins to sailing on leg m ), C m s ( k ) (the ship is still sailing on leg m ), B m b ( k ) (the ship begins to berth at port j ), and C j b ( k ) (the ship is still berthing at port j ); and f B represents the Boolean function, which consists of logical rules in the FSM module. Unlike the system plant, which includes all the ships in the fleet, the prediction model only considers a single ship; therefore, it involves fewer binary variables compared with the system plant and, thus, its computation capacity is much lower. This decentralized arrangement increases the computing speed of a local controller and enables real-time optimization of the ships’ speed (Negenborn and Maestre, 2014; Zheng et al., 2017; Wang et al., 2020).

Formula (33) is the optimization model in the DMPC controller, which considers minimizing tracking errors, fuel consumption, and carbon emissions as objectives. Q 1 , Q 2 , and Q 3 are the weights of the three performance indicators ( Q 1 = 4 , Q 2 = Q 3 = 1 ). The control variable is the sailing speed of each ship. We executed (33) using a receding horizon method to compensate for the uncertainties in the disturbances; moreover, binary variables are involved, such as the binary state variables and events in (30)-(32). Therefore, (33) is a mixed integer programming (MIP) problem.

In the receding horizon scheme, the prediction models predict future states in the prediction horizon after observing the current states; then, they calculate the tracking errors, resulting in every local controller computing the MIP optimization problem in (33) to obtain the control sequence in the control horizon. At each step, only the first element of the control sequence is implemented and the rest are ignored. Then, a new control sequence is obtained in the next step. We repeatedly execute this operation to realize online rolling optimization (Algorithm 2).

subject to

for i = 1 ,   … ,   n v ; j = 1 ,   … , n p :

initial states x ( 0 ) , n j ( 0 ) , s m ( 0 ) , b j ( 0 ) ;

predictive model dynamics (29), (30), (31), (32);

EEXI limitations (22);

x i , j , v i , j ≥ 0 .

Various types of liner shipping routes exist, such as end-to-end, pendulum, and circular services (Notteboom, 2004). We choose the classical circular service for our case study, which is commonly used in practical liner shipping route design. Figure 3A shows the actual shipping route layout, and Table 1 displays the route’s parameters. The route can be abstracted as a concept of the circular voyage shown in Figure 3B, which means the ship visits a fixed port rotation on a pre-published schedule repeatedly. This route is served by four ships, and we assume that the fleet consists of four sister ships with the same technical conditions. Our calculations and solutions to improve the ships’ EEXI are shown in Section 3.1.3. All the data on shipping routes and sample ships are provided by CU Lines. The simulation period is the time for the fleet to accomplish two circle voyages, which is 1050 steps. The sampling time is 1h; furthermore, both the prediction and the control horizon are 10 steps. Table 2 shows the data from customer booking platform containing the due quantities of containers to be handled.

Table 3 displays the anticipated port handling efficiency used by shipping lines for tactical planning. The port handling efficiency is not the actual operation rate of the container terminal equipment; instead, it is an average handling efficiency rate, which is obtained by dividing the total number of containers handled by the total time that ships stay at port. However, in practice, port handling efficiency can fluctuate because of unpredictable factors, such as extreme weather, which makes loading and unloading operations difficult; port congestion, which increases the waiting time for ships to berth; or strikes, which result in low handling efficiency. Figure 4 shows the actual port handling efficiency during the simulation period, which we input into the optimization model as disturbances.

We designed two reference trajectories and compared their performances: (a) Reference 1 ignores the berthing process. We set only E T A i , j w as relay points to calculate the position–time coordinates. Algorithm 3 in Appendix A shows the detailed method to generate Reference 1. (b) Reference 2 considers the berthing process and is generated by Algorithm 1. Figure 5 illustrates the two reference trajectories.

Table 4 shows the ship’s parameters, which we used to calculate EEXI.

The relationship between main engine power and speed can be described in a power–speed curve. The power–speed curve should be converted from the service draft to the EEXI draft. Figure 6 shows the transferred power–speed curve. The calculation equation is as follows:

where v r e f is the speed at EEXI draft, v s e r is the speed at service draft, Δ s e r v i c e is the displacement at service draft, and Δ E E X I is the displacement at EEXI draft. For container ships, the service draft is the summer draft; furthermore, the EEXI draft is 70% summer draft.

The required EEXI/reference line is based on the below formula:

where a = 174.22 , c = 0.201 , and b = D .

The attained EEXI is 17.50 g / t · nm while the required EEXI with a 5% margin is 17.38 g / t · nm ; thus, the EPL has to be employed to limit the main engine power. We calculated the MCR after limitations as follows:

The limited MCR is 10787.9 kW, and the maximum speed under this limited MCR is 18 knots.