Enhanced MILP Approach for Long-Term Multi-Vessel Maritime Inventory Routing with Application to Antarctic Logistics

Abstract

The maritime inventory routing problem (MIRP) integrates vessel routing and inventory management over a planning horizon to optimize logistical operations in marine environments. While existing models predominantly address short-term planning with single vessels, this research advances the field by presenting a tightened mixed-integer linear programming (MILP) model designed for long-term planning with multiple vessels. The proposed model leverages an improved mathematical formulation and state-of-the-art optimization solvers to enhance computational performance. To demonstrate its applicability, the model was evaluated using benchmark instances from the literature and new instances derived from the logistics of Chilean scientific bases in Antarctica, a challenging and underexplored maritime environment. The results show computational time reductions of up to 98% for small to medium-sized instances, achieved through the incorporation of valid inequalities into the model and the use of advanced hardware and solvers. For larger instances, optimal or near-optimal solutions were achieved within one hour for a planning horizon of 60 time units, with optimality gaps below 24.7% for a 120-time-unit horizon. These findings highlight the potential of the model to support decision-making in complex maritime logistics scenarios, extending its application to long-term, multi-vessel operations in remote and environmentally sensitive regions. The proposed framework provides a valuable tool for enhancing the sustainability and efficiency of maritime logistics systems.

1. Introduction

The coordination of the route of cargo vessels, considering the levels of inventory of products in the ports along the route, presents a complex management task. This situation leads to the formulation of the Maritime Inventory Routing Problem (MIRP), which represents a challenging optimization problem that can be solved through mathematical modeling. The MIRP is critical for different types of industries, such as oil and gas [1], chemicals [2], bulk commodities [3], and waste recovery [4]. The MIRP generally aims to optimize transportation and inventory costs while reducing environmental impact by minimizing fuel consumption, emissions, and inventory maintenance costs. The MIRP involves various factors, such as demand and supply locations, production and consumption rates, vessel capacities, or sailing times. When addressing the combined issue of distribution and inventory management, it is possible to reduce costs by 10% and the number of vehicles employed [5].

Researchers have extensively studied Inventory Routing Problems (IRPs) in the literature, proposing various formulations and solution methods for nearly four decades. These formulations vary in planning horizon, policies, objectives, and decisions [6]. Alternative IRP approaches include partially satisfying a customer’s demand by dropping cargo to a nearby second customer [7] or considering time-variable demands [8]. The MIRP is a particular case of IRP in which the transportation of goods is restricted to maritime vessels. A classical formulation of the MIRP considers discrete time and multi-vessel, with valid inequalities to tighten the search space [9]. Determining an optimal solution for large-scale instances of the problem, particularly those with extended planning horizons, multiple vessels, and numerous ports, remains a significant computational challenge.

Despite extensive research on IRPs and their maritime counterparts (MIRP), a significant research gap exists in addressing large-scale, multi-vessel scenarios over extended planning horizons. Most existing studies focus on short-term planning horizons (e.g., 60 time units or shorter) or consider only a single vessel, which limits their applicability to real-world situations where coordination of multiple vessels over longer periods is essential. Additionally, achieving optimal solutions for larger instances remains a computational challenge due to the increased complexity associated with additional vessels and extended time frames.

To address this gap, we propose an improved and tightened MILP model for the Multi-Vessel Maritime Inventory Routing Problem (MV-MIRP). Our model extends existing formulations [10,11,12,13,14] by incorporating multiple heterogeneous vessels and longer planning horizons, capturing the complexities and dynamics of real-world maritime logistics operations. By introducing valid inequalities and leveraging advanced computational tools, our model enhances computational efficiency and scalability, enabling it to solve large-size instances that were previously intractable. To evaluate the practicality of the proposed model, we generated problem instances using real-world data from Chilean Antarctic bases. These data include the geographical coordinates of seven bases, derived from Google Earth, and logistical requirements based on the daily meal demands of the scientific personnel stationed at these locations.

Antarctic logistics presents a unique and challenging context for maritime inventory routing due to the extreme environment, long distances between bases, and the need for efficient and sustainable resource allocation. The Chilean scientific bases serve as an ideal case study to validate the applicability of our model to real-world scenarios requiring multi-vessel coordination and long-term planning.

The key contributions of this paper are as follows:

• A novel MILP model for the MIRP: We introduce a new mixed-integer linear programming (MILP) model that addresses the Multi-Vessel Maritime Inventory Routing Problem (MV-MIRP), incorporating a fleet of heterogeneous vessels, multiple depots, and a continuous time horizon.

• Tightened formulation: Our model incorporates valid inequalities that improve computational efficiency, enabling it to solve large-scale instances that previously resulted in unsolved gaps and reducing the optimality gap in other challenging cases.

• Comparative computational experiments: We conducted experiments to evaluate the performance of our model using benchmark data sets from the existing literature, specifically discrete-time formulations.

• Original real-world instances: We generate new problem instances based on real-world data from Chilean scientific bases in Antarctica, adding a layer of practical application to our theoretical model.

• Evaluation of new instances: Our model solves several generated instances to optimality, while others show reduced gaps, demonstrating its capability to address real-world complexities and offering valuable insights into its scalability and efficiency.

The remainder of this paper is structured as follows. Section 2 reviews the relevant literature, Section 3 offers an in-depth description of the problem and explains the formulations used, Section 4 presents the findings, and Section 5 presents a discussion about scalability, uncertainty management, and sustainability. Section 6 draws conclusions and outlines potential directions for future research.

2. Literature Review

The MIRP is a complex challenge faced by shipping companies and logistics providers involved in maritime transport. This problem involves intricate decision-making processes due to various factors, such as multiple products, time windows, vessel capacities, port constraints, inventory holding costs, and uncertain demand. The MIRP is a computationally complex and multidimensional problem that requires the expertise of various fields, such as operations research, transportation, logistics, and supply chain management. Efficient logistics of shipping inventory directly impacts the cost-effectiveness of global trade and logistics operations. By optimizing inventory allocation, vessel routes, and schedules, companies can minimize transportation costs, reduce inventory holding costs, and improve overall supply chain performance. For a more in-depth view, see [15,16].

The MIRP can be formulated considering a continuous or discrete time horizon. In models that use continuous time, the entire planning period is treated in a continuous way. Decision variables can be made for any time between zero and the end of the planning period, and constraints, such as inventory limits, are applied at all times. However, the parameters must remain constant over time. In contrast, the discrete-time formulation subdivides the planning horizon into intervals. Continuous-time formulations are more appropriate for applications with constant production and consumption rates. In contrast, discrete-time formulations are better suited for applications where production and consumption rates vary over time [17]. In practice, the choice of formulation depends on the specific characteristics of the problem considered. For a given problem instance, the best algorithmic performance occurs with a tight continuous-time formulation [18].

Regardless of the formulation, the MIRP often must address end-of-horizon difficulties, which can be addressed by defining a parameter to limit the difference between initial and final inventories in ports [10] or by producing a repeatable cyclic solution [14]. However, these approaches have only been tested in instances with up to six ports and five vessels, which did not reach optimality even after two hours of computational time. Optimality was achieved only in instances with up to four ports and three vessels. Stochasticity can be modeled using a robust optimization approach and a Benders decomposition algorithm [19], which has achieved optimality in instances with up to eight ports and six vessels, with the largest instances with sixteen ports and seven vessels. Furthermore, different techniques have been compared to model uncertainty in the MIRP, including deterministic models with safety stocks, stochastic programming models, robust optimization models, and models based on conditional value-at-risk measures [20]. These techniques have been tested in instances with up to six ports and five vessels, with deterministic and stochastic approaches achieving the least computational time. In turn, robust optimization and conditional value-at-risk approaches obtained better results, but with a higher computational cost. The critical parameters in the formulations are each port’s consumption and production rates, which can be treated as decision variables to reduce costs over the entire time horizon. This approach achieved an average cost reduction of 11.4% in all instances with up to six ports and five vessels [11]. Other critical parameters of the model are the speed and load of the vessels, which influence fuel consumption. Taking these parameters into account, the costs of sailing can be reduced to 38% [12]. Furthermore, ref. [21] aimed to minimize fuel consumption by modeling speed as a continuous decision variable; however, they ended up with a nonlinear formulation. They solved instances with up to ten ports and six vessels using a metaheuristic algorithm.

In discrete-time formulations, ref. [22] developed a practical approach with a flexible modeling framework for the more general IRP, and successfully applied it to a MIRP. Furthermore, ref. [18] developed discrete- and continuous-time formulations for the MIRP, and ref. [23] tested several matheuristics approaches to solving the MIRP, which outperformed MILP formulations solved with commercial solvers. Some works use alternative time approaches, such as [24], which developed a hybrid continuous-discrete-time formulation. Ref. [25] proposes a formulation for the short-term planning of LNG deliveries, solving instances with up to 62 vessels using a metaheuristic approach. In summary, although novel and innovative modeling techniques have been devised, achieving optimality remains a challenge for larger instances, particularly when compared to other problems such as the vehicle routing problem and the traveling salesman problem.

