An iterative matheuristic for the inventory routing problem

(1)

Computers & Operations Research 131 (2021) 105262

Available online 8 March 2021

An iterative matheuristic for the inventory routing problem

Simen T. Vadseth

^*

, Henrik Andersson , Magnus Stålhane

Department of Industrial Economics and Technology Management, Norwegian University of Science and Technology, Alfred Getz veg 3, 7491 Trondheim, Norway

A R T I C L E I N F O Keywords:

Transportation Inventory routing Matheuristic

A B S T R A C T

The paper considers the inventory routing problem with the Maximum Level replenishment policy. Here, the supplier is in charge of replenishing goods to a number of customers and can decide when, and in what order, these customers should be visited over a defined time period. The goal is to minimize transportation costs and inventory holding costs at both the supplier and the customers. We present a matheuristic that uses a giant tour and simple operators to heuristically create routes that are used in a path-flow formulation. The proposed method iterates between solving a path-flow model with a small set of routes and updating the route set based on the optimal solution from the previous iteration. Computational results on known benchmark instances show that it outperforms state-of-the-art exact methods and heuristics on larger and more difficult instances. It finds the best-known solution on 179 out of 240 larger multi-vehicle benchmark instances, where 178 of them are strictly improving upon the previously best-known solution, and does so in considerably shorter time compared with other methods. In addition, when tested on another set of benchmark instances consisting of 798 smaller instances, the matheuristic finds the optimal solution in 44.7% of the 642 instances with known optima and has an average gap of 1.75% on the others. It also improves the best-known solution of 14 out of 156 open instances.

1. Introduction

With an increasing amount of online retailing and goods delivery, research on efficient routing and inventory control is highly relevant and applicable. Better goods transportation and inventory control can lead to significant savings for companies and potentially lower greenhouse gas emissions. The focus of this paper is on the standard inventory routing problem (IRP) with the Maximum Level replenishment policy, a well- studied problem described in great detail by several authors including Archetti et al. (2017), Coelho et al. (2014), Adulyasak et al. (2014), Desaulniers et al. (2016) and Archetti et al. (2017). The IRP is a part of a business practice called vendor-managed inventory (VMI) where the supplier decides how much quantity of goods it should deliver to each of its customers and when to do so. The entire supply chain may benefit as a result of this comprehensive planning as it can lead to better routing and inventory control, while ensuring that the customers’ storage limits are respected. However, the IRP has proven to be a very challenging and computationally hard problem to solve.

Several papers have been written on the standard IRP, and there exist two relevant surveys from the last decade. The first one by Andersson et al. (2010) focused on different applications of the IRP, while the second one by Coelho et al. (2014) studied the methodological aspects.

Integrating inventory management and vehicle routing in the scientific literature started with the paper of Bell et al. (1983), and the first exact method on the standard IRP itself was proposed by Archetti et al. (2007).

The authors used a branch-and-cut algorithm to solve the single-vehicle IRP and solved instances up to 50 customers with three time periods and up to 30 customers with six time periods to optimality. They also introduced one of the two sets of benchmark instances that most researchers have worked on since. This set of instances consists of 798 small instances for up to five vehicles. The multi-vehicle version of the standard IRP was solved exactly by Coelho and Laporte (2013) and Adulyasak et al. (2014) with branch-and-cut algorithms, while Desaul- niers et al. (2016) used a branch-cut-and-price algorithm. Avella et al.

(2018) defined a new generic family of valid inequalities for the IRP and solved the problem with a branch-and-cut algorithm using two specific subclasses of the proposed valid inequalities. This method tightened the duality gap of several of the small benchmark instances.

Two additional branch-and-cut algorithms were developed by Guimar˜aes et al. (2020) and Manousakis et al. (2020). The former employed two techniques to find improved primal solutions during the branch-and-cut search. The authors were able to close the duality gap for several instances and found new best-known solutions for 129 instances in the set of small benchmark instances. The latter proposed a two-

* Corresponding author.

E-mail address: [email protected] (S.T. Vadseth).

Contents lists available at ScienceDirect

Computers and Operations Research

journal homepage: www.elsevier.com/locate/cor

https://doi.org/10.1016/j.cor.2021.105262

Received 20 October 2020; Received in revised form 3 February 2021; Accepted 19 February 2021

(2)

commodity flow formulation for the problem and was with the help of a good starting heuristic able to improve the best-known solution of 139 instances in a set of 300 large benchmark instances introduced by Archetti et al. (2012).

None of the exact methods described above have been able to solve larger multi-vehicle instances to optimality. Further, the duality gap can become very large when the size of the instances increases, and often no feasible solution is found within a reasonable time frame when an exact method is used. There is, consequently, a need for good heuristics for the IRP, and several have been developed over the years. Coelho et al.

