An ALNS-based matheuristic algorithm for a multi-product many-to-many maritime inventory routing problem

In this paper, we propose an adaptive large neighborhood search-based matheuristic algorithm to solve a multi-product many-to-many maritime inventory routing problem. The problem addresses a short sea inventory routing problem that aims to find the best route and distribution plan for multiple products with a heterogeneous fleet of vessels through a network including several producers and customers. Each port can be visited a given number of times during the planning horizon, and the stock level for each product should lie within the predefined bound limits. The problem was introduced by Hemmati et al. (Eur J Oper Res 252:775–788, 2016). They developed a mixed integer programming formulation and proposed a matheuristic algorithm to solve the problem. Although their proposed algorithm worked well in terms of running time, it suffers from disregarding a part of the solution space. In this study, we propose a new matheuristic algorithm to find better solutions by exploring the entire solution space for the same problem. In our solution methodology, we split the variables into routing and non-routing variables. Then in an iterative process, we determine the values of the routing variables with an adaptive large neighborhood search algorithm, and we pass them as input to a penalized model which is a relaxed and modified version of the mathematical model introduced in Hemmati et al. (2016). The information from solving the penalized model, including the values of the non-routing variables, is then passed to the adaptive large neighborhood search algorithm for the next iteration. Several problem-dependent operators are defined. The operators use the information they get from the penalized model and focus on decreasing the penalty values. Computational results show up to 26% improvement in the quality of the solutions for the group of instances with a large feasible solution space. We get the optimal value for the remaining instances matched with the reported results.