Valid inequalities are linear constraints that are added to an MILP formulation to reduce the size of the feasible region and improve the quality of the linear relaxation [26]. They are often used to solve complex optimization problems that arise in maritime transportation, such as the MIRP. Ref. [9] proposes a discrete-time formulation for the MIRP and derives several classes of valid inequalities based on flow, inventory, and capacity constraints, improving the performance of the branch-and-cut algorithm for large-scale instances. Ref. [27] introduces a new formulation for the MIRP based on customer delivery patterns, solving it with a branch-price-and-cut algorithm using valid inequalities based on loading operations and timing cuts. Their algorithm produce lower bounds that are stronger than those of the default solving algorithm on a commercial solver. Ref. [28] proposes a discrete-time formulation for a crude oil distribution problem using cascading knapsack inequalities, which are linear constraints obtained by reformulating inventory balance constraints. This type of inequality has a structure that many commercial solvers can exploit, improving the quality of the solution.

Despite numerous studies proposing models and solution methods for the Inventory Routing Problem (IRP), scalability remains a critical challenge. Researchers have employed decomposition techniques such as two-phase decomposition methods [29], Benders decomposition [30,31], and Lagrangian decomposition [32] to address this issue. Other approaches include leveraging parallel computing resources [33] and developing novel algorithms like the Three-Front Parallel Branch-and-Cut Algorithm (3FP-B&C) [34]. Additionally, various metaheuristic algorithms—such as Simulated Annealing and Iterated Local Search [35], Genetic Algorithms (GAs) [36], and Modified Adaptive GAs [37]—have been employed. Matheuristic methods that combine mathematical programming and heuristics have also been proposed [38].

Although these methods have shown effectiveness in finding high-quality solutions within acceptable computational times, they often sacrifice the optimality guarantees inherent in exact methods. Our work contributes to this field by presenting an exact optimization model capable of solving medium-sized instances efficiently. Recognizing the scalability limitations of exact approaches, we aim to build upon this foundation by incorporating strategies to enhance computational performance for larger instances in future research.

A comparison of the main characteristics of the existing literature on the MIRP is presented in Table 1. The table includes the following aspects: the nature of the time horizon (discrete (D) or continuous (C)), the length of the time horizon, the type of algorithm used (exact (E) or metaheuristic (M)), and the number of considered vessels. The first column lists the authors of each work, whereas the second and third columns indicate whether the time horizon is discrete or continuous. The fourth column indicates time horizons of 60 time units or shorter, while the fifth column indicates time horizons longer than 60 time units. The sixth and seventh columns classify the algorithm used as exact or metaheuristic. The eighth and ninth columns report the number of vessels involved in the problem. Table 1 shows that only two studies have considered time horizons longer than 60 time units: Agra et al. (2022) [10], who used a single vessel and a metaheuristic method to solve instances of up to 240 time units, and Papageorgiou et al. (2018) [23], who used a metaheuristic method and fictitious markets in each port to solve instances of up to 360 time units. This work fills a gap in the literature by proposing an exact approach to solving multi-vessel problems with time horizons exceeding 60 time units.

3. Mathematical Model

In this section, we consider a new approach to address the MV-MIRP. We adopt a multi-vessel approach, capturing the complexities and dynamics in real-world shipping industries, such as fleet optimization and scheduling. With this approach, our objective is to provide a more comprehensive framework for solving the MIRP, improving the decision support for maritime logistic practitioners.

The MV-MIRP considers a set $N=\{1,\dots ,n\}$ of ports, where each port can be visited a maximum of ${\overline{m}}_i$ times. Consequently, we define a visit by the pair $(i,m)$, with $i\in N$ and $m\in \{1,\dots ,{\overline{m}}_i\}$. The port i either produces ($J_i=1$) or consumes ($J_i=-1)$ inventory at a constant rate $R_i$. Let $G=(V,A)$ be a directed graph where the set V contains an origin node $\left\{O\right\}$ and all nodes formed by the ports and their respective visits. The set of arcs A comprises all possible combinations of nodes. However, an arc $\left(\right(i,m),(j,n\left)\right)\in A$ is considered valid only if no consecutive visits occur at the same port i, satisfying the condition $i\ne j$. The port visits are executed by a fleet of heterogeneous vessels, denoted $K=\{1,\dots ,k\}$, each with capacity $C_k$. The trip of these vessels starts at the node of origin $\left\{O\right\}$ and can end at any node on the graph. Finally, we define the parameters and variables used in our model to facilitate further analysis.

3.1. Parameters

• $J_i$: 1 if port i is a supplier, and −1 if port i is a consumer;

• ${\overline{m}}_i$: Maximum possible visits to port i;

• $T_i^Q$: Rate of time/(unit of cargo) loaded or unloaded in port i;

• $C_{ij}^k$: Variable cost for traveling from port i to port j by vessel k;

• $C_0^k$: Fixed cost for operating vessel k;

• $C_y^k$: Fixed cost per visit to port y by vessel k;

• ${\overline{S}}_i$: Maximum inventory level in port i;

• ${\underline{\mathit{S}}}_i$: Minimum inventory level in port i;

• $S_i^0$: Initial stock level in port i;

• $R_i$: Production or consumption rate in port i;

• $Q_k^0$: Initial cargo in vessel k on the origin O;

• $C_k$: Maximum capacity of vessel k;

• $T_{ij}^S$: Travel time from port i to port j;

• $T_i^0$: Travel time from the origin O to port i;

• T: Planning horizon.

Although time and cost units are arbitrary, it is crucial to maintain consistency throughout the entire model. For instance, if the planning horizon T is measured in days, then the travel times $T_i^0$ and $T_{ij}^S$ must also be measured in days, and the rate $T_i^Q$ must be measured in days per unit of cargo. A similar approach is required for costs $C_{ij}^k$, $C_0^k$, and $C_y^k$.

3.2. Binary Variables

• $x_{imjnk}$: 1 if vessel k travels directly from the node $(i,m)$ to the node $(j,n)$, 0 otherwise;

• $x_{imk}^0$: 1 if vessel k travels directly from the origin to the node $(i,m)$, 0 otherwise;

• $y_{imk}$: 1 if vessel k makes the $m^{th}$ visit to port i, 0 otherwise;

• $z_{imk}$: 1 if vessel k ends its journey at the node $(i,m)$, 0 otherwise.

3.3. Continuous Variables

• $q_{imk}$: Quantity collected or delivered to the node $(i,m)$ by vessel k;

• $f_{imjnk}$: Quantity transported from the node $(i,m)$ to the node $(j,n)$ by vessel k;

• $f_{imk}^0$: Quantity transported from the origin to the node $(i,m)$ by vessel k;

• $f_{imk}^D$: Quantity transported from the node $(i,m)$ to the destination by vessel k;

• $t_{im}$: Start time of operation at the node $(i,m)$;

• $s_{im}$: Stock level at the beginning of the operation at the node $(i,m)$.

We present an example in which two vessels navigate the graph; thus, $K=\{k_1,k_2\}$. Figure 1 illustrates the spatial distribution of ports 1, 2, and 3, forming the network over which the vessels route. Understanding the layout of these ports is essential to model routing decisions and travel times in the MIRP. Figure 2 displays an expanded graph that includes the origin node, the destination node, and the nodes representing each pair of port and visit, demonstrating the scheduling and routing complexities inherent in coordinating multiple vessels over multiple visits in the MV-MIRP. We consider that ports 1 and 2 can each be visited a maximum of two times—thus, ${\overline{m}}_1={\overline{m}}_2=2$—while port 3 can be visited up to three times, which means ${\overline{m}}_3=3$.

The solution to this example unfolds as follows: vessel 1 departs from the origin to the node $(1,1)$; thus, $x_{111}^0=1$. Then, it travels directly from the node $(1,1)$ to the node $(3,1)$, indicated by $x_{11311}=1$. The remainder of the voyage of vessel 1 is as follows: $(3,1)\to (2,1)\to (3,2)\to (1,2)$. Since $(1,2)$ is the last node visited by vessel 1 before proceeding to the destination, we have $z_{121}=1$.

Vessel 2 departs from the origin to the node $(3,3)$, resulting in $x_{332}^0=1$. Then, it travels directly to the node $(2,2)$; thus, $x_{33222}=1$. Since the node $(2,2)$ is the last node visited by vessel 2 before heading to the destination, we have $z_{222}=1$. To improve the understandability of Figure 2, the edges used by vessel 2 are colored blue.

Regarding the cargo transported, in this example, vessel 1 carries 10 units of cargo from the origin ($f_{111}^0=10$), while vessel 2 carries 0 units ($f_{332}^0=0$). In the example, vessel 1 unloads 5 units of cargo at the node $(1,1)$—therefore, $q_{111}=5$—and transports the remaining 5 units to the node $(3,1)$, which is represented by $f_{11311}=5$. To avoid cluttering Figure 2, we proceed directly to the variable $f^D$ for vessel 1. In this case, vessel 1 transports 1 unit of cargo to the destination; thus, $f_{121}^D=1$. For vessel 2, the variables function analogously: $f_{332}^0=0$, $q_{332}=6$, $f_{33222}=6$, $q_{222}=6$, and $f_{222}^D=0$.

Next, we present the mathematical programming model.

3.4. Objective Function

(1) $minz=\sum _{k\in K}\sum _{\left(\right(i,m),(j,n\left)\right)\in A}{C}_{ij}^k{x}_{imjnk}+\sum _{k\in K}\sum _{(i,m)\in V}{C}_0^k{x}_{imk}^0+\sum _{k\in K}\sum _{(i,m)\in V}{C}_y^k{y}_{imk}$

The objective function (1) minimizes total costs, considering travel costs and fixed costs per used vessel, and the cost per visit.

To enhance readability, the constraints are grouped into thematic subsections that align with key operational aspects of the MV-MIRP. Each group addresses a specific set of requirements:

• Routing constraints (Section 3.5) ensure vessel connectivity and visit feasibility. This group primarily involves the variables x, y, and z.

• Pickup and delivery constraints (Section 3.6) maintain cargo balance and respect vessel capacities. This group primarily involves the variables f and q, along with x, y, and z.

• Time constraints (Section 3.7) regulate service times at ports and travel durations. This group primarily involves the variables t and x.

• Inventory constraints (Section 3.8) manage stock levels and production/consumption at ports. This group primarily involves the variables s, q, and t.

3.5. Routing Constraints

$\begin{array}{}\mathrm{(2)}& \hfill {\displaystyle \sum _{(i,m)\in V}}{x}_{imk}^0& \le 1\hfill & \hfill & \forall k\in K\hfill \mathrm{(3)}& \hfill y_{imk}-{\displaystyle \sum _{(j,n)\in V|j\ne i}}{x}_{jnimk}-x_{imk}^0& =0\hfill & \hfill & \forall (i,m)\in V,k\in K\hfill \mathrm{(4)}& \hfill y_{imk}-{\displaystyle \sum _{(j,n)\in V|j\ne i}}{x}_{imjnk}-z_{imk}& =0\hfill & \hfill & \forall (i,m)\in V,k\in K\hfill \mathrm{(5)}& \hfill {\displaystyle \sum _{k\in K}}{y}_{imk}& \le 1\hfill & \hfill & \forall (i,m)\in V\hfill \mathrm{(6)}& \hfill {\displaystyle \sum _{k\in K}}{y}_{i(m-1)k}-{\displaystyle \sum _{k\in K}}{y}_{imk}& \ge 0\hfill & \hfill & \forall (i,m)\in V|m\ge 2\hfill \mathrm{(7)}& \hfill x_{imk}^0,y_{imk},z_{imk}& \in \{0,1\}\hfill & \hfill & \forall (i,m)\in V,k\in K\hfill \mathrm{(8)}& \hfill x_{imjnk}& \in \{0,1\}\hfill & \hfill & \forall \left(\right(i,m),(j,n\left)\right)\in A,k\in K\hfill \end{array}$

Constraints (2) states that all vessels are located first at the origin. Constraints (3) and (4) are vessel connectivity constraints. Constraints (5) ensures that each node $(i,m)$ can be visited by at most one vessel. Constraints (6) ensures the feasibility of consecutive visits to port i (that is, node $(1,2)$ cannot be visited before node $(1,1)$). Constraints (7) and (8) are variable definitions.

3.6. Pickup and Delivery Constraints

(9) $\begin{array}{c}\sum _{k\in K}{f}_{imk}^0+\sum _{k\in K}\sum _{(j,n)\in V|j\ne i}{f}_{jnimk}+J_i\sum _{k\in K}{q}_{imk}\hfill \\ \hfill =\sum _{k\in K}\sum _{(j,n)\in V|j\ne i}{f}_{imjnk}+\sum _{k\in K}{f}_{imk}^D\forall (i,m)\in V\end{array}$

$\begin{array}{}\mathrm{(10)}& \hfill f_{imk}^0& =Q_k^0{x}_{imk}^0\hfill & \hfill & \forall (i,m)\in V,k\in K\hfill \mathrm{(11)}& \hfill f_{imjnk}& \le C_k{x}_{imjnk}\hfill & \hfill & \forall \left(\right(i,m),(j,n\left)\right)\in A,k\in K\hfill \mathrm{(12)}& \hfill f_{imk}^D& \le C_k{z}_{imk}\hfill & \hfill & \forall (i,m)\in V,k\in K\hfill \mathrm{(13)}& \hfill q_{imk}& \le C_k{y}_{imk}\hfill & \hfill & \forall (i,m)\in V,k\in K\hfill \mathrm{(14)}& \hfill f_{imk}^0,f_{imk}^D,q_{imk}& \ge 0\hfill & \hfill & \forall (i,m)\in V,k\in K\hfill \mathrm{(15)}& \hfill f_{imjnk}& \ge 0\hfill & \hfill & \forall \left(\right(i,m),(j,n\left)\right)\in A,k\in K\hfill \end{array}$

Constraints (9) ensures the balance of the cargo at each node $(i,m)$. Constraints (10) states that the cargo transported from the origin node to the node $(i,m)$ is equal to the initial cargo on vessel k. Constraints (11) and (12) ensure that the cargo transported cannot exceed the capacity of vessel k. Constraints (13) imposes upper limits on the pickup or delivery of cargo from each node $(i,m)$. Constraints (14) and (15) ensure that the variables are nonnegative.

3.7. Time Constraints

$\begin{array}{}\mathrm{(16)}& \hfill t_{im}+T_i^Q{q}_{imk}+(T+T_{ij}^S)x_{imjnk}& \le T+t_{jn}\hfill & \hfill & \forall \left(\right(i,m),(j,n\left)\right)\in A,k\in K\hfill \mathrm{(17)}& \hfill T_i^0{\displaystyle \sum _{k\in K}}{x}_{imk}^0& \le t_{im}\hfill & \hfill & \forall (i,m)\in V\hfill \mathrm{(18)}& \hfill 0\le t_{im}& \le T\hfill & \hfill & \forall (i,m)\in V\hfill \end{array}$

Constraints (16) gives a lower limit on the start time of service at node $(j,n)$ by vessel k, after sailing directly from the node $(i,m)$. Constraints (17) gives a lower bound to the start time of service on the node $(i,m)$ by vessel k. Constraints (18) provides a minimum and maximum value for the t variables.

3.8. Inventory Constraints

$\begin{array}{}\mathrm{(19)}& \hfill S_i^0+J_i{R}_i{t}_{i1}& =s_{i1}\hfill & \hfill & \forall i\in N\hfill \mathrm{(20)}& \hfill s_{i(m-1)}-{\displaystyle \sum _{k\in K}}{J}_i{q}_{i(m-1)k}+J_i{R}_i(t_{im}-t_{i(m-1)})& =s_{im}\hfill & \hfill & \forall (i,m)\in V|m\ge 2\hfill \mathrm{(21)}& \hfill s_{im}+{\displaystyle \sum _{k\in K}}{q}_{imk}-{\displaystyle \sum _{k\in K}}{R}_i{T}_i^Q{q}_{imk}& \le {\overline{S}}_i\hfill & \hfill & \forall (i,m)\in V|J_i=-1\hfill \mathrm{(22)}& \hfill s_{im}-{\displaystyle \sum _{k\in K}}{q}_{imk}+{\displaystyle \sum _{k\in K}}{R}_i{T}_i^Q{q}_{imk}& \ge {{\displaystyle \underline{\mathit{S}}}}_i\hfill & \hfill & \forall (i,m)\in V|J_i=1\hfill \mathrm{(23)}& \hfill s_{i\overline{m}}+{\displaystyle \sum _{k\in K}}{q}_{i\overline{m}k}-R_i(T-t_{i\overline{m}})& \ge {{\displaystyle \underline{\mathit{S}}}}_i\hfill & \hfill & \forall i\in N|J_i=-1\hfill \mathrm{(24)}& \hfill s_{i\overline{m}}-{\displaystyle \sum _{k\in K}}{q}_{i\overline{m}k}+R_i(T-t_{i\overline{m}})& \le {\overline{S}}_i\hfill & \hfill & \forall i\in N|J_i=1\hfill \mathrm{(25)}& \hfill {\underline{\mathit{S}}}_i\le s_{im}& \le {\overline{S}}_i\hfill & \hfill & \forall (i,m)\in V\hfill \end{array}$

Constraints (19) specifies that the inventory level at the first visit to port i is determined by the initial inventory at port i ($S_i^0$) and the inventory produced (for $J_i=1$) or consumed (for $J_i=-1$) between the beginning of the time horizon and the time of the first visit. Constraints (20) ensures that the inventory level at port i during visit m depends on the inventory level from the previous visit $m-1$. This relationship accounts for the amount of cargo collected (for $J_i=1$) or delivered (for $J_i=-1$) during the visit $m-1$, and the inventory produced or consumed at port i between visits $m-1$ and m. Constraints (21) and (22) ensure that inventory levels remain within their prescribed lower and upper limits at each port at the end of each visit. These constraints account for the cargo collected or delivered and the inventory produced or consumed during the handling of the cargo $q_{imk}$. Constraints (23) and (24) ensure that the inventory levels during the final visits remain within the allowed limits. This takes into account the cargo handled during the last visit and the inventory produced or consumed between the last visit and the end of the time horizon. Constraints (25) provides an upper and lower bound for the inventory level at each port at the beginning of each visit.

The constraints presented ensure a comprehensive formulation of the MV-MIRP by rigorously addressing its multifaceted requirements. Although the detailed structure of the constraints reflects the complexity of the problem, their thematic grouping provides a systematic framework for understanding their purpose. This detailed formulation is essential to achieve optimal solutions and accurately capture the intricacies of real-world maritime logistics.

3.9. Valid Inequalities

To improve the computational efficiency of our MV-MIRP model, we incorporate four sets of valid inequalities. These inequalities are designed to tighten the formulation by reducing the feasible region of the linear relaxation without excluding any feasible integer solutions. In doing so, they improve the quality of the lower bounds and accelerate the convergence of the solver. The valid inequalities included are the following: (i) minimum number of visits to a port (Section 3.9.1), (ii) logic inequalities (Section 3.9.2), and (iii) symmetry-breaking inequalities (Section 3.9.3).