(2012b) developed an adaptive large neighborhood search (ALNS) heuristic for the IRP with transshipment, which the authors also tested on the single-vehicle IRP. An ALNS for the multi-vehicle version has been designed by Adulyasak et al. (2014), who solved instances with two and three vehicles.

The complexity of the IRP has led most researchers to integrate mathematical programming techniques into their heuristics. These types of heuristics, regardless of the type of problems they are applied to, have come to be known as matheuristics. A definition of what a matheuristic is was proposed by Boschetti et al. (2009) as: “Matheuristics are heuristic algorithms made by the interoperation of metaheuristics and mathematical programming techniques”. This is also the definition used by Archetti and Speranza (2014) in their survey on matheuristics for routing problems. The authors classified matheuristics in three cate- gories: Decomposition approaches, improvement heuristics and branch- and-price/column generation-based approaches. They also show that matheuristics are used to a large extent on problems similar to the IRP, e.

g. the vehicle routing problem (VRP), the production routing problem and the location routing problem.

To the authors’ knowledge, the first matheuristic used to solve the standard IRP was developed by Archetti et al. (2012). It is a hybrid heuristic that combines tabu search with the solution of mixed integer programs (MIP) and can be classified as an improvement heuristic. The authors studied the single-vehicle case and, as previously mentioned, released the second set of benchmark instances that most researchers use today. The same authors extended the method for the multi-vehicle version in Archetti et al. (2017). It is both a decomposition and an improvement search, combining a tabu search heuristic with solving MIPs. The solutions of 92% of the larger multi-vehicles instances were improved. However, many of these solutions were further improved by Chitsaz et al. (2019) with their three-phase decomposition matheuristic which relies on the iterative solution of different subproblems. Although designed for the assembly routing problem, the algorithm was able to find new best-known solutions for 194 out of 300 large instances for the IRP. Another matheuristic was proposed by Alvarez et al. (2020). They developed a hybrid heuristic, combining an iterative local search metaheuristic and two mathematical programming components, to solve the IRP with perishable products. The authors also tested the algorithm on the standard IRP and improved the best-known solution for a few smaller instances. An additional matheuristic for the IRP was presented by Diniz et al. (2020). Here, the authors combined an iterative local search with a randomized variable neighborhood descent, and were able to find and improve the best-known solutions of several small benchmark instances.

As seen in the paragraphs above, there are several exact methods and heuristics designed for, and applied to, the IRP. Smaller instances of the IRP can be solved to optimality by exact methods, while heuristics outperform exact methods on larger instances. Even though there exists heuristics for the IRP, there is still room for improvement, especially when it comes to larger instances. Many of the proposed solution methods cannot be applied to the largest instances due to the computational complexity, and those that can, suffer from long computing times.

The purpose of this paper is to present a new matheuristic to solve large instances of the IRP in shorter computing times. The matheuristic iteratively solves an exact mathematical model with a limited number of

routes. The set of routes used in the first iteration is generated from a giant tour, and is modified by different operators between each iteration. It is tested on known benchmark instances for the IRP and finds new best solutions on 178 out of 240 instances with multiple vehicles, and a further seven new best-known solutions for the single-vehicle case.

These improved solutions have also been found in only a small fraction of the time spent by other heuristics in the literature.

The remainder of the paper is organized as follows. In Section 2, the standard IRP is defined and presented mathematically, while our matheuristic is presented in detail in Section 3. Our computational results are reported in Section 4 and concluding remarks are presented in Section 5.

2. Problem definition and formulation

The inventory routing problem concerns the repeated distribution of goods from a supplier to a set of customers over a given planning horizon. We formulate this problem on a graph G(N,A)where N is a set of nodes N = {0,1,…,N} consisting of N customers and a supplier denoted 0. We also introduce N^′= {1,…,N}as the set of customers.

The set of arcs A defines movements between each pair of nodes. The problem is defined over a time horizon T = {0,1,…,T}and we also introduce the set of planning time periods T^′ = {1,…,T}.

In each time period, V vehicles with capacity Q can be used to deliver the goods. There is a driving cost Cij, associated with each arc (i,j) ∈A. Customer i has a known demand for the commodity, Rit, in each time period t and a maximum, Ui, and minimum, Li, inventory capacity. The supplier produces R0t units of the commodity at the beginning of each time period t. Both the supplier and the customers have an inventory holding cost C^H_i per unit of commodity at the end of each time period.

Each customer can only be visited once per time period. The problem consists of minimizing the transportation and inventory costs of the entire supply chain while making sure that no stock-outs occur.

2.1. A path-flow formulation

