A Two-Stage Matheuristic Algorithm for Classical Inventory Routing Problem

1 Introduction

The rise of e-business brings great opportunities for the supply chain industry, along with challenges due to the growing demands for higher throughput and lower latency. Fortunately, the developments of Internet of Things and statistical learning makes it possible to collect sufficient data from every link in the value chain and provide constructive analysis and predictions. By utilizing these data, the related companies are able to reduce their operation expenses, increases the profits, improve the efficiency, and bring better user experiences. Thus, more and more companies start to pay attentions to systematic optimization for the whole supply chain.

∗Corresponding author

Traditionally, suppliers or retailers manage their own warehouses and the inventory replenishment by themselves. They work out specific orders according to their production and consumption, and submit them to the logistics companies and the latter ones just do the deliveries. Due to the uncertainty of the call-in orders, it is hard for the logistics companies to make reasonable plans, which usually affects the operation efficiency and introduces unnecessary costs. In order to tackle the drawbacks brought by asymmetric information and achieve global optimization, vendor managed inventory (VMI) comes into play. Distinguished from the call-in mode, vendors (logistics companies) provide integrated management for all their customers (suppliers and retailers) in VMI mode. In other words, vendors are in charge of monitoring their clients' inventory levels, predicting product consumption, scheduling replenishment, and planning delivery routes. It is believed that the integration of information and resources will bring benefits for both vendors and customers.

Inventory routing problem (IRP) is one of the most important challenges to apply the VMI mode in supply chain management. It combines two classical optimization problems, which are inventory problem and vehicle routing problem, respectively. Regarding the inventory part, it aims to minimize product backlog. Due to the risks for storing too many products that the commodities may expire, evaporate, be stolen, or encounter other accidents, there will be inventory holding costs proportional to the inventory level at each customer. As for the routing part, it requires to reduce the traveling costs which are usually measured by the total distance of the delivery routes. When the inventory management and the transportation optimization is jointly considered in the IRP, the entire planning horizon is divided into multiple independent periods, we need to decide the amount of product to replenish for each customer in each period, along with the delivery route in each period, so that no customer runs out of stock and the total cost for backlog and transportation is minimized. These two factors often contradict to each other, because less backlog means less quantity per delivery, which means more frequent deliveries if the consumption holds, which increases the traveling costs, and vice versa. Therefore, it is natural to consider the inventory management and the vehicle routing together, in order to balance and optimize the total costs.

As a challenging NP-hard problem [Coelho et al., 2013] with valuable applications, IRP has been drawing attentions from both industrial and academic communities. Commercial software such as Aspen Fleet Optimizer was developed and widely used in the petroleum industry [Becraft et al., 2012]. Air Liquide proposed a complex inventory routing problem in its healthcare business as the topic of the ROADEF/EUROChallenge 2016 1 . It is also one of the topics of the 12th DIMACS Implementation Challenge 2 . Apart from the industrial applications, many researchers investigated the IRP and its variants in recent decades. There are mainly four characteristics that distinguishes different versions of IRPs, which are fleet composition, planning horizon, demand mode, and replenishment policy [Coelho et al., 2013], respectively. The single-vehicle multi-period deterministicdemand inventory routing problem under maximum-level policy proposed in Bertazzi et al. [2002] is one of the most representative and classical versions [Archetti et al., 2007; Archetti et al., 2012; Coelho et al., 2012; Alvarez et al., 2018; Chitsaz et al., 2019]. Based on this well-known IRP, several variants concerning homogeneous fleet [Coelho and Laporte, 2013; Adulyasak et al., 2013; Coelho and Laporte, 2014; Archetti et al., 2017] or multiple depots [Bertazzi et al., 2019] have been widely studied. Regarding the solution methods, several exact algorithms, heuristics, and hybrid algorithms were proposed for solving the IRP. The exact algorithms such as branch-and-cut algorithms [Archetti et al., 2007; Coelho and Laporte, 2013; Avella et al., 2017] have made great progress on solving small IRP instances. Desaulniers et al. [2015] presented a branch-cut-and-price algorithm for the IRP, which utilizes the column generation approach to calculate the lower bound, and employs the branch-and-cut algorithm to tighten the bound. In addition, several metaheuristic algorithms were proposed to tackle the IRP. Aksen et al. [2014] defined multiple neighborhood structures and integrated them into an adaptive large neighborhood search algorithm to solve the IRP. Park et al. [2016] presented a genetic algorithm inspired by the CPLEX optimizer. Iterated local search and simulated annealing were also used for tackling the IRP [Alvarez et al., 2018], and they can obtain satisfactory solutions in relatively short time.