3.9.1. Minimum Number of Visits to Port

The first set of valid inequalities ensures that the minimum required number of visits to each port is accounted for, based on the net demand over the planning horizon. This approach prevents the solver from considering solutions that do not meet the essential supply or demand requirements of each port, thereby eliminating infeasible or suboptimal solutions early in the process.

For their formulation, we calculate the net demand $Q_i^N$ for each port i over the time horizon T. This net demand represents the total quantity that must be transported to or from port i to satisfy its production or consumption needs while respecting inventory limits.

For supplier ports ($J_i=1$), we use Equation (26):

(26)$Q_i^N=max\{T{R}_i+S_i^0-{\overline{S}}_i,0\}\forall i\in N,|J_i=1$

For consumer ports ($J_i=-1$), we use Equation (27):

(27)$Q_i^N=max\{T{R}_i-S_i^0,0\}\forall i\in N,|J_i=-1$

Here, $T{R}_i$ represents the total production or consumption over the planning horizon, $S_i^0$ is the initial inventory level, and ${\overline{S}}_i$ is the maximum inventory level at port i.

Next, with Equation (28), we estimate the minimum number of visits ${\mu }_i$ required for each port i by considering the largest vessel capacity available:

(28)${\mu }_i=\lceil Q_i^N/max{C}_k\rceil \forall i\in N,k\in K$

To enforce these minimum visits, we introduce the following constraints (29) and (30):

$\begin{array}{}\mathrm{(29)}& \hfill {\displaystyle \sum _{k\in K}}{y}_{i{\mu }_{i}k}& =1\hfill & \hfill & \forall i\in N\hfill \mathrm{(30)}& \hfill {\displaystyle \sum _{n<m}}{z}_{ink}+y_{imk}& \le 1\hfill & \hfill & \forall (i,m)\in V|m>{\mu }_i,k\in K\hfill \end{array}$

Constraint (29) ensures that a visit ${\mu }_i$ is made by a vessel $k\in K$. Constraint (30) ensures that if the vessel $k\in K$ visits port i for the $m^{th}$ time, it cannot have made its last visit before in the same port.

These inequalities are based on the capacity and demand balance principles commonly used in inventory routing problems. By enforcing the minimum number of visits required to meet net demand, we eliminate infeasible or suboptimal solutions from the feasible region. This approach is supported by previous research that demonstrates the effectiveness of such constraints in strengthening optimization models [9].

3.9.2. Logic Inequalities

The logic inequalities capture the dependency between a vessel’s departure from the origin and its ability to make port visits. They ensure that a vessel cannot perform any operations unless it has first departed from the origin, which is in line with real-world operational logic.

(31)$\sum _{(i,m)\in V}{y}_{imk}\le M\sum _{(i,m)\in V}{x}_{imk}^0\forall k\in K$

Here, M represents a sufficiently large constant, such as the total number of nodes in V. For instance, if there are five ports and the maximum allowable number of visits per port is four, then V would contain 21 nodes, including the origin node, and M would be assigned the value 21.

• If vessel k does not depart from the origin $({\sum }_{(i,m)\in V}{x}_{imk}^0=0)$, then ${\sum }_{(i,m)\in V}{y}_{imk}\le 0$, implying $y_{imk}=0$ for all $(i,m)$. This ensures that no visits are scheduled for a vessel that remains at the origin.

• If vessel k leaves $({\sum }_{(i,m)\in V}{x}_{imk}^0=1)$, the constraints allow for any number of visits up to M.

By adding logical constraint (31), we eliminate infeasible scenarios where a vessel makes port visits without departing from the origin. This refinement reduces the solution space and assists the solver in converging more quickly.

Logic constraints are a standard technique in integer programming to model dependencies and conditional relationships. They have been shown to effectively reduce computational complexity by pruning infeasible or illogical solutions [39].

3.9.3. Symmetry Breaking Inequalities

Symmetry in optimization models can lead to redundant exploration of equivalent solutions, increasing computational effort without providing additional insights. Symmetry-breaking inequalities reduce this redundancy by imposing an order on homogeneous vessels.

Constraint (32) stipulates that vessel k can only leave the origin if vessel $(k-1)$ has already departed. It effectively orders the departures of identical vessels, breaking symmetry.

(32)$\sum _{(i,m)\in V}{x}_{imk}^0\le \sum _{(i,m)\in V}{x}_{im,k-1}^0\forall k\in S|k>1$

where $S\subseteq K$ is the subset of homogeneous vessels.

By reducing the number of symmetric solutions, we decrease the computational burden on the solver. These constraints help prevent the exploration of multiple permutations of vessel departures that are equivalent in terms of the objective function.

Symmetry-breaking constraints are widely recognized for improving the efficiency of solving integer programming models. They have been successfully applied in various contexts to reduce solution times without sacrificing optimality [40].

3.10. Case Study: Antarctic Logistics

To evaluate the practicality of the proposed model, we created a set of problem instances based on real-world data from Chilean Antarctic scientific bases. This case study reflects the logistical challenges of supplying research stations in a remote and environmentally sensitive region.

• 1.. Data Sources: The geographical coordinates of seven Chilean Antarctic bases were obtained using Google Earth, providing accurate distance matrices for modeling travel times and costs. Chile maintains 10 scientific bases and 2 shelters in Antarctica, and we selected 7 specific bases for this study to ensure a representative and diverse set of scenarios.

  The selection was not random; the chosen bases are strategically distributed across key locations, including the South Shetland Islands (e.g., Base Escudero) and continental Antarctica (e.g., Base Yelcho). This distribution allows the model to account for both short-distance travels within the islands and long-distance travels between the islands and continental Antarctica, reflecting realistic logistical challenges. The locations of the selected bases are shown in Figure 3. These bases are Base Yelcho, Base Presidente Eduardo Frei Montalva, Base Escudero, Risopatrón Base, Captain Arturo Prat Base, Doctor Guillermo Mann Base, and Base General Bernardo O’Higgins.

  Logistical requirements, including cargo demand, were estimated based on the daily meal needs of the scientific personnel at each base. Each time period in the model represents one day, aligning with the planning horizon of a typical scientific mission in Antarctica (approximately four months for a 120-day horizon).

• 2.. Instance Generation: The data were used to define the parameters of the problem, such as the number of ports (Antarctic bases), the capacities of the vessels, the production and consumption rates, and the levels of inventory. The generated instances reflect realistic logistical scenarios that capture the unique operational conditions of Antarctic supply chains.

• 3.. Model Application: The Antarctic instances were solved using the proposed MILP model to optimize multi-vessel routing and inventory management over long planning horizons. The results, including computational performance and optimality gaps, are presented in Section 4, which compares the performance of weak and strong formulations.

• 4.. Assumptions: For simplicity, the analysis assumes constant production and consumption rates at the bases, as well as fixed vessel capacities. These assumptions allow the model to focus on routing and inventory decisions while maintaining computational tractability.

Figure 3

Locations of the seven Chilean Antartic scientific bases considered in this work.

[Figure omitted. See PDF]

4. Results

We present the experimental results of our mathematical model proposed for the MV-MIRP. Our experimental setup consists of two phases. In the first phase, we compare the elapsed times of our model with those of the literature [10], as shown in Table 2, Table 3 and Table 4. Although we had access to key data from the instances provided by Agra et al. (2022)—including the distance matrix, ship speeds, and supply and demand at each port—shared with us upon request, we did not have access to the exact parameters of the objective function. As a result, we were unable to compute an initial lower bound value or a gap in this phase. In the second phase, we apply our model to our original instances, whose results are shown in Table 5, Table 6, Table 7 and Table 8. The results show that our method can solve most large instances optimally and in a computationally efficient time.

We implemented our approach with Gurobi 10 optimization software on a Windows 10 machine with a 12-core and 24-thread Intel Xeon E5-2696v2 at 2.5 GHz CPU, and 16 GB of RAM at 1866 MHz. We selected the default tolerance level of ${10}^{-6}$ and set a maximum solver time of 3600 s. We recorded the elapsed time taken by Gurobi 10 to find an optimal solution or prove infeasibility on each MV-MIRP instance, providing insights into our methodology’s computational complexity and efficiency. By following this setup, we expected to deliver a high-quality analysis of the performance of our formulation to solve problem instances of different sizes.

4.1. Model Performance Comparison

Table 2 provides an overview of the time elapsed for each instance of both studies, over a time horizon of $T=60$. The first column of the table indicates the name of each instance. The following two columns show the number of loading ports and the number of unloading ports. We compare our work with three different methods utilized in [10]: (i) branch-and-cut with a weak formulation (BCw), (ii) branch-and-cut with a tightened formulation (BCs), and (iii) a custom algorithm developed by the authors (ALG1). The results are given in the following three columns. The results of this work (TW) are presented in the last column.

Table 2

Elapsed time comparison between branch-and-cut weak formulation (BCw), branch-and-cut tightened formulation (BCs), custom algorithm (ALG1), and this work (TW) for a time horizon $T=60$.

Instance	Loading Ports	Unloading Ports	Elapsed Time [s]
			BCw	BCs	ALG1	TW