The proposed path-flow formulation requires some additional nota- tion. The set R contains all routes. A route is a Hamiltonian cycle through a subset of the nodes including the supplier. Introducing Aijr as 1 if route r traverses arc (i,j), and 0 otherwise, the cost of route r can be defined as C^T_r =∑

(i,j)∈ACijAijr. The variable λrt is 1 if route r is used by a vehicle in time period t, and 0 otherwise. The amount of commodity delivered at node i in time period t is denoted qit and the inventory level at node i at the end of time period t is denoted sit. The inventory at node i at the beginning of the planning horizon is represented by Ii0. Finally, let lijt be the flow of commodity on arc (i,j)in time period t. With this no- tation, the model can be formulated as follows:

min∑

i∈N

∑

t∈T^′

C^H_isit+∑

r∈R

∑

t∈T^′

C^T_rλrt (1)

si0=Ii0 i∈N (2)

s0t− s0(t−1)− R0t+

∑

i∈N^′

qit=0 t∈T ^′ (3)

sit− si(t−1)+Rit− qit=0, i∈N^′,t∈T ^′ (4)

Li⩽sit⩽Ui i∈N,t∈T (5)

si(t−1)+qit⩽Ui i∈N^′,t∈T^′ (6)

∑

j∈N

ljit− qit− ∑

j∈N

lijt=0 i∈N^′,t∈T

′ (7)

(3)

lijt− Q∑

r∈R

Aijrλrt⩽0 (i,j) ∈A,t∈T^′ (8)

∑

r∈R

∑

j∈N

Aijrλrt⩽1 i∈N^′,t∈T^′ (9)

∑

r∈R

λrt⩽V t∈T^′ (10)

λrt∈ {0,1} r∈R,t∈T^′ (11)

qit⩾0 i∈N^′,t∈T^′ (12)

lijt⩾0 (i,j) ∈A,t∈T^′ (13)

The objective function (1) minimizes the transportation and inventory holding costs over the entire planning horizon, while constraints (2) set the starting inventory level at each node. The inventory balance at the supplier and the customers are taken care of by constraints (3) and (4). Moreover, the upper and lower limits on the inventory level at each node are handled by constraints (5) and (6). Constraints (7) ensure that the flow of goods out of a node is equal to what comes in except for the amount that is delivered. Constraints (8) make sure that the flow on an arc does not exceed the vehicle capacity, while constraints (9) state that two routes that visit the same node are not used in the same time period.

With constraints (10) we make sure that the maximum number of available vehicles is not exceeded. Moreover, constraints (11) state that a route is either used or not while constraints (12) and (13) impose non-

negativity for the other variables.

2.2. Valid inequalities

The following valid inequalities, first introduced by Coelho and Laporte (2014) for the IRP, have been adapted for the path-flow formulation and added to the model.

∑

r∈R

∑^t²

t^′=t1

Airλ_rt^′⩾⌈

∑_t

2

t^′=t1R_it^′− Ui

min{Q,Ui} ⌉ i∈N^′,t1,t2∈T ^′, t2⩾t1 (14)

∑

r∈R

∑^t²

t^′=t1

Airλ_rt^′⩾

∑_t₂

t^′=t1R_it^′− si(t1−1)

min{Q,Ui,∑_t

2

t^′=t1R_it^′} i∈N^′,t1,t2∈T^′, t2⩾t1 (15) Constraints (14) state that if the sum of the demands in a node over time periods t1 to t2 is greater than the inventory limit, then there must be at least one visit to this node in the interval. Constraints (15) state the same, but here the actual inventory at the node at the start of time period t1 is taken into account instead. Rounding this up is not possible since the expression then becomes non-linear. The inequalities are also strengthened by adding the total demand for node i to the denominator.

3. Matheuristic

The formulation given in Section 2.1 is a valid and complete formulation for the IRP given that the set of routes, R, includes every feasible route in the graph G. However, the number of feasible routes

Fig. 1.The matheuristic first creates a set of routes using a giant tour and a split algorithm. These routes are used to solve a modified version of the path-flow model.

The model is from there on solved iteratively where the routes of the current best solution are used to create a new set of routes.

(4)

grows exponentially with the number of customers and thus becomes so large that, even for small instances, it is not possible to generate all of them within a reasonable amount of time. However, the number of routes used in any feasible solution is bounded by V⋅T. Thus, it is possible to obtain feasible (and optimal) solutions to the problem by replacing R with a small set of routes ̂R, given that |̂R|is in the same magnitude as V⋅T.

The matheuristic presented in this paper exploits this fact by iteratively generating a set of promising routes and then solving the path- flow model presented in Section 2.1 using this subset. The outline of our heuristic solution method is illustrated in Fig. 1. The method starts by generating a giant tour, which is split into routes in different ways to generate an initial set of promising routes. Then, for a number of iterations, the method alternates between solving the path-flow model and a method that updates the set of routes based on the optimal solution to the path-flow model. In the following, we go through the details of each part of the matheuristic. Section 3.1 describes how the initial set of routes is generated. In Section 3.2, we explain how the set of routes is modified based on the previous solution to the path-flow model.

