A Vessel Schedule Recovery Problem at the Liner Shipping Route with Emission Control Areas

: Liner shipping is a vital component of the world trade. Liner shipping companies usually operate ﬁxed routes and announce their schedules. However, disruptions in sea and / or at ports a ﬀ ect the planned vessel schedules. Moreover, some liner shipping routes pass through the areas, designated by the International Maritime Organization (IMO) as emission control areas (ECAs). IMO imposes restrictions on the type of fuel that can be used by vessels within ECAs. The vessel schedule recovery problem becomes more complex when disruptions occur at such liner shipping routes, as liner shipping companies must comply with the IMO regulations. This study presents a novel mixed-integer nonlinear mathematical model for the green vessel schedule recovery problem, which considers two recovery strategies, including vessel sailing speed adjustment and port skipping. The objective aims to minimize the total proﬁt loss, endured by a given liner shipping company due to disruptions in the planned operations. The nonlinear model is linearized and solved using CPLEX. A number of computational experiments are conducted for the liner shipping route, passing through ECAs. Important managerial insights reveal that the proposed methodology can assist liner shipping companies with e ﬃ cient vessel schedule recovery, minimize the monetary losses due to disruptions in vessel schedules, and improve energy e ﬃ ciency as well as environmental sustainability. experiments. M.K. assisted with encoding the developed solution methodology within the MATLAB environment, performing the computational experiments, and interpretation of the obtained results.


Background
Maritime transportation is a critical component of the global trade. According to the United Nations Conference on Trade and Development (UNCTAD), more than 80% of the global trade tonnage and 70% of the global trade value are carried by oceangoing vessels around the world [1]. Over the past three decades, liner shipping, which involves transporting containerized cargoes for various customers, has become a significant part of maritime transportation [1]. The latter tendency has been predicted to be sustained. In liner shipping, a vessel transports containerized cargo to a set of ports based on a fixed schedule. Usually, liner shipping companies form a round trip and call at several ports during a voyage (the term "voyage" refers to a single round-trip). Moreover, their schedules and transit routes are published. According to the World Shipping Council [2], there are approximately 400 liner services around the world, which adopt a regular service frequency (mostly weekly). For a weekly service frequency, each port of call at the liner shipping route is visited the same day every week. Generally, the service network is designed by the liner shipping company at the strategic level of decision-making over a fixed time period [3]. Other important components of liner shipping operations, including environmental regulations, enforced by IMO, and proposes two vessel schedule recovery strategies. The proposed nonlinear mathematical model is linearized and solved to optimality using CPLEX.
A number of computational experiments are conducted for the Asia-North Europe LL5 route, which is served by OOCL liner shipping company and passes through ECAs with the SO x emission control. Important managerial insights, which would be of interest to liner shipping companies, are revealed using the proposed mathematical model and the developed solution approach. The proposed methodology (i.e., the mathematical model and the solution approach) will assist liner shipping companies with efficient vessel schedule recovery, minimizing monetary losses endured as a result of disruptions in vessel schedules, improving energy efficiency and environmental sustainability throughout the transportation process. The latter aspects are considered as critical [12][13][14][15].
The rest of the manuscript is organized in the following manner. The Section 2 presents a review of the state-of-the-art with a focus on vessel schedule recovery and green vessel scheduling. The Section 3 provides a detailed description of the green vessel schedule recovery problem, investigated in this study, while the Section 4 presents a mixed-integer nonlinear model for the problem. The Section 5 presents the solution methodology, adopted in this study, while the Section 6 describes in detail the computational experiments that were conducted in this study and showcases the important managerial insights that were revealed. The Section 7 summarizes the entire research efforts and highlights potential future research directions.

Review of the Relevant Literature
The studies, collected from the state-of-the-art, were classified into two main categories: (a) vessel schedule recovery (studies that focus on disruptions in liner shipping and suggest various vessel schedule recovery strategies); and (b) green vessel scheduling (studies that propose various mathematical models and solution approaches, which address the environmental issues in vessel routing and scheduling). Both study categories are presented in the following sections of the manuscript.

Vessel Schedule Recovery
Various uncertainties, associated with the liner shipping operations, have received attention in the literature; however, only a few studies examined the vessel schedule recovery problem in liner shipping. Paul and Maloni [16] presented a nonlinear programming model for disruptions at seaports due to disasters, minimizing the total operational cost. The results from the numerical experiments, conducted for the North American container port network, indicated that single and multiple port shutdown scenarios incurred additional costs. Brouer et al. [4] conducted a pioneering research on vessel schedule recovery. The authors studied a mixed-integer vessel schedule recovery problem of NP-hard complexity to evaluate the effects of disruptions. The considered vessel schedule recovery strategies included: (i) swap port calls; (ii) omit a port of call; and (iii) adjust vessel sailing speeds. It was found that a 58% cost reduction was achieved when the liner shipping company decided to omit a port of call or swap ports of call, aiming to ensure timely cargo delivery.
Lee et al. [17] proposed a model that incorporated port time uncertainty and slow steaming in determining service reliability, fuel consumption, and emissions. The study underlined that delays at ports might require liner shipping companies to speed up in order to recover the vessel schedules. The decision to implement the latter strategy could increase the fuel consumption and the associated emissions. Li et al. [9] formulated the vessel schedule recovery problem and considered the following recovery options: (a) speeding up only; (b) port skipping only; and (c) port swapping. The results showed that if the delay was not too large, the option of speeding up was preferred for the vessel schedule recovery, as the vessel might not need to sail at the maximum speed to recover the schedule. Furthermore, port swapping and port skipping were found to be better alternatives with significant cost savings if the delay was large. It was also found that additional buffer time could yield cost savings for the recovered vessel schedules. Wang and Meng [18] and Qi and Song [19] highlighted that short segments of a given liner shipping route require more buffer time due to weak flexibility. Li et al. [10] focused on a real-time schedule recovery problem for a liner shipping service. Two types of uncertainties in liner shipping operations were defined: (1) regular uncertainties; and (2) disruptive events. The results from numerical experiments indicated that, for scenarios without earliest handling time constraints at ports, skipping the disruption port might allow the liner shipping company effectively recovering the vessel schedule. Cheraghchi et al. [11] formulated the vessel schedule recovery problem as a bi-objective optimization problem, minimizing the total monetary losses and the total delay due to disruptions. The authors focused on a speed-based vessel schedule recovery. A total of six different multi-objective algorithms were applied to solve the problem. The results from computational experiments showed that the NSGA-II algorithm performed better than other multi-objective metaheuristics used.

Green Vessel Scheduling
Although vessel scheduling and routing has received a lot of attention from researchers [20][21][22][23][24][25][26][27][28], only a very few studies have addressed the vessel scheduling problem, taking into account the environmental concerns [29]. For example, Chang and Wang [30] conducted a comprehensive assessment of the approaches used in commercial maritime transportation to determine the optimal sailing speed of vessels. The results demonstrated that speed reduction was the most efficient for the cases with high unit fuel costs and low charter rates. Kontovas [31] presented a generalized formulation for the green vessel routing and scheduling problem and discussed the need for adopting more accurate models in order to estimate the vessel emissions and the fuel consumption. Psaraftis and Kontovas [32] investigated the factors that influence the choice of vessel sailing speed. Some of the critical factors in determining the vessel sailing speed included the unit fuel cost, state of the market, freight rates, inventory costs, and dependency of the fuel consumption on payload.
Fagerholt et al. [33] evaluated the impacts of environmental regulations in ECAs on sailing paths, sailing speed, and operational costs in maritime shipping. It was found that liner shipping companies might choose to sail along longer voyage legs in order to avoid sailing at a lower speed through ECAs. Fagerholt and Psaraftis [34] presented two speed optimization models for the vessels sailing through ECAs. The numerical computations revealed that the best alternative for the vessels would be sailing through the shortest route within the ECA in order to minimize the total fuel cost and maximize the total profit. Mansouri et al. [29] examined the use of multi-objective decision support tools for enhancing the environmental sustainability of maritime transportation. The authors indicated that simulation and optimization are promising approaches to address the emission reduction challenges in maritime shipping. The importance of modeling emissions in multi-objective vessel scheduling was also discussed by Dulebenets [35]. Song et al. [36] formulated a stochastic multi-objective vessel scheduling problem with uncertain port times, which included the following objectives: (i) minimization of the annual vessel operational costs; (ii) minimization of the average schedule unreliability; and (iii) minimization of the annual total CO 2 emissions. The numerical experiments indicated that about 5-10% reduction in CO 2 emissions could be achieved when either the operating cost is minimized or the schedule unreliability is minimized.
Dulebenets [37] extended the work, conducted by Dulebenets et al. [38], and assessed the benefits and the drawbacks of introducing restrictions on the emissions, produced by oceangoing vessels within ECAs. The study compared the existing IMO regulations regarding the fuel sulfur content while sailing through ECAs with an alternative environmental policy, which imposed restrictions on the amount of emissions, produced within ECAs, along with the existing IMO regulations. The results from computational experiments demonstrated that the emission restrictions reduced the SO 2 emissions by ≈40% but increased the total route service cost. Dulebenets [39] also focused on the green vessel scheduling problem, considering the liner shipping routes with ECAs. It was found that imposing the transit time constraints along with the emission restrictions could lead to significant changes in the vessel schedules and incur additional route service costs. Dulebenets [40] designed a mathematical model for green vessel scheduling, considering CO 2 emission costs of oceangoing vessels in sea and at Energies 2019, 12, 2380 5 of 28 ports of call due to container handling. The numerical experiments showed how changes in the CO 2 tax could affect vessel schedules.