A1	1	3	31.6	21.7	20.4	1.4
A2	1	3	35.2	21.7	17.3	1.6
B1	1	2	11.5	7.7	8.0	1.4
B2	1	2	13.3	9.2	9.0	2.0
C1	2	2	20.5	19.9	29.1	1.4
C2	2	2	20.4	14.0	16.5	3.6
D1	2	3	89.8	32.5	24.3	2.6
D2	2	3	195.0	37.2	37.0	5.2
E1	2	3	56.7	30.2	59.4	0.8
E2	2	3	166.0	43.0	72.2	0.8
F1	2	2	31.4	19.0	23.8	0.3
F2	2	2	27.8	16.0	17.6	1.0
G1	3	3	20.5	11.5	13.0	3.9
G2	3	3	21.0	13.7	13.2	2.0

The custom algorithm (ALG1) is a two-level branching approach to solve the MIRP to optimality. The main idea is to divide the problem into two subproblems based on the total number of visits by the vessel. The first subproblem has a constraint that limits the number of visits to be less than or equal to half of a given upper bound, while the second subproblem has a constraint that requires the number of visits to be greater than or equal to half of the upper bound plus one. The first subproblem is solved using branch-and-cut considering a strong formulation with valid inequalities (BCs). The best feasible solution found is used as a cut-off value for the second subproblem. The algorithm returns the best solution among the two subproblems. Since each subproblem has a time limit of 3600 s, the entire algorithm has a maximum time limit of 7200 s.

We present the results for a longer planning horizon ($T=120$) in Table 3. This table shows the name of the instance, the number of loading ports (LP), the number of unloading ports (UP), and the solution time for each method. If the optimal solution is not found within 3600 [s], we report the % gap between the best value of the upper bound (UB) and the best value of the lower bound (LB), where % gap = ((UB − LB)/UB) 100. We label instances that were proven to be unfeasible within the time limit as (UNF). We use a hyphen to denote inconclusive instances, which means that the model does not reach optimality, find a gap, or prove infeasibility. We do not use any tightening constraints for these model tests. For the time horizon $T=60$, we obtain optimal solutions for all instances, with significantly shorter solving times than those reported in the literature. For the $T=120$ time horizon, our approach outperforms previous studies by solving nine instances to optimality, while only five were previously solved. Moreover, we observe a gap in two instances, whereas they observed a gap in three. Four instances were unresolved in the previous study and zero in ours. The shorter running times achieved with Gurobi and a high-speed processor highlight the importance of using efficient computational tools.

Table 3

Performance comparison between branch-and-cut weak formulation (BCw), branch-and-cut tightened formulation (BCs), custom algorithm (ALG1), and this work (TW) for a time horizon $T=120$.

Instance	LP	UP	Elapsed Time [s]				% Gap
			BCw	BCs	ALG1	TW	BCw	BCs	ALG1	TW
A1	1	3	3600.0	3600.0	3922.5	460.8	-	-	-	0.0
A2	1	3	3600.0	3600.0	5113.1	883.5	-	-	19.0	0.0
B1	1	2	3600.0	2870.9	3619.4	1722.4	-	UNF	-	UNF
B2	1	2	195.0	19.5	21.0	9.8	0.0	0.0	0.0	0.0
C1	2	2	204.0	104.2	115.6	10.6	0.0	0.0	0.0	0.0
C2	2	2	2907.0	242.2	320.9	47.5	0.0	0.0	0.0	0.0
D1	2	3	3600.0	3600.0	3736.1	3600.0	-	18.9	-	8.1
D2	2	3	3600.0	3600.0	3698.9	3600.0	-	-	-	21.8
E1	2	3	3600.0	2795.0	5566.7	80.2	-	0.0	0.0	0.0
E2	2	3	3600.0	3600.0	7200.0	2.8	-	18.9	21.5	0.0
F1	2	2	3600.0	180.2	210.7	1.2	-	0.0	0.0	0.0
F2	2	2	3600.0	3.9	3.4	14.2	-	UNF	UNF	UNF
G1	3	3	3600.0	617.8	198.2	18.8	-	UNF	UNF	UNF
G2	3	3	3600.0	3600.0	6451.4	6.9	-	-	0.8	0.0

We compare the performance of three formulations of our work in the same instances: (i) a weak formulation without valid inequalities (WF), (ii) a strong formulation with valid inequalities (29) and (30), which consider the minimum number of visits to a port (SA); and (iii) a strong formulation with logic and symmetry-breaking valid inequalities (31) and (32) (SO). Table 4 shows the results. The first column shows the name of the instance. The next three columns show the elapsed time in seconds for each formulation. The SO formulation is the fastest on four instances, and the SA formulation is the fastest on instances E1 to G2. The WF formulation is faster only in instances A1 and C1. Columns 5 to 7 show the nodes explored by each formulation. Usually, fewer explored nodes mean less computational time, but this is not true for instances B1 and C1. The SA formulation shows the fewer nodes explored in two instances, while the SO formulation does so in four. Columns 8 to 10 show the initial lower bound of the linear relaxation for each formulation. The SA formulation always improves the initial lower bound of the WF formulation, while the SO formulation does so only in instances B1 and B2.

Table 4

Performance comparison between formulations WF, SA, and SO for time horizon $T=60$.

Instance	Elapsed Time [s]			Explored Nodes			Initial Lower Bound
	WF	SA	SO	WF	SA	SO	WF	SA	SO
A1	1.4	3.7	4.1	1179	4426	3871	601,510	683,252	601,510
A2	1.6	1.2	1.2	2145	1243	1015	412,552	553,423	412,552
B1	1.4	1.4	1.4	2030	793	1176	791,184	800,783	796,769
B2	2.0	1.6	1.1	2342	2461	973	859,426	961,331	865,491
C1	1.4	1.5	1.8	215	292	185	413,037	426,614	413,037
C2	3.6	1.7	2.7	2614	373	1611	509,115	524,363	509,115
D1	2.6	1.4	2.0	96	1	1	449,030	644,628	449,030
D2	5.2	5.0	3.7	2040	2506	959	580,720	713,767	580,720
E1	0.8	0.5	0.9	1	1	1	231,233	433,827	231,233
E2	0.8	0.6	1.2	1	1	1	8451	352,682	8451
F1	0.3	0.2	0.2	1	1	1	79,196	156,100	79,196
F2	1.0	0.2	1.2	1	1	1	200,225	390,100	200,225
G1	3.9	1.3	4.8	1	1	1	140,548	720,611	140,548
G2	2.0	1.1	2.6	1	1	1	120,296	583,900	120,296

4.2. Model Performance with New Instances

To evaluate the practical applicability of our MILP formulation, we create a set of new MV-MIRP instances using real-world data from Chilean scientific bases in Antarctica. These bases require regular deliveries of supplies to support their research activities, making them ideal candidates for our case study.

We use Google Earth to obtain the geographical coordinates of seven Antarctic Chilean bases. The coordinates allow us to calculate accurate distance matrices between the bases, which serve as ports in our model. The cargo at each base is estimated on the basis of the number of meals required to sustain the scientific personnel. We assume an average number of personnel at each base and calculate daily meal requirements, which is translated into cargo demand in our model. This approach reflects the logistical needs of supplying food and essentials to the bases. We consider each time period as one day. The longer planning horizon of 120 time periods corresponds to the typical duration of a scientific mission in Antarctica, which spans approximately four months.

We report the experimental results for the new set of problem instances on a time horizon of $T=60$ in Table 5. The name of the instance, the number of loading and unloading ports, and the number of vessels are given in columns 1 to 4, respectively. Columns 5 to 7 show the running times for each formulation until optimality or the time limit is reached. Columns 8 to 10 show the percentage gap between the best value and the lower bounds for each formulation, or a hyphen if no gap was achieved.

The computational performance of the solver varies depending on the characteristics of the instances. We notice that larger instances (from 1 to 6) generally require more time, nodes, and iterations than smaller ones. For example, while instance 1 is solved optimally in less than a second, instance 6 only reaches a gap after an hour of computation. However, instance 6b shows that the size of the instance is not the only factor affecting the solver efficiency, as it obtains an optimal solution in less than 5 s for $T=60$. Therefore, other aspects, such as production/consumption rates, vessel capacities, and inventory levels at the bases, also play a significant role in the solver performance.

Table 5

Elapsed time and percentage gap for MV-MIRP on original instances for $T=60$.

Instance	LP	UP	V	Elapsed Time [s]			% Gap
				WF	SA	SO	WF	SA	SO
1	1	3	1	0.4	0.2	0.6	0.0	0.0	0.0
1b	1	3	1	0.2	0.2	0.5	0.0	0.0	0.0
2	1	2	2	0.8	0.0	0.6	0.0	0.0	0.0
2b	1	2	2	1.8	0.4	1.6	0.0	0.0	0.0
3	2	2	2	7.1	6.7	0.6	0.0	0.0	0.0
3b	2	2	2	6.9	5.9	0.9	0.0	0.0	0.0
4	2	3	2	585.0	2452.6	1.9	0.0	0.0	0.0
4b	2	3	2	894.1	377.7	1.7	0.0	0.0	0.0
5	2	3	3	2049.2	3600.0	374.4	0.0	15.6	0.0
5b	2	3	3	2615.0	801.7	318.9	0.0	0.0	0.0
6	2	4	2	3600.0	28.1	3600.0	24.7	0.0	24.9
6b	2	4	2	124.6	224.4	9.4	0.0	0.0	0.0
7	2	5	3	1146.2	726.5	389.6	0.0	0.0	0.0
7b	2	5	3	3600.0	956.3	133.2	5.9	0.0	0.0