3.1. Generating an initial set of routes

The first step of the heuristic is to generate an initial set of promising routes. The details of this step are outlined in Algorithm 1. First, a giant tour (GT) with the minimal total distance on the graph G = (N^′,A)is created by solving a travelling salesman problem using the function solveTSP(N^′,A). In the implementation, we obtain the giant tour by using the TSP-solver released by Helsgaun (2009), which is considered to be the fastest implementation of the Lin-Kernighan algorithm (Lin and Kernighan, 1973).

After the giant tour is created, it is split into segments, and routes are created by inserting copies of the supplier node at the start and end of each segment. To do this, we utilize the split algorithm proposed by Vidal (2016) for the capacitated vehicle routing problem (CVRP). The split algorithm aims to partition a giant tour solution into separate routes. It does this by solving a shortest path problem on an acyclic graph ̂G = (N,̂A) where ̂A includes one arc (i,j)with cost CÂ=

C0(i+1)+∑ ij

k=i+1,…,j−1Ck(k+1)+Cj0 for any feasible route visiting customers i+1 to j. Here, the node order follows the ordering of the giant tour where 1 represents the first node in the giant tour. Traditionally, a version of Bellman’s equation has been used to solve the problem, but the aforementioned paper introduces a more efficient labeling algorithm with additional dominance rules. By using the open-source code released by the author, we are able to solve the splitting problem with a limited fleet size in O(N⋅V).

However, to utilize the split algorithm on the obtained giant tour, the giant tour must be turned into a sequence of nodes by choosing a starting node (i^Start). In addition, since, unlike the CVRP, the IRP does not have pre-determined demands at each customer, the algorithm must create an array d of customer demands. Thus, the function Split(GT,i^Start,d) returns the set of routes obtained by applying the split algorithm to a given giant tour, start node, and demand combination.

To create a diverse set of routes in the initial set of routes ̂R, the function Split(GT,i^Start,d)is called n times with different input combi- nations. In a given iteration, the algorithm first selects the start node using the getStartNode(j)function, which returns the j-th closest node to the supplier. Since the first node in the giant tour is the first node on the first route, it is likely that this node should lie close to the supplier, and the number j is thus varied between 1 and 5. To create customer demand d[i]at node i in a given iteration, the algorithm multiplies the upper inventory limit Ui with a percentage P. This percentage is reduced for each iteration of the algorithm by a given value D. Thus, we get a varied set of routes as splits of the giant tour using large customer demands give short routes, and splits using small customer demands give longer

routes.

Algorithm 1: Generating initial set of routes 1. initialize d – array of customer demands 2. initialize P – percentage

3. ̂R =∅

4. GT=solveTSP(N^′,A) 5. for n iterations do

6. i^Start =getStartNode((nmod5) +1) 7. for i∈N^′do

8. d[i] =P⋅Ui

9. end for

10. ̂R =̂R∪Split(GT,i^Start,d); 11. P=P−D

12. end for 13. return ̂R

3.2. VRP heuristic and operators Algorithm 2: Updating route set

1. input: λ^*_rt,r∈̂R,t∈T 2. input: q^*_it,i∈N^′,t∈T 3. ̂R = {r∈R:∃t∈T,λ^*_rt=1} 4. for t∈T do

5. ̂R =̂R⋃

solveVRP(q^*_it) 6. end for

7. ̂N =∅ 8. ̂R=̂R⋃

RemoveNodes(̂R,̂N)(Algorithm 3) 9. ̂R=̂R⋃

InsertNodes(̂R,̂N)(Algorithm 4) 10. for r∈̂R do

11. r=solveTSP(r) 12. end for 13. return ̂R

After a solution from the path-flow model is obtained, the route set is updated with the aim of improving the solution of the path-flow model in the next iteration. A pseudo-code outlining how the route set is updated is given in Algorithm 2. First, the algorithm solves a VRP heuristically for each time period (lines 4–6), where the optimal quantities from the last solution of the path-flow model, q^*_it, are used to define the quantities delivered to each customer. This gives (near) optimal routing for these quantities, which may be an improvement on the initial route set. The new routes are added to the set of routes, ̂R, which is added to the path-flow model in the next iteration together with the routes given from the last solution. In this paper, the genetic algorithm proposed by Vidal et al. (2012) has been used as the VRP-heuristic. The authors of this paper have used the open-source implementation which is described in Vidal (2020).