Introduction
The Inventory Routing Problem (IRP) addresses two main logistic activities, inventory management, and vehicle routing, and is a well-studied problem in the context of vendor-managed inventory.In the IRP, the goal is to manage the flow of products across different producers and consumers to satisfy customer demands and minimize transportation costs.In IRPs, the decisions are to determine the quantity of products required to be handled at each location, the time of visiting each location, and the vehicle routes to serve different customers (Coelho et al. 2013).Various types of inventory routing problems have been studied in the literature for different modes of transportation.Two surveys by Andersson et al. (2010) and Coelho et al. (2013) have explored different aspects and variants of IRPs.
The problem we are addressing is a Short Sea Inventory Routing Problem (SSIRP), which was introduced by Hemmati et al. (2016).The problem is a multi-product, many-to-many Maritime Inventory Routing Problem (MIRP) in a continuous time framework.In this study, we propose a new method to find better solutions for the same problem.We use the same realistic-sized instances and compare the quality of the solutions with the results of the Hybrid Cargo Generating and Routing algorithm (HCGR) proposed by Hemmati et al. (2016).
Our focus here is to review different solution approaches that have been applied to similar maritime inventory routing problems.According to the literature, different mathematical models were developed to formulate MIRPs, and various commercial solvers were applied to solve the problem.(Al-Khayyal and Hwang 2007;Diz et al. 2017;Agra et al. 2017), and (Misra et al. 2020).Moreover, different exact algorithms, including but not limited to branch-and-bound (Christiansen 1999) and (Agra et al. 2013), branch-and-price (Hewitt et al. 2013), and branchprice-and-cut (Engineer et al. 2012) have been studied in the past.
Among the relevant literature, Christiansen (1999) solved a single-product MIRP with time windows by applying a branch-and-bound algorithm and using the Dantzig-Wolf decomposition approach.They decomposed the overall problem into ship routing and inventory management sub-problems.Later, Agra et al. (2013) applied a branch-and-bound approach to solve a single-product MIRP.They proposed two discrete time formulations with valid inequalities.Moreover, Engineer et al. (2012) implemented a branch-price-and-cut algorithm to solve a practical-sized single-product MIRP in an oil company.They used different cuts, including capacity cuts and other special cuts targeting fractional solutions, to find the optimal solution.Furthermore, Hewitt et al. (2013) proposed a branchand-price algorithm to solve a MIRP and mentioned that their results are near-tooptimal.They reduced the required time for finding good solutions by developing local search schemes.
However, in addition to the aforementioned exact methods, various heuristics, metaheuristics, and matheuristics were developed to solve MIRPs.For example, Christiansen et al. (2011) studied a multi-product MIRP in the cement industry where they proposed a constructive heuristic embedded in a Genetic Algorithm to solve realistic-sized instances within a reasonable time.Moreover, Song and Furman (2013) presented a Large Neighborhood Search (LNS) through a simple algorithmic framework to solve MIRP.Later, Agra et al. (2016) developed a local search heuristic for the stochastic MIRP.Recently, Friske and Buriol (2020) developed a solution approach including two metaheuristics: a multi-start algorithm and a large neighborhood search.They used a commercial solver to solve the reduced mixed integer problem with the solutions obtained by LNS.In addition, Friske et al. (2022) studied the use of Relax-and-Fix and Fix-and-Optimize metaheuristics over two discrete-time formulations.They evaluated the contribution of different components of the formulations to the performance of the algorithm.
As an extension of the LNS algorithm, Adaptive Large Neighborhood Search (ALNS) is known as an efficient heuristic to deal with vehicle routing problems.(Ropke and Pisinger 2006;Pisinger and Ropke 2007;Ribeiro and Laporte 2012;François et al. 2019;Yu et al. 2020;Chen et al. 2021).The effectiveness of ALNS in maritime routing problems has also been shown by different researchers.(Lianes et al. 2021;Brekkå et al. 2022).
Several published studies are relevant to our problem setting, which are discussed in more detail as follows.Christiansen (1999) formulated a MIRP in which a heterogeneous fleet of ships was planned to transport a single product through a network including producers and customers.Each port could be visited a given number of times during the planning horizon.They assumed a time window limitation for starting service at each port and the possibility of serving multiple ships at the same time.The problem was solved using branch-and-bound and applying the Dantzig-Wolf decomposition approach.In each instance, a maximum of 5 vessels, 16 ports, and 36 days as the planning horizon were considered.Later, Al-Khayyal and Hwang (2007) formulated the same problem of transporting several products on a many-tomany distribution network.The model assumed dedicated compartments for each product in each ship without time window limitations.The small-size instances were solved using a commercial solver.They mentioned the impact of the number of port visits on the solution time and the need for particular algorithms to take advantage of the structure of the model.They considered up to 3 products, 4 vessels, 4 ports, and 10 days for the planning horizon.
A few years later, Siswanto et al. (2011) proposed a set of heuristics with four groups of rules to solve the same problem with undedicated compartments.The vessels had different numbers of compartments that were not dedicated to a specific product.They provided high-quality solutions for the set of instances with at most 2 products, 3 vessels, 4 ports, and 15 days as the planning horizon.Following that, Agra et al. (2014) proposed a tightened model with valid inequalities for the problem described by Al-Khayyal and Hwang (2007).Given the continuous-time 44 Page 4 of 23 formulation, they presented three heuristics: rolling horizon, feasibility pump, and local branching.They proposed heuristics to solve real instances with a maximum of 4 products, 2 vessels, 7 ports, and 15 days for the planning horizon.According to their results, the quality of solutions improved by combining all three heuristics.
Later, Hemmati et al. (2016) proposed the HCGR algorithm to solve the multiproduct SSIRP in a continuous time framework.A heterogeneous fleet of vessels was assumed to distribute multiple products within a many-to-many distribution network.They assumed different compatibility between vessels, ports, and products.There was no limitation on the compartment and no possibility of serving multiple vessels at the same time in their problem setting.They converted the IRP into a routing and scheduling problem to solve realistic-size instances.They considered up to 3 products, 10 vessels, 16 ports, and 60 days for the planning horizon.First, they solved two mathematical models, including transportation and time window models, to find the number of products picked up and delivered from/to each port and the corresponding time scheduling.In this step, the pair of producer and customer and the time of pickup and delivery operations became fixed based on the direct transportation and operational costs.In the next step, based on the best scheduling plan, the routing problem was solved using an ALNS.The algorithm disregarded some parts of the solution space because of making a fixed pair of producers and customers.The proposed algorithm worked well in terms of running time.
Following that, Agra et al. (2017) presented discrete-time and continuous-time formulations and different valid inequalities for the problem introduced by Al- Khayyal and Hwang (2007).They defined different production/consumption rates in the discrete-time formulation for various periods.The rates were fixed during the whole planning horizon in the continuous-time formulation.They compared different proposed formulations in terms of their size, running time, and integrality gap for both discrete and continuous time.Using a commercial solver, they reported the computational results for real test instances.They found the optimal solutions for the cases with a maximum of 4 products, 2 vessels, 7 ports, and a fifteen-day planning horizon.
As mentioned earlier, the SSIRP proposed by Hemmati et al. (2016) was to find the best route and distribution plan for multiple products through a network with several producers and customers.They studied two groups of instances with different compatibility between vessels, ports, and products.The size of the instances was considerably larger compared to the literature.However, their proposed solution method could not find high-quality solutions for the instances where all the vessels were fully compatible with the ports and the products.In this paper, we propose a new solution approach to find better solutions for the problem introduced by Hemmati et al. (2016).To this end, we develop an ALNS-based Matheuristic algorithm (ALNSM) in which an ALNS algorithm determines routing variables.Moreover, the non-routing variables are taken care of by a Mixed Integer Programming (MIP) solver, given the value of routing variables.We compare our results on the test instances in Hemmati et al. (2016) and analyze the quality of the results accordingly.
The rest of the paper is organized as follows: problem characteristics are explained in Sect.2, the proposed ALNSM algorithm is described in Sect.3, and computational results are reported in Sect.4, followed by the conclusion in Sect. 5.