We present additional results on node exploration and the initial lower bound for each instance in Table 6. The table lists the name of the instance, the number of loading ports (LP), the number of unloading ports (UP), and the number of vessels (V) in the first four columns. The next three columns show the nodes explored by Gurobi with different formulations, and the last three columns show the initial lower bound obtained with each formulation. The SA formulation reduces node exploration in two instances, but improves the initial lower bound in all instances. The SO formulation reduces node exploration in 12 instances, but improves the initial lower bound only in 4 instances.

Table 6

Explored nodes and initial lower bound for MV-MIRP on original instances for $T=60$.

Instance	LP	UP	V	Explored Nodes			Initial Lower Bound
				WF	SA	SO	WF	SA	SO
1	1	3	1	1	1	1	261,558	315,411	261,558
1b	1	3	1	1	1	1	285,973	337,822	285,973
2	1	2	2	436	1	1	374,223	649,000	374,223
2b	1	2	2	2331	1	1242	687,550	711,333	687,550
3	2	2	2	18,185	37,813	1	26,570	88,900	30,231
3b	2	2	2	14,968	26,280	1	26,787	88,900	31,109
4	2	3	2	828,851	363,8473	1	31,542	49,100	34,828
4b	2	3	2	1,042,948	610,822	1	14,438	34,423	17,384
5	2	3	3	1,630,249	3,390,406	363,247	86,074	202,829	86,074
5b	2	3	3	2,105,200	585,315	333,075	116,571	214,369	116,571
6	2	4	2	1,317,710	23,196	1,396,386	165,538	600,939	165,538
6b	2	4	2	175,889	312,753	12,237	229,373	300,622	229,373
7	2	5	3	731,809	448,182	302,965	156,077	424,932	156,077
7b	2	5	3	1,453,616	749,996	74,234	99,454	196,153	99,454

The results in Table 7 show that increasing the time horizon from $T=60$ to $T=120$ has a significant impact on the computational performance of the model. Most instances require much longer running times to reach optimality or terminate without finding an optimal solution in 3600 s. Only instance 3 has a negligible difference in runtime, suggesting that the complexity and difficulty of the problem vary depending on the instance characteristics. Comparison with Table 5 confirms that all instances are more difficult to solve with a longer time horizon and that only five instances achieve optimality in this setting. These findings highlight the trade-off between the accuracy and tractability of the model and the need for further research on efficient solution methods for the MIRP.

Table 7

Elapsed time and percentage gap for MV-MIRP on original instances for $T=120$.

Instance	LP	UP	V	Elapsed Time [s]			% Gap
				WF	SA	SO	WF	SA	SO
1	1	3	1	3600.0	3600.0	3600.0	INC	INC	INC
1b	1	3	1	3600.0	3600.0	3600.0	INC	INC	INC
2	1	2	2	1775.4	1266.8	1185.7	0.0	0.0	0.0
2b	1	2	2	42.0	114.3	27.5	0.0	0.0	0.0
3	2	2	2	11.5	8.7	7.2	0.0	0.0	0.0
3b	2	2	2	15.2	13.4	3.1	0.0	0.0	0.0
4	2	3	2	3600.0	3600.0	3154.0	17.8	14.6	0.0
4b	2	3	2	3600.0	3600.0	3600.0	13.8	12.7	3.3
5	2	3	3	3600.0	3600.0	3600.0	7.8	13.5	7.8
5b	2	3	3	3600.0	309.9	3600.0	21.5	0.0	29.1
6	2	4	2	3600.0	419.1	3600.0	6.1	0.0	7.2
6b	2	4	2	3600.0	3600.0	3600.0	4.5	1.6	4.6
7	2	5	3	3600.0	3600.0	3600.0	INC	INC	INC
7b	2	5	3	3600.0	3600.0	3600.0	20.0	24.7	22.1

Table 8 shows the number of nodes explored and the initial lower bound values for the original instances with a time horizon of $T=120$. The table has the same columns as Table 6. The results suggest that longer time horizons increase the difficulty of solving the instances, as the number of nodes explored increases significantly compared to Table 6. The SA formulation outperforms the other formulations in seven instances, by exploring fewer nodes and enhancing all initial lower bound values. The SO formulation also explores fewer nodes in seven instances, but only enhances the initial lower bound value in three instances.

Table 8

Explored nodes and initial lower bound for MV-MIRP on original instances for $T=120$.

Instance	LP	UP	V	Explored Nodes			Initial Lower Bound
				WF	SA	SO	WF	SA	SO
1	1	3	1	4,435,944	8,900,858	1,890,043	3,036,711	3,043,122	3,036,711
1b	1	3	1	7,715,299	5,524,034	8,524,051	2,955,044	2,956,662	2,955,044
2	1	2	2	6,755,676	5,945,934	1,323,876	2,289,638	2,525,500	2,289,638
2b	1	2	2	77,532	381,348	43,824	2,899,500	2,963,133	2,899,500
3	2	2	2	7399	5566	1892	164,670	182,612	196,086
3b	2	2	2	17,407	13,489	693	164,887	182,804	196,964
4	2	3	2	2,646,762	1,832,312	2,603,886	279,596	390,409	287,353
4b	2	3	2	2,523,325	2,131,800	2,927,427	289,798	358,631	289,798
5	2	3	3	2,205,410	2,646,163	1,756,032	651,129	830,380	651,129
5b	2	3	3	1,996,743	282,437	1,417,788	772,523	1,230,566	772,523
6	2	4	2	1,915,121	299,373	1,664,805	1,382,327	1,946,887	1,382,327
6b	2	4	2	2,129,875	1,680,105	1,950,939	1,438,538	1,942,467	1,438,538
7	2	5	3	931,465	1,181,111	982,733	1,782,577	2,326,201	1,782,577
7b	2	5	3	716,540	688,821	760,832	1,196,418	1,468,096	1,196,418

5. Discussion and Future Research

In this section, we delve into three critical aspects of the proposed MILP model and its applicability to the MV-MIRP. Section 5.1 addresses the scalability of the model, exploring its computational performance as the complexity of the problem increases, and discussing potential strategies to improve the scalability. Section 5.2 considers the incorporation of uncertainty into the problem, emphasizing the challenges posed by dynamic and stochastic elements in real-world maritime logistics and suggesting future directions for extending the model. Finally, Section 5.3 discusses environmental considerations, highlighting the potential to integrate sustainability objectives, such as reducing emissions and optimizing fuel efficiency, into the model. Together, these subsections provide a comprehensive discussion on improving the practical relevance, robustness and alignment of the model with modern logistical and environmental demands.

5.1. Scalability

The scalability of our MILP model is a critical factor for its practical applicability in real-world maritime logistics. As the number of vessels, ports, and time units increases, the computational complexity of the model grows significantly. Specifically, the number of binary variables related to routing decisions and scheduling increases combinatorially with the addition of more vessels and ports. This growth leads to longer computational times and a greater demand for memory resources.

Our computational experiments, presented in Table 7 and Table 8, demonstrate that instances with larger sizes or longer planning horizons (e.g., $T=120$) require substantially more computational time, and in some cases, the solver could not find optimal solutions within the allotted time. For example, instances 7 and 7b showed increased difficulty as the number of ports and vessels increased, indicating that scalability is a limiting factor of our current approach.

To address scalability concerns and provide a broader perspective on our model’s performance relative to other solution techniques, we compared our results with other works that use metaheuristics to solve related problems. Ref. [35] addressed the IRP using Simulated Annealing, solving instances with up to 50 customers, five potential vehicles, and six time periods. For the longest time horizon of six time periods, they achieved a best optimality gap of 0.66% for instance H6 with one potential vehicle. Ref. [37] tackled a Multi-Product Multi-Vehicle Inventory Routing Problem using a Modified Adaptive Genetic Algorithm, solving instances with up to 30 customers, five vehicles, and seven time periods. They reported an average “closeness” value of 29.05, a performance metric they used to evaluate their results, but did not reach optimality in any instance. Ref. [36] solved instances of a location-inventory-routing problem with up to 50 customers and five time periods using a Genetic Algorithm; however, the number of vehicles was not explicitly stated.

Although our work encounters gaps of up to 24.9% on the largest tested instance 7b, we consider instances with up to 120 time periods, which is over 17 times longer than the longest time horizon considered by [37]. This demonstrates that our exact optimization approach can handle significantly longer planning horizons than those addressed by metaheuristic methods, albeit with a trade-off in computational performance and solution optimality for larger instances.

Although metaheuristic approaches have shown promise in addressing large-scale problems with acceptable computational times, they often sacrifice optimality guarantees. Our exact MILP model provides optimal solutions for small to medium-sized instances and longer planning horizons, making it a valuable tool for strategic decision-making where accuracy is paramount. Future research should focus on integrating advanced decomposition techniques, parallel computing, and hybrid methods to address scalability challenges and extend the model’s applicability to even larger instances. In addition to tackling computational complexity, addressing the inherent uncertainty present in real-world maritime logistics is crucial to enhancing the robustness of the model.

5.2. Uncertainty Management

An important avenue for future research is the incorporation of dynamic and stochastic elements into the MV-MIRP. In real-world maritime logistics, parameters such as demand in ports, production and consumption rates, vessel capacities, and operational conditions are rarely constant and often subject to uncertainty. Factors such as time-varying demand, fluctuating supply rates, and unforeseen events such as vessel breakdowns or adverse weather conditions can significantly impact routing and inventory decisions.