Recently, more and more hybrid algorithms demonstrate their effectiveness in solving this intractable problem. These algorithms usually decompose complex problems into several sub-problems, and utilizes different kinds of methods for solving these reduced problems in certain hierarchy. The matheuristics, which integrate metaheuristic algorithms and mathematical programming techniques, is one of the most common and successful hybrid algorithms. For instance, Archetti et al. [2012] proposed a hybrid heuristic for the IRP, which uses the tabu search to find good solutions, and refines these solutions by optimizing two mixed-integer programming (MIP) models. Coelho et al. [2012] presented an algorithm based on adaptive large neighborhood search framework. It adopts a linear programming (LP) model to determine delivery quantities for each neighboring solution, and tunes the best solution found by another MIP model in a specified frequency. Campbell and Savelsbergh [2004] decomposed the IRP into two phases. It first constructs highlevel base plan using an integer programming model, then utilizes an insertion heuristic to generate the complete schedule. In addition, Cordeau et al. [2015] solved the IRP by a multi-phase hybrid algorithm. It determines replenishment plans by a Lagrangian-based method in the first phase, generates vehicle routes by a constructive heuristic in the second phase, and reoptimizes the solution via MIP.

1 http://www.roadef.org/challenge/2016/en/sujet.php

2 http://dimacs.rutgers.edu/events/details?eID=1090

This paper proposes a two-stage matheuristic (TSMH) algorithm which consists of two optimization stages in different granularities for the classical IRP. The TSMH algorithm begins with a coarse-grained structural optimization which determines the initial replenishment schedule. In the second stage, the TSMH algorithm adopts a solution-based tabu search procedure to bring sophisticated refinement to the delivery routes, timing, and quantities. In both stages, the proposed algorithm employs both metaheuristics and mathematical programming to tackle the routing and inventory subproblems, respectively. Tested on two sets of totally 220 widely used instances of the IRP, the proposed TSMH algorithm improves the best known results on 40 out of 60 large-scale instances and matches the records on most small instances. Moreover, comprehensive tests and comparisons demonstrate that the combination of the two stages and corresponding techniques are essential for the proposed algorithm.

2 Classical Inventory Routing Problem

The classical inventory routing problem is to deliver products to a number of clients in multiple periods. It aims to minimize the total inventory holding costs and vehicle traveling costs, while preventing any customer from running out of stock in each period. In order to accomplish this goal, a solution method needs to determine when, where, and how many products at each delivery, as well as the sequence of visits to each client in each period, as illustrated in Figure 1.

More precisely, the IRP problem is defined on an undirected complete graph G = (V,E). V = {0, 1, ..., n} is the set of vertices, where vertex 0 stands for the depot and V ′ = V - {0} represents the client set. E = {(i, j) : i, j ∈ V, i ̸= j} denotes the edge set, where a corresponding traveling cost Cij is associated with each edge (i, j). There are totally p periods in the planning horizon, and the capacity of the vehicle is Q. The inventory holding cost and maximum inventory level for each depot or client i ∈ V are Hi and Ui, respectively. Within each period t ∈ T = {1, ..., p}, the depot produces rt0 units of products and client i ∈ V ′ consumes rt i units of products. To keep it simple and clear, we assume that the production and delivery always happen before loading and consumption at each vertex, respectively. This simplification follows the same convention in the previous works [Archetti et al., 2007; Archetti et al., 2012; Coelho et al., 2012]. We define variable xti to be the quantity of products delivered to vertex i ∈ V in period t ∈ T. Obviously, xt0 is always non-positive since the vehicle can only load instead of delivering products at depot. We use Boolean variable yt ij to denote if edge (i, j) exists in the delivery route in period t ∈ T . Let variable I t i be the inventory level of vertex i ∈ V at the end of period t ∈ T ∪ {0}, except that the initial inventory level I 0 i is already known. Then, the classical IRP can be formulated as the following MIP model.

Objective (1) minimizes the sum of the inventory costs and routing costs. Constraints (2) ensure that each vertex can be visited at most once in each period, and its in-degree equals to the out-degree. Constraints (3) are for subtour elimination where a subtour is a cycle which does not start from the depot [Dantzig et al., 1954]. Constraints (4) indicate that there can be positive delivery quantity for each client in each period iff the client is visited, and the delivery quantity must not exceed the capacity of the client. Constraints (5) require that in each period, the quantity loaded at the depot must match the total quantity delivered to all clients. Constraints (6) guarantee that the vehicle will never be overloaded. Constraints (7) are the inventory conservation equalities. Constraints (8) and (9) impose the capacity of each vertex. Constraints (10) forbid the shortage at customers. Constraints (11) present the domains of the decision variables.