Problem description
The problem introduced by Hemmati et al. (2016) addresses the distribution management of different products within a many-to-many distribution structure, including several producers and customers.It is assumed that a heterogeneous fleet of vessels is planned to transport products through the network.The vessels differ in terms of capacities, compatibility with ports and products, and transportation cost and time.There is a limitation on the number of visits for each port during the planning horizon.Furthermore, the stock level of each product should lie within the predefined minimum and maximum levels.The initial stock level and production or consumption rates for each product at each port are given.The quantity of the products to be handled at each port is limited, and the handling operation time at each port for each product is predefined.The problem is to minimize transportation and operational costs, including fixed and variable handling costs.
The following mathematical model is presented by Hemmati et al. (2016).However, for the sake of a better explanation of our solution approach, their mathematical formulation is given as follows: Ports are indexed by i and j, and vessels, products, and ports visit numbers are represented by v, k, and m/n, respectively.The m th visit of port i is denoted by (i, m), and direct travel from port visit (i, m) to port visit (j, n) is represented by (i, m, j, n).N, V, and K are defined as the set of ports, the set of vessels, and the set of products, respectively.Moreover, the set of ports that vessel v can visit and the set of products k that can be transported with vessel v are denoted by N v and K v .Also, the set of all possible port visits, the set of possible port visits for vessel v, and the set of all possible direct travels for vessel v are represented by S A , S A v , and S X v , respectively.

Mathematical formulation
Parameters and variables of the model are listed below: Parameters Minimum and maximum amount of product k that is allowed to be han- dled at port i Fixed operational cost for product k at port i

Variables
x imjnv : 1 if there is a direct movement from port visit (i, m) to port visit (j, n) with vessel v and 0 otherwise.
x O imv : 1 if there is a movement from the vessel v's origin to port visit (i, m) and 0 otherwise.
1 if the operation of vessel v ends at port visit (i, m) and 0 otherwise.
1 if vessel v is not used and 0 otherwise.The mathematical model is defined as follows: Objective function: Subject to Routing constraints: (1) 44 Page 8 of 23 Loading and unloading constraints Time constraints