Extending our deterministic model to accommodate these uncertainties would enhance its practical applicability and robustness. One potential approach is to develop a stochastic programming version of the MV-MIRP, where uncertain parameters are modeled as random variables with known probability distributions. This would allow for the optimization of expected costs or service levels while considering the variability in demand and supply. Alternatively, robust optimization techniques could be employed to find solutions that perform well under the worst-case scenarios of uncertainty, providing a safeguard against extreme events.

Incorporating dynamic elements, such as the time-varying demand and supply rates, would require adapting the model to handle parameters that change over time within the planning horizon. This could involve segmenting the time horizon into smaller intervals where the parameters are assumed to be constant or developing continuous-time models that can capture dynamic behaviors more effectively.

Implementing these extensions poses several challenges, including increased computational complexity and the need for advanced solution methodologies. Techniques such as scenario decomposition, sampling methods, or hybrid approaches combining exact and heuristic algorithms might be necessary to solve the stochastic MV-MIRP efficiently.

5.3. Environmental Considerations

An important aspect of modern maritime logistics is the consideration of environmental sustainability, particularly in reducing emissions and optimizing fuel efficiency. While our current model focuses on minimizing total costs—which inherently accounts for fuel consumption through variable travel costs—it does not explicitly incorporate environmental objectives or constraints. Integrating environmental considerations into the MV-MIRP could involve introducing emission factors associated with fuel consumption, vessel speed, and cargo load, as suggested by [12].

For example, the model could be extended to minimize greenhouse gas emissions by incorporating emission coefficients into the objective function or by adding constraints that limit the total emissions over the planning horizon. Furthermore, optimizing vessel speeds for fuel efficiency, as explored by [21], could be integrated into the model by treating speed as a decision variable and modeling its nonlinear relationship with fuel consumption and emissions.

Incorporating environmental sustainability into the MV-MIRP would enable practitioners to optimize both operational efficiency and ecological impact, aligning with growing global regulatory demands for greener shipping practices. Future work could integrate emission factors, fuel efficiency, and vessel speeds into the model, contributing to a reduction in carbon footprints while maintaining cost-effectiveness. By combining sustainability objectives with scalability and uncertainty management, our model can evolve into a comprehensive tool that balances economic, operational, and environmental considerations, thus addressing the multifaceted challenges of modern maritime logistics.

6. Conclusions

Our research extends the best-performance continuous-time model for the MIRP, introducing several novel features. Unlike the base model, which allows only the use of a single vessel, our model enables the routing of a fleet of heterogeneous vessels. Furthermore, our model incorporates an optimization component for the size of the fleet, selecting the combination of vessels that minimizes the objective value. To enhance the efficiency of the model, we applied tightening constraints based on valid inequalities from the existing literature and also proposed two new valid inequalities. Using these cuts, we improved the initial lower bounds and the efficiency of the model. To evaluate the effectiveness of our extended model, we created a new set of instances based on real-world data on Chilean scientific bases in Antarctica, obtained from Google Earth. These instances reflect the logistical challenges faced by Antarctic bases, such as long distances, and the need for efficient resource allocation to maintain an appropriate inventory level of meals at all times. Through these enhancements and the use of real-world data, our model represents a significant advance in addressing the MIRP. These advancements provide maritime logistics practitioners with a more powerful tool to optimize fleet routing and inventory management, reducing both operational costs and environmental impact.

To evaluate the performance of our model, we tested it on reference instances from the literature and obtained computational times lower than those reported in previous studies. Since the exact parameters of the reformulated instances for one vessel are not available, we cannot make a direct comparison of the results. However, we can demonstrate that computational times can be significantly reduced by using more advanced hardware and a commercial solver, such as Gurobi.

We ran the model with our original instances, obtaining results that met our expectations based on the results of the instances in the literature, which are presented in Table 5, Table 6, Table 7 and Table 8. Using the weak formulation, we could solve all but one instance to optimality within a time horizon of 60 time units, while the remaining instance had a small gap. However, when we extended the time horizon to 120 time units, the weak formulation became less effective and could only solve four instances to optimality. The strong formulations (SA and SO) performed better and solved six and five instances to optimality, respectively. In total, we solved seven instances to optimality using SA or SO. For the remaining instances, we either obtained a gap (five instances) or could not determine the feasibility or optimality (three instances). We also observed that the strong formulations reduced the computational time or the gap for most instances after 3600 s, except for instance 2b, where the SA formulation increased the elapsed time.

The authors acknowledge that the computational complexity of the instance, which includes the number of ports and vessels, cannot be determined solely by its size. We observe that instance 6b had a faster running time considering a time horizon of 60 time units, compared with smaller instances, such as 3 and 3b. Therefore, the smallest instances (1 and 1b) exhibited the highest complexity with a time horizon of 120 time units, making it challenging to achieve a gap or demonstrate infeasibility.

We evaluated the performance of MV-MIRP, a method for solving the multi-vessel MIRP with long time horizons and constant port rates, using a continuous-time formulation. Our results showed that MV-MIRP was effective, fast, and reliable. This method has theoretical and practical implications for researchers and practitioners in the shipping industry, as it can handle complex problems with fewer computational and timely resources. It also opens new avenues for future research on stochasticity, dynamic effects, climate prediction, etc.

In addition to the contributions presented in this paper, future research could further explore the integration of stochastic elements and environmental sustainability, as discussed. These enhancements would extend the model’s applicability to real-world dynamic scenarios, making it even more robust for maritime logistics challenges.

Author Contributions

All authors contributed to the conception and design of the study. Conceptualization, formal analysis, investigation, and methodology were performed by D.C.-L., L.P. and V.P.; data curation was performed by D.C.-L. and V.P.; software, validation, and visualization were performed by D.C.-L.; the first draft of the manuscript was written by D.C.-L.; and the review and editing processes were performed by L.P. and V.P. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

Article Processing Charges were supported by Universidad Católica de la Santísima Concepción, grants Proyectos Ingeniería 2030 (ING222010004) and InES Ciencia Abierta (INCA210005).

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

IRP	Inventory Routing Problem
MILP	Mixed-Integer Linear Programming
MIRP	Maritime Inventory Routing Problem
MV-MIRP	Multi-Vessel Maritime Inventory Routing Problem

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. A geographical representation of the three ports in the MV-MIRP example.

Figure 2. An expanded graph representation of the MV-MIRP example with two vessels, including the origin and destination nodes.

Table 1

Literature review comparison.

Reference	Time Horizon		TH Length		Algorithm		Vessels
	D	C	≤60	>60	E	M	1	>1
[1]	✓		✓		✓			✓
[9]	✓	✓	✓		✓			✓
[18]	✓	✓	✓		✓			✓
[10]		✓		✓	✓		✓
[11]		✓	✓		✓			✓
[7]	✓		✓		✓			✓
[21]	✓		✓			✓		✓
[12]		✓	✓		✓			✓
[23]	✓			✓		✓		✓
[27]	✓		✓			✓		✓
[13]		✓	✓			✓		✓
[8]	✓		✓		✓			✓
[22]	✓		✓			✓		✓
[3]	✓		✓		✓			✓
[14]		✓	✓		✓			✓
This work		✓		✓	✓			✓

References