Next, nodes from each of these routes are removed using the method RemoveNodes(̂R,̂N)(line 8). The resulting routes are added to the set R̂. In addition, every node that is removed from a route is added to the set ̂N. These nodes are then re-inserted into different routes in the method InsertNodes (line 9). The resulting routes from these insertions are included in the set ̂R. The two methods are described in detail in Algorithms 3 and 4.

Once this process is completed, all routes to be used in the path-flow model in the next iteration have been generated. However, as a final step, the TSP-heuristic, solveTSP(r), is used on each route r to (possibly) reduce the distance driven (line 10–12). This is done because even though a node is removed from or inserted into a route in the cheapest possible way there is no guarantee that the resulting route is optimal.

Finally, the route set ̂R is returned and the path-flow model can be re- solved. Since the new set of routes, R̂, includes the optimal routes from the previous iteration, we know that this route set provides a solution that is at least as good as in the previous iteration. This is ensured by

(5)

warm starting the solver with the solution of the previous iteration.

Algorithm 3: RemoveNodes 1. input: ̂R,N ̂

2. R^new=∅ 3. for r∈̂R do

4. for k∈K⧹{RandomRemoval}do 5. r^′ =N^k(r)

6. ̂N =̂N⋃ (N̂_r⧹N̂_r^′) 7. R^new=R^new⋃

{r^′} 8. end for

9. for cnt=1,…,Y do 10. k =Random Removal 11. r^′ =N^k(r)

12. N̂=̂N⋃ (N̂_r⧹N̂_r^′) 13. R^new=R^new⋃

{r^′} 14. end for

15. end for 16. return R^new

The method RemoveNodes is described in Algorithm 3 and consists of using a set, K, of operators on each route in the set ̂R. These operators, represented by the functions N^k(r), remove one node from a route in different ways before returning the resulting route. The set of operators K consists of Cheapest Removal, Least Served Removal, Most Served Removal and Random Removal. All operators have time complexity O(n), except Random Removal that has O(1).

– Cheapest removal: Removes the node from the route that gives the largest reduction in route cost, if removed.

– Least served removal: The node in the route that received the smallest quantity of goods in the last solution of the path-flow model is removed from the route.

– Most served removal: The node in the route that received the largest quantity of goods in the last solution of the path-flow model is removed from the route.

– Random removal: A randomly selected node is removed from the route.

An element of stochasticity is added to the algorithm in the form of Random Removal so that it is highly unlikely that two consecutive iter- ations of the path-flow model optimize over the same set of routes ̂R. Unlike the other operators, Random Removal is run several times. The number of times Random Removal is used per route is a parameter denoted Y.

All nodes removed from a route are added to the set N ̂(line 6) and

̂N_ris the subset of nodes visited on route r. The routes generated using the removal operators are added to the set R^newwhich thereafter is included in the set ̂R.

Algorithm 4: InsertNodes 1. input: ̂R,N ̂ 2. R^new2=∅ 3. R¹,R²⊆̂R 4. for n∈̂N do 5. R^new2=R^new2⋃

{I(R¹,n)}⋃ {I(R²,n)}

6. end for 7. for r∈R¹do 8. R^new2=R^new2⋃

{C(r)}

9. end for 10. return R^new2

The method InsertNodes is described in Algorithm 4 and consists of inserting the nodes in N ̂back into different routes and hence creating new routes for the path-flow model. Here, the increase in route cost when inserting a node into a route decides which route it is inserted into.

This is described as function I(R^*,n), which receives a set of routes, R^*, and a node n. The function calculates the cheapest position to insert the node into each route. The node is inserted into the route with the lowest marginal cost and the resulting route is returned by the function. I(R^*, n)is run twice - once for the routes in the current optimal solution of the path-flow model and the routes generated by the VRP-heuristic (described by the set R¹), and once for the routes generated by RemoveNodes (described by the set R²). By partitioning the route set ̂R into two subsets, we ensure to create routes that are longer than our original ones and that nodes are inserted into routes that have been shortened.

InsertNodes also includes an operator named Closest Insertion which is described as function C(r)in the algorithm. Closest Insertion takes a route and finds the closest neighbor to each of its nodes. All the closest neighbors are added to the route given that they are not already included. The resulting route is added to R^new2which thereafter is included in the set ̂R.

4. Computational study

To evaluate the proposed matheuristic, we have tested it on known benchmark instances for the IRP found in the literature. Section 4.1 introduces the benchmark instances used. In Section 4.2, implementation issues and parameter testing is discussed. The computational results are presented in Section 4.3, while the effect of the improvement phase is investigated in Section 4.4.

4.1. Benchmark instances and solution methods