Inventory constraints
The objective function (1) aims to minimize the total transportation and operational costs.Equations (2) determine whether vessel v is used or not.In (3), the previous location of each port visit is defined, which can be either the origin of vessel v or another port visit of that vessel (j, n).Equations ( 4) state that each vessel can either end its route at a port visit (i, m) or continue to another port visit (j, n).Constraints (5) ensure that each port visit (i, m) is carried out by one of the available vessels, and equations ( 6) ensure that the port visits are conducted successively.Constraints ( 7)-( 10) define binary variables while constraints ( 11) and (12) determine the onboard quantity of each product after each port visit.Equations ( 13) describe the capacity constraint for the vessels.In ( 14), the number of products allowed to be loaded/ unloaded is limited.Equations ( 15) and ( 16) define the possibility of loading and (20) unloading operations with the required port visits.Constraints ( 17) and ( 18) introduce loading variables.Constraints (19) determine the end-time operation at port visit (i, m) based on the number of products loaded/unloaded.Equations (20) state that the operation at port i cannot start before ending its previous port visit.Constraints ( 21) and ( 22) relate the beginning operation time at each port visit with the transportation time from the last port visit and its prior ending operation time.Equations ( 23) define time variables.
Constraints (24) determine the stock level of each product at each port at its first visit.Equations ( 25) and ( 26) define the inventory level of each product at the beginning and end of each port visit.Constraints ( 27) and ( 28) ensure that each product's stock level lies within the predefined maximum and minimum levels at each port visit and at the end of the planning horizon.Constraints (29) define stock variables.

ALNS-based matheuristic algorithm
We propose a matheuristic algorithm to solve a multi-product many-to-many inventory routing problem.In our algorithm, we take advantage of the power of metaheuristics to search the discrete domain and the capabilities of MIP solvers to deal with continuous variables.
In the proposed algorithm, called Adaptive Large Neighborhood Search-based Matheuristic (ALNSM), all variables in the mathematical model are divided into two main groups: routing and non-routing variables.The routing variables are taken care of by a metaheuristic algorithm known as Adaptive Large Neighborhood Search.Following that, the non-routing variables are determined by the MIP solver, given the values of the routing variables.
Considering the difficulty of generating feasible solutions for the problem, we define a penalized model where the inventory constraints are relaxed, and instead, penalty functions for the violated constraints are added to the objective function.As finding a feasible solution for the penalized model is still difficult, we initiate our algorithm with a sub-model that offers a more relaxed version of the penalized model.This sub-model, explained in Sect.3.2, is called only once at the beginning of the algorithm to generate a feasible initial solution.
In each iteration of the ALNSM algorithm, a new solution is generated using a group of designed operators to determine routing variables.To this end, we design different operators for the penalized and the original problems.Once the search reaches the solution space of the original problem, i.e., the penalty variables become zero, it continues using operators designed solely for the original problem.
Next, the routing variables are passed to the MIP solver to calculate the objective function, determine the values of the non-routing and penalty variables, and check the feasibility of the solution.The ALNSM algorithm proceeds with searching the solution space until it meets the stop-criterion.In order to escape from local optima, an escape-algorithm is applied which is described in Sect.3.4.Ultimately, the bestfound solution is reported.The main steps of the algorithm are summarized in Algorithm 1.The proposed penalized model and initial solution are described in Sects.3.1 and 3.2, respectively, followed by the ALNS algorithm in Sect.3.3.