1. Agra, A.; Christiansen, M.; Delgado, A. Mixed integer formulations for a short sea fuel oil distribution problem. Transp. Sci.; 2013; 47, pp. 108-124. [DOI: https://dx.doi.org/10.1287/trsc.1120.0416]

2. Prause, F.; Prause, G.; Philipp, R. Inventory Routing for Ammonia Supply in German Ports. Energies; 2022; 15, 6485. [DOI: https://dx.doi.org/10.3390/en15176485]

3. Yusuf, F.K.; Ridwan, A.Y.; Pambudi, H.K. Maritime Inventory Routing Problem: Application on Discharge the Load of the Ship in Cement Companies to Minimize the Total Transportation Cost. IOP Conf. Ser. Mater. Sci. Eng.; 2020; 982, 012056. [DOI: https://dx.doi.org/10.1088/1757-899X/982/1/012056]

4. Montagné, R.; Gamache, M.; Gendreau, M. A shortest path-based algorithm for the inventory routing problem of waste vegetable oil collection. J. Oper. Res. Soc.; 2019; 70, pp. 986-997. [DOI: https://dx.doi.org/10.1080/01605682.2018.1476801]

5. Archetti, C.; Speranza, M.G. The inventory routing problem: The value of integration. Int. Trans. Oper. Res.; 2016; 23, pp. 393-407. [DOI: https://dx.doi.org/10.1111/itor.12226]

6. Bertazzi, L.; Speranza, M.G. Inventory routing problems: An introduction. EURO J. Transp. Logist.; 2012; 1, pp. 307-326. [DOI: https://dx.doi.org/10.1007/s13676-012-0016-7]

7. Baller, A.C.; Dabia, S.; Desaulniers, G.; Dullaert, W.E. The Inventory Routing Problem with Demand Moves. Oper. Res. Forum; 2021; 2, 6. [DOI: https://dx.doi.org/10.1007/s43069-020-00042-z]

8. Skålnes, J.; Andersson, H.; Desaulniers, G.; Stålhane, M. An improved formulation for the inventory routing problem with time-varying demands. Eur. J. Oper. Res.; 2022; 302, pp. 1189-1201. [DOI: https://dx.doi.org/10.1016/j.ejor.2022.02.011]

9. Agra, A.; Andersson, H.; Christiansen, M.; Wolsey, L. A maritime inventory routing problem: Discrete time formulations and valid inequalities. Networks; 2013; 62, pp. 297-314. [DOI: https://dx.doi.org/10.1002/net.21518]

10. Agra, A.; Christiansen, M.; Wolsey, L. Improved models for a single vehicle continuous-time inventory routing problem with pickups and deliveries. Eur. J. Oper. Res.; 2022; 297, pp. 164-179. [DOI: https://dx.doi.org/10.1016/j.ejor.2021.04.027]

11. Algendi, A.; Urrutia, S.; Hvattum, L.M. Optimizing production levels in maritime inventory routing with load-dependent speed optimization. Flex. Serv. Manuf. J.; 2023; 35, pp. 111-141. [DOI: https://dx.doi.org/10.1007/s10696-022-09460-z]

12. Eide, L.; Ardal, G.C.H.; Evsikova, N.; Hvattum, L.M.; Urrutia, S. Load-dependent speed optimization in maritime inventory routing. Comput. Oper. Res.; 2020; 123, 105051. [DOI: https://dx.doi.org/10.1016/j.cor.2020.105051]

13. Shaabani, H.; Hoff, A.; Hvattum, L.M.; Laporte, G. A matheuristic for the multi-product maritime inventory routing problem. Comput. Oper. Res.; 2023; 154, 106214. [DOI: https://dx.doi.org/10.1016/j.cor.2023.106214]

14. Zojaji, A.; Soltaniani, K.; Hvattum, L.M.; Urrutia, S. Cyclic solutions to a maritime inventory routing problem. Marit. Transp. Res.; 2022; 3, 100074. [DOI: https://dx.doi.org/10.1016/j.martra.2022.100074]

15. Fagerholt, K.; Hvattum, L.M.; Papageorgiou, D.J.; Urrutia, S. Maritime inventory routing: Recent trends and future directions. Int. Trans. Oper. Res.; 2023; 30, pp. 3013-3056. [DOI: https://dx.doi.org/10.1111/itor.13313]

16. Christiansen, M.; Fagerholt, K. Ship Routing and Scheduling in Industrial and Tramp Shipping. Vehicle Routing; MOS-SIAM Series on Optimization; SIAM: University City, PA, USA, 2014; Chapter 13 pp. 381-408. [DOI: https://dx.doi.org/10.1137/1.9781611973594.ch13]

17. Andersson, H.; Hoff, A.; Christiansen, M.; Hasle, G.; Løkketangen, A. Industrial aspects and literature survey: Combined inventory management and routing. Comput. Oper. Res.; 2010; 37, pp. 1515-1536. [DOI: https://dx.doi.org/10.1016/j.cor.2009.11.009]

18. Agra, A.; Christiansen, M.; Delgado, A. Discrete time and continuous time formulations for a short sea inventory routing problem. Optim. Eng.; 2017; 18, pp. 269-297. [DOI: https://dx.doi.org/10.1007/s11081-016-9319-0]

19. Liu, B.; Zhang, Q.; Yuan, Z. Two-stage distributionally robust optimization for maritime inventory routing. Comput. Chem. Eng.; 2021; 149, 107307. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2021.107307]

20. Rodrigues, F.; Agra, A.; Christiansen, M.; Hvattum, L.M.; Requejo, C. Comparing techniques for modelling uncertainty in a maritime inventory routing problem. Eur. J. Oper. Res.; 2019; 277, pp. 831-845. [DOI: https://dx.doi.org/10.1016/j.ejor.2019.03.015]

21. De, A.; Mamanduru, V.K.R.; Gunasekaran, A.; Subramanian, N.; Tiwari, M.K. Composite particle algorithm for sustainable integrated dynamic ship routing and scheduling optimization. Comput. Ind. Eng.; 2016; 96, pp. 201-215. [DOI: https://dx.doi.org/10.1016/j.cie.2016.04.002]

22. Song, J.H.; Furman, K.C. A maritime inventory routing problem: Practical approach. Comput. Oper. Res.; 2013; 40, pp. 657-665. [DOI: https://dx.doi.org/10.1016/j.cor.2010.10.031]

23. Papageorgiou, D.J.; Cheon, M.S.; Harwood, S.; Trespalacios, F.; Nemhauser, G.L. Recent progress using matheuristics for strategic maritime inventory routing. Intell. Syst. Ref. Libr.; 2018; 131, pp. 59-94. [DOI: https://dx.doi.org/10.1007/978-3-319-61801-2_3]

24. Misra, S.; Kapadi, M.; Gudi, R.D. Hybrid Time-Based Framework for Maritime Inventory Routing Problem. Ind. Eng. Chem. Res.; 2020; 59, pp. 20394-20409. [DOI: https://dx.doi.org/10.1021/acs.iecr.0c03186]

25. Msakni, M.K.; Haouari, M. Short-term planning of liquefied natural gas deliveries. Transp. Res. Part C Emerg. Technol.; 2018; 90, pp. 393-410. [DOI: https://dx.doi.org/10.1016/j.trc.2018.03.013]

26. Wolsey, L. Strong Valid Inequalities. Integer Programming; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2020; Chapter 9 pp. 167-194. [DOI: https://dx.doi.org/10.1002/9781119606475.ch9]

27. Rakke, J.G.; Andersson, H.; Christiansen, M.; Desaulniers, G. A New Formulation Based on Customer Delivery Patterns for a Maritime Inventory Routing Problem. Transp. Sci.; 2015; 49, pp. 384-401. [DOI: https://dx.doi.org/10.1287/trsc.2013.0503]

28. Rocha, R.; Grossmann, I.E.; de Aragão, M.V.S.P. Cascading Knapsack Inequalities: Reformulation of a crude oil distribution problem. Ann. Oper. Res.; 2013; 203, pp. 217-234. [DOI: https://dx.doi.org/10.1007/s10479-011-0857-8]

29. Wang, S.; Chu, F. A Decomposition-Based Heuristic Method for Inventory Routing Problem. IEEE Trans. Intell. Transp. Syst.; 2022; 23, pp. 18352-18360. [DOI: https://dx.doi.org/10.1109/TITS.2022.3170569]

30. Alkaabneh, F.; Diabat, A.H.; Gao, H.O. Benders decomposition for the inventory vehicle routing problem with perishable products and environmental costs. Comput. Oper. Res.; 2020; 113, 104751. [DOI: https://dx.doi.org/10.1016/j.cor.2019.07.009]

31. Lefever, W.; Touzout, F.A.; Hadj-Hamou, K.; Aghezzaf, E. Benders’ decomposition for robust travel time-constrained inventory routing problem. Int. J. Prod. Res.; 2021; 59, pp. 342-366. [DOI: https://dx.doi.org/10.1080/00207543.2019.1695167]

32. Cordeau, J.; Laganà, D.; Musmanno, R.; Vocaturo, F. A decomposition-based heuristic for the multiple-product inventory-routing problem. Comput. Oper. Res.; 2015; 55, pp. 153-166. [DOI: https://dx.doi.org/10.1016/j.cor.2014.06.007]

33. Adulyasak, Y.; Cordeau, J.; Jans, R. Formulations and Branch-and-Cut Algorithms for Multivehicle Production and Inventory Routing Problems. INFORMS J. Comput.; 2012; 26, pp. 103-120. [DOI: https://dx.doi.org/10.1287/ijoc.2013.0550]

34. Schenekemberg, C.M.; Guimarães, T.A.; Chaves, A.A.; Coelho, L.C. A Three-Front Parallel Branch-and-Cut Algorithm for Production and Inventory Routing Problems. Transp. Sci.; 2023; 58, pp. 687-707. [DOI: https://dx.doi.org/10.1287/trsc.2022.0261]

35. Alvarez, A.; Munari, P.; Morabito, R. Iterated local search and simulated annealing algorithms for the inventory routing problem. Int. Trans. Oper. Res.; 2018; 25, pp. 1785-1809. [DOI: https://dx.doi.org/10.1111/itor.12547]

36. Hiassat, A.; Diabat, A.H.; Rahwan, I. A genetic algorithm approach for location-inventory-routing problem with perishable products. J. Manuf. Syst.; 2017; 42, pp. 93-103. [DOI: https://dx.doi.org/10.1016/j.jmsy.2016.10.004]

37. Mahjoob, M.; Fazeli, S.S.; Milanlouei, S.; Tavassoli, L.S.; Mirmozaffari, M. A modified adaptive genetic algorithm for multi-product multi-period inventory routing problem. Sustain. Oper. Comput.; 2021; 3, pp. 1-9. [DOI: https://dx.doi.org/10.1016/j.susoc.2021.08.002]

38. Solyalı, O.; Süral, H. An Effective Matheuristic for the Multivehicle Inventory Routing Problem. Transp. Sci.; 2022; 56, pp. 1044-1057. [DOI: https://dx.doi.org/10.1287/trsc.2021.1123]

39. Williams, H.P. Model Building in Mathematical Programming; Wiley: Hoboken, NJ, USA, 2013.

40. Margot, F. Symmetry in Integer Linear Programming. 50 Years of Integer Programming 1958–2008: From the Early Years to the State-of-the-Art; Jünger, M.; Liebling, T.M.; Naddef, D.; Nemhauser, G.L.; Pulleyblank, W.R.; Reinelt, G.; Rinaldi, G.; Wolsey, L.A. Springer: Berlin/Heidelberg, Germany, 2010; pp. 647-686. [DOI: https://dx.doi.org/10.1007/978-3-540-68279-0_17]