Stability metrics for a maritime inventory routing problem under sailing time uncertainty

Shaabani, Homayoun; Hvattum, Lars Magnus; Laporte, Gilbert; Hoff, Arild Article Stability metrics for a maritime inventory routing problem under sailing time uncertainty EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Shaabani, Homayoun; Hvattum, Lars Magnus; Laporte, Gilbert; Hoff, Arild (2024) : Stability metrics for a maritime inventory routing problem under sailing time uncertainty, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 13, Iss. 1, pp. 1-15, https://doi.org/10.1016/j.ejtl.2024.100146 This Version is available at: https://hdl.handle.net/10419/325216 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Stability metrics for a maritime inventory routing problem under sailing time uncertainty Homayoun Shaabani a,* , Lars Magnus Hvattum a , Gilbert Laporte a,b , Arild Hoff a a Faculty of Logistics, Molde University College, PO Box 2110, NO 6402, Molde, Norway b HEC Montr´ eal, Montr´ eal, H3T 2A7, Canada ARTICLE INFO Keywords: Reoptimization Uncertainty Stability metrics Maritime inventory routing ABSTRACT We study a multi-product maritime inventory routing problem (MIRP) with sailing time uncertainty. We explicitly consider the replanning that happens after uncertainty is revealed. The objective is to determine the stability of the adjusted plans after the occurrence of an uncertain event and to evaluate the effect of incorporating different stability metrics in the rescheduling process. Five stability metrics are introduced, and mathematical formulations of the MIRP incorporating each metric are presented. A reoptimization framework is then used to analyze the impact of each stability metric. Calculations are performed using 360 instances. The main result is that adjustments to the original plan occur at no additional cost almost 50% of the time. If decision makers want a more stable plan, they should accept a 5% cost deterioration, resulting in 20% more stable solutions. 1. Introduction In 2022, more than 80% of the volume of goods in international trade was carried by maritime transport, corresponding to 12.03 billion tons. It is expected that the volume of maritime trade will grow by more than 2% annually between 2024 and 2028 (UNCTAD, 2023), and therefore optimized maritime transportation is of great importance. We study the maritime inventory routing problem (MIRP) which is a particular maritime transportation planning problem (Papageorgiou et al., 2014). The MIRP is a variant of the inventory routing problem (IRP) in a maritime context. The IRP integrates inventory management decisions with routing decisions under a vendor-managed inventory (VMI) system, where the supplier is responsible for determining the delivery schedule for a given customer, the delivery quantity for that customer, and the assignment of customers to vehicle routes. In the MIRP there are five time elements, shown in Fig. 1 along with six event points labelled from “a" to “f". The routing is between points “a" and “b" and after “f". Although the vessel is stationed at a port between points “b" and “f", the temporal status of this interval affects the routing after “f" to reach “b" at the next port. There are always several uncertain parameters in maritime trade. According to UNCTAD (2021), supply chain disruptions, changes in globalization patterns, transportation costs, port congestion, and pandemics are the main uncertain elements. In the MIRP, there are several problem features that could be influenced by uncertainty, some of which are listed below. •The sailing time can be uncertain due to reasons such as bad weather conditions (Rodrigues and Agra, 2022), mechanical failure of vessels (Rodrigues and Agra, 2022), or the ice conditions in the Arctic region (Choi et al., 2015). •The waiting time can be uncertain due to port congestion (Agra et al., 2015). •The port delay time can be uncertain for reasons such as strikes and equipment failure at ports (Christiansen and Nygreen, 2005). •Demand, which is the main feature with uncertainty in inland IRPs (Touzout et al., 2021), can also be uncertain in a maritime setting (Cheng and Duran, 2004;Soroush and Al-Yakoob, 2018). Previous research on the MIRP under uncertainty has mostly focused on sailing time as an uncertain parameter. Papageorgiou et al. (2014) stated that the sailing time is one of the primary features influenced by uncertainty in maritime applications. Accordingly, sailing time is considered as the only source of uncertainty in the current study. * Corresponding author. E-mail addresses: [email protected] (H. Shaabani), [email protected] (L.M. Hvattum), [email protected] (G. Laporte), Arild.Hoff@himolde. no (A. Hoff). Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.sciencedirect.com/journal/euro-journal-on-transportation-and-logistics https://doi.org/10.1016/j.ejtl.2024.100146 Received 2 January 2023; Received in revised form 14 June 2024; Accepted 21 September 2024 EURO Journal on Transportation and Logistics 13 (2024) 100146 Available online 25 September 2024 2192-4376/© 2024 The Authors. Published by Elsevier B.V. on behalf of Association of European Operational Research Societies (EURO). This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). Three different approaches for dealing with uncertainty in an optimization problem have been proposed in the literature (Rodrigues and Agra (2022);De Maio et al. (2021);Aytug et al. (2005)), defined in Table 1. In a MIRP with reactive approaches, frequent adjustments to the original plan can lead to inefficiencies from the perspective of the port planner (Liu et al., 2017), since these adjustments can trigger a series of changes in subsequent decisions such as staff scheduling and container storage (Xu et al., 2012). It is therefore important to know how stable the adjusted plans are, i.e., how large the deviation from the original plan is when reactive actions are applied. To incorporate the frequency information, Cui et al. (2022) introduced a distributionally robust optimization aimed at ensuring the robustness of inventory replenishment and routing decisions against the impact of distributional ambiguity. Hence, the research question of the current study is: how to measure the stability of solutions to a MIRP under the uncertainty of sailing time? We introduce stability metrics that examine the sequence of routes, which port visits are made, and the quantities loaded and unloaded in each visit. The only paper having used a reactive approach for the MIRP is by Dong et al. (2018) who solved an uncertain MIRP using a mixed integer linear programming model and then reviewed the information revealed in each period. Whenever the solution obtained after considering this information is infeasible, a reoptimization is performed for the entire planning horizon, ignoring the amount of deviation from the original plan. After the reoptimization, the horizon is rolled forward and the procedure is iterated until the end of the horizon. Touzout et al. (2021) attempted to measure the stability of solutions to the IRP under uncertain demand using reoptimization models. They stated that their method could be extended to other applications and proposed to consider other sources of uncertainty. In this regard, the current study aims to introduce a reoptimization framework for the MIRP under sailing time uncertainty in which stability metrics are introduced. The main contribution of this paper is threefold. 1. Unlike Dong et al. (2018) and Touzout et al. (2021), who performed reoptimization at specific time intervals, in the current study reoptimization can occur at any point in time. This allows us to respond to uncertainties whenever they are revealed. This is explained in Section 5. 2. Relevant stability metrics for the MIRP are identified, and mathematical formulations for each of these metrics are proposed. This is explained in Section 6. 3. Each of the formulations is tested by performing computational experiments to determine the impact of each stability metric. This is explained in Section 7. The remainder of the paper is organized as follows. Section 2reviews the literature on MIRPs with uncertainties and classifies the papers according to uncertain parameters, approaches, and models. Section 3 provides a description of the problem. Mathematical notations are explained in Section 4. Section 5is devoted to the reoptimization framework. Stability metrics are introduced in Section 6, followed by their analysis in Section 7, where numerical results and findings are presented. Finally, Section 8provides concluding remarks. 2. Literature review The most recent review of the MIRP was presented by Papageorgiou et al. (2014), who studied a deterministic single-product MIRP. The authors stated that robustness is a challenge for MIRP and therefore recommended the development of approaches that can deal with uncertainty. Ksciuk et al. (2022) provided a review of uncertainty in maritime ship routing and scheduling, examining uncertainty in eight different problems, including the MIRP. The authors mentioned that in the MIRP, there are no fixed pickup and delivery port pairs, and no predetermined number of port calls. Therefore, they concluded that this makes the MIRP a challenging problem even without uncertainty. The current section focuses on MIRPs with uncertainty. The reviewed papers are summarized in Table 2, which indicates the uncertain parameters, the approaches used to deal with uncertainty, and the employed model. The remainder of this section first introduces some of the commonly used modelling techniques in an uncertain environment and then discusses each of the approaches used to deal with uncertainty. Some of the commonly used techniques for modelling of optimization problems under uncertainty are the following. •Stochastic programming is a modeling framework for optimization problems under uncertainty (Klein Haneveld et al., 2020), in which the uncertain parameters are assumed to follow known (or partially known) probability distributions (Rodrigues and Agra, 2022). •Recourse models are a class of models in stochastic programming, including two-stage and multistage models. When the true value of an uncertain parameter is observed, corrective actions can be taken in recourse models (Klein Haneveld et al., 2020). •Chance-constrained programming, introduced by Charnes and Cooper (1959), provides a tool for solving optimization problems Fig. 1. Five time elements. Table 1 Three approaches for dealing with uncertainty. Approaches Which decisions are made before uncertainty is revealed Which decisions are made after uncertainty is revealed Notes Considering uncertain information explicitly Considering deterministic parameters Proactive All decisions –No adjustment These approaches are better suited for problems with low uncertainty and where the original plan can be maintained without any adjustment. Reactive –An initial plan All the decisions are recourse actions These approaches are better suited for problems with high uncertainty. Mixed An initial plan –Some of the decisions are recourse actions This is known as a priori optimization, a concept introduced by Bertsimas et al. (1990) H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 2 under uncertainty. This method optimizes the problem in such a way that the constraints are satisfied with a given probability. The minimum required reliability should be set by the decision maker to a value between zero and one. If this value is zero, the decision maker is extremely risk seeking, and if it is one, it indicates an extremely conservative attitude (risk averse). •Robust optimization accounts for uncertainty sets where the probability distribution is unknown or does not exist (Rodrigues and Agra, 2022). The decision maker constructs a solution that is feasible for each realization of uncertainty in the given set (Bertsimas et al., 2011). In other words, it optimizes the problem based on the worst possible outcome within the uncertainty set. Unlike stochastic programming, where the uncertain parameter is assumed to be a random variable that follows a known (or partially known) probability distribution, the uncertainty model in robust optimization is usually deterministic and set-based (Bertsimas et al., 2011). Therefore, the uncertain parameters can take any value within the uncertainty set. •Reoptimization can be used to deal with uncertainty or in situations where the planning horizon is shorter than the horizon of the actual problem (Dong et al., 2018). •The integration of simulation and optimization can be used to deal with uncertainty. Zhou et al. (2021) studied different types of integration approaches in maritime logistics. Touzout et al. (2021) expressed that a priori approaches proactively address uncertainties by formulating robust replenishment plans. In terms of proactive approaches, Cheng and Duran (2004) considered a decision support system that uses a simulation model and an optimization model. The simulation model represents the inventory and transportation system, and the optimization model is formulated as a discrete-time Markov decision process that deals with the uncertainty of sailing time and of demand. A deterministic model with penalty costs is used as a proactive approach in two studies. First, Christiansen and Nygreen (2005) considered the uncertainty of sailing time and waiting time for a single-product MIRP. They applied soft inventory constraints, where levels should lie within a certain interval, and introduced lower and upper alarm intervals with artificial penalty costs to increase the robustness of their model. Second, Rakke et al. (2011) introduced a deterministic model with penalty costs for deviating from long-term customer contracts, maximizing revenue based on the spot market price and the quantity of sales in that market. In another proactive approach by Soroush and Al-Yakoob (2018) for a single-product MIRP, demand was assumed to be a normally distributed random variable, and penalties for understocking or overstocking were considered. The authors proposed a stochastic optimization model with linear constraints and a convex objective function. They used DICOPT as a commercial solver to solve the problem. Zhang et al. (2018) employed time windows to model sailing time uncertainty for a single-product MIRP. They defined flexible solutions as those that can accommodate unplanned disruptions by adjusting routing solutions where delivery dates and total delivery quantities cannot be changed. Furthermore, a Lagrangian heuristic was implemented to find flexible solutions using soft constraints, and a simulator was introduced that generates a disruption in each simulation run to evaluate the robustness of the solutions. Diz et al. (2019) considered the uncertainty of the total time vessels spend in ports due to delays in vessel operations for a single-product MIRP. They developed a robust optimization scheme using more vessels to protect the solution against delays. The risk of infeasibility was quantified for different levels of robustness and Gurobi used to solve the problem. Regarding the mixed approaches, three studies have considered recourse models that take into account routing, the quantities to be loaded and unloaded, the order of port visits in the first stage, as well as Table 2 Summary of MIRP papers with uncertainty. Year Authors Uncertain parameters Approach Models P a R b M c 2004 Cheng &Duran •Sailing times •Demand ✓  Simulation and optimization 2005 Christiansen & Nygreen •Sailing times •Port delay times ✓  Deterministic model with penalty cost for inventory violation 2011 Rakke et al. •Spot market price ✓  Deterministic model with penalty cost for deviation from the customer long-term contracts 2015 Agra et al. •Sailing times •Waiting times   ✓Stochastic programming 2016 Agra et al. •Sailing times   ✓Stochastic programming 2018 Soroush &Al-Yakoob •Demand ✓  Chance-constrained programming 2018 Agra et al. •Sailing times •Port delay times   ✓Robust optimization 2018 Cho et al. •Sailing times   ✓Stochastic programming 2018 Zhang et al. •Sailing times ✓  Stochastic programming 2018 Dong et al. •Vessel availability •Trip delays •Pick-up window information •Consumption and production rates ✓Reoptimization 2019 Diz et al. •Waiting times •Port delay times ✓  Robust optimization 2019 Rodrigues et al. •Sailing times ✓✓•Deterministic models with inventory buffers •Robust optimization •Stochastic programming •Conditional value-at-risk 2021 Liu et al. •Sailing times •Waiting times   ✓Two-stage distributionally robust optimization 2023 Nikolaisen et al. •Departure times •Sailing times   ✓Optimization and simulation 2024 Current study •Sailing times ✓Reoptimization including stability metrics a P: Proactive approaches. b R: Reactive approaches. c M: Mixed approaches. H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 3 visit times to ports and inventory decisions in the second stage, which can be adjusted to the scenario. The first study, by Agra et al. (2015), introduced a two-stage stochastic programming model with recourse and solved this model using a decomposition algorithm in which optimality cuts are added dynamically. The second study, by Agra et al. (2016), used a model similar to that of the previous study, but solved it with a combination of a commercial solver and local search heuristics. In the third study, by Agra et al. (2018), robust optimization was used and a decomposition algorithm was suggested. Also, an iterated local search heuristic was introduced to improve the decomposition algorithm. Several techniques to handle uncertainty in MIRPs were compared by Rodrigues et al. (2019), who considered uncertain sailing times for a single-product MIRP and employed different models and algorithms to handle uncertainty. They discovered that three methods provide a good trade-off between the amount and probability of inventory limit violations and routing costs. These methods are 1) deterministic modeling with inventory buffers, 2) stochastic programming with high penalties for inventory bounds violations, and 3) a hybrid algorithm that solves a deterministic approach with inventory buffers derived from a conditional value-at-risk approach. Another mixed methodology applied to MIRP under uncertainty comes from Cho et al. (2018), who proposed a two-stage stochastic programming model in which production inventory schedule decisions are made in the first stage and the production rate is adjusted for each scenario in the second stage. Liu et al. (2021) applied a two-stage distributionally robust optimization algorithm in which the routing decisions are made in the first stage, while decisions regarding quantities to be loaded and unloaded, visit time to ports, and inventory levels are made in the second stage after observing uncertainties. A reactive methodology was applied once in the context of MIRP by Dong et al. (2018). The authors developed stochastic simulations to account for several sources of uncertainty, presented in Table 2, and an algorithm that integrates reoptimization and stochastic simulation results. They reoptimized the model at a specified frequency, typically once per day. At each stage, the parameters are updated as uncertainties are observed, and the optimization problem is solved. This procedure is repeated for each day of the time horizon of the current problem. 3. Problem description The MIRP considers the transportation of products between multiple ports while meeting inventory requirements. Different ports produce and consume multiple products at a given production and consumption rate. Initial inventories, minimum inventory levels, and maximum inventory levels are specified for each port. A heterogeneous fleet of vessels with a given capacity, a fixed speed, and a daily operating cost is given. The position of a vessel at the beginning of the planning horizon is referred to as its origin, which can be a port or any location at sea. Sailing times from the origin to each port and between each pair of ports are determined based on the given distance and speed of the vessel. The sailing costs are also derived from the sailing time multiplied by the daily cost of a vessel. The maximum unloading quantities are determined by the consumption ports based on the vessel capacity and the maximum inventory of the port. The maximum number of visits to each port is predetermined. The holding cost and penalty cost for each product in each port are known. The objective of the problem is to minimize the sum of three components: sailing costs, inventory holding costs, penalty costs for backlogs and overstocks. The sailing times are assumed to be subject to uncertainty due to weather conditions. Although a planning horizon is specified, the uncertainty in sailing times may cause the planning horizon to be exceeded. The problem is solved under deterministic conditions and whenever the uncertainty is revealed, reoptimization is performed. 4. Mathematical notations This section explains some of the most frequently used notations throughout this paper, whereas the complete list of notations can be found in Appendix A. The problem consists of some ports represented by i,jand h, and each port can be visited at most mtimes. There is a set of products denoted by Kand a set of vessels denoted by V. We define a network in which the nodes are represented by (i,m), denoting the visit mto port i. The vessels movement from node (i,m)to node (j,n)are represented by (i,m,j,n). The set of possible port arrivals (i,m)is defined as SAand the set of port arrivals that may be made by vessel vis defined as SA v. The set of all possible vessel movements (i,m,j,n)is defined as SX and the set of all possible moves for vessel vis defined as SX v. The binary variable oimvk is one if and only if product kis loaded onto or unloaded from vessel vat the port visit (i,m). The amount of product k loaded onto or unloaded from vessel vat port visit (i,m)is denoted by qimvk. The amount of product kthat vessel vtransports from port visit (i,m)to port visit (j,n)is denoted by fimjnvk. Let simk represent the inventory level of product kat the start of port visit (i,m)and sE imk represent the inventory level of product kat the end of port visit (i,m). The sailing of vessel vfrom port arrival (i,m)directly to port arrival (j,n)is denoted by ximjnv, sailing of vessel vfrom its initial position to port arrival (i,m)is denoted by xO imv, the port visit (i,m)is denoted by yim, the visit to port iby vessel vat port arrival (i,m)is denoted by wimv. Let tim be the start time for port arrival (i,m)and tE im be the end time for port arrival (i,m). Fig. 2 depicts the route of one vessel as an example of this network, where Ovis the origin of v. 5. Reoptimization framework The mathematical formulation of the MIRP given in Appendix A is the same as that of Shaabani et al. (2023). It is used as the basic model of the reoptimization framework. The solution of this deterministic model is considered as the initial plan for the reoptimization. Periodic reoptimization considers the problem periodically at fixed time intervals. Continuous reoptimization, on the other hand, solves the problem throughout the day and whenever data change; a procedure collects the information up to that point and then starts the reoptimization (Pillac et al., 2013). In the current study, unlike Dong et al. (2018) and Touzout et al. (2021), a continuous-time model is used that reacts to uncertainties as soon as they appear, hence continuous reoptimization is performed. In this context, TUis defined as the time at which an uncertain event occurs. Therefore, the nominal sailing times are used until TUand then the sailing times are changed due to the uncertain event. Since TUis an uncertain event, it can occur at any time Fig. 2. Example of the network. H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 4 in [0,T]. We assume that there is only one TUvalue in each planning horizon. Whenever the TUvalue is revealed, the following changes are made to the basic model to prepare the model for reoptimization. •The set SAin the basic model, is replaced with SB, excluding port arrivals visited before TU. In the same manner, the set SA vis replaced with SB v, the set SXis replaced with SY, the set SX vis replaced with SY v. •The solution to the deterministic problem is extracted and defined as the data set for the reoptimization model. The sailing of vessel vfrom port arrival (i,m)directly to port arrival (j,n)is denoted by Ximjnv, the sailing of vessel vfrom its initial position to port arrival (i,m)is denoted by XO imv, the visit to port arrival (i,m)is denoted by Yim, the visit of port iby vessel vat port arrival (i,m)is denoted by Wimv, and the amount of product kloaded onto or unloaded from vessel vat port visit (i,m)is denoted by Qimvk. •Due to the occurrence of the uncertain event, two new time constraints are defined: tim ≥TU(i,m)∈SB(1) tE im ≥TU(i,m)∈SB.(2) If a vessel was on route from (i,m)to (j,n)when time hit TU, then the vessel is forced to visit (j,n), but the planned arrival time may be affected by the updated sailing times. •The initial inventory levels are updated when the new problem starts after TU. If Jik = − 1, the amount of inventory consumed up to TUis subtracted from the initial inventory, and if Jik =1, the amount of inventory produced up to TUis added to the initial inventory. •The values of the decision variables visited before TUare fixed. These decision variables are as follows: ximjnv,xO imv,oimvk,qimvk,fimjnvk,simk, sE imk,tE im. •Constraints (A28) and (A29) are deleted because the uncertainty of the sailing time is considered, which may lead to exceeding the planning horizon. Now the modified model is ready, and we call it “Model 0”, which represents the metric “cost”. The reoptimized solution represents the sailing costs and the port operation costs, plus penalty costs for backlogs and overstocks without stability metrics. The sailing costs and port operation costs are called C*. Therefore, the reoptimized solution may differ from the initial solution. In this context, stability metrics are introduced in the next section to reduce this discrepancy. 6. Stability metrics In this section we present five stability metrics for the MIRP and new constraints added to the model are then given for each metric. The objective function of the mathematical model for each of the stability metrics consists of two parts: first it minimizes the violation of the metrics, second it minimizes the penalty cost for backlogs and overstocks. Based on the values of consumption and production rates, vessel capacity, and minimum and maximum inventories, the sizes of the different elements in the objective function are such that we implicitly prioritize the first part of the objective function before the second part. 6.1. Sequence preservation The sequence preservation metric, called SP, means that the sequence of the reoptimized solution should not differ significantly from that of the original solution (Dettenbach and Ubber, 2015). Applications of the sequence preservation metric mostly belong to routing and scheduling problems (Touzout et al., 2021). In the MIRP, traveling times are typically much longer than in an inland IRP, and because the uncertain event can occur at any time, changing the sequence and rerouting may be costly. If the sequence of shipments has changed, more lifting operations are required at the new port in order to reach the unscheduled product unloads, resulting in higher costs. Another case where sequence preservation is critical in maritime transport occurs on transshipment routes where another vessel is waiting at a port of transshipment. The mathematical formulation of the SP metric contains two new binary variables. The binary variable zSP imjnv is defined to indicate whether or not there is a sequence change, and zSPO imv is a binary variable equal to one if and only if there is a change in the first visit made by the vessel. Therefore, two new constraints are defined as follows: zSP imjnv =Ximjnv −ximjnvv∈V,(i,m,j,n) ∈ SY v(3) zSPO imv =XO imv −xO imvv∈V,(i,m)∈SB v.(4) Constraints (3) and (4) are nonlinear but can be linearized into constraints (6) to (9) as shown by Touzout et al. (2021). Constraints (3) count a sequence change when an arc from (i,m)to (j,n)is visited by vessel vin the original solution but not in the reoptimized solution, and vice versa. Constraints Eq. (4) count a sequence change if vessel vsails directly from its initial position to port arrival (i,m)in the original solution but does not in the reoptimized solution, and vice versa. The mathematical formulation for the SP metric is as follows: Model 1: Reoptimization based on sequence preservation (SP) metric Minimize∑v∈V∑(i,m,j,n)∈SX vzSP imjnv +∑v∈V∑(i,m)∈SA vzSPO imv + ∑(i,m)∈SA∑k∈Kv|Jik=− 1CP ik(rimk +rE imk)+∑i∈N∑k∈Kv|Jik =− 1CP ikrT ik + ∑i∈N∑k∈Kv|Jik=1CPP ik rPT ik (5) subject to   (1) and (2)   (A2) to (A27) and (A30) to (A45)   zSP imjnv ≥Ximjnv −ximjnv v∈V,(i,m,j,n) ∈ SY v(6) zSP imjnv ≥ximjnv −Ximjnv v∈V,(i,m,j,n) ∈ SY v(7) zSPO imv ≥XO imv −xO imv v∈V,(i,m) ∈ SB v(8) zSPO imv ≥xO imv −XO imv v∈V,(i,m) ∈ SB v(9) zSP imjnv ∈ {0,1}v∈V,(i,m,j,n) ∈ SY v(10) zSPO imv ∈ {0,1}v∈V,(i,m) ∈ SB v.(11) 6.2. Sequence preservation with vessel replacement The sequence preservation with vessel replacement metric, called SPV, is similar to the SP metric with one difference. Touzout et al. (2021) proposed the SP metric where the sequence of an original solution must be preserved in the reoptimized solution; in comparison, the SPV metric measures the preservation of the sequence of the original solution for the reoptimized solution, even if the vessels are replaced with others. For example, an arc from (i,m)to (j,n)that was traversed by vessel 1 in the original solution can be traversed by vessel 3 in the reoptimized solution. A sequence here is meant to be the sequence of deliveries to ports regardless of vessel number. With respect to the SP metric, the SPV metric places greater emphasis on maintaining the delivery sequence, ensuring the order of deliveries to ports is preserved even with a new vessel, which is more stringent compared to the SP metric. This metric aims to enhance operational efficiency and logistical coordination at ports. Ports typically allocate schedules and resources based on the expected sequence of arrivals. Disrupting this sequence can lead to delays and increased waiting times for vessels. Additionally, cargo handling and distribution are often planned according to a specific sequence. By maintaining this sequence, the right resources and personnel can be ensured to be available at the right time and place. Besides, the SPV metric may make a correct measurement when the fleet is homogeneous, whereas SP metric may make a correct measurement when the fleet is heterogeneous. The SPV implicitly considers the arrival time because preserving the H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 5 sequence can also preserve the visiting times of the ports. This is shown in Fig. 3 which also shows the difference of SP and SPV. Fig. 3 shows an example for a part of a representation of a solution of the original plan (3a) and two examples of the reoptimized plans (3b,3c). The values of SP and SPV for each example are calculated and presented in Figure (3). In Figure (3b), the reoptimized solution maintains the same vessel as the original plan but changes the sequence of visits. In contrast, Figure (3c) shows a reoptimized solution that employs different vessels compared to the original plan while preserving the original sequence. When comparing the solutions in Figure (3b) and Figure (3c), the altered sequence in Figure (3b) can result in significant changes to the arrival times of most port visits, whereas the solution in Figure (3c) preserves the original sequence and, consequently, the same arrival times as the original plan. The mathematical formulation of the SPV metric uses two new binary variables. The binary variable zSPV imjn is defined to indicate whether or not there is sequence change with vessel replacement, and zSPVO im is a binary variable equal to one if and only if there is sequence change with vessel replacement for the initial position of the vessels. Therefore, two new constraints are defined as follows: zSPV imjn =∑ v∈V Ximjnv −∑ v∈V ximjnv (i,m,j,n) ∈ SY(12) zSPVO im =∑ v∈V XO imv −∑ v∈V xO imv (i,m) ∈ SB.(13) Like constraints (3) and (4), constraints (12) and (13) can be linearized, as is the case in constraints (15) to (18). Constraints (12) count a change of sequence with vessel replacement if an arc from (i,m)to (j,n) was visited in the original solution but not in the reoptimized solution and vice versa. Constraints (13) count a change of sequence with vessel replacement if a vessel sails directly from its initial position to port arrival (i,m)in the original solution but not in the reoptimized solution Table 3 Probability distributions for sailing time. Probability distributions Sailing time ×1×1.5×2 P=1 0.85 0.10 0.05 P=2 0.50 0.30 0.20 P=3 0.15 0.45 0.40 Fig. 3. Example for the calculation of the values for the SP and SPV. H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 6 and vice versa. The mathematical formulation for the SPV metric is as follows: Model 2: Reoptimization based on sequence preservation with vessel replacement (SPV) metric Minimize∑(i,m,j,n)∈SXzSPV imjn +∑(i,m)∈SAzSPVO im +∑(i,m)∈SA∑k∈Kv|Jik=− 1CP ik(rimk + rE imk)+∑i∈N∑k∈Kv|Jik =− 1CP ikrT ik +∑i∈N∑k∈Kv|Jik =1CPP ik rPT ik (14) subject to   (1) and (2)   (A2) to (A27) and (A30) to (A45)   zSPV imjn ≥∑v∈VXimjnv −∑v∈Vximjnv (i,m,j,n) ∈ SY(15) zSPV imjn ≥∑v∈Vximjnv −∑v∈VXimjnv (i,m,j,n) ∈ SY(16) zSPVO im ≥∑v∈VXO imv −∑v∈VxO imv (i,m) ∈ SB(17) zSPVO im ≥∑v∈VxO imv −∑v∈VXO imv (i,m) ∈ SB(18) zSPV imjn ∈ {0,1} (i,m,j,n) ∈ SY(19) zSPVO im ∈ {0,1} (i,m) ∈ SB.(20) 6.3. Visit deviation The visit deviation metric, called VD, compares port visits in the reoptimized solution to those in the original solution and counts visit violations which should be minimized (Touzout et al., 2021). This metric does not consider the vessel number. As an example, port iat visit mthat is visited by vessel 1 in the original solution, may be visited by vessel 3 in the reoptimized solution. This metric is helpful in situations where it is very costly to miss a scheduled visit. If a planned port visit is omitted in the reoptimized solution, this may result in wasted time and resources, and if a new port visit occurs in the reoptimized solution that was not planned in the original solution, this may result in higher operating costs due to the unavailability of resources at a port. The mathematical formulation of the VD metric has a new binary variable, zVD im which is defined to denote whether or not there is a visit deviation for port arrival (i,m). Thus, the new constraint is defined as follows: zVD im = |Yim −yim| (i,m) ∈ SB.(21) Like constraints (3) and (4), constraint (21) can be linearized, resulting in constraints (23) and (24). Constraints (21) count a visit violation if port arrival (i,m)was visited in the original solution and not in the reoptimized solution, and vice versa. The mathematical formulation for the VD metric is as follows: Model 3: Reoptimization based on visit deviation (VD) metric Minimize∑(i,m)∈SAzVD im +∑(i,m)∈SA∑k∈Kv|Jik =− 1CP ik(rimk +rE imk)+ ∑i∈N∑k∈Kv|Jik =− 1CP ikrT ik +∑i∈N∑k∈Kv|Jik=1CPP ik rPT ik (22) subject to   (1) and (2)   (A2) to (A27) and (A30) to (A45)   zVD im ≥Yim −yim (i,m) ∈ SB(23) zVD im ≥yim −Yim (i,m) ∈ SB(24) zVD im ∈ {0,1} (i,m) ∈ SB.(25) 6.4. Visit deviation without vessel replacement The visit deviation without vessel replacement metric, called VDV, computes the number of ports that are not visited in the reoptimized solution but are visited in the original solution with a certain vessel number. This metric differs with previous metric in terms of vessel number. For example, if W221 =1 and w223 =1 then VDV metric is equal to 2 but VD metric is equal to 0 because Y22 =1 and y22 =1. We have introduced the new VDV metric, which differs from the VD metric described by Touzout et al. (2021). Since these authors consider a single-product problem, all the vehicles transport the same product. Therefore, if the same customer is visited in the reoptimized solution compared to the original solution, then it is less important by which vehicle it has been visited. However, since the current study considers a multi-product problem, it could be important to use the same vessel in the reoptimized solution as in the original solution because each vessel can carry a different product mix. When determining which vessel will be assigned to a particular cargo, it is important to plan for any additional equipment or services needed for port operations, such as pilot, tugboats, and port services. These arrangements can be made in advance. However, the fleet of vessels is heterogeneous, and if a cargo is rescheduled and assigned to a different vessel, it may necessitate the procurement of additional equipment or services. In such cases, the existing arrangements must be cancelled, and new ones must be established. Making these changes, even if possible, can be a labour-intensive process (Fagerholt et al., 2009). The mathematical formulation of the VDV metric contains a new binary variable. Let zVDV idenote whether there is a visit deviation without vessel replacement for port ior not. Hence, the new constraint is defined as follows: zVDV imv = |Wimv −wimv|v∈V,(i,m) ∈ SB v.(26) Constraints (26) count a violation of visit without vessel replacement if the number of visits to port iin the original solution is not the same as in the reoptimized solution. Like constraints (3) and (4), constraints (26) can be linearized as shown in constraints (28) and (29). The mathematical formulation for the VDV metric is as follows: Model 4: Reoptimization based on visit deviation without vessel replacement (VDV) metric Minimize∑(i,m)∈SA v∑v∈VzVDV imv +∑(i,m)∈SA∑k∈Kv|Jik=− 1CP ik(rimk + rE imk)+∑i∈N∑k∈Kv|Jik =− 1CP ikrT ik +∑i∈N∑k∈Kv|Jik =1CPP ik rPT ik (27) subject to   (1) and (2)   (A2) to (A27) and (A30) to (A45)   zVDV imv ≥Wimv −wimv v∈V,(i,m) ∈ SB v(28) zVDV imv ≥wimv −Wimv v∈V,(i,m) ∈ SB v(29) zVDV imv ∈ {0,1}v∈V,(i,m) ∈ SB v.(30) 6.5. Quantity deviation The quantity deviation metric, called QD, calculates the difference between the quantity of product loaded onto or unloaded from a vessel at a port in the original solution and the loaded or unloaded quantity in the reoptimized solution. Minimizing this metric leads to fewer planning issues (Touzout et al., 2021) and it is significant because it is the only metric that addresses the inventory component of the MIRP. The mathematical formulation of the QD metric contains a new variable. Let zQD imvk be the difference in the quantity of product kloaded onto or unloaded from vessel vupon arrival at port (i,m)in the original solution and the reoptimized solution. Thus, the new constraint is defined as follows: zQD imvk = |Qimvk −qimvk|v∈V,(i,m) ∈ SB v,k∈Kv:Jik ∕= 0.(31) Constraints (31) calculate the difference between the loaded or unloaded quantity in the original solution and the reoptimized solution. Like constraints (3) and (4), constraints (31) can be linearized as in constraints (33) and (34). The mathematical formulation for the QD metric is as follows: Model 5: Reoptimization based on quantity deviation (QD) metric Minimize∑(i,m)∈SA v∑v∈V∑k∈Kv Jik∕=0zQD imvk +∑(i,m)∈SA∑k∈Kv|Jik=− 1CP ik(rimk + rE imk)+∑i∈N∑k∈Kv|Jik =− 1CP ikrT ik +∑i∈N∑k∈Kv|Jik =1CPP ik rPT ik (32) subject to   (1) and (2)   (A2) to (A27) and (A30) to (A45)   zQD imvk ≥Qimvk −qimvk v∈V,(i,m) ∈ SB v,k∈Kv:Jik ∕= 0(33) zQD imvk ≥qimvk −Qimvk v∈V,(i,m) ∈ SB v,k∈Kv:Jik ∕= 0(34) zQD imvk ≥0v∈V,(i,m) ∈ SB v,k∈Kv:Jik ∕= 0.(35) H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 7 7. Analysis of the stability metrics In this section, each stability metric is analyzed. First, the problem instances are presented in Section 7.1. The evaluation procedure is explained in Section 7.2 and numerical results and findings are presented in Section 7.3. 7.1. Problem instances There are 360 instances in this paper. Instances are divided into three groups called A, B and C. Group A consists of instances with one product, three vessels and eight ports; Group B consists of instances with two products, four vessels and 16 ports; and group C is similar to group A but with four products. In group A and B, each port is limited to at most one product, but in group C, there are sometimes more than one product in a given port. The input parameters of instances are derived from base cases I1 to I10 of Shaabani et al. (2023) where each case differs by the initial inventory of product kat port i, the maximum inventory of product kat port i, and the demand rate of port ifor product k. Two different time horizons (T=30,60 days)are considered for the problem. The instances are available upon request. The sailing times are assumed to be subject to uncertainty due to weather conditions. As was done by Agra et al. (2015), in the current study we introduce two possible changes in sailing times, considering that sailing times can also remain unchanged. In the first change, sailing times are increased to 1.5 times the original value, and in the second change, they are increased to twice the original value. Based on any of these changes, the sailing times for each port may change. Since uncertainty may affect an area at sea, sailing times are selected based on the combination of arrival and departure ports. Therefore, for example, if an event occurs in the area of port 1, all sailing times from port 1 to other ports and from other ports to port 1 are affected. Table 3 shows the three probability distributions considered, where each column represents one of the possible changes that can occur for sailing times and each row represents a probability distribution and gives the probability for each of the three possible changes. If all sailing times are assigned to the first column, then no uncertainty has occurred. However, since the current problem examines uncertainty in sailing times, this is not considered, and another change is created based on the same probability distribution until at least one sailing time is assigned to the second or third column. Besides, four scenarios are generated for each probability distribution. Table 4 shows the characteristics of the problem instances where the total number of instances is 360 (3 groups ×5 base cases × 2 time horizons ×3 probability distributions ×4 scenarios). Let g be a group of instances and G= {A,B,C}be the set of groups. Also, let ube a base case and Ua set of base cases, then U={I1,…,I5,g∈ {A,B} I6,…,I10,g∈ {C}. According to Shaabani et al. (2023), the structure of instances can make the MIRP difficult to solve; therefore, these instances are selected such that the optimal solutions can be found in less than 21,600 s, thereby enabling a fair analysis of the stability metrics. 7.2. Evaluation procedure A set of metrics which represents the ʹʹcos tʹʹ metric and five introduced stability metrics in Section 6is defined as Θ= {θ0,…,θ5}where θ0=cost,θ1=SP,θ2=SPV,θ3=VD,θ4=VDV,θ5=QD. All instances are solved directly by CPLEX 20.1. First the deterministic model given in Section 5is solved according to Algorithm 1. Then, the evaluation procedure for an instance is given in Algorithm 2 and matrix Rrepresents the structure of the outcomes for all metrics for an instance and is shown in Table 5. The details of the evaluation procedure are as follows. After solving the deterministic model the modifications to the basic model explained in Section 5are applied. Now the modified model can be solved, minimizing the ʹʹcos tʹʹ metric, and the reoptimized solution is obtained, which includes the sailing costs and port operation costs, plus penalty costs for backlogs and overstocks, but only the sailing costs and port operation costs, which is called C*, is reported. Since all models include penalty cost terms and the values obtained for these terms are mostly identical, the solutions for all models are reported without the value for Table 4 Characteristics of the problem instances. Groups gNumber of products |K| Number of vessels |V|Number of ports |N| Base cases U Time horizon T Probability distributions PScenarios S A1 3 8 {I1,…,I5} {30,60} {1,2,3} {1,2,3,4} B2 4 16 {I1,…,I5} C4 3 8 {I6,…,I10} Table 5 Matrix Rshowing the structure of outcomes for an instance. Cost SP SPV VD VDV QD Cost C*SP*(C*)SPV*(C*)VD*(C*)VDV*(C*)QD*(C*) SP C*(SP*)SP*SPV*(SP*)VD*(SP*)VDV*(SP*)QD*(SP*) SPV C*(SPV*)SP*(SPV*)SPV*VD*(SPV*)VDV*(SPV*)QD*(SPV*) VD C*(VD*)SP*(VD*)SPV*(VD*)VD*VDV*(VD*)QD*(VD*) VDV C*(VDV*)SP*(VDV*)SPV*(VDV*)VD*(VDV*)VDV*QD*(VDV*) QD C*(QD*)SP*(QD*)SPV*(QD*)VD*(QD*)VDV*(QD*)QD* Table 6 Additional constraints. Additional constraints indicator Added constraints α 0∑v∈V∑(i,m,j,n)∈SX vCT ijvximjnv +∑v∈V∑(i,m)∈SA vCTO oiv xO imv + ∑v∈V∑(i,m)∈SA v∑k∈Kv|Jik∕=0CO ikoimvk =C* α 1Constraints (6)to (11), and ∑v∈V∑(i,m,j,n)∈SX vzSP imjnv + ∑v∈V∑(i,m)∈SA vzSPO imv =SP* α 2Constraints (15)to (20), and ∑(i,m,j,n)∈SXzSPV imjn +∑(i,m)∈SAzSPVO im = SPV* α 3Constraints (23)to (25), and ∑(i,m)∈SAzVD im =VD* α 4Constraints (28)to (30), and ∑(i,m)∈SA v∑v∈VzVDV imv =VDV* α 5Constraints (33)to (35), and ∑(i,m)∈SA v∑v∈V∑k∈Kv Jik∕=0zQD imvk =QD* H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 8 Diz, G.S., dos, S., Hamacher, S., Oliveira, F., 2019. A robust optimization model for the maritime inventory routing problem. Flex. Serv. Manuf. J. 31, 675–701. https://doi. org/10.1007/s10696-018-9327-9. Dong, Y., Maravelias, C.T., Jerome, N.F., 2018. Reoptimization framework and policy analysis for maritime inventory routing under uncertainty. Optim. Eng. 19 (4), 937–976. https://doi.org/10.1007/s11081-018-9383-8. Fagerholt, K., Korsvik, J.E., Løkketangen, A., 2009. Ship routing scheduling with persistence and distance objecives. Lect. Notes Econ. Math. Syst. 89–107. https:// doi.org/10.1007/978-3-540-92944-4. Klein Haneveld, W.K., van der Vlerk, M.H., Romeijnders, W., 2020. Stochastic Programming. Springer International Publishing, Cham. https://doi.org/10.1007/ 978-3-030-29219-5. Ksciuk, J., Kuhlemann, S., Tierney, K., Koberstein, A., 2022. Uncertainty in maritime ship routing and scheduling: a Literature review. Eur. J. Oper. Res. https://doi.org/ 10.1016/j.ejor.2022.08.006. Liu, B., Zhang, Q., Yuan, Z., 2021. Two-stage distributionally robust optimization for maritime inventory routing. Comput. Chem. Eng. 149. https://doi.org/10.1016/j. compchemeng.2021.107307. Article 107307. Liu, C., Xiang, X., Zheng, L., 2017. Two decision models for berth allocation problem under uncertainty considering service level. Flex. Serv. Manuf. J. 29 (3–4), 312–344. https://doi.org/10.1007/s10696-017-9295-5. Nikolaisen, J.B., Vågen, S.S., Schütz, P., 2023. Solving a maritime inventory routing problem under uncertainty using optimization and simulation. Comput. Manag. Sci. 20 (1). https://doi.org/10.1007/s10287-023-00459-x. Papageorgiou, D.J., Nemhauser, G.L., Sokol, J., Cheon, M.S., Keha, A.B., 2014. MIRPLib - a library of maritime inventory routing problem instances: survey, core model, and benchmark results. Eur. J. Oper. Res. 235 (2), 350–366. https://doi.org/10.1016/j. ejor.2013.12.013. Pillac, V., Gendreau, M., Gu´ eret, C., Medaglia, A.L., 2013. A review of dynamic vehicle routing problems. Eur. J. Oper. Res. 225 (1), 1–11. https://doi.org/10.1016/j. ejor.2012.08.015. Rakke, J.G., Stålhane, M., Moe, C.R., Christiansen, M., Andersson, H., Fagerholt, K., Norstad, I., 2011. A rolling horizon heuristic for creating a liquefied natural gas annual delivery program. Transport. Res. C Emerg. Technol. 19 (5), 896–911. https://doi.org/10.1016/j.trc.2010.09.006. Rodrigues, F., Agra, A., 2022. Berth allocation and quay crane assignment/scheduling problem under uncertainty: a survey. Eur. J. Oper. Res. 303 (2), 501–524. https:// doi.org/10.1016/j.ejor.2021.12.040. Rodrigues, F., Agra, A., Christiansen, M., Hvattum, L.M., Requejo, C., 2019. Comparing techniques for modelling uncertainty in a maritime inventory routing problem. Eur. J. Oper. Res. 277 (3), 831–845. https://doi.org/10.1016/j.ejor.2019.03.015. Shaabani, H., Hoff, A., Hvattum, L.M., Laporte, G., 2023. A matheuristic for the multiproduct maritime inventory routing problem. Comput. Oper. Res. 154. https://doi. org/10.1016/j.cor.2023.106214. Article 106214. Soroush, H.M., Al-Yakoob, S.M., 2018. A maritime scheduling transportation-inventory problem with normally distributed demands and fully loaded/unloaded vessels. Appl. Math. Model. 53, 540–566. https://doi.org/10.1016/j.apm.2017.08.015. Touzout, F.A., Ladier, A.-L., Hadj-Hamou, K., 2021. Modelling and comparison of stability metrics for a re-optimisation approach of the inventory routing problem under demand uncertainty. EURO J. Transport. Logist. 10. https://doi.org/10.1016/ j.ejtl.2021.100050. UNCTAD, 2021. Review of Maritime Transport. United Nations, Geneva. Retrieved from. https://unctad.org/webflyer/review-maritime-transport-2021. UNCTAD, 2023. Review of Maritime Transport. United Nations, Geneva. Retrieved from. https://unctad.org/meeting/launch-review-maritime-transport-2023. Xu, Y., Chen, Q., Quan, X., 2012. Robust berth scheduling with uncertain vessel delay and handling time. Ann. Oper. Res. 192 (1), 123–140. https://doi.org/10.1007/ s10479-010-0820-0. Zhang, C., Nemhauser, G.L., Sokol, J., Cheon, M.S., Keha, A., 2018. Flexible solutions to maritime inventory routing problems with delivery time windows. Comput. Oper. Res. 89, 153–162. https://doi.org/10.1016/j.cor.2017.08.011. Zhou, C., Ma, N., Cao, X., Lee, L.H., Chew, E.P., 2021. Classification and literature review on the integration of simulation and optimization in maritime logistics studies. IISE Transactions 53 (10), 1157–1176. https://doi.org/10.1080/ 24725854.2020.1856981. H. Shaabani et al. EURO Journal on Transportation and Logistics 13 (2024) 100146 15