Penalized model
Taking the complexity of the problem into account for finding a feasible solution, a set of constraints in the mathematical model, i.e., inventory constraints, are relaxed, and the corresponding penalties for the violated constraints are added to the objective function.We define the stock penalty variables as follows: g imk : The amount of stock shortage for product k at port visit (i, m) at the begin- ning of handling operation compared to the lower bound for the stock level, S ik .
g E imk : The amount of stock shortage for product k at port visit (i, m) at the end "E" of handling operation compared to the lower bound for the stock level, S ik .
e imk : The excess amount of stock for product k at port visit (i, m) at the beginning of handling operation compared to the upper bound for the stock level, S ik .
e E imk : The excess amount of stock for product k at port visit (i, m) at the end "E" of handling operation compared to the upper bound for the stock level, S ik .
g T ik : The amount of stock shortage for product k at port i at the end of planning horizon "T" compared to the lower bound for the stock level at the end of planning horizon, S T ik .
e T ik : The excess amount of stock for product k at port i at the end of planning horizon "T" compared to the upper bound for the stock level at the end of planning horizon, S T ik .
With the new variables defined above, the inventory constraints (27, 28) are relaxed as follows: Moreover, we define "P" as a sufficiently big number, and we update the objective function by adding the following penalty function: The main challenge here is to minimize the penalty function and find a feasible solution for the original problem.To this end, we design a set of operators to focus on this purpose described in Sect.3.3.2.

Initial solution
Although the stock constraints in the penalized model are relaxed, the routing, loading/unloading, and time constraints should be satisfied to start with a feasible initial solution. (30) To make it easier to find a feasible initial solution, we start the algorithm with a solution in which each port is visited only once, and a set of constraints are satisfied.To this end, the sub-model includes routing constraints (2)-( 5) and ( 7)-( 10), loading and unloading constraints ( 11)-( 18), and the following constraints: Constraints ( 33) and (34) ensure that each port is visited only once.Equations ( 35)-( 38) along with the other mentioned constraints determine one feasible route for each used vessel.

Adaptive large neighborhood search algorithm
The core of our proposed ALNS-based matheuristic algorithm is based on an adaptive large neighborhood search framework inspired by Ropke and Pisinger (2006).In this framework, as an efficient heuristic for vehicle routing problems, different operators are defined to create new solutions.In the selection procedure, initially, all the operators have the same chance to be selected to generate new solutions.After a number of iterations called a segment, the probability of choosing different operators is updated based on their performance in the previous iterations.Therefore, the more efficient operators have a higher chance of being selected in the next iterations.The algorithm selects the operators using a roulette wheel selection principle.
Algorithm 2 presents the framework of the ALNS algorithm.The acceptance method that we used in the following framework (lines 7-12) is based on the wellknown simulated annealing algorithm's acceptance criteria (Kirkpatrick et al. 1983). (33) 44 Page 14 of 23

representation
The solution is represented with a two-dimensional vector defined as follows: in the first dimension, the sequence of port visits for different vessels is defined, and in the second dimension, the corresponding visit counter for each port is specified.For example, suppose that the number of ports is 6, the number of vessels is 3, and the maximum number of visits for each port is 2. Therefore, one possible solution can be represented as follows: 1 3 0 1 2 6 0 5 3 4 2 1 0 1 1 1 0 1 2 1 The solution depicts that the first vessel makes the second port visit of port 1, and the first port visit of port 3. The second vessel makes the first port visits of ports 1, 2, and 6, and the third vessel makes the first port visit of port 5, the second port visit of port 3, and the first port visit of port 4. In this solution representation, vessels and their routes are separated by (0, 0).