Literature Summary and Contributions
A review of the literature revealed that the vessel schedule recovery problem in liner shipping has not received a lot of attention from researchers. A number of strategies have been proposed in the literature for the vessel schedule recovery. The most common strategies include: (a) vessel sailing speed adjustment; (b) port skipping; and (c) port swapping. Some studies have developed mathematical models, which incorporated one or more vessel schedule recovery strategies. On the other hand, green vessel scheduling has received growing attention from researchers. A number of vessel scheduling models, which account for the environmental issues, energy efficiency, environmental sustainability, and emissions produced, have been proposed in the past. From a detailed review of the literature, none of the existing studies on vessel schedule recovery have taken explicitly into account the environmental issues and modeled the liner shipping routes, which pass through ECAs and where additional regulations are imposed by IMO. Considering the importance of vessel schedule recovery and growing attention of the community to the environmental issues, this study aims to fill the gaps in the state-of-the-art by making the following contributions:

1.
Present a novel mixed-integer nonlinear programming model for the green vessel schedule recovery problem, which can be applied to the liner shipping routes, passing through ECAs; 2.
Explicitly account for the regulations, imposed by IMO at voyage legs within ECAs; 3.
Consider two common strategies for the vessel schedule recovery (i.e., vessel sailing speed adjustment and port skipping); 4.
Propose an exact solution methodology for solving the formulated mathematical model; 5.
Reveal a set of managerial insights, which will be of interest to liner shipping companies, using the proposed mathematical model and the developed solution approach.

Problem Description
A detailed description of the vessel schedule recovery problem for the liner shipping route, passing through ECAs, and the main modeling assumptions are presented in this section of the manuscript.

Liner Shipping Route Description
A liner shipping route, which connects several ports of call, can be either served by one liner shipping company or multiple liner shipping companies that form an alliance. Throughout this study, the liner shipping route, served by one liner shipping company, will be considered. Let P = 1, . . . , m 1 be the set of ports, which will be visited by vessels of the considered liner shipping company. The vessels, which are dispatched for service of the liner shipping route, sail from port p to p + 1 along voyage leg p. The sequence of ports that have to be visited is assumed to be known. Generally, each port of the liner shipping route is called once. However, in certain cases, a port could be visited more than once during the voyage. Liner shipping companies typically determine the sequence of port visits (which is also referred to as "port rotation" or "port network" in the liner shipping literature) and the port service frequency at the strategic and tactical levels, respectively [3]. The Ports of New York, Norfolk, and Savannah are visited once, while the Ports of Le Havre and Antwerp are visited twice. In order to model multiple visits to the same ports of call during the voyage, two dummy nodes are included in the port rotation graph to represent additional visits to the same ports. The port rotation graph, which will be used for modeling of the liner shipping route, is presented in Figure 2. Two dummy nodes "Le Havre*"and "Antwerp*" are added and placed after the node for