3 Two-Stage Matheuristic Algorithm

The IRP integrates both inventory management and vehicle routing, thus involves both continuous and discrete optimizations which are usually solved by different methodologies. As a hybrid algorithm which combines the advantages of mathematical programming and metaheuristics, the matheuristic algorithm is a promising approach for such kind of problems. Meanwhile, due to the complexity of the IRP, it is very challenging to optimize the whole problem directly. Therefore, the proposed TSMH algorithm consists of two stages whose search granularities vary from high-level structure to low-level details. Specifically, the first stage generates the initial solution by row generation approach with additional repairing procedure (Section 3.1). The second stage adopts a solution-based tabu search procedure to implement the finegrained optimization (Section 3.2).

In order to find out a proper overall structure of the delivery schedule, the proposed algorithm relaxes the MIP model for the IRP to obtain an easier subproblem which optimizes the timing and the quantity of the deliveries while considering the routing as accurate as possible. Specifically, due to the exponential complexity of the subtour elimination constraints (3), relaxing themwill greatly reduce the intractability of the induced subproblem. However, once a solution for the relaxed model is found, the subtour elimination constraints will be lazily added to claim its infeasibility. Although the relaxed solutions may not form connected tours, we can decide whether a vertex i is visited or not in period t by examining Z t i = ∑j∈V,j̸=i yt ij . Furthermore, in each period, when the vertices to visit (black nodes in Figure 1) are fixed, repairing the route is equivalent to solving the traveling salesman problem (TSP) on the subgraph composed of the visited vertices. The TSMH algorithm will repeat the above solvetighten-repair procedure until the time limit is reached.

When the quality of the best found solution gradually converges in the first stage, the TSMH algorithmmoves on to the second stage. This stage employs a solution-based iterated tabu search algorithm started from the best solution found so far to carry out a fine-grained optimization. This tabu search algorithm perturbs the current best solution when it fails to improve the best solution found for a number of consecutive iterations. Finally, the TSMH stops after the time limit is reached and returns the best solution found so far.

3.1 Stage 1: Structural Optimization by Relax-and-Repair Method

We consider a relaxed IRP (RIRP) model which consists of Eqs. (1)-(2) and (4)-(11) in the first stage. As described in Algorithm 1, the TSMH solves the RIRP model to obtain relaxed solutions (line 3). In these solutions, the inventoryrelated constraints are respected but there may be subtours disconnected to the depot. Hence, for each period, we drop the clients with zero delivery quantity or Z t i = 0, and solve the TSP on the remaining subgraph using Lin-Kernighan heuristic [Lin and Kernighan, 1973] (line 4). This repairing procedure is able to efficiently produce feasible solutions while preserving the overall structure of the relaxed solutions. If the cost of the repaired solution is less than the best solu-

Algorithm 1 StructuralOptimization

Input: Instance I;

Output: Best found solution Sbest

• 1: Sbest ←∅
• 2: repeat
• 3: Solve RIRP model to get a relaxed solution Srel
• 4: Repair the tours in Srel to get a feasible solution Scur
• 5: if f(Scur) < f(Sbest) then
• 6: Sbest ← Scur
• 7: end if
• 8: Add a constraint which eliminate the minimum subtour (if there exists) to RIRP model
• 9: until time limit is reached or no subtour in Srel

tion found so far, the best solution will be updated (lines 5-7). Next, one of the violated subtour elimination constraints involving the least number of edges will be added to the RIRP model as lazy constraints (line 11). After adding the lazy constraints, TSMH resumes the optimization and iteratively executes the above procedure until the time limit is reached (line 2) or a feasible solution is found (lines 8-10). This stage is able to produce good initial solutions, and it may directly obtain the optimal solutions on some small instances.

3.2 Stage 2: Detailed Refinement by Solution-based Tabu Search

This stage adopts a solution-based iterated tabu search [Lai et al., 2018] to perform a fine-grained search starting from the best solution obtained in the previous stage. Similar to typical tabu search algorithms, the proposed solution-based tabu search iteratively evaluates the neighborhood of the current solution and performs the best non-tabu neighborhood move. If it is unable to improve the best solution found for a given number of iterations, a perturbation procedure is launched. The above procedures are executed repeatedly until the time limit is reached. However, unlike the traditional attributebased tabu search, the proposed TSMH algorithm employs a solution-based tabu strategy which records the complete information of each evaluated solution. If we regard the limited characteristics in the attribute-based tabu list as directions or landmarks in the solution space, the ones in the solutionbased tabu list will be the exact coordinates. Next, we will present comprehensive descriptions for the essential ingredients of the solution-based tabu search subroutine.