Operators
To find the best routing plan for the fleet of vessels, we design different operators to focus on the sequence of visits performed by vessels.The main challenge is moving from the penalized model's feasible solution space to the original feasible region.To achieve this, in the first part of the algorithm, we design several operators to focus on decreasing the penalty function and moving toward the feasible region of the original problem.Once the algorithm reaches the original feasible solution space, the other operators are applied to find the best solution to the original problem.The following three principles are used in all designed operators: 1.Each used vessel should start its route by visiting one producer and end it with one customer.2. If a given port is producing only one specific product and there is only one customer for that product, both the producer and the corresponding customer must be visited by the same vessel.3. Changing the position of ports in a solution may result in rearranging port visit counters accordingly.
Among the following operators, the first two are specifically designed to decrease the penalty function, while the rest of them are applied to search the solution space for both penalized and original problems.
• Stock-Balanced: The Stock-Balanced operator removes one port visit with a positive stock penalty value and inserts it in a possible place in a route.Posig imk and g E imk indicate that the corresponding vessel arrives late at the port visit (i, m).Moreover, positive g T ik denotes that the last visit to port i is performed early so that the stock level of product k at the end of the planning horizon is less than the minimum required level.Similarly, the positive values for e imk /e E imk and e T ik indicate a high stock level of product k at the port visit (i, m), and port i at the end of the planning horizon, respectively.The Stock-Balanced operator removes one of the port visits with a positive stock penalty value.Ports with higher penalty values have more chances to be removed.If g imk , g E imk , e imk , and e E imk are positive, the removed port visit will not be inserted at any place in a route of its allocated vessel later than its present point.Instead, it can be done earlier or assigned to another compatible vessel.However, positive g T ik and e T ik represent that the last visit to port i should be done later so that the stock level of product k at port i remains in a range of predefined bound limits at the end of the planning horizon.In this case, the last port visit to port i should not be done earlier than its present time.Therefore, either it can be done later or moved to another compatible vessel.
• Random-Worst-Remove-Random-Insert (RWRR): Compared to the Stock-Balanced operator, the RWRR operator selects more than one port visit with positive penalty values and inserts them in random possible places in a route.The maximum number of port visits to remove from each route is limited such that each route includes at least one producer and one customer.Among port visits with positive values, the number of possible port visits to remove is determined randomly, defining at most four.Port visits with higher penalty values get more chances to be selected and removed.Next, selected port visits with their connected ports (based on the second principle described above) are removed and inserted in a random possible place in a solution, considering the compatibility of ports and vessels.Since more than one port visit is selected to be removed, the insertion rule in the Stock-Balanced operator is not followed by this operator.Changing the position of several port visits may affect the status of other removed port visits.
44 Page 16 of 23 • Random-Remove-Random-Insert (RRR): Unlike the RWRR, which is applied to the penalized model, RRR is used for both penalized and original problems.Herein, at most, four port visits are selected randomly to remove.Ensuring that at least one producer and one customer are assigned to each vessel, selected port visits are removed and inserted in a random possible place in a route.• Random-Remove-First-Improvement-Insert (RRFI): Here, one port visit in a solution is selected randomly, ensuring that it does not belong to a vessel with only one producer and one customer.Next, the removed port visit is inserted in a place where the first improvement is found, i.e., a point in a solution in which the new solution is still feasible, and the objective value is less than its present value.To calculate the objective value and check the feasibility of the solution, the solver is called for each step of insertion.If there is no better possible place for insertion, the removed port visit remains in its current place in a solution.• Swap: To diversify the search, we change the position of two port visits in a solution with a swap operator.It is necessary to make sure that selected ports are compatible with their newly assigned vessels, the logic of starting a route with one producer and ending it with one customer being satisfied, and each used vessel contains at least one producer and one customer after the swap operation.• Route-Add: To generate solutions with fewer idle vessels, one producer and one customer are selected and successively inserted in one possible place in a solution.It is required to make sure that a common compatible vessel is used to carry a product for a producer and its customer.Therefore, it is possible to use more vessels for the transportation plan, which may create more and better feasible solutions.• Vessel-Change: Compared to the Route-Add operator, the Vessel-Change operator generates new solutions where one used vessel is discarded, and its allocated route is assigned to the other compatible vessels.This operator helps to escape from local optima when it is necessary to reduce the number of used vessels.• Visit-Add: In case we are allowed to visit each port more than once, the Visit-Add operator is applied.The minimum required number of visits for each port is calculated based on the initial inventory of each product at each port, the production and consumption rates of each product at each port, the stock bounds at each port for each product, and the maximum amount of each product allowed to be handled at each port.It is necessary to ensure that each port is visited at least to its minimum required number of visits.In the Visit-Add operator, extra port visits are added to the current solution to check if more port visits could lead to a better solution at a lower cost.• Visit-Remove: Similar to the Visit-Add operator, the minimum required number of visits for each port is calculated.In this operator, the extra port visits compared to its minimum required number of visits (if any) are removed from the current solution to check the possibility of making a better new feasible solution.• Vessel-Swap: In this operator, the routes of two vessels with common characteristics but different transportation costs are changed.The Vessel-Swap operator helps to create better feasible solutions where compatible vessels are strictly limited, and the transportation costs for compatible vessels are considerably different.
• Farthest-Remove-Nearest-Insert (FRN): The FRN operator is designed to generate solutions with shorter distances among port visits in one single route.First, one port visit with the longest distance from its subsequent port visit is selected randomly.Port visits are arranged based on their distance in descending order to avoid choosing the same port visit in successive iterations.Port visits with longer distances have more chances to be selected to remove.Finally, the removed port visit is inserted into a solution next to its nearest possible place.

Escape-algorithm
Although the ALNSM algorithm benefits from diversifying components, including the simulated annealing and randomness inside the operators, we also exploit the escape-algorithm to escape from local optima.In Algorithm 1, if there is no improvement after a number of iterations (m), we activate the escape-algorithm.In this iterative algorithm, two operators are randomly selected to diversify the search.
The new solution in each iteration is accepted regardless of its quality.The algorithm stops once a better solution than the best-found solution is found.Otherwise, the algorithm continues until it meets the stopping criterion, which is reaching a certain number of iterations.

Computational result
A set of realistic-sized instances, as defined by Hemmati et al. (2016), is used to evaluate the performance of the ALNSM algorithm.The instances are categorized into two groups.In the first group, instances are deliberately designed to allow for problem segmentation into smaller sub-problems, enabling the reporting of optimal solutions.Conversely, the instances in the second group, indexed as 22-28 in Table 1, cannot be partitioned into smaller sub-problems, thus preventing the reporting of optimal solutions.The instances are characterized based on the number of producers and customers, the number of vessels, and the number of products for a fixed planning horizon.Furthermore, the test instances are divided into two groups based on their compatibility: with full compatibility and with limited compatibility between vessels, ports, and products.It is assumed that there is limited compatibility between vessels, ports, and products for all instances except for those indexed 22-28.In the group of instances with limited compatibility, we can only visit each port and carry each product with a few compatible vessels, while in the cases with full compatibility, all vessels, ports, and products are compatible.
Gurobi is used as the MIP solver in our proposed algorithm.In the ALNS algorithm, the number of iterations in each segment is 50.Moreover, in the escape-algorithm, m is defined as 100, and the stop-criterion is set at 20 iterations.The two selected operators in the escape-algorithm include swap and RRR.The ALNSM algorithm stops after 6000 iterations.
Table 1 illustrates the computational results achieved by the ALNSM and its comparison with the HCGR presented by Hemmati et al. (2016).The first four columns of Table 1 specify the instance index, the number of ports, the number of vessels, and the number of products, respectively.The next column shows the best-found solutions by the HCGR.The following three columns illustrate the average objective values, the best-found solution, and the average running time (in seconds) achieved by running the ALNSM 10 times on each instance.The last column shows the gap between the ALNSM and the HCGR.Bold values are reported as the optimal solutions by Hemmati et al. (2016).
As the results show, our proposed algorithm finds the optimal solution for all instances with limited compatibility between vessels, ports, and products.However, for the remaining cases, our proposed algorithm outperforms the HCGR.Specifically, for instances with full compatibility, the gap between the results of the ALNSM and the HCGR is improved by up to 26%, respectively.
The ALNSM and the HCGR algorithms are tackling the problem differently.The HCGR algorithm is initiated by specifying the minimum required quantity to be loaded/unloaded from each port.This is calculated based on the given production/ consumption rates, initial stocks, and the minimum and maximum bounds of stock levels.Then, the best amount of products to be loaded/unloaded from each port and the relevant time windows are determined.Here, the HCGR algorithm specifies the fixed pair of pickup and delivery ports based on their direct transportation costs and operational costs.Next, the best route of fixed pickup and delivery ports is determined using adaptive large neighborhood search.The main drawback of the HCGR algorithm is disregarding some parts of the feasible solution space and considering only the immediate impact of decisions.For example, consider a network with six ports, where ports 1 and 5 act as producers for one product, and ports 2, 3, 4, and 6 are customers.Assume that the minimum required quantities to be loaded from ports 1 and 5 are 660 and 220 units, respectively, and each customer demands a minimum of 220 units.The HCGR algorithm yields the following fixed pickup and delivery pairs: (1, 2) with 440 units of products to be loaded from port 1 and unloaded at port 2 (1, 3) with 220 units of products to be loaded from port 1 and unloaded at port 3 (5, 4) with 220 units of products to be loaded from port 5 and unloaded at port 4 (5, 6) with 220 units of products to be loaded from port 5 and unloaded at port 6.The above solution is obtained due to the consideration of direct transportation costs between the paired pickup and delivery ports.However, a feasible and more improved solution is possible by keeping products onboard the vessels and unloading them later.Therefore, a better feasible solution is the route of 1-2-3-4 with the loading of 660 units at port 1, unloading 220 units at ports 2, 3, and 4, and the route of 5-6 with the loading of 220 units at port 5, and unloading 220 units at port 6.Our proposed algorithm can achieve this by determining different load/unload quantities with no limitation in searching the entire solution space.In the ALNSM, routes are defined using an adaptive large neighborhood search algorithm, and a MIP solver is utilized to obtain optimal values of non-routing variables, subject to the specified routing variables.This process continues until the stop-criterion is met, reporting the best-found feasible solution.
Two groups of instances are analyzed separately to highlight the effect of problem characteristics on the running time.Figures 1 and 2 depict the relationship between the running time and problem characteristics for the first and second groups of instances, respectively.
In this context, the problem characteristic is defined as the multiplication of the number of ports, the number of vessels, and the number of products.As the results show, the required time to solve different instances is affected by the problem characteristics.Moreover, according to Fig. 1, within the group of instances with identical problem characteristics, the number of active producers and customers affects the running time.Active producers and customers handle different product types in their stock.Increasing the number of inventories in a problem leads to a corresponding increase in the running time.