Two sets of benchmark instances have been used. The first set of instances are those created by Archetti et al. (2007) for the IRP with a single vehicle. It consists of 100 instances with three time periods ranging from 5 to 50 customers, where one-half of the instances has high inventory costs while the other half has low. For the rest of this paper, we call this subset of instances SV-S3 (single-vehicle - small 3). There are also 60 instances with six time periods ranging from 5 to 30 customers.

Again, one half of them has high inventory costs while the other half has low. These are called SV-S6. The instances were modified for the multi- vehicle version by Adulyasak et al. (2014) and Coelho et al. (2012a) by dividing the original vehicle capacity by the number of vehicles available and rounding to the nearest integer. Up to five vehicles are considered. These are called MV-S3 and MV-S6, where M stands for multiple. The second set of benchmark instances was created by Archetti et al. (2012). The instances are larger and more challenging, consisting of 100 instances with 50 customers, 100 instances with 100 customers and 100 instances with 200 customers. They are ranging from one to five vehicles and one half of them has high inventory costs, while the other half has low. We call the single-vehicle version SV-L6 and the multi- vehicle version MV-L6 where L stands for large. All the instances have six time periods.

To evaluate the proposed matheuristic, our results are compared with the results presented in other papers that have used the same instances. Some of these methods are exact and some are heuristics. Most of the exact methods have focused on solving the smaller instances, while the newer heuristics have also focused on finding good solutions to the larger ones. As mentioned previously, our main focus has been on improving the solution quality of the larger instances. However, the heuristic is also tested on the smaller ones to evaluate its performance on instances with known optima. The methods and results presented in the literature span several years and have all used different CPUs and soft- ware. Hence, making a fair comparison of their time consumption is hard. In Table 1, we summarize the other solution methods and their most important features. A complete overview of which instances each paper has solved can be found in Table 2.

(6)

4.2. Algorithmic implementation and parameters

The parameter values used in this computational study are a result of preliminary testing, where different parameter configurations of the matheuristic were tested and compared. The final parameter values selected are presented in Table 3. There are big differences in the number of customers, vehicles, and time periods between the various instances and different parameter settings are thus more suitable for different sets of instances. However, to avoid overfitting of the parameters to the benchmark instances, we have chosen to keep the parameters, except for the time limit, the same for all sets of instances.

The number of iterations, n, of the split algorithm is set as a result of the initial value of P which is multiplied by Ui, and the percentage, D, which P is decreased by in each iteration. The preliminary testing showed that having P larger than 80% usually results in the split algorithm splitting the giant tour into routes visiting a single customer while

having P less than 40% gives a single route visiting all customers. Hence, having P outside of this range rarely results in additional routes.

Decreasing P by 2.5% per iteration is, in most cases, sufficient to alter the routes between each iteration, and having P starting at 80% and ending at 40% gives n =16.

Sometimes, commercial solvers can spend a large amount of time to close the duality gap of a MIP without improving the primal solution.

Since the path-flow model is solved with a small set of routes, the dual bound is not meaningful, only the primal solution is. Therefore, a maximum time limit on the time spent solving the path-flow formulation in each iteration of the matheuristic is set. The parameter, 1. Time limit, refers to the time limit put on the first run of the path-flow model, while Time limit refers to the consecutive runs as described in Section 3. The number of iterations and the number of times the Random Removal operator is used in RemoveNodes are both set based on the preliminary testing.

Moreover, there are stochastic elements in the matheuristic, and hence two different runs on the same instance will not necessarily lead to the same solution. As a consequence, every instance has been run ten times in our computational study.

4.3. Computational results

In this section, the results of the proposed matheuristic are compared with the results of the benchmark solution methods. First, a comparison between the results obtained by our matheuristic and the best-known solutions for the smaller benchmark instances is made in Section 4.3.1. Since the majority of these have been solved to optimality, a comparison can give valuable insight into how close to optimality the obtained solutions are. In Section 4.3.2, we present the results for the larger instances to see how good the proposed matheuristic is compared with existing heuristics and exact methods. The detailed computational results, and solutions, can be found at http://axiomresearchproject.

com/publications/.

Table 1

Benchmark solution methods. We present the solution approach, running plat- form, number of threads and standard MIP solver. Note: Sol: Solution approach, E: Exact, H: Heuristic/Metaheuristic, M: Matheuristic, Def: Default.

Reference Name Sol CPU #Threads Solver

Archetti et al.

(2007) A-BC E Pentium IV 2.8

GHz Def Cplex 9.0

Coelho and

Laporte (2013) CL-

BC E Xeon 2.66 GHz 6 Cplex 12.3

Desaulniers et al.

(2016) D-

BPC E Core i7-2600 3.4

GHz 1 Cplex 12.2

Avella et al.

(2018) AV-

BC E Core i7-2620,

2.70 GHz 1 Xpress 7.6