Neighborhood Structure and Evaluation

We define four neighborhood structures to optimize the delivery schedule. Unlike the classical ones for routing problems, these neighborhoods manipulate the timing of visits to the clients (Z t i ) instead of reordering the routes. At each iteration of the tabu search, the current solution S is replaced with its neighboring solution by performing the best move M selected from all the following neighborhoods.

• Addition operation Ma(t, i). Visit unvisited client i in period t. There are O(pn) addition operations.
• Removal operationMr(t, i). Cancel the existing visit to client i in period t. There are O(pn) removal operations.
• Move operation Mm(t1, t2, i) = Ma(t1, i) +Mr(t2, i). Visit client i in period t1 if it is not visited, and simultaneously cancel the existing visit to the same client in period t2. There are O(p2n) move operations.
• Swap operation Ms(ti, tj, i, j) = Mm(tj, ti, i) + Mm(ti, tj, j). Exchange the visits to a pair of clients in a pair of periods where Z ti i = Z tj j = 1,Z tj i = Z ti j = 0. There are O(p2n2 ) swap operations.

The neighborhood evaluation is the bottleneck of the tabu search algorithm, since it evaluates a lot of moves but performs only one of them at each iteration. In order to improve the search efficiency, our TSMH algorithm adopts a filtering technique and only applies exact evaluations on promising moves. The proposed filtering technique evaluates the entire neighborhood in an incremental and approximate way, and then it abandons the moves with poor approximate objective improvements. Specifically, the approximate objective function only considers the routing costs, and the objective improvements ∆M of move M are calculated as Eqs. (12).

Eqs. (12) are incremental evaluations which assume that only the predecessors i1 and successors i2 of the added or removed visits will be affected in the original tour. For example, when removing client i from tour (..., i1, i, i2, ...) in period t, it assumes that the optimal tour for the resulting vertex combination is (..., i1, i2, ...). Based on Eqs. (12), we can sort the neighborhood moves by ∆M values in an ascending order, and only take the top-m elite moves into further consideration. Thus, the filtering technique can efficiently identify the non-promising moves, at the cost of possible inaccuracy.

Next, the proposed TSMH figures out the exact holding costs of the elite solutions. Since the visiting states are the only band to connect the inventory and routing subproblems as Eqs. (4) impose, the holding and routing costs are independent to each other once the visiting states Z t i are fixed, which is true at this moment. So, we just leave out the routing part and minimize the inventory holding costs fh(S) by solving the LP model composed of (5)-(11) and (13)-(14). Note that we treat the move as the worst if its LP model is infeasible.

Finally, we make the move with the smallest fh(S ⊕M)+ ∆M value and re-optimize the routes in each affected period using Lin-Kernighan heuristic. As discussed in Section 3.1, the route optimization in period t is equivalent to solving the TSP in a subgraph each of whose vertex i satisfies Z t i = 1.

Solution-based Tabu Strategy

In order to jump out of the local optima during the neighborhood search, we integrate a solution-based tabu strategy into the TSMH algorithm. Comparing to its classical attribute-based counterpart, the solution-based tabu strategy records more complete information of each evaluated solution. Therefore, we can omit the tabu tenure parameter and prohibit visiting the evaluated solutions forever. This feature overcomes the drawbacks of the attribute-based tabu search that inappropriate parameter settings may falsely forbid some promising moves or accept moving back to visited solutions. For efficiency consideration, we only store the hash values of the evaluated solutions instead of the entire solutions. We define a hash function h : S → {0, 1, 2, ...,K - 1} to map the solution vector S into non-negative integers ranging from 0 to K - 1. Then the tabu list is implemented by a Boolean vector B containing K entries, where Bk is true indicates that the solution with the hash value h(S) = k is in the tabu state. Due to the possibility of hash collision, we employ multiple hash functions hj and corresponding tabu lists Bj for j = 1, 2, ..., J to improve the robustness of the solutionbased tabu strategy. In the proposed TSMH algorithm, we take J = 3, and a solution S is in tabu state iff ∧Jj=1 Bjhj (S) is true. Correspondingly, we need to set the above three entries to true for each evaluated solution. We define hash function hj as follows, where γj is a parameter.

Perturbation Operator

When the tabu search fails to improve the current best solution for consecutive α iterations, TSMH utilizes a perturbation operator to escape from the local optima. The perturbation starts from the best solution found so far, and successively performs πa addition, πr removal, and πm move operations. Each operation is selected from all neighborhood moves in the corresponding neighborhood from a uniform distribution. In addition, an operation that neutralizes the effect of another one will not be selected. If the resulting solution is infeasible or in tabu state, the perturbation will restart until a proper solution is obtained.