Vessel Service at Ports
The arrival time of vessels at each port along the liner shipping route is planned by the liner shipping company at the tactical level. Let , ∈ (hours) denote the planned arrival time of vessels at port . Generally, the marine container terminal (MCT) operator at each port of call and the liner shipping company coordinate the arrival time of each vessel. The planned arrival time at each port falls within the time window (TW), which was negotiated between the MCT operator and the liner shipping company. The duration of the arrival TW at port is defined based on the start of TW at port ( , ∈ -hours) and the end of TW at port ( , ∈ -hours). Some of the factors, which may affect the planned arrival TW duration at each port of call, include the availability of MCT berthing positions, the availability of handling equipment, the vessel arrival frequency, and others [35,[39][40][41]. The vessels, which arrive at a port of call before the start of negotiated arrival TW, are moored at a dedicated waiting area before service. The planned waiting time of vessels at port + 1 ( +1 , ∈ -hours) can be calculated based on the planned vessel departure time from port ( , ∈ -hours), the planned sailing time between ports and + 1 ( , ∈ -hours), and the start of TW at port + 1 using the following relationships: Note that estimation of the vessel waiting time at the first port of call (i.e., 1 ) also includes the product of the agreed port service frequency ( -hours) and the total number of vessels deployed for service of ports along the liner shipping route ( -vessels), which represents the planned total turnaround time of vessels at the given liner shipping route (see Section 3.4 of the manuscript for more details). Introduction of the latter term is required in order to capture a round voyage of vessels and return at the first port of call [35,39,40].
Furthermore, throughout the vessel schedule design at the tactical level, the liner shipping company negotiates the handling rate with the MCT operator at each port of call of the given liner shipping route. It is assumed that the MCT operator at port offers a set of handling rates = {1, … , 2 }, ∈ to the liner shipping company for service of the arriving vessels. Based on the selected handling rate, the MCT operator has to provide a specific handling productivity ( ℎ , ∈ , ℎ ∈ , measured in twenty-foot equivalent units per hour-TEUs/hour) throughout the vessel service at a given port of call. The container demand ( , ∈ -TEUs) is assumed to be known. The planned handling time of vessels at port ( ℎ , ∈ -hours) can be calculated based on the

Vessel Service at Ports
The arrival time of vessels at each port along the liner shipping route is planned by the liner shipping company at the tactical level. Let , ∈ (hours) denote the planned arrival time of vessels at port . Generally, the marine container terminal (MCT) operator at each port of call and the liner shipping company coordinate the arrival time of each vessel. The planned arrival time at each port falls within the time window (TW), which was negotiated between the MCT operator and the liner shipping company. The duration of the arrival TW at port is defined based on the start of TW at port ( , ∈ -hours) and the end of TW at port ( , ∈ -hours). Some of the factors, which may affect the planned arrival TW duration at each port of call, include the availability of MCT berthing positions, the availability of handling equipment, the vessel arrival frequency, and others [35,[39][40][41]. The vessels, which arrive at a port of call before the start of negotiated arrival TW, are moored at a dedicated waiting area before service. The planned waiting time of vessels at port + 1 ( +1 , ∈ -hours) can be calculated based on the planned vessel departure time from port ( , ∈ -hours), the planned sailing time between ports and + 1 ( , ∈ -hours), and the start of TW at port + 1 using the following relationships: Note that estimation of the vessel waiting time at the first port of call (i.e., 1 ) also includes the product of the agreed port service frequency ( -hours) and the total number of vessels deployed for service of ports along the liner shipping route ( -vessels), which represents the planned total turnaround time of vessels at the given liner shipping route (see Section 3.4 of the manuscript for more details). Introduction of the latter term is required in order to capture a round voyage of vessels and return at the first port of call [35,39,40].
Furthermore, throughout the vessel schedule design at the tactical level, the liner shipping company negotiates the handling rate with the MCT operator at each port of call of the given liner shipping route. It is assumed that the MCT operator at port offers a set of handling rates = {1, … , 2 }, ∈ to the liner shipping company for service of the arriving vessels. Based on the selected handling rate, the MCT operator has to provide a specific handling productivity ( ℎ , ∈ , ℎ ∈ , measured in twenty-foot equivalent units per hour-TEUs/hour) throughout the vessel service at a given port of call. The container demand ( , ∈ -TEUs) is assumed to be known. The planned handling time of vessels at port ( ℎ , ∈ -hours) can be calculated based on the

Vessel Service at Ports
The arrival time of vessels at each port along the liner shipping route is planned by the liner shipping company at the tactical level. Let τ arr p , p ∈ P (hours) denote the planned arrival time of vessels at port p. Generally, the marine container terminal (MCT) operator at each port of call and the liner shipping company coordinate the arrival time of each vessel. The planned arrival time at each port falls within the time window (TW), which was negotiated between the MCT operator and the liner shipping company. The duration of the arrival TW at port p is defined based on the start of TW at port p (τ st p , p ∈ P-hours) and the end of TW at port p (τ f t p , p ∈ P-hours). Some of the factors, which may affect the planned arrival TW duration at each port of call, include the availability of MCT berthing positions, the availability of handling equipment, the vessel arrival frequency, and others [35,[39][40][41]. The vessels, which arrive at a port of call before the start of negotiated arrival TW, are moored at a dedicated waiting area before service. The planned waiting time of vessels at port p + 1 (τ wait p+1 , p ∈ P-hours) can be calculated based on the planned vessel departure time from port p (τ dep p , p ∈ P-hours), the planned sailing time between ports p and p + 1 (τ sail p , p ∈ P-hours), and the start of TW at port p + 1 using the following relationships: Note that estimation of the vessel waiting time at the first port of call (i.e., τ wait 1 ) also includes the product of the agreed port service frequency (ϕ-hours) and the total number of vessels deployed for service of ports along the liner shipping route (V-vessels), which represents the planned total turnaround time of vessels at the given liner shipping route (see Section 3.4 of the manuscript for more details). Introduction of the latter term is required in order to capture a round voyage of vessels and return at the first port of call [35,39,40].
Furthermore, throughout the vessel schedule design at the tactical level, the liner shipping company negotiates the handling rate with the MCT operator at each port of call of the given liner shipping route. It is assumed that the MCT operator at port p offers a set of handling rates H p = 1, . . . , m 2 p , p ∈ P to the liner shipping company for service of the arriving vessels. Based on the selected handling rate, the MCT operator has to provide a specific handling productivity (pr ph , p ∈ P, h ∈ H p , measured in twenty-foot equivalent units per hour-TEUs/hour) throughout the vessel service at a given port of call. The container demand (d port p , p ∈ P-TEUs) is assumed to be known. The planned handling time of vessels at port p (τ hand p , p ∈ P-hours) can be calculated based on the container demand at that port and the requested handling productivity using the following relationship: where: x ph , p ∈ P, h ∈ H p -is the handling rate selection parameter (=1 if handling rate h is requested at port p; = 0 otherwise). The handling cost (c port ph , p ∈ P, h ∈ H p -USD/TEU) will be imposed to the liner shipping company for requesting handling rate h at port p to serve the vessels. Note that disruptions at MCTs (e.g., disruptions due to labor strikes, natural disasters resulting in port infrastructure damages, port congestion) may further cause an increase in the planned waiting and handling times of vessels.
The actual waiting and handling times of vessels will be further referred to as τ wait p , p ∈ (hours) and τ hand p , p ∈ P (hours), respectively. Furthermore, disruptions at preceding ports of call may cause deviations from the planned arrival time of vessels at a given port of call. The actual vessel arrival time will be denoted as τ arr p , p ∈ P (hours). The vessel arrival delay (τ del p , p ∈ P-hours) at port p can be calculated based on the actual and planned vessel arrival times using the following relationship: Delays in the planned arrival time of vessels at ports of call are not desirable, as they may negatively affect both liner shipping and MCT operations. Therefore, an additional delayed vessel arrival time cost (c del p , p ∈ P-USD/hour) will be imposed to the liner shipping company for any deviations from the planned arrival time of vessels at ports of call.

Fuel Consumption
One of the common assumptions used for the vessel schedule design at the tactical level is a homogeneous nature of the vessel fleet, which is deployed to serve the ports of the given liner shipping route [35,39,40,42,43]. Homogeneous vessels have the same major technical specifications (e.g., structure of the main and auxiliary vessel engines, capacity of the main and auxiliary vessel engines, maximum possible sailing speed). However, the assumption of homogeneous vessels may not be fully accurate in real-life liner shipping operations for certain routes, as the vessels of exactly the same types may have differences in their technical characteristics due to vessel repairs that were conducted in the past, age, utilization, and other factors. Furthermore, a large number of different factors may impact the fuel consumption, including but not limited to sailing speed, payload, vessel geometric characteristics, and weather conditions. However, the vessel sailing speed has been identified as the key factor that directly influences the fuel consumption [6,18,35]. Based on the available liner shipping literature [6], the planned fuel consumption per nautical mile (nmi) at voyage leg p ( f p , p ∈ P-tons/nmi) can be calculated using the following relationship: where: α, γ-are the fuel consumption function coefficients; s p , p ∈ P-is the planned sailing speed adopted for vessels at voyage leg p (knots). The fuel consumption coefficients (α, γ) are typically estimated based on the historical data, which contain the information regarding the fuel consumption and sailing speed for the vessel fleet, serving a given liner shipping route. The planned sailing time at voyage leg p can be calculated based on the length of that voyage leg (d leg p , p ∈ P-nmi) and the planned sailing speed using the following relationship: Note that disruptions in sea (e.g., inclement weather, mechanical failure of vessel engines, inexperience/errors of the vessel crew) may further cause a decrease in the planned sailing speed of vessels. Hence, the actual sailing speed at a given voyage leg (s p , p ∈ P-knots) would be lower than the planned sailing speed: s p < s p , p ∈ P. A reduction of the vessel sailing speed will decrease the fuel consumption and may result in the actual fuel cost savings. However, reduction in the vessel sailing speed will increase the actual sailing time at a given voyage leg (τ sail p , p ∈ P-hours) and cause a delayed arrival at the next port of call.

Port Service Frequency
Throughout the vessel schedule design at the tactical level, the liner shipping company needs to ensure that the agreed port service frequency is met at ports of the liner shipping route, which can be achieved using the following relationship [39,40,42,44]: The planned total vessel turnaround time at the liner shipping route (i.e., the time taken by a vessel to sail from the port of origin, visit other ports of the liner shipping route, and return to the port of origin) is estimated based on the right-hand side of Equation (7). The planned total vessel turnaround time is computed based on the following three components: (a) the planned total sailing time of vessels; (b) the planned total waiting time of vessels at ports of call; and (c) the planned total port handling time. However, due to disruptions at ports and in sea, the actual total vessel turnaround time (VTT-hours) can deviate from the planned total vessel turnaround time (VTT-hours). The actual total vessel turnaround time can be calculated using the following relationship:

Vessel Sailing Speed Selection
The planned vessel sailing speed (s p , p ∈ P-knots) is determined by the liner shipping company at each voyage leg of a given liner shipping route at the vessel scheduling stage [42,45]. The sailing speed of vessels at each voyage leg along the liner shipping route is generally affected by a wide range of factors. According to Wang et al. [45], the lower bound on the vessel sailing speed (s min -knots) is selected to minimize the wear of the vessel's engine. The upper bound on the vessel sailing speed (s max -knots) is generally determined by the vessel's engine capacity [7]. As discussed earlier, disruptions in sea may significantly affect the planned vessel sailing speed. The latter will further result in fluctuations of the fuel consumption of vessels and cause late arrivals at the consecutive ports of call. However, the vessel sailing speed adjustment may serve as an efficient mean for recovery of the vessel schedules from disruptive events (as will be discussed more in Section 3.8 of the manuscript).

Container Inventory Cost
Throughout the vessel schedule design, it is assumed that the liner shipping company knows the number of containers (i.e., TEUs) to be transported at each voyage leg of the liner shipping route. Based on the liner shipping literature, the planned container inventory cost (CIC-USD) can be estimated based on the unit cost of shipping a TEU (c inv -USD/TEU/hour), the total number of TEUs transported Energies 2019, 12, 2380 9 of 28 at voyage legs (d sea p , p ∈ P-TEUs), and the planned total vessel sailing time at voyage legs of the liner shipping route [39,40,42] using the following relationship: Disruptions in sea and/or at ports will affect the planned vessel sailing speed at voyage legs. Specifically, at the voyage legs that experience disruptions, the planned vessel sailing time is expected to increase due to reduction in sailing speed (s p < s p , p ∈ P). A vessel schedule recovery can be achieved by increasing the vessel sailing speed at the consecutive voyage legs. The latter action will reduce the planned sailing time. The actual container inventory cost (CIC-USD) can be estimated by replacing the planned vessel sailing time in Equation (9) with the actual sailing time of vessels at voyage legs (τ sail p , p ∈ P-hours) as follows:

Existing International Maritime Organization Regulations on Emissions
In order to improve energy efficiency (i.e., cost-effective utilization of fuel) and environmental sustainability of liner shipping, this study considers the existing IMO environmental policies under "Sulphur oxides (SO x ) and Particulate Matter (PM)-Regulation 14", which limits the emissions of SO x and PM by vessels in the areas that are designated as ECAs [8]. Specifically, regulation 14 requires that the sulfur content in fuel, used by the vessels that sail through ECAs, must not exceed 0.10% m/m [8].
In order to comply with the existing IMO regulations, this study assumes that the liner shipping company will use two types of fuel (namely, marine gas oil (MGO) and heavy fuel oil (HFO)) at the considered liner shipping route. It is assumed that the liner shipping company will use a more costly MGO with the sulfur content of 0.10% only when sailing through ECAs (SC p = 0.10% ∀p ∈ P * , where P * -is a set of voyage legs, passing through ECAs; SC p , p ∈ P-is the sulfur content of fuel used at voyage leg p), while HFO with the sulfur content of 3.50% will be used at other voyage legs of the liner shipping route (SC p = 3.50% ∀p ∈ P 0 , where P 0 -is a set of voyage legs, passing outside ECAs). However, without loss of generality, the proposed methodology will be still applicable after 1 January 2020 (when the vessels that are sailing outside ECAs would be required to use fuel with the sulfur content of 0.50%) by setting an appropriate value for the unit fuel cost.
Note that along with fuel switching the liner shipping company may consider other alternatives to meet the ECA sulfur regulations, which include, but not limited to, installation of scrubbers and utilization of liquefied natural gas (LNG) [34]. In order to comply with the IMO regulations on nitrogen oxide (NO x ) pollutants, produced by vessels in certain ECAs (e.g., the North American area and the United States Caribbean Sea), this study assumes that the diesel engines of vessels, deployed by the liner shipping company, satisfy the Tier III limits. The latter requirement can be relaxed for the liner shipping routes, which pass through ECAs in the Baltic Sea and the North Sea (which have the SO x emission control only).

Vessel Schedule Disruptions and Recovery
Disruptions may occur in sea and/or at ports of call. Disruptions at ports may be caused by labor strikes, natural disasters that result in port infrastructure damages, port congestion, and other events. On the other hand, disruptions in sea may occur due to inclement weather, mechanical failure of vessel engines, inexperience/errors of the vessel crew, and other events. Let y port p = 1, p ∈ P if a disruption occurred at port p (=0 otherwise). Let y sea p = 1, p ∈ P if a disruption occurred at voyage leg p (=0 otherwise). Note that both y port p , p ∈ P and y sea p , p ∈ P are treated as parameters in this study, as the considered vessel schedule recovery problem for the liner shipping route, passing through ECAs, is an operational-level decision problem, and the liner shipping company already knows where Energies 2019, 12, 2380 10 of 28 disruption(s) occurred. The liner shipping company also has the information regarding the expected duration of a disruptive event. The expected duration of a disruptive event at ports will be further referred to as δ port p , p ∈ P (hours), while the expected vessel sailing speed change due to a disruptive event in sea will be denoted as δ sea p , p ∈ P (knots). In order to recover the vessel schedule due to disruptions at the liner shipping route, passing through ECAs, this study considers the following recovery strategies: (1) vessel sailing speed adjustment; and (2) port skipping. The liner shipping company may decide to increase the vessel sailing speed at voyage legs of the liner shipping route in order to recover the vessel schedule and compensate for the delays, incurred as a result of disruptions in sea and/or at ports of call. Denote ∆ sea p , p ∈ P (knots) as the vessel sailing speed adjustment at voyage leg p of the liner shipping route. As discussed earlier, increasing vessel sailing speed will further increase the fuel consumption by the vessels, deployed for service of the given liner shipping route, which will increase the actual total fuel cost. It is assumed that the liner shipping company will not be able to adjust the vessel sailing speed at the voyage leg that experienced a disruptive event. If a disruptive event substantially affected the liner shipping operations at a certain voyage leg of the liner shipping route, even selection of the maximum vessel sailing speed at the consecutive voyage legs may not be sufficient to fully recover the vessel schedule. Figure 3 presents an example of the vessel schedule recovery strategy based on the vessel sailing speed adjustment. The horizontal and vertical axes of the time-space network represent the time within the planning horizon (in days) and the port position, respectively. A vessel, scheduled to visit three ports (the Ports of Le Havre, New York, and Norfolk), experienced a disruption after departure from the Port of Le Havre, which resulted in a late arrival at the Port of New York. In order to recover the schedule and maintain the planned arrival time at the Port of Norfolk, the liner shipping company selected a higher sailing speed at the voyage leg between the Ports of New York and Norfolk.
Energies 2019, 12, x FOR PEER REVIEW 10 of 28 adjustment; and (2) port skipping. The liner shipping company may decide to increase the vessel sailing speed at voyage legs of the liner shipping route in order to recover the vessel schedule and compensate for the delays, incurred as a result of disruptions in sea and/or at ports of call. Denote ∆ , ∈ (knots) as the vessel sailing speed adjustment at voyage leg of the liner shipping route. As discussed earlier, increasing vessel sailing speed will further increase the fuel consumption by the vessels, deployed for service of the given liner shipping route, which will increase the actual total fuel cost. It is assumed that the liner shipping company will not be able to adjust the vessel sailing speed at the voyage leg that experienced a disruptive event. If a disruptive event substantially affected the liner shipping operations at a certain voyage leg of the liner shipping route, even selection of the maximum vessel sailing speed at the consecutive voyage legs may not be sufficient to fully recover the vessel schedule. In order to recover the schedule, the liner shipping company may decide not to skip port and endure the delay. Let , ∈ be the port skipping decision variable (=1 if port is skipped; =0 otherwise). However, if duration of a disruptive event at a given port of call is expected to be significant, the liner shipping company will have to skip that port in order to reduce the additional liner shipping route service costs due to a disruptive event. Figure 4 presents an example of port skipping due to a disruptive event, where the vessel, scheduled to visit the Port of Le Havre, the Port of New York, the Port of Norfolk, the Port of Savannah, and the Port of Antwerp, respectively, skips the Port of Norfolk due to a disruption and visits only the Port of Le Havre, the Port of New York, the Port of Savannah, and the Port of Antwerp. If the liner shipping company decides to skip the port, which experiences a disruptive event, the sequence of visits to the consecutive ports of the port rotation is assumed to remain unchanged. In order to recover the schedule, the liner shipping company may decide not to skip port p and endure the delay. Let x skip p , p ∈ P be the port skipping decision variable (=1 if port p is skipped; =0 otherwise). However, if duration of a disruptive event at a given port of call is expected to be significant, the liner shipping company will have to skip that port in order to reduce the additional liner shipping route service costs due to a disruptive event. Figure 4 presents an example of port skipping due to a disruptive event, where the vessel, scheduled to visit the Port of Le Havre, the Port of New York, the Port of Norfolk, the Port of Savannah, and the Port of Antwerp, respectively, skips the Port of Norfolk due to a disruption and visits only the Port of Le Havre, the Port of New York, the Port of Savannah, and the Port of Antwerp. If the liner shipping company decides to skip the port, which experiences a disruptive event, the sequence of visits to the consecutive ports of the port rotation is assumed to remain unchanged.  If the liner shipping company decides to skip a given port of call due to a disruptive event, it will not be able to unload the import containers and load the export containers at that particular port. A port skipping cost ( , ∈ -USD) is imposed to the liner shipping company for skipping the port affected by a disruption. The penalty may be imposed for every single skip or over a certain time period (e.g., quarterly or annually). The latter aspect is generally negotiated between the MCT operators and the liner shipping company in the corresponding contractual agreements. The port skipping cost is introduced to compensate the MCT operator for the unmet demand, reserved handling equipment, as well as reserved storage yard space.
The vessel schedule recovery strategies considered (i.e., vessel sailing speed adjustment and port skipping) can be implemented by the liner shipping company to recover the vessel schedule from disruptions in sea and/or at ports of call. Depending on the effects of a disruptive event, the liner shipping company may be willing to use both strategies (if they allow reducing the additional liner shipping route service costs due to that disruptive event). However, if the cost of implementing the aforementioned vessel schedule recovery strategies is substantial, the liner shipping company may decide to endure the effects of a disruptive event without application of any vessel schedule recovery actions.

Decisions
The decision problem, examined in this study, can be classified as the operational-level problem in liner shipping and will be referred to as the green vessel schedule recovery problem. The major decisions that need to be addressed in this problem by the liner shipping company are as follows: (1) identify the voyage legs of the liner shipping route, where the vessel sailing speed adjustment will be a favorable alternative to compensate for the delays, caused by a disruptive event; (2) determine whether energy efficient vessel schedule recovery strategies would be favorable to compensate for the delays, caused by a disruptive event (i.e., a limited increase in the vessel sailing speed will not significantly increase the actual fuel consumption and the associated cost, but may not be sufficient to recover the delays); (3) determine whether the vessel sailing speed adjustment would be a favorable alternative at the voyage legs, passing through ECAs; (4) skip a port of call due to a disruption or wait until the port recovers from that disruption and is able to provide service of vessels; (5) select effective vessel schedule recovery strategies, aiming to reduce the vessel arrival delays at ports of the liner shipping route. The overall objective of the liner shipping company is to effectively recover the vessel schedule and minimize the total profit loss as a result of disruptions in sea and/or at ports of call. If the liner shipping company decides to skip a given port of call due to a disruptive event, it will not be able to unload the import containers and load the export containers at that particular port.
A port skipping cost (c skip p , p ∈ P-USD) is imposed to the liner shipping company for skipping the port affected by a disruption. The penalty may be imposed for every single skip or over a certain time period (e.g., quarterly or annually). The latter aspect is generally negotiated between the MCT operators and the liner shipping company in the corresponding contractual agreements. The port skipping cost is introduced to compensate the MCT operator for the unmet demand, reserved handling equipment, as well as reserved storage yard space.
The vessel schedule recovery strategies considered (i.e., vessel sailing speed adjustment and port skipping) can be implemented by the liner shipping company to recover the vessel schedule from disruptions in sea and/or at ports of call. Depending on the effects of a disruptive event, the liner shipping company may be willing to use both strategies (if they allow reducing the additional liner shipping route service costs due to that disruptive event). However, if the cost of implementing the aforementioned vessel schedule recovery strategies is substantial, the liner shipping company may decide to endure the effects of a disruptive event without application of any vessel schedule recovery actions.

Decisions
The decision problem, examined in this study, can be classified as the operational-level problem in liner shipping and will be referred to as the green vessel schedule recovery problem. The major decisions that need to be addressed in this problem by the liner shipping company are as follows: (1) identify the voyage legs of the liner shipping route, where the vessel sailing speed adjustment will be a favorable alternative to compensate for the delays, caused by a disruptive event; (2) determine whether energy efficient vessel schedule recovery strategies would be favorable to compensate for the delays, caused by a disruptive event (i.e., a limited increase in the vessel sailing speed will not significantly increase the actual fuel consumption and the associated cost, but may not be sufficient to recover the delays); (3) determine whether the vessel sailing speed adjustment would be a favorable alternative at the voyage legs, passing through ECAs; (4) skip a port of call due to a disruption or wait until the port recovers from that disruption and is able to provide service of vessels; (5) select effective vessel schedule recovery strategies, aiming to reduce the vessel arrival delays at ports of the liner shipping route. The overall objective of the liner shipping company is to effectively recover the vessel schedule and minimize the total profit loss as a result of disruptions in sea and/or at ports of call.

Mathematical Model
This section of the manuscript presents a mixed-integer nonlinear programming model for the green vessel schedule recovery problem (the model will be referred to as GVSRP), which takes into account the existing IMO regulations at the voyage legs that pass through ECAs. The GVSRP mathematical model guarantees that low-sulfur fuel is used when sailing through ECAs. The nomenclature, which will be adopted throughout the manuscript, is also presented in this section.
Nomenclature Sets P = 1, . . . , m 1 set of ports to be visited (ports) H p = 1, . . . , m 2 p , p ∈ P set of available handling rates at port p (handling rates)

Decision Variables
∆ sea p ∈ R ∀p ∈ P vessel sailing speed adjustment at voyage leg p (knots) Auxiliary Variables s p ∈ R + ∀p ∈ P actual sailing speed of vessels at voyage leg p (knots) τ arr p ∈ R + ∀p ∈ P actual arrival time of vessels at port p (hours) τ wait p ∈ R + ∀p ∈ P actual waiting time of vessels at port p (hours) τ hand p ∈ R + ∀p ∈ P actual handling time of vessels at port p (hours) τ dep p ∈ R + ∀p ∈ P actual departure time of vessels from port p (hours) τ sail p ∈ R + ∀p ∈ P actual sailing time of vessels at voyage leg p that connects ports p and p + 1 (hours) f p ∈ R + ∀p ∈ P actual fuel consumption at voyage leg p (tons/nmi) τ del p ∈ R + ∀p ∈ P vessel arrival delay at port p (hours) VTT ∈ R + actual total vessel turnaround time (hours) REV ∈ R + actual total revenue to be generated by the liner shipping company (USD) PHC ∈ R + actual total port handling cost (USD) LAC ∈ R + actual total late vessel arrival cost (USD) FCC ∈ R + actual total fuel cost (USD) CIC ∈ R + actual total container inventory cost (USD) TP ∈ R + actual total profit to be generated by the liner shipping company (USD) Parameters m 1 ∈ N number of ports to be visited (ports) m 2 p ∈ N ∀p ∈ P number of available handling rates at port p (handling rates) τ arr p ∈ R + ∀p ∈ P planned arrival time of vessels at port p (hours) τ st p ∈ R + ∀p ∈ P start of TW at port (hours) τ f t p ∈ R + ∀p ∈ P end of TW at port p (hours) ϕ ∈ N agreed port service frequency (hours) V ∈ N total number of vessels deployed for service of ports along the liner shipping route (vessels) d port p ∈ R + ∀p ∈ P container demand at port p (TEUs) pr ph ∈ R + ∀p ∈ P, h ∈ H p handling productivity at port p under handling rate h (TEUs/hour) x ph ∈ B ∀p ∈ P, h ∈ H p =1 if handling rate h is requested at port p (=0 otherwise) d leg p ∈ R + ∀p ∈ P length of voyage leg p (nmi) α, γ ∈ R + fuel consumption coefficients s p ∈ R + ∀p ∈ P planned vessel sailing speed at port p (knots) s min ∈ R + minimum vessel sailing speed (knots) s max ∈ R + maximum vessel sailing speed (knots) d sea p ∈ R + ∀p ∈ P number of containers transported at voyage leg p (TEUs) y port p ∈ B ∀p ∈ P =1 if a disruptive event occurred at port p (=0 otherwise) y sea p ∈ B ∀p ∈ P =1 if a disruptive event occurred at voyage leg p (=0 otherwise) δ port p ∈ R + ∀p ∈ P expected duration of a disruptive event at port p (hours) δ sea p ∈ R ∀p ∈ P expected vessel sailing speed change due to a disruptive event at voyage leg p (knots) c cargo p ∈ R + ∀p ∈ P freight rate for shipping a TEU from port p to port p + 1 (USD/TEU) c port ph ∈ R + ∀p ∈ P, h ∈ H p handling cost at port p under handling rate h (USD/TEU) c del p ∈ R + ∀p ∈ P delayed vessel arrival cost at port p (USD/hour) c f p ∈ R + ∀p ∈ P unit fuel cost when sailing at voyage leg p (USD/ton) c oper ∈ R + vessel operational cost (USD/hour) c inv ∈ R + unit container inventory cost (USD/TEU/hour) c skip p ∈ R + ∀p ∈ P cost for skipping port p due to a disruptive event (USD) VOC ∈ R + planned total vessel operational cost (USD) TP ∈ R + planned total profit to be generated by the liner shipping company (USD) Green Vessel Schedule Recovery Problem (GVSRP): Subject to: s p ≥ s min + δ sea p y sea p ∀p ∈ P (14) τ arr p+1 = τ dep p + τ sail p ∀p ∈ P, p < m 1 (17) In the proposed GVSRP mathematical model, the objective function (11) aims to minimize the total profit loss as a result of disruptions in sea and/or at ports of call for the liner shipping route, passing through ECAs. Constraint set (12) computes the vessel sailing speed for the recovered vessel schedule at each voyage leg of the given liner shipping route, passing through ECAs. It is assumed that the vessel sailing speed adjustment strategy cannot be used for the vessel schedule recovery by the liner shipping company at the voyage legs, where a disruption occurs in sea (i.e., the vessel sailing speed cannot be increased at the voyage legs, experiencing disruptions, to reduce the delays). Constraint sets (13) and (14) (24) guarantees that the liner shipping company can skip a port of the given liner shipping route, passing through ECAs, only if that particular port experienced a disruption. Constraint set (25) calculates the actual total vessel turnaround time. Constraint sets (26)-(32) estimate the individual cost components of the objective function (11) for the GVSRP mathematical model, which include the following: (i) the actual total revenue; (ii) the actual total port handling cost; (iii) the actual total late vessel arrival cost; (iv) the actual total fuel cost; (v) the actual total container inventory cost; (vi) the total vessel operational cost; and (vii) the actual total profit for the recovered vessel schedule. Note that the unit fuel cost c f p , p ∈ P (USD/ton) varies when the vessel is sailing at voyage legs within and outside ECAs. Based on the existing IMO regulations, the liner shipping company is mandated to use a more expensive low-sulfur fuel when sailing through ECAs (c f p = c MGO ∀p ∈ P * , where c MGO -is the unit cost of MGO). On the other hand, liner shipping company will switch to cheaper HFO fuel at the voyage legs outside ECAs (c f p = c HFO ∀p ∈ P 0 , where c HFO p -is the unit cost of HFO).

Solution Methodology
The proposed GVSRP mathematical model is nonlinear due to constraint set (15), which was used for estimating the actual vessel sailing time at each voyage leg, and constraint set (16), which was used for estimating the actual fuel consumption at each voyage leg. In order to linearize constraint set (15), the actual vessel sailing speed s p , p ∈ P was replaced with its reciprocal v p = 1 s p ∀p ∈ P. Similarly, the planned vessel sailing speed s p , p ∈ P will be replaced with its reciprocal . v p = 1 s p ∀p ∈ P. Also, the expected vessel sailing speed change due to a disruptive event at a given voyage leg (δ sea p , p ∈ P-measured in knots) should be adjusted accordingly. Denote . δ sea p , p ∈ P as the expected vessel sailing speed change due to a disruptive event at a given voyage, measured in knots −1 (i.e., converted to the reciprocal of the vessel sailing speed change for consistency). Denote . ∆ sea p , p ∈ P as the vessel sailing speed adjustment at a given voyage, measured in knots −1 .
Let the fuel consumption function, which was computed based on the reciprocal of the vessel sailing speed v p , p ∈ P, be represented by FC p , p ∈ (tons/nmi). Throughout this study, the piecewise linear static secant approximation will be adopted to linearize the GVSRP mathematical model. For static piecewise linear secant approximations, the number of secant lines is predefined [40,45,46]. The piecewise function FC pr , p ∈ P, r ∈ K approximates the nonlinear function FC p , p ∈ P using a defined number of linear segments, where K = {1, . . . , r} represents a set of linear segments in the piecewise function. Let b pk = 1 if segment k is selected to approximate the fuel consumption function at voyage leg p of the liner shipping route (=0 otherwise). Denote st k , k ∈ K and ed k , k ∈ K as the vessel speed reciprocal values at the start of linear segment k and at the end of linear segment k, respectively. Let SL k , k ∈ K and IN k , k ∈ K be the slope of linear segment k and the intercept of linear segment k, respectively. Denote M 1 and M 2 as sufficiently large positive numbers. The mixed-integer nonlinear GVSRP mathematical model can be further reduced to a mixed-integer linear mathematical model (the model will be referred to as GVSRPL) as follows.
Linearized Green Vessel Schedule Recovery Problem (GVSRPL): Subject to: Constraint sets (17)- (28) and (30) In the proposed GVSRPL mathematical model, the objective function (33) aims to minimize the total profit loss as a result of disruptions in sea and/or at ports of call for the liner shipping route, passing through ECAs. Constraint set (34) computes the reciprocal of the vessel sailing speed for the recovered vessel schedule at each voyage leg of the given liner shipping route, passing through ECAs. Constraint sets (35) and (36) define the limits for the reciprocal of the vessel sailing speed at each voyage leg of the liner shipping route, passing through ECAs. Constraint set (37) calculates the vessel sailing time based on the reciprocal of the vessel sailing speed for the recovered vessel schedule at each voyage leg of the given liner shipping route, passing through ECAs. Constraint set (38) ensures that only one linear segment will be selected to approximate the fuel consumption function for the recovered vessel schedule at each voyage leg of the liner shipping route, passing through ECAs. Constraint sets (39) and (40) define the ranges for the reciprocal of the vessel sailing speed, when a given linear segment k is selected to approximate the fuel consumption function, at each voyage leg of the liner shipping route, passing through ECAs. Constraint set (41) estimates the fuel consumption based on the reciprocal of the vessel sailing speed for the recovered vessel schedule at each voyage leg of the given liner shipping route, passing through ECAs. Constraint set (42) estimates the actual total fuel cost.
Based on the preliminary numerical experiments, it was found that the GVSRPL mathematical model can be solved to the global optimality using CPLEX for the realistic-size liner shipping routes in a reasonable computational time. Therefore, development of the approximate solution approaches is not required for the GVSRPL mathematical model, and CPLEX will be further used as a solution approach. Note that increasing the number of linear segments in the piecewise approximation for the fuel consumption function typically enhances its accuracy but may cause an increase in the computational time, required to solve the GVSRPL mathematical model (due to increasing number of variables in the GVSRPL mathematical model). The latter tradeoff will be further investigated throughout the numerical experiments (see Appendix A that accompanies the manuscript for more details).

Computational Experiments
This section of the manuscript presents the numerical experiments, which were performed to exhibit how the proposed GVSRPL mathematical model could be used for decision making by liner shipping companies for the liner shipping routes, passing through ECAs, in case of disruptive events. All the numerical experiments in this study were performed using a Dell Inspiron AMD FX-9800P processor with 16 GB of RAM. The MATLAB 2016a software [47] was used to develop a procedure for generating the piecewise secant approximations, which were further applied to linearize the fuel consumption function. The General Algebraic Modeling System [48] was used to encode the GVSRPL mathematical model and solve it to the global optimality with CPLEX. The following aspects are further discussed in this section: (1) input data generation; and (2) managerial insights.

Input Data Generation
This study considered the Asia-North Europe LL5 liner shipping route, passing through ECAs. The liner shipping route is served by OOCL liner shipping company [41] and is presented in Figure 5. Note that Figure 5 was developed based on the data, provided by OOCL [41]. The route connects Asia, Arabian Sea, Red Sea, North Africa, Malta, and North Europe. As illustrated in Figure 5, the Asia-North Europe LL5 liner shipping route passes through the North Sea and the English Channel, which are designated as ECAs with the SO x emission control. The port rotation for the Asia-North Europe LL5 liner shipping route includes 14 ports of call, where the Port of Le Havre (France) is visited twice. Each one of the ports is scheduled to be visited by vessels on a weekly basis. The distances between the consecutive ports, measured in nautical miles, were obtained from the world seaports catalogue [49] and are presented in parenthesis: 1 In order to conduct the computational experiments for this study, the data from the liner shipping literature and the MCT operations literature [6][7][8][9][10][11][37][38][39][40][50][51][52][53][54][55] were used to generate the parameter values for the GVSRPL mathematical model. The adopted parameter values are presented in Table 1. The start of TW at each port of the port rotation was estimated based on a mathematical relationship between the start of TW at the preceding port, the upper and lower bounds of the vessel sailing speed, and the length of a voyage leg between the consecutive ports as ∀ ∈ (hours). The notation will be further adopted for the values that are drawn from a set of uniformly distributed pseudorandom numbers within a specified range. The planned handling time of vessels at port was calculated as follows: ℎ = ∑ ( ℎ ) ℎ∈ ℎ ∀ ∈ (hours). The negotiated handling rate (i.e., the value of parameter ℎ , ∈ , ℎ ∈ ) was selected randomly from the available handling rates at each port. The planned arrival time of vessels at port + 1 was set based on the following relationship: +1 = + ℎ + ∀ ∈ (hours).  In order to conduct the computational experiments for this study, the data from the liner shipping literature and the MCT operations literature [6][7][8][9][10][11][37][38][39][40][50][51][52][53][54][55] were used to generate the parameter values for the GVSRPL mathematical model. The adopted parameter values are presented in Table 1. The start of TW at each port of the port rotation was estimated based on a mathematical relationship between the start of TW at the preceding port, the upper and lower bounds of the vessel sailing speed, and the length of a voyage leg between the consecutive ports as follows: τ st p+1 = τ st p + d leg p U[s min ; s max ] ∀p ∈ P (hours). The notation U will be further adopted for the values that are drawn from a set of uniformly distributed pseudorandom numbers within a specified range. The planned handling time of vessels at port p was calculated as follows: τ hand p = h∈H p d port p pr ph x ph ∀p ∈ P (hours). The negotiated handling rate (i.e., the value of parameter x ph , p ∈ P, h ∈ H p ) was selected randomly from the available handling rates at each port. The planned arrival time of vessels at port p + 1 was set based on the following relationship: τ arr p+1 = τ arr p + τ hand p + d leg p s p ∀p ∈ P (hours). Different cases of disruption occurrence in sea and at ports of call will be modeled throughout the numerical experiments for the considered liner shipping route (see Section 6.2 of the manuscript for more details). Note that the unit cost of HFO fuel was originally assumed to be c HFO = 200 USD/ton, while the unit cost of MGO fuel was set to c MGO = 500 USD/ton. However, an additional sensitivity analysis will be conducted for the unit cost of HFO and the unit cost of MGO (details will be presented in the following sections of the manuscript). The port skipping inconvenience coefficient was originally set to µ = 1.10. However, an additional sensitivity analysis will be conducted for the port skipping inconvenience coefficient (details will be presented in the following sections of the manuscript). (1) Case 1: no disruptions in sea and at ports This is an ideal condition, as there are no disruptions in sea, which may significantly affect the planned vessel sailing speed at voyage legs of the considered liner shipping route. Also, there are no disruptions at ports of call, which may cause delays in vessel service. The vessels visit each port of the port rotation as planned.
(2) Case 2: disruptions in sea and at ports within ECAs There is a strike at the Port of Antwerp (Belgium) with an expected duration of δ port 4 = 50 hours. Also, there is a disruption in sea due to inclement weather at voyage leg "4", connecting the Port of Antwerp (Belgium) and the Port of Le Havre (France), which is expected to reduce the planned vessel sailing speed at voyage leg "4" by δ sea 4 = −4.0 knots. Note that the Port of Antwerp (Belgium) and the Port of Le Havre (France) are located within ECAs, where the liner shipping company must use low-sulfur MGO fuel only. Note that the ports, connected by voyage legs "1" to "4", are within the areas that are designated as ECAs, while the other ports are outside the ECAs. Table 2. The considered unit fuel cost values for the generated problem instances.

Instance
Unit HFO Cost (USD/Ton) Unit MGO Cost (USD /Ton)   1  200  500  2  250  550  3  300  600  4  350  650  5  400  700  6  450  750  7  500  800  8  550  850  9  600  900  10 650 950 Sensitivity of the GVSRPL mathematical model to changes in the unit fuel costs (both HFO and MGO) is further presented in this section of the manuscript with a primary emphasis on the following two aspects: (a) port skipping decisions; and (b) vessel sailing speed and sailing time decisions.

Port Skipping Decisions
The port skipping decisions were retrieved for each one of the considered problem instances and each one of the three disruption cases using the GVSRPL mathematical model. Throughout the numerical experiments, it was found that port skipping decisions were not affected with changes in the unit cost of HFO and MGO. The primary factor, influencing the port skipping decisions, was determined to be the disruption type. The port skipping decisions for all the considered disruption cases are presented in Figure 6. In disruption case 1 (no disruptions occur in sea and at ports within and outside ECAs), it was observed that the vessels did not skip any port of call-see Figure 6a. Hence, the liner shipping company tended to adhere to the planned vessel schedule even after increasing the unit cost of HFO and MGO. On the other hand, in disruption case 2 (disruptions occur in sea and at ports within ECAs), skipping the Port of Antwerp (Belgium) was found to be the most favorable decision for the liner shipping company in order to recover the vessel schedule-see Figure 6b. As for disruption case 3 (disruptions occur in sea and at ports within and outside ECAs), skipping the Port of Antwerp (Belgium) and the Port of Le Havre (France) was found to be the most favorable decision for the liner shipping company in order to recover the vessel schedule-see Figure 6c. Therefore, the results from conducted numerical experiments demonstrate that an increasing number of disruptive events along the given liner shipping route caused significant changes in the vessel schedule and required the liner shipping company to implement more radical recovery strategies (i.e., port skipping rather than vessel sailing speed adjustment). Moreover, increasing unit fuel cost further imposes limitations on implementation of the vessel sailing speed adjustment strategy, since increasing unit fuel cost and increasing vessel sailing speed are expected to significantly increase the actual total fuel cost. The latter will reduce the actual total profit to be generated by the liner shipping company.

Vessel Sailing Speed and Sailing Time Decisions
The average vessel sailing speed adjustment was calculated for each one of the generated problem instances and disruption cases 2 and 3 using the GVSRPL mathematical model, and the results are presented in Figure 7. Furthermore, the average vessel sailing speed was estimated for each one of the considered problem instances and disruption cases 1, 2, and 3 using the GVSRPL mathematical model, and the results are presented in Figure 8. Note that the vessel sailing speed adjustment was not computed for disruption case 1, as no disruptions in sea and at ports were modeled for the latter case. Although changes in the average vessel sailing speed were observed after increasing the unit cost of HFO and MGO even for disruption case 1 (see Figure 8), the vessel sailing speed adjustment was not used as the vessel schedule recovery strategy. Along with the average vessel sailing speed adjustment and the average vessel sailing speed, the total vessel sailing time was retrieved for each one of the generated problem instances and disruption cases 1, 2, and 3 using the GVSRPL mathematical model, and the results are presented in Figure 9. Therefore, the results from conducted numerical experiments demonstrate that an increasing number of disruptive events along the given liner shipping route caused significant changes in the vessel schedule and required the liner shipping company to implement more radical recovery strategies (i.e., port skipping rather than vessel sailing speed adjustment). Moreover, increasing unit fuel cost further imposes limitations on implementation of the vessel sailing speed adjustment strategy, since increasing unit fuel cost and increasing vessel sailing speed are expected to significantly increase the actual total fuel cost. The latter will reduce the actual total profit to be generated by the liner shipping company.

Vessel Sailing Speed and Sailing Time Decisions
The average vessel sailing speed adjustment was calculated for each one of the generated problem instances and disruption cases 2 and 3 using the GVSRPL mathematical model, and the results are presented in Figure 7. Furthermore, the average vessel sailing speed was estimated for each one of the considered problem instances and disruption cases 1, 2, and 3 using the GVSRPL mathematical model, and the results are presented in Figure 8. Note that the vessel sailing speed adjustment was not computed for disruption case 1, as no disruptions in sea and at ports were modeled for the latter case. Although changes in the average vessel sailing speed were observed after increasing the unit cost of HFO and MGO even for disruption case 1 (see Figure 8), the vessel sailing speed adjustment was not used as the vessel schedule recovery strategy. Along with the average vessel sailing speed adjustment and the average vessel sailing speed, the total vessel sailing time was retrieved for each one of the generated problem instances and disruption cases 1, 2, and 3 using the GVSRPL mathematical model, and the results are presented in Figure 9.
The numerical experiments demonstrate that the average vessel sailing speed was generally lower for disruption cases 2 and 3 as compared to disruption case 1, which further increased the total vessel sailing time. The latter pattern can be supported by the fact that disruptive events at certain voyage legs caused a substantial reduction in the vessel sailing speed (i.e., the vessel sailing speed values under disruption cases 2 and 3). Even application of the vessel sailing speed adjustment strategy did not allow the liner shipping company reaching the planned vessel sailing speed (i.e., the vessel sailing speed values under disruption case 1). Furthermore, an increase in the unit cost of HFO and MGO typically reduced the average vessel sailing speed for all the considered disruption cases. Such reduction can be justified the fact that the liner shipping company aimed to decrease the actual fuel consumption and the associated cost. Moreover, an increase in the unit cost of HFO and MGO also resulted in reduction of the average vessel sailing speed adjustment for the recovered vessel schedules under disruption cases 2 and 3.       However, some fluctuations in the average vessel sailing speed can be noticed from one problem instance to another, which can be explained by the fact that the arithmetic average vessel sailing speed does not consider the length of voyage legs. On the other hand, more clear patterns were recorded for the total vessel sailing time, as it accounts not only for the vessel sailing speed but also for the length of voyage legs: = ∀ ∈ (hours). Higher total vessel sailing time values were generally obtained for the case, when disruptions occur in sea and at ports within and outside ECAs (i.e., disruption case 3), as compared to the case, when disruptions occur in sea and at ports within ECAs only (i.e., disruption case 2). Hence, an increasing number of disruptive events along the given liner shipping route caused significant changes in the vessel schedule and increased the total vessel sailing time. The scope of numerical experiments also included a detailed assessment of changing port skipping cost effects on the recovered vessel schedules. More information regarding the latter analysis is provided in Appendix B that accompanies this manuscript.

Concluding Remarks
Liner shipping has maintained a steady growth over the years. However, disruptions in sea and at ports of call significantly affect the planned vessel schedules and result in negative externalities for liner shipping companies, including monetary losses. Effective vessel schedule recovery is critical to reduce monetary losses. The vessel schedule recovery problem becomes more complex when disruptions occur at the liner shipping routes, passing through emission control areas (ECAs), as liner shipping companies must comply with the established International Maritime Organization (IMO) regulations. This study presented a novel mixed-integer nonlinear mathematical model for the green vessel schedule recovery problem at the liner shipping route, which passes through ECAs and where additional regulations are imposed by IMO on the fuel sulfur content. The objective of the proposed model aimed to minimize the total profit loss, endured by a given liner shipping company due to disruptions in the planned operations. Two of the most common vessel recovery strategies were considered, including vessel sailing speed adjustment and port skipping. The nonlinear model was linearized and solved to the global optimality by applying CPLEX.
A number of computational experiments were conducted for the Asia-North Europe LL5 liner shipping route, served by OOCL liner shipping company and passing through ECAs. It was found that an increasing number of disruptive events along the given liner shipping route generally resulted in significant vessel schedule changes and required the liner shipping company to implement more radical recovery strategies (i.e., port skipping became preferential over vessel sailing speed adjustment). Furthermore, increasing unit fuel cost imposed limitations on implementation of the vessel sailing speed adjustment strategy, since increasing unit fuel cost and generally obtained for the case, when disruptions occur in sea and at ports within and outside ECAs (i.e., disruption case 3), as compared to the case, when disruptions occur in sea and at ports within ECAs only (i.e., disruption case 2). Hence, an increasing number of disruptive events along the given liner shipping route caused significant changes in the vessel schedule and increased the total vessel sailing time. The scope of numerical experiments also included a detailed assessment of changing port skipping cost effects on the recovered vessel schedules. More information regarding the latter analysis is provided in Appendix B that accompanies this manuscript.

Concluding Remarks
Liner shipping has maintained a steady growth over the years. However, disruptions in sea and at ports of call significantly affect the planned vessel schedules and result in negative externalities for liner shipping companies, including monetary losses. Effective vessel schedule recovery is critical to reduce monetary losses. The vessel schedule recovery problem becomes more complex when disruptions occur at the liner shipping routes, passing through emission control areas (ECAs), as liner shipping companies must comply with the established International Maritime Organization (IMO) regulations. This study presented a novel mixed-integer nonlinear mathematical model for the green vessel schedule recovery problem at the liner shipping route, which passes through ECAs and where additional regulations are imposed by IMO on the fuel sulfur content. The objective of the proposed model aimed to minimize the total profit loss, endured by a given liner shipping company due to disruptions in the planned operations. Two of the most common vessel recovery strategies were considered, including vessel sailing speed adjustment and port skipping. The nonlinear model was linearized and solved to the global optimality by applying CPLEX.
A number of computational experiments were conducted for the Asia-North Europe LL5 liner shipping route, served by OOCL liner shipping company and passing through ECAs. It was found that an increasing number of disruptive events along the given liner shipping route generally resulted in significant vessel schedule changes and required the liner shipping company to implement more radical recovery strategies (i.e., port skipping became preferential over vessel sailing speed adjustment). Furthermore, increasing unit fuel cost imposed limitations on implementation of the vessel sailing speed adjustment strategy, since increasing unit fuel cost and increasing vessel sailing speed would substantially increase the actual total fuel cost. The findings also indicate that the port skipping decisions could be significantly affected with the port skipping cost for the cases, when disruptions occur in sea and at ports within and outside ECAs. The computational experiments showcase that the proposed mathematical model and the developed solution methodology can assist liner shipping companies with efficient vessel schedule recovery, minimize the monetary losses due to disruptions in vessel schedules, and improve energy efficiency as well as environmental sustainability.
This study can be extended in several dimensions that include, but are not limited to the following: (1) apply the proposed methodology for the liner shipping routes, passing through ECAs with SO x , NO x , and PM emission control; (2) modeling other types of emissions (e.g., nitrogen oxide-NO x , carbon monoxide-CO, non-methane volatile organic compounds-VO); (3) evaluating potential re-routing options for the vessels that sail through ECAs; (4) modeling stochastic disruptive events in sea and ports of call, where the expected vessel sailing speed change in sea and the expected disruption duration at ports of call are uncertain; (5) consider alternative vessel recovery strategies in the presented mathematical model (e.g., port skipping with cargo diversion to other ports, port swapping, handling rate adjustment); (6) collect the operational data from liner shipping companies in order to derive a more comprehensive fuel consumption function, which accounts for sailing speed, payload, and vessel geometric characteristics; and (7) evaluation of alternative methods for linearizing the presented mathematical model (e.g., discretization, outer approximations, second order cone programming).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
A set of computational experiments were conducted to investigate the relationship between the number of linear segments in the piecewise function, accuracy of the fuel consumption function approximation, and the computational time, required to solve the GVSRPL mathematical model. A total of 19 scenarios were developed by varying the number of linear segments in the piecewise function from 2 to 100 (see Table A1). The examples of piecewise approximations for the scenarios with 1 segment, 2 segments, 5 segments, and 10 segments are illustrated in Figure A1  . It can be observed that the accuracy of approximation for FC(v) generally increases with increasing number of linear segments. The GVSRPL mathematical model was solved for each scenario with a specified number of linear segments. It was assumed that no disruptions occurred in sea and at ports throughout the analysis at this stage. The results, retrieved from the numerical experiments (see Table A1), include: (i) the scenario ID; (ii) the number of linear segments-|K|; (iii) the GVSRPL objective function value-Z; (iv) the value of nonlinear objective function at the solution, which was provided by GVSRPL-Z * ; (v) the objective gap ∇ = Z * −Z Z * ; and (vi) the mean value of CPU time over 10 replications. A substantial increase in the GVSRPL computational time was observed when the number of linear segments was greater than 30 without any significant change in the objective gap. Therefore, a piecewise function with 30 linear segments will be further adopted throughout this study for the approximation of the fuel consumption function and analysis of the managerial insights.  The GVSRPL mathematical model was solved for each scenario with a specified number of linear segments. It was assumed that no disruptions occurred in sea and at ports throughout the analysis at this stage. The results, retrieved from the numerical experiments (see Table A1), include: (i) the scenario ID; (ii) the number of linear segments-| |; (iii) the GVSRPL objective function value-; (iv) the value of nonlinear objective function at the solution, which was provided by GVSRPL- * ; (v) the objective gap = | * − * |; and (vi) the mean value of CPU time over 10 replications. A substantial increase in the GVSRPL computational time was observed when the number of linear segments was greater than 30 without any significant change in the objective gap. Therefore, a The port skipping decisions for all the generated problem instances of disruption case 2 and the Port of Antwerp (Belgium) are presented in Figure A2. It was observed that the Port of Antwerp (Belgium) was skipped by the liner shipping company for all the considered problem instances, when the port skipping inconvenience coefficient was set to µ = [1.00]-see Figure A2a. Such finding can be supported by the fact that skipping the Port of Antwerp (Belgium) was the most favorable decision for the liner shipping company in order to recover the vessel schedule (i.e., the port skipping cost was lower as compared to the other additional route service costs, which could be endured as a result of disruptions).
piecewise function with 30 linear segments will be further adopted throughout this study for the approximation of the fuel consumption function and analysis of the managerial insights.

Appendix B
This section of the manuscript discusses the effects of changing port skipping cost on the recovered vessel schedules. Throughout the numerical experiments, the value of the port skipping inconvenience coefficient ( ) was altered from = [1.00] to = [10.00] with an increment of [1.00]. The GVSRPL mathematical model was solved for all the considered problem instances and disruption cases 2 and 3. The results from computational experiments are presented next for disruption cases 2 and 3. The port skipping decisions for all the generated problem instances of disruption case 2 and the Port of Antwerp (Belgium) are presented in Figure A2. It was observed that the Port of Antwerp (Belgium) was skipped by the liner shipping company for all the considered problem instances, when the port skipping inconvenience coefficient was set to = [1.00]-see Figure A2a. Such finding can be supported by the fact that skipping the Port of Antwerp (Belgium) was the most favorable decision for the liner shipping company in order to recover the vessel schedule (i.e., the port skipping cost was lower as compared to the other additional route service costs, which could be endured as a result of disruptions). On the other hand, an increase in the port skipping inconvenience coefficient to = [2.00] could significantly increase the port skipping cost to be endured by the liner shipping company. Therefore, the liner shipping company did not skip the Port of Antwerp (Belgium) for the majority of the considered problem instances and the port skipping inconvenience coefficient of = [2.00]see Figure A2b. The Port of Antwerp (Belgium) was skipped only for problem instance 10 with the highest unit cost of HFO and the highest unit cost of MGO. The latter can be explained by the fact that the liner shipping company was not able to make substantial changes in the vessel sailing speed due to increasing unit cost of both fuel types, and port skipping was the most favorable vessel schedule recovery strategy for problem instance 10. Based on the conducted numerical experiments, it can be concluded that the port skipping decisions could be significantly affected with the port skipping cost for the cases, when disruptions occur in sea and at ports within ECAs. On the other hand, an increase in the port skipping inconvenience coefficient to µ = [2.00] could significantly increase the port skipping cost to be endured by the liner shipping company. Therefore, the liner shipping company did not skip the Port of Antwerp (Belgium) for the majority of the considered problem instances and the port skipping inconvenience coefficient of µ = [2.00]-see Figure A2b. The Port of Antwerp (Belgium) was skipped only for problem instance 10 with the highest unit cost of HFO and the highest unit cost of MGO. The latter can be explained by the fact that the liner shipping company was not able to make substantial changes in the vessel sailing speed due to increasing unit cost of both fuel types, and port skipping was the most favorable vessel schedule recovery strategy for problem instance 10. Based on the conducted numerical experiments, it can be concluded that the port skipping decisions could be significantly affected with the port skipping cost for the cases, when disruptions occur in sea and at ports within ECAs. The port skipping decisions for all the generated problem instances of disruption case 3 are presented in Figure A3 for the following ports: (1) the Port of Antwerp (Belgium); (2) the Port of Le Havre (France); and (3) the Port of Marsaxlokk (Malta). It was observed that the Port of Antwerp (Belgium) and the Port of Le Havre (France), which are both located within ECAs, were skipped by the liner shipping company for all the considered problem instances, when the port skipping inconvenience coefficient was set to µ = [1.0] ÷ [7.0]-see Figure A3a. Such finding can be supported by the fact that skipping the Port of Antwerp (Belgium) and the Port of Le Havre (France) was the most favorable decision for the liner shipping company in order to recover the vessel schedule (i.e., the port skipping costs were lower as compared to the other additional route service costs, which could be endured as a result of disruptions).
On the other hand, an increase in the port skipping inconvenience coefficient to µ = [8.0] ÷ [10.0] could significantly increase the port skipping cost to be endured by the liner shipping company. Therefore, the liner shipping company did not skip the Port of Antwerp (Belgium) for all the considered problem instances and the port skipping inconvenience coefficient of µ = [8.0] ÷ [10.0]-see Figure A3b. Changes in the unit cost of HFO and MGO from one problem instance to another did not influence the port skipping decisions. However, the Port of Le Havre (France) was still skipped by the liner shipping company even after increasing the port skipping inconvenience coefficient to µ = [8.0] ÷ [10.0] in order to effectively recover the vessel schedule. The port skipping decisions for all the generated problem instances of disruption case 3 are presented in Figure A3 for the following ports: (1) the Port of Antwerp (Belgium); (2) the Port of Le Havre (France); and (3) the Port of Marsaxlokk (Malta). It was observed that the Port of Antwerp (Belgium) and the Port of Le Havre (France), which are both located within ECAs, were skipped by the liner shipping company for all the considered problem instances, when the port skipping inconvenience coefficient was set to = [1.0] ÷ [7.0] -see Figure A3a. Such finding can be supported by the fact that skipping the Port of Antwerp (Belgium) and the Port of Le Havre (France) was the most favorable decision for the liner shipping company in order to recover the vessel schedule (i.e., the port skipping costs were lower as compared to the other additional route service costs, which could be endured as a result of disruptions). On the other hand, an increase in the port skipping inconvenience coefficient to = [8.0] ÷ [10.0] could significantly increase the port skipping cost to be endured by the liner shipping company. Therefore, the liner shipping company did not skip the Port of Antwerp (Belgium) for all the considered problem instances and the port skipping inconvenience coefficient of = [8.0] ÷ [10.0]-see Figure A3b. Changes in the unit cost of HFO and MGO from one problem instance to another did not influence the port skipping decisions. However, the Port of Le Havre (France) was still skipped by the liner shipping company even after increasing the port skipping inconvenience coefficient to = [8.0] ÷ [10.0] in order to effectively recover the vessel schedule.