Guimar˜aes et al.

(2020) G-BC E Xeon E5-2630 v2

2.60 GHz 6 Gurobi 8.1

Manousakis et al.

(2020) M-BC E Intel Core i7- 7700 CPU 3.60 GHz

8 Gurobi 8.1

Coelho et al.

(2012b) CL-H H Intel T7700, 2.4

GHz Def –

Adulyasak et al.

(2014) AB-H H 2.10 GHz Duo

CPU PC Def Cplex 12.3

Archetti et al.

(2012) AR-

H1 M Intel Dual Core

1.86 GHz Def Cplex 10.1

Archetti et al.

(2017) AR-

H2 M Xeon W3680,

3.33 GHz 8 Cplex 12.5

Chitsaz et al.

(2019) C-H M Xeon X5650 2.67

GHz 1 Cplex 12.6

Alvarez et al.,

2020 AL-H M Xeon X5650 2.67

GHz 1 Cplex 12.8

Diniz et al.

(2020) D-H M Intel Core i7-

8700 K 3.7 GHz 1 LEMON

library This paper V-H M Xeon Gold 6144

3.5 GHz 1 Gurobi 9.0

Table 2

Number of instances solved by each solution method. Note: V: number of vehicles.

Name A-BC CL-BC D-BPC AV-BC G-BC M-BC CL-H AB-H AR-H1 AR-H2 C-H AL-H D-H V-H

Set V Size E E E E E E H H M M M M M M

SV-S3 1 100 100 100 – – 100 – 100 – 100 – 100 100 100 100

SV-S6 1 60 60 60 – – 60 – 60 – 60 – 60 60 60 60

MV-S3 2 100 – 100 100 10 100 100 – 100 – 100 100 100 100 100

3 100 – 100 100 10 100 100 – 100 – 100 100 100 100 100

4 100 – 100 100 10 100 100 – – – 100 100 100 100 100

5 100 – 100 100 10 100 100 – – – 100 100 100 100 100

MV-S6 2 60 – 60 60 40 60 60 – 50 – 60 60 60 60 60

3 60 – 60 60 40 60 60 – 50 – 60 60 60 60 60

4 60 – 60 60 40 60 60 – – – 60 60 60 60 60

5 58 – 58 58 38 58 58 – – – 58 58 58 58 58

SV-L6 1 60 – 60 – – 60 60 – – – 60 – – 60

MV-L6 2 60 – 40 – – 60 40 – – – 60 60 – – 60

3 60 – 40 – – 60 40 – – – 60 60 – – 60

4 60 – – – – 60 40 – – – 60 60 – – 60

5 60 – – – – 60 40 – – – 60 60 – – 60

Table 3 Parameter values

Parameter Value

Split algorithm: Start Percentage, P 80%

Split algorithm: Decr. Percentage, D 2.5%

Split algorithm: Iterations, n 16

Path Flow: 1. Time Limit 2 s ⋅ (number of customers) Path Flow: Time Limit 1.5 s ⋅ (number of customers)

Path Flow: Iterations, maxIt 5

RemoveNode: number of Random Removal, Y 3

(7)

4.3.1. Small benchmark instances

The results of running the proposed matheuristic (denoted V-H) 10 times on each of the 798 small test instances are compared with existing results from the literature. All the computational results are summarized in Tables 4–6. In Table 4, the average gap of each solution method is presented. The gap of a solution method on an instance is calculated as the method’s objective value minus the best known objective value (not considering the objective value of V-H) divided by the best known objective value. The average per instance set is then calculated. For V-H, two values are given. V-H Best is the average gap of the best solution obtained over all instances in each set. For V-H Average, the average gap over the ten runs has been calculated for each instance, and then the average over these for each instance set is reported.

It is clear that the most recent exact methods, G-BC and M-BC, outperform all the heuristics in terms of solution gap. This is natural since 642 of the 798 small test instances have been solved to optimality.

In fact, the only subsets of instances where the majority of instances still have not been solved to proven optimality are MV-S6 with four and five vehicles. V-H Best and V-H Average have the smallest solution gaps for these instances among all heuristics. Overall, the results of V-H are competitive compared with the other heuristics and V-H Best has an average gap of 0.99% over all instances. D-H is the only heuristic that has a better average gap with 0.37% over all instances. However, there is a significant performance difference between the three and six time period instances for V-H. V-H Best has an average gap of 0.40% for the multi-vehicle instances with six time periods which is lower than all the other heuristics.