Conclusion
In this study, we proposed an ALNS-based matheuristic algorithm to solve the multi-product many-to-many inventory routing problem introduced by Hemmati et al. (2016).The proposed ALNSM follows an iterative process where routing variables are determined using an adaptive large neighborhood search algorithm, while a MIP solver is employed to determine the optimal values of stock and time variables based on the obtained routing values.Comparing the performance of the ALNSM with the HCGR algorithm reveals that the ALNSM significantly outperforms the HCGR on different sets of realistic instances.Particularly, our algorithm excels in improving solution quality when there is full compatibility between vessels, ports, and products.In such instances, the feasible solution space is larger, necessitating more exploration to find better solutions.The ALNSM's approach of emphasizing finding optimal values for non-routing variables for each specific vessel routs allows for a more comprehensive search compared to the HCGR, where some areas of the solution space are not even explored.Our results demonstrate notable improvements of up to 26% in solution quality for this group of instances.Furthermore, for the other group of instances with limited compatibility, the ALNSM achieves the optimal solutions reported by Hemmati et al. (2016).
Funding Open access funding provided by University of Bergen (incl Haukeland University Hospital).

Declarations
Conflict of interest None of the authors have any financial or non-financial interests that are directly or indirectly related to this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.

Fig. 1
Fig. 1 With full compatibility between vessels, ports, and products meets m th visit of port i and 0 otherwise y im :1 if m th visit of port i is done and 0 otherwise.
o imvk : 1 if product k is loaded/unloaded at port visit (i, m) with vessel v and 0 otherwise l imvk : Onboard quantity of product k on vessel v after port visit (i, m) q imvk : Loaded/unloaded amount of product k at port visit (i, m) with vessel v imk : Level of inventory for product k at the start of handling operation at port visit (i, m)s E imk :Level of inventory of product k at the end of handling operation at port visit (i, m)

Table 1
Computational results for ALNSM and HCGR