4 Computational Experiment and Analysis

4.1 Experimental Protocols

The proposed TSMH algorithm is evaluated on two sets of most commonly used datasets, which are 160 small instances proposed by Archetti et al. [2007] and 60 large ones proposed by Archetti et al. [2012], respectively. In the small-scale dataset, there are 100 three-period instances with 5, 10, ..., 50 clients and 60 six-period ones whose numbers of clients are 5, 10, ..., 30. The large instances consist of 6 periods and 50, 100, 200 clients. Each configuration can be further divided into two groups according to the distributions of their unit inventory costs, which are [0.01, 0.05] (low cost) and [0.1, 0.5] (high cost), respectively. In addition, the traveling cost between each pair of vertices i, j is calculated by rounding the Euclidean distance between them.

Our TSMH algorithm is developed in C++ and tested on Intel Xeon E5-2698 v3 2.30GHz CPU. The MIP models are solved by Gurobi 8.0 [Gurobi Optimization, 2018], and the TSP subproblems are solved by LKH3 [Helsgaun, 2017]. We perform 20 independent runs for each instance under a onehour time limit and record the best results. We also extend the time limit to 3 hours to investigate the potential of the proposed TSMH algorithm. The time limit for the first stage is 1/6 of the total time limit. The values of other parameters for TSMH are listed in Table 1.

4.2 Computational Results

In order to justify the effectiveness and efficiency of the proposed algorithm, we compare TSMH with the best known exact algorithm CL-BC [Coelho and Laporte, 2013; Coelho and Laporte, 2014] and heuristic algorithm HAIR [Archetti et al., 2012]. The computational platforms for CL-BC and HAIR are Xeon 2.66GHz and Intel Dual Core 1.86GHz CPU, respectively. By assuming the hardware performance is proportional to the CPU frequency, we will scale down the computational time of HAIR to 80% for a fair comparison.

The overall computational results of each group of instances are shown in Table 2. Column #best indicates the number of instances where the corresponding algorithm can match the best known objective values. Column Gap reports the average gap between the best result obtained by each algorithm and the best known one. Column Time gives the normalized average computational time in seconds for obtaining the reported results. From Table 2 we can observe that, TSMH is a robust algorithm whose average gap on all 220 instances are 7 and 150 times smaller than HAIR and CLBC. Moreover, the proposed algorithm is highly scalable that it obtains most best known solutions on the two sets of large instances. In addition, TSMH is also very efficient as the average computational time is 3 times shorter than CL-BC.

Tables 3 and 4 present the detailed results and comparisons on the 60 large instances. Column Instance gives the reformatted name of each large instance according to their characteristics. For example, p6c2n50.4 indicates the fourth instance where there are 6 periods and 50 customers and the scale of the unit inventory costs are 10-2 (low cost). The next five columns report the objective values obtained by each algorithm under the specified time limits. Column Time presents the computational time taken by TSMH to obtain its best results. The numbers in bold indicate the best known results. As we can see, although CL-BC obtains the best results on all 20 instances with 50 clients under a four-hour time limit, TSMH shows its advantage on larger instances. As a result, TSMH obtains the best known results on 40 instances, while CL-BC only found 26 ones. Regarding the one-hour time limit, TSMH outperforms HAIR in terms of both the number of best solutions and the average objective values. In sum, the proposed TSMH algorithm is competitive comparing to the state-of-the-art algorithms in the literature.

4.3 Effectiveness of Combining Two Stages

To demonstrate the importance of each stage, we conduct extensive experiments to compare the complete TSMH algorithm with its two simplified versions TSMH-1 and TSMH-2 under one-hour time limit. TSMH-1 only includes the structural optimization stage. TSMH-2 generates the initial solution by solving the RIRP model without adding lazy constraints, and incorporates the detailed refinement stage based on the tabu search algorithm. Figure 2 plots the distribution of the solution qualities obtained by these three versions. We can observe that it fails to match any best known result if there is only a single stage. In addition, the solution-based tabu search obtains relatively better results comparing to the relax-and-repair method. These results highlight the merit of the hybridization of different methodologies.

5 Conclusion

We proposed a two-stage matheuristic algorithm which includes a structural optimization stage and a detailed refinement stage. Tested on 220 most commonly used instances, the TSMH algorithm improves the best known results on 40 out of 60 large instances. In addition, extensive experiments show that the combination of the two stages is essential to the performance of TSMH.

Acknowledgments

Research partially supported by Alkaid lab of Huawei Cloud.