In Table 5, the number of best-known solutions found by each solution method is presented. V-H Best has the same meaning as above, but its column has an extra number written within parentheses. This shows the number of best-known solutions that have a strictly lower objective value than the previously best-known solution, i.e. the number of new best-known solutions found. V-H Best finds more best-known solutions than any of the other heuristics, except D-H, with 306 solutions, of which 14 are new solutions. The column V-H Worst shows the number of instances where all the ten runs of V-H have found the best-known solution. With 215 best-known solutions, it is still competitive compared with the other heuristics. As in Table 4, V-H performs significantly better on the six time period instances, compared with those that have three time periods. For the six time period instances, V-H Best finds the best-

known solution for 59 out of 180 instances for the three, four and five vehicle instances. This is significantly higher compared with the other heuristics.

Table 6 presents the average computing time for each solution method in seconds on each set of instances. V-H Average represents the average computing time over the ten runs, while V-H Max is the average of the runs with the maximal computing times for each instances. The results show that V-H’s computing time is similar to C-H and D-H, and it is substantially shorter than the other solution methods. The best exact solution method on the MV-S6 instances, M-BC, spends on average 4,210 s with an average gap of 0.02%, while V-H spends, on average, 138 s with an average gap of 0.75%. Thus, the results indicate that V-H produces close to optimal solutions within a small fraction of the computing time used by exact methods.

4.3.2. Large benchmark instances

The results of running the proposed matheuristic (denoted V-H) ten times on each of the 300 large test instances are compared with existing results from the literature. The results are summarized in Tables 7–9. In Table 7, average gaps are presented. Here, V-H Best and V-H Average have the same meaning as in Section 4.3.1. Both G-BC and M-BC as well as AR-H1 produce lower gaps for the single-vehicle instances. However, the results for V-H Best, an average gap of 1.65%, are in fact quite good taken into account that in this subset of instances all 50 customer instances are solved to optimality and most 100 customers instances have a duality gap of less than one percent for its best-known solution. On the multi-vehicle instances, V-H reports lower average gaps than all the other solution methods except for two vehicles. V-H Best has an average gap of − 0.54% for MV-L6, while M-BC, which has the second-lowest average gaps, has an average gap of 0.28%. However, M-BC has not been tested on the 200 customers instances which are the largest and most challenging ones. In fact, the authors themselves state that they are not solved, because their exact solution method provides no insightful results within 2 h on these instances. C-H, which has been tested on all instances, has in comparison a gap of 1.24%.

In Table 8, the number of best-known solutions for each solution method is presented. V-H Best and V-H Worst have the same meaning as in Section 4.3.1 and the extra number shows the number of strictly improving best-known solutions. V-H Best has the best-known solution for 179 out of 240 multi-vehicle instances and naturally outperforms the

Table 4

Average gaps for the different solution methods. The best average gap for each subset of instances is highlighted.

Set V # Instances A- BC CL-

BC D-

BPC AV-

BC G-

BC M-

BC CL-

H AB-

H AR-

H1 AR-

H2 C-H AL-

H D-H V-H

Best V-H

Average

SV-S3 1 100 0 0 – – 0 – 0.44 – 0 – 1.93 1.03 0.02 0.16 0.16

SV-S6 1 60 0 0 – – 0 – 0.49 – 0.14 – 1.09 2.68 0.15 0.77 1.15

MV- S3 2 100 – 0.01 12.32 1.13 0 0.01 – 8.48 – 0.14 1.79 2.58 0.04 1.17 1.32

3 100 – 0.77 7.92 2.90 0.03 0.05 – 9.14 – 0.38 1.30 2.70 0.11 2.07 2.38

4 100 – 3.27 6.73 7.29 0.30 0.11 – – – 1.07 1.58 2.84 0.51 1.61 2.03

5 100 – 5.86 5.92 11.24 0.85 0.06 – – – 1.71 2.08 2.89 0.65 1.47 1.85

MV-S3:

Average Gap

400 – 2.48 8.22 5.64 0.30 0.06 – 8.81 – 0.83 1.69 2.75 0.33 1.58 1.90

MV- S6 2 60 – 1.43 33.11 1.38 0 0.03 – 4.05 – 0.35 3.98 3,01 0.12 0.44 0.73

3 60 – 1.84 36.80 3.46 0.17 0.04 – 3.71 – 1.74 4.76 2.68 0.47 0.43 0.81

4 60 – 3.75 36.41 4.91 0.63 0.02 – – – 3.24 5.61 2.59 0.85 0.41 0.75

5 58 – 4.89 34.39 6.53 1.05 0 – – – 4.00 6.80 2.52 1.19 0.31 0.70

MV-S6:

Average Gap

238 – 2.98 35.18 4.07 0.46 0.02 – 3.88 – 2.33 5.28 2.70 0.66 0.40 0.75

Total:

Average Gap

798 0 2.03 18.33 4.32 0.29 0.04 0.46 7.20 0.05 1.38 2.74 2.52 0.37 0.99 1.28