Heuristics for vehicle routing problems: Sequence or set optimization?
\SetKwContinuecontinue \SetAlCapSkip0.5em \DontPrintSemicolon\SetKwInputKwDataInput \SetKwInputKwResultOutput \SetKwForForforendfor \SetCommentStyalgcommentfont
Heuristics for vehicle routing problems: Sequence or set optimization?
Túlio A. M. Toffolo, Thibaut Vidal, Tony Wauters
KU Leuven, Department of Computer Science, CODeS & imec  Belgium
Federal University of Ouro Preto, Department of Computing  Brazil
Pontifícia Universidade Católica do Rio de Janeiro, Computer Science Department  Brazil
Abstract. We investigate a structural decomposition for the capacitated vehicle routing problem (CVRP) based on vehicletocustomer “assignment” and visits “sequencing” decision variables. We show that an heuristic search focused on assignment decisions with a systematic optimal choice of sequences (using Concorde TSP solver) during each move evaluation is promising but requires a prohibitive computational effort. We therefore introduce an intermediate search space, based on the dynamic programming procedure of Balas & Simonetti, which finds a good compromise between intensification and computational efficiency. A variety of speedup techniques are proposed for a fast exploration: neighborhood reductions, dynamic move filters, memory structures, and concatenation techniques. Finally, a tunneling strategy is designed to reshape the search space as the algorithm progresses.
The combination of these techniques within a classical local search, as well as in the unified hybrid genetic search (UHGS) leads to significant improvements of solution accuracy. New best solutions are found for surprisingly small instances with as few as 256 customers. These solutions had not been attained up to now with classic neighborhoods. Overall, this research permits to better evaluate the respective impact of sequence and assignment optimization, proposes new ways of combining the optimization of these two decision sets, and opens promising research perspectices for the CVRP and its variants.
Keywords. Decisionset decompositions, Metaheuristics, Dynamic programming, Integer programming, Large neighborhood search, Vehicle routing problem
1 Introduction
The Capacitated Vehicle Routing Problem (CVRP) is classically described as a combination of a Traveling Salesman Problem (TSP) with an additional capacity constraint which lends a Bin Packing (BP) substructure to the problem (Toth and Vigo 2014). It can be seen as a Set Packing (SP) problem in which the cost of each set corresponds to the distance of the associated optimal TSP tour (Balinski and Quandt 1964). These problem representations emphasize the two decision sets at play: customertovehicle Assignments, and Sequencing choices for each route (Vidal et al. 2013b), a duality that has left longstanding impressions in the literature, from the early developments of routefirst clustersecond (Bodin and Berman 1979, Beasley 1983) and clusterfirst routesecond constructive methods (Fisher and Jaikumar 1981), all the way to the setcoveringbased exact methods and matheuristics which are currently gaining popularity.
Examining the recent progress on metaheuristics for the CVRP, little has changed in recent years concerning intraroute neighborhood search: Relocate, Swap and 2opt neighborhoods and their immediate generalizations are employed, and these neighborhoods alone are sufficient to guarantee that most solutions resulting from a local search contain TSP–optimal routes. This is generally because classical CVRP instances involve short routes with up to 15 or 20 visits. For such small problems, even simple neighborhood search methods for the TSP tend to produce optimal tours.
Based on this observation, a larger effort dedicated to TSP tour optimization, as a standalone neighborhood, is unlikely to result in further improvements. For this reason, it is very uncommon to observe the use of larger intraroute neighborhoods (e.g., 3Opt or beyond) in recent stateoftheart metaheuristics. Nevertheless, does this mean that Sequencing optimization should be abandoned in favor of more extensive search concerning Assignment choices? Certainly not. Indeed, even if local minima exhibit optimal TSP tours, interroute moves frequently lead to TSPsuboptimal tours which are rejected due to their higher cost, but would be accepted otherwise if the tours were optimized. Such solution improvements would then not arise from separate Assignment or Sequencing optimizations, but from a careful combination of both.
Figure 1 schematically represents the solution set of the CVRP, whose decision variables are split into Sequencing decisions (axis) and Assignment decisions (axis). The axis also represents the solutions in terms of their Assignment decisions solely, ignoring Sequencing choices. These partial solutions can be viewed as a projection (Geoffrion 1970) of the original solutions on the space defined by a single decision subset (Assignment). Moreover, from a solution represented in terms of Assignment decisions, it is possible to find the best associated complete solution by solving each TSP associated with the routes.
With this picture in mind, it is tempting to conduct a search in the space rather than . After all, the size of is exponentially smaller than that of , the solutions of are in average of better quality, and the average size of a path in is smaller, such that fewer LS iterations are expected for convergence. However, the obvious drawback is that each move evaluation in requires solving one or several small TSPs to optimality, leading to a significant computational effort. Still, note that considerable progress has been made in the past 30 years with regard to the efficient solution of TSPs, and small problems with approximately 20 customers are solvable in a few milliseconds. Based on these observations, this article takes a fresh look at heuristic searches for the CVRP to answer two essential questions about the search space :
 [nosep]

Is it practical and worthwhile to search in the space rather than ?

If searching in requires an excessive effort, can we define a search space which maintains most of the key properties of but can be more efficiently explored ?
As will be demonstrated in Section 4, our experiments led us to answer the first question negatively: even with nontrivial memory and speedup techniques (hashtables and move filters) the computational overhead related to the exact solution of TSPs during each move evaluation, for a complete search in , does not appear worth the gain in terms of solution quality.
By contrast, our answer to the second question is positive. Rather than requiring a complete exact solution of each TSP, the dynamic programming approach of Balas and Simonetti (2001), hereafter referred to as B&S, can be employed to perform a restricted route optimization during move evaluations. Given a range parameter and an initial tour, the B&S algorithm finds, in operations, the vertex sequence with minimum cost such that no vertex is displaced by more than positions. This allows us to define a search space such that and . Moreover, even for a fixed , we propose tunneling techniques that exploit the memory of past solutions to dynamically reshape the search space, in such a way that converges towards as the search progresses.
To evaluate experimentally the potential of the new search spaces, we conduct experiments with a simple multistart local search (MSLS), and with the unified hybrid genetic search (UHGS) of Vidal et al. (2012, 2014a). The use of for appears to lead to solutions of higher quality on the new instances from Uchoa et al. (2017). New best solutions were also found for surprisingly small instances with as few as 242 or 256 customers. These solutions had not been attained up to now with classic neighborhoods. Overall, this research allows to better evaluate the respective impact of Sequencing and Assignment optimization, proposing new ways to combine the optimization of these two decision sets, and leading to new stateoftheart algorithms for the CVRP.
2 Related Literature
This section reviews some key milestones concerning the management of Sequencing and Assignment decisions in vehicle routing heuristics, as well as decisionset decompositions.
Sequencing and Assignment decisions were first optimized separately in early constructive heuristics, giving rise to different families of methods. Routefirst clustersecond algorithms (Bodin and Berman 1979, Beasley 1983) first produce a giant TSP tour, before subsequently assigning consecutive visits into separate trips to produce a complete solution. In clusterfirst routesecond methods (Fisher and Jaikumar 1981), a clustering algorithm is employed to group customer visits into clusters, followed by TSP optimizations. Finally, petal algorithms (Foster and Ryan 1976, Renaud et al. 1996) are based on an apriori generation of candidate routes (petals), followed by the solution of a set packing or covering problem.
In the development of local search and metaheuristic algorithms which ensued in the 1990s and thereafter, Assignment and Sequencing optimizations began to be better integrated. Classical moves such as Relocate, Swap, 2opt* and their close variants allow to optimize both decision subsets. The associated neighborhood search methods form the basis of the vast majority of stateoftheart algorithms. Petal algorithms have withstood the test of time, and highquality routes are now extracted from local minima of metaheuristics instead of being enumerated in advance (see, e.g., Muter et al. 2010, Subramanian et al. 2013).
The variety of vehicle routing problem variants has also triggered studies concerning problem decompositions. Vidal et al. (2013b) established a review of the classical variants and their associated constraints, objectives, and decision sets, called attributes. The attributes were classified in relation to their impact on Sequencing, Assignment decisions, and Route evaluations in heuristics, leading to a structural problem decomposition which serves as a basis for the Unified Hybrid Genetic Search (UHGS) algorithm of Vidal et al. (2014a) and allows to produce stateoftheart results for dozens of VRP variants. Many problem attributes come jointly with new decision subsets, e.g., when optimizing vehicle routing with packing, timing or scheduling constraints (Goel and Vidal 2014, Pollaris et al. 2015, Vidal et al. 2015b), visit choices (Vidal et al. 2014b) or servicemode choices (Vidal et al. 2015a, Vidal 2017).
Decisionset decompositions are employed throughout many of the aforementioned papers to perform a search in the space of Sequencing and Assignment choices and optimally determine the remaining decision variables during each route and move evaluation. In this paper, the decisionset decomposition does not result from supplementary problem attributes, but is instead used to define exponentialsize polynomiallysearchable neighborhoods and transform the search space. Exponentialsize neighborhoods have a long history in the combinatorial optimization literature (Deineko and Woeginger 2000, Ahuja et al. 2002, Bompadre 2012). Most of these neighborhoods are based on shortest path or matching subproblems, as well as specific graph and distance matrix structures with which some NPhard problems become tractable (consider, as examples, Halin graphs or Monge matrices). As a rule of thumb, larger neighborhoods and faster search procedures are generally desirable. There are, however, theoretical limitations to the size of polynomiallysearchable neighborhoods. Gutin and Yeo (2003) proved that, for the TSP, no neighborhood of cardinality at least for a given constant can be searched unless NP P/poly.
The neighborhood of Balas and Simonetti (2001) is an exponentialsize neighborhood for the TSP. Given an incumbent tour represented as a permutation and a value , it contains all permutations such that fulfills and for all such that . In other words, if precedes by more than positions in , then precedes . Setting allows to fix the origin location (e.g., depot). This neighborhood contains solutions, and can be explored in using dynamic programming. This is a linear time complexity when is constant, and a polynomial complexity when . Balas and Simonetti (2001) performed extensive experiments, and demonstrated that this dynamic programming procedure can be used as a standalone neighborhood to improve high quality local minima of the TSP and its immediate variants. Later, Irnich (2008), Gschwind and Drexl (2016), and Hintsch and Irnich (2017) employed this neighborhood to solve arcrouting problems with possible cluster constraints, and dialaride problems. One common characteristic of these studies is that they employed B&S as a standalone neighborhood for route improvements. Only in one conference presentation (Irnich 2013), the possibility of using the B&S neighborhood in combination with some classical CVRP moves has been highlighted, but the performance of such an approach remains largely unexplored.
We seek to go one step further. Rather than applying this tour optimization procedure as a standalone optimization technique or in combination with a single classical neighborhood, we investigate its systematic use in combination with every move of a classical CVRP local search. As discussed in the following, the methodological implications of such a redefinition of the search space are noteworthy.
3 Proposed Methodology
We will describe the methodology as a local search on indirect solution representations, using a decoder. This algorithm can be readily extended into a wide range of vehicle routing metaheuristics, e.g., tabu search, iterated local search, or hybrid genetic algorithm (Gendreau and Potvin 2010, Laporte et al. 2014). There is no widely accepted term, in the current heuristic literature, for referring to the elements which represent such indirect solutions. The evolutionary literature usually refers to a genotype to denote solution encodings (and phenotype for the solutions themselves), whereas the localsearch based metaheuristic literature refers to incomplete or indirect solutions (which are converted into complete solutions via a decoder function). Incomplete lets us think that the representation is necessarily a subset of a complete solution, and unnecessarily restricts the application scope. To circumvent this issue, we henceforth employ the term primitive solutions. We first recall some basic definitions related to neighborhood search and indirect solution representations, then proceed with an analysis of alternative search spaces and the description of the proposed local search algorithm.
Definition 1 (Primitive solutions and search space).
We consider a combinatorial optimization problem of the form , where is the solution space, and is an objective function to minimize. Let be the set of primitive solutions, and let the decoder be an injective application that transforms any into a complete solution . A neighborhood is defined as a mapping that associates with each primitive solution a set of neighbors . The graph induced by and is referred to as the search space.
Definition 2 (TSP–optimal tour).
A tour is TSP–optimal if there exists no other permutation of its visits such that with a shorter total distance.
Definition 3 (B–optimal tour).
A tour is B–optimal if there exists no other permutation of its visits with a shorter total distance such that and for all with . The parameter is the range of the B&S neighborhood.
3.1 A Choice of Search Space
In this section, we examine the search spaces associated with the set of all solutions (), of those with TSP–optimal tours (), and of those with B–optimal tours () and discuss their relative merits.
Search Space . Classical local search methods for the CVRP do not distinguish between primitive and complete solutions. In the usual search space , solution sets and are equal and the decoder is the identity function. The neighborhood is based on the definition of one or several classes of moves. A move is a local modification that can be applied to a primitive solution to generate a neighbor . For each search space considered in this paper, we will eventually refer to several classes of moves, but to a single neighborhood only, which corresponds to the union of all primitive solutions attainable from via one single move. Classical moves for the CVRP are based on relocations and exchanges of a bounded number of vertices, or replacements of a bounded number of edges. Most common neighborhoods have a quadratic cardinality ( for all ). We refer to Vidal et al. (2013b) for a comprehensive survey on classical local searches for the CVRP.
Figure 2 represents the search space associated with Relocate moves only, for a small asymmetric CVRP instance with three customers.
There are 13 possible solutions for this problem, each represented by a set of ordered customer visits. Solution ‘[1,2,3]’, for example, employs one vehicle to visit customers 1, 2 and 3, while solution ‘[1][2][3]’ employs three vehicles, one per customer.
Each solution is represented by a node, positioned on the axis according to its quality (the more to the right, the better a solution is). The set of outgoing arcs of each solution points towards its neighbors.
Moreover, solutions with identical customerroute assignments are grouped within dashed areas.
Note that, for this instance size, it is always possible to reach the optimum solution from any starting point in two successive moves.
In a local search that explores the neighborhood in random order and applies an improving move as soon as it is found, the worst case corresponds to five moves (when the initial solution is ‘[2][1,3]’ or ‘[1][2][3]’).
Search Space . As discussed earlier in this work, CVRP solutions can also be represented in terms of their Assignment decisions, excluding the Sequencing decisions in the representation and delegating the choices of the best visit sequences to the decoder. With such a paradigm, one can define a local search in the space of primitive solutions, where each represents a partition of the customer set into subsets whose sums of demands do not exceed the vehicle capacity. The decoder is based on an exact TSP solver, responsible for generating the best visit sequence originating and finishing at the depot for each subset of customers. In this sense, the image contains exclusively solutions with TSP–optimal tours.
The neighborhood used to explore the search space can remain similar to classical CVRP neighborhoods, based on relocations or exchanges of customers between subsets, or involve other families of moves specialized for partition problems. Figure 3 represents the resulting search space with simple Relocate moves. Only TSP–optimal tours are explored and therefore the size of the search space reduces down to six primitive solutions. The other solutions and their connections are represented in light gray. Note, in our small example, that now at most three successive improving moves may be applied to attain the optimum from ‘[1][2][3]’.
Search space is smaller than , and our computational experiments (Section 4) demonstrate that a search in this space indeed leads to solutions with higher quality.
However, each move evaluation in this space requires executing an algorithm with exponential worstcase time complexity, a TSP solver, in order to decode each primitive solution in the neighborhood for cost evaluation.
Although research on the TSP has culminated in very efficient algorithms over the past thirty years, thousands (or millions) of small TSP instances should be solved during a local search in , and thus the total computational effort dedicated towards decoding can grow prohibitively large.
Moreover, bad behavior in a single case (e.g., due to an unusual long route with many customers or bad branching decisions) can be sufficient, without any other safeguard, to stall the entire algorithm.
Search Space . To circumvent the aforementioned issues, we study an alternative search space in which the set of primitive solutions is a subset of the complete solutions (with their Assignment and Sequencing decisions), but where the decoder is nontrivial, and consists of applying B&S multiple times to each route with a fixed range ( value) until the tours become B–optimal. With these assumptions, the image contains exclusively complete solutions with B–optimal tours. As such, the application of B&S can be viewed as a postoptimization step during classical CVRP move evaluations, opening the way for additional solution improvements. A careful analysis of the resulting search space gives even more significance to this approach, due to three properties:
Property 1.
From an initial solution containing a B–optimal tour, a local search in the space explores only B–optimal tours.
Property 2.
For a fixed range , each move evaluation and subsequent solution decoding is done in polynomial time as a function of and the number of applications of B&S.
Property 3.
The search space is such that and , with being the number of customers.
These three properties are all fundamental for the methodology that follows. Property 1 demonstrates how space contains fewer solutions than , and that the overall quality of these solutions tends to be higher (since nonBoptimal tours are filtered out). Moreover, Property 2 gives some computational time guarantees: even if the computational effort grows quickly with the range , the effort of the decoder is guaranteed to remain stable for all when is constant, eliminating the possibility of a computational effort peak for specific TSP instances. Finally, Property 3 demonstrates how balances the effort dedicated towards the optimization of the Assignment and Sequencing decision sets, and establishes as an intermediate search space generalizing and .
Figure 4 illustrates the search space for the same example as previous figures with . It is an intermediate between the spaces depicted by Figures 2 and 3, which correspond to and , respectively. Note how the solution is now a neighbor of and , as highlighted by the two dotted arcs.
Discussions and Choice. In light of these observations, we have conducted computational experiments on search space , using the dynamic programming algorithm of Balas and Simonetti (2001) to decode each solution, as well as on search space using the TSP solver Concorde (Applegate et al. 2006). Despite several speedup techniques (Section 3.2), the search in space remained inefficient throughout our current experiments, especially for instances with a large number of customers per route. We therefore decided to focus on search space , and devised several speedup techniques to enable its efficient exploration.
3.2 Efficient Local Search
To efficiently explore space , we developed a local search algorithm which exploits static neighborhood reductions, dynamic move filters, efficient memory structures and concatenation techniques. Most of these techniques seek to limit the search effort in . This resulting method, displayed in Algorithm 3.2, can be easily integrated in stateoftheart metaheuristics for vehicle routing problems.
[p] \KwDataAn initial complete solution , an evaluation threshold and a granularity threshold
is a local minimum \tcpEnumerating moves  candidate lists based on vertex proximity \Foreach move involving a vertex pair ,
The move modifies up to two routes of . Let be the sum of the costs of these two routes, and let and be the new routes in .
First, filter infeasible moves with respect to capacity constraints in :
or \Continue.
Second, consider the cost of the classical CVRP move to filter nonpromising solutions in :
\Continue.
Third, decode the routes and to evaluate the move in :
each route with \tcpCompute hash key in and check memory in :
If not in memory, repeatedly apply the dynamic programming algorithm of Balas and Simonetti (2001) until the route becomes : \If
Filter nonimproving moves:
\Continue.
At this stage, apply since it is an improving move in : Set ; Replace the routes (,) by (,) in \Return
Neighborhood Reductions. First of all, as in the majority of recent local search based metaheuristics for the CVRP, the neighborhood of each incumbent solution is limited to moves that involve close vertices (Algorithm 1, Line 3). In particular, we use the classical intraroute and interroute Relocate and Swap moves, for single vertices or generalized to pairs of consecutive vertices, as well as the 2Opt and 2Opt* moves. The resulting neighborhood contains a quadratic number of moves. As detailed in Vidal et al. (2013a), and in a similar way as Johnson and McGeoch (1997) and Toth and Vigo (2003), the search can be restricted to a subset of these moves that reconnect at least one vertex with a vertex belonging to the closest vertices of . The neighborhood size becomes , enabling a significant speedup for largescale problem instances.
Dynamic move filters. To further restrict the search to promising moves, each move ’s feasibility is evaluated in , in terms of capacity constraints, being discarded if it leads to an infeasible solution. The total cost of the solution generated by prior to its optimization by the B&S decoder is evaluated subsequently. This cost represents an upper bound for the final cost of the move in after the application of the decoder. The move evaluation is pursued only if the solution cost has increased by a factor or less due to its application, that is, only if Condition (1) is satisfied. Otherwise, the move is discarded.
(1) 
Parameter plays an important role in defining how many moves are evaluated. The higher the value of , the less pruning is induced by Equation (1). Contrastingly, when only immediately improving neighbors are evaluated. Defining a good value for is nontrivial, given that it is an instancedependent parameter. Since a fixed value would not suit instances with different sizes and characteristics, we suggest to use an adaptive parameter. The principle consists in adjusting to ensure a target range for the fraction of filtered moves. After each 1,000 move evaluations, the fraction of filtered moves is collected and whenever it falls outside of the desired range, is updated. If this fraction is too large, then is increased by a multiplicative factor . Conversely, if is insufficient, then the parameter is decreased:
(2) 
Global memory. The B&S algorithm, used as a decoder, requires a computational effort which grows linearly with the route size and exponentially with parameter . It is thus essential to restrain the use of this procedure to a strict minimum and avoid decoding twice the same route over the course of the search. To that end, we rely on a global memory to store the routes that have been decoded, avoiding recalculations (Lines 12–14). We use a hashtable for this task, since it allows queries given the hash key associated with a route.
Two important aspects should be discussed. First, since the available memory space is finite, some strategy is necessary to limit the memory size in case of an excessive space consumption. To that end, we define an upper bound on the number of routes stored in the memory and eliminate half of the entries, those used with less frequency, whenever this limit is attained.
The second aspect to be discussed concerns the effort spent querying the memory. It is well known that local searches for the CVRP evaluate millions of moves, and that constanttime move evaluations are essential for a good performance. Capacity checks (Line 5) and simple distance computations (Line 7) can be easily achieved in using incremental move evaluations or concatenation techniques (Vidal et al. 2014a). Moreover, querying the memory for a given route supposes the availability of a hash index which characterizes the associated sequence of visits, but a direct approach that sweeps through the route to compute this index already takes time. To avoid this bottleneck, we employ specific hash functions and calculation techniques based, again, on concatenations in . These concepts are discussed in the following section.
3.3 ConstantTime Evaluations
We use the concatenation strategy of Vidal et al. (2014a, 2015b) to perform efficient cost and loadfeasibility evaluations. This strategy exploits the fact that any route obtained from a classical move on an incumbent solution corresponds to a recombination of a bounded number of (customer and depot) visit sequences of . As such, the new routes can be expressed as a concatenation of sequences . We also extend this approach to enable computations of hash keys.
To efficiently evaluate the cost, load, and compute the hash keys, we perform a preliminary preprocessing on the subsequences of consecutive visits which compose the solution . Four quantities are calculated: the total demand of a sequence , its distance , and its hash keys and . For a sequence containing a single visit with demand , , , and , where is a prime number. Moreover, Equations (3–5) extend these quantities, by induction, for any sequence of visits expressed as the concatenation of two sequences and . In these equations, expresses the distance between visits and .
(3)  
(4)  
(5)  
(6) 
As in Vidal et al. (2014a), Equations (3–6) are first employed iteratively, in lexicographic order, to obtain information concerning all sequences during the preprocessing phase. Afterwards, the same equations are used for move evaluations. Since any route obtained from a classical move corresponds to the concatenation of a bounded number of sequences, it is possible to obtain the associated load, distance, and hash keys by applying these equations a limited number of times. Then, the information on subsequences is updated every time an improving move is applied, a rare occurrence in comparison to the number of moves evaluated.
The two hash functions defined in Equations (5–6) are employed together as a means of reducing chances of two distinct sequences having identical hashes. The function is a multiplicative hash which depends on the visit permutation (Knuth 1973). Note that, when implementing such a function, the values must be precomputed and bounded (taking the rest of the integer division by a large number) to prevent overflow during multiplication. The second function is an additive hash which only depends on the set of visited customers, and not on the visit sequence. These hash functions are easily recognized when reformulated as follows:
(7)  
(8) 
These functions fit well our purposes due to their inductive definition based on the concatenation operation. They are employed, along with the route distance and its number of visits, to verify a correct match in the memory in without a complete route comparison in . To minimize the risk of two routes having identical hashes, we duplicated these hash functions with different values for , leading to four hash values overall. In the first case, is set to the smallest prime number greater than the number of customers. In the second case, is set to (multiplier of Kernighan and Ritchie 1988). Despite this strategy, a tiny chance of false positives remains. However, no false positive was registered within our computational experiments considering multiple runs on 100 different instances.
3.4 Reshaping the search space – Tunneling strategy
Until now, the purpose of the global memory has been focused on saving computational effort by avoiding duplicate calls to the B&S algorithm. Yet, as shown in the following, this structure can be exploited to a larger extent to promote the discovery of good solutions.
Consider two input routes and indexed in the order of their appearance in the local search, and representing the same set of customer visits. The route is first saved in the global memory along with its associated B–optimal route of cost . Subsequently, the LS considers without finding a match in the memory, triggering a new execution of the B&S algorithm and leading to a cost . if , then the algorithm has failed to recognize that a better TSP tour has been found in prior search for the same customers.
To improve the behavior of the algorithm in such situations, we introduce a guidance mechanism called tunneling. Guidance techniques are a set of strategies which analyze and exploit the search history to direct the search towards promising or unexplored regions of the search space (Crainic and Toulouse 2008). Our method works as follows: every route issued from a nonfiltered LS move and absent from the memory is decoded by the B&S algorithm; yet, rather than directly returning the output of B&S, the algorithm finds and returns the best known permutation of the visits for this customer set, found over previous B&S executions. The goal of this strategy is to intensify the search around known highquality tours without jeopardizing the discovery of better route configurations. It can be efficiently implemented with a refinement of the hashtablebased memory structure, by grouping the routes into different buckets according to their visit set, and using the additive hash function of Equation (6) for queries. Figure 5 summarizes this process.
This tunneling strategy has a significant impact on the search space. Initially, as the search starts, the algorithm explores the space . Then, as the search progresses, the memory starts to be filled, and the algorithm reintroduces more and more frequently the best known routes in its solutions. In a hypothetical situation where all feasible routes have already been memorized (hypothetical due to the needed exponential memory size), the TSP–optimal routes would be systematically returned, and the algorithm behaves as if it was searching in . With this limit case in mind, the tunneling strategy contributes to reshape the space into as the search progresses. Moreover, note that this strategy remains fully relevant for a classical neighborhood search, even without B&S decoder.
Consider the same example as depicted in Section 3.1 (Figures 2, 3 and 4). Suppose that the route has been identified in the past search along with its associated B–optimal tour and that the tunneling strategy is employed. Figure 6 illustrates the resulting search space: all solutions that include ’’ in their neighborhood now point towards solution ’’ instead, which has identical customertovehicle Assignments but a lower cost due to better Sequencing decisions. With only one route in memory, the resulting search space already becomes equivalent to .
4 Computational experiments
We conducted extensive computational analyses to measure the benefits of a search in and , the impact of the tunneling strategy and move evaluation filters. For these tests, we considered the simple local search (LS) described in Section 3, as well as a more advanced metaheuristic, the UHGS of Vidal et al. (2012, 2014a), which was adapted by replacing its native local search by the proposed method. This extension of the method will be referred as UHGSBS. Except this modification of the local search, all other parameters and procedures of UHGSBS remain the same as in the original article, and the termination criterion is set to consecutive iterations without improvement.
A single local search requires only a limited computational effort. Thus, the first analyses (Section 4.1) on the performance of the LS in and could be done for multiple values of the range parameter . It was also possible to evaluate the impact of without the interference of dynamic move filters and tunneling techniques. Since UHGSBS performs a more extensive search with multiple local search runs, our analyses with this method (Section 4.2) are focused on a smaller set of values for the range parameter , in the presence of dynamic move filters.
All algorithms were implemented in C++ and executed on a single thread of an Intel(R) Xeon(R) E52680v3 CPU. A RAM limit of 8GB was imposed for each run. To solve the TSP problems when considering the space , we used the Concorde solver (Applegate et al. 2006). To conduct the experiments with the UHGSBS, we used the code base made available at https://github.com/vidalthi/HGSCARP, from Vidal (2017).
4.1 Preliminary experiments with a simple local search
In a first experiment, we tested the local search of Section 3 on spaces and for , to observe the growth of its computational time as a function of the parameter , and identify a range of values over which the approach remains practical. As an initial solution, we used the result of the savings algorithm of Clarke and Wright (1964). We set in order to observe the results without the interference of move filters.
We considered the 100 recent benchmark instances of Uchoa et al. (2017), as these instances remain highly challenging for metaheuristics and cover a larger variety of instance size and characteristics: demand and customer distribution, depot location, and route length. For each instance, we ran the method 20 times with different random seeds. The results are summarized in Figure 7, in the form of boxplots. The leftmost graph represents the percentage gap in terms of solution quality, relative to that of the best known solution (BKS) collected from Uchoa et al. (2017): , where is the solution value of the method and is the BKS value. The rightmost graph represents the CPU time of the method, using a logarithmic scale.
As illustrated by these experiments, the search in space visibly leads to solutions of higher average quality, albeit in a CPU time largely greater than that of a classical local search in . The solution quality of a search in the space consistently increases as grows, along with the needed CPU time. When , the algorithm behaves as a classical local search in . When is large, the method becomes more similar to a search in . A difference of solution quality can still be noticed between and , due to some instances containing up to 25 deliveries per route.
Each value of the range parameter establishes a tradeoff between computational effort and solution quality. The computational time of the local search does not exceed two seconds when exploring with , but it culminates to 1000 seconds on the largest instances when exploring . Clearly, the additional CPU time required by Concorde is not compatible with the repeated application of local searches within stateoftheart metaheuristic algorithms, such that we should concentrate our efforts on the exploration of the space with moderate values of .
Moreover, a careful analysis of these results on subsets of instances with a different average number of customers per route gives additional insights. This analysis is reported in Figure 8, considering the 20 instances with smallest average route cardinality, in the range , and the 20 instances with largest average route cardinality, in .
As illustrated in Figure 8, the benefits of a search in or are small for instances with a small number of customer visits per route. In particular, all runs with lead to a similar solution quality and CPU time. This is due to the fact that B&S does an exact TSP optimization when is greater or equal to the route cardinality. In contrast, the benefits in terms of solution quality are larger on instances with a high number of customer visits per vehicle. We observe a significant improvement of the solutions when varies from to , from an average gap of down to . Subsequently, as increases beyond , the rate of improvement is smaller. Increasing up to the maximum route size would still be beneficial, but impracticable in terms of CPU time.
4.2 Experiments with UHGSBS – range parameter, move filters and tunneling
As viewed in the previous section, a local search in the space can lead to solutions of better quality than a search in , at the expense of a higher computational effort. Still, even if solution improvements were observed for simple local searches, it is an open question whether the inclusion of these extended search procedures into stateoftheart metaheuristics translates back into significant quality improvements.
This section now analyses the performance of the UHGSBS metaheuristic, considering in combination with dynamic move filters. The move filters are designed to eliminate a large fraction of complete move evaluations, such that a larger computational effort can be spent in the evaluation of each remaining move without significantly impacting the overall CPU time of the method.
After preliminary experiments, we observed that setting and establishes a good compromise between solution quality and computational effort. This configuration, without tunneling strategy, was used as baseline for the experiments of this section. We then modified each parameter and design choice, in turn, to investigate its impact. This experimentation was restricted to a subset of 40 instances with a number of customer visits , as these instances require limited CPU time and remain challenging for stateoftheart metaheuristics. For each instance, ten runs have been conducted with different seeds.
Impact of the range parameter. Figure 9 compares the performance, in terms of average percentage gap and CPU time, of the classical UHGS with that of its extended version searching in , for . As in previous figures, the results are presented in the form of boxplots.
As expected, the gaps obtained with the variants of UHGSBS are much smaller than those of simple local searches, due to the better exploration capabilities of the method. Average gaps range from when exploring the classical search space , to when exploring . As highlighted by pairwise Wilcoxon tests, statistically significant differences exist between the results of , and (pvalues ). A decreasing returns effect can also be observed; the difference of quality between the solutions obtained by and is larger than between by , which is turn is larger than between by , and so on.
The configuration , in particular, achieves a good tradeoff between quality and search effort. As demonstrated by the outliers with negative gap in the figure, this configuration has led to new best known solutions for surprisingly small instances with 256 and 294 customers, on which hundreds of test runs had been conducted in past. This is an indication that the search in has the potential to lead to structurally different solutions, which were not attained with more classical searches. The main challenge therefore is to harness this ability without sacrificing too much computational effort.
Impact of the dynamic move filter. Figure 10 investigates the impact of different target intervals (desired quantity of filtered moves – Section 3.2) for the dynamic move filters. The range parameter remains fixed to . It also indicates the results obtained when filtering all nonimproving moves (), which is equivalent to using B&S only as a postoptimization procedure, after the discovery of each improving move.
These experiments demonstrate that move filters have a large incidence on the solution quality. These filters, are however, essential to maintain a low computational effort. In particular, filtering all nonimproving moves prior to the evaluation of B&S () leads to an average gap of 0.26%, compared to 0.21% when setting as a target for the dynamic move filter and evaluating the nonfiltered moves in combination with B&S.
This validates an important hypothesis explored in this article: many moves that are usually discarded in regular local search methods and metaheuristics can lead to improved solutions when applied in combination with a route optimization procedure. Naturally, this capability goes along with an increased CPU time. Nonetheless, by an adequate calibration of the move filters, the total CPU time dedicated to B&S can be restricted enough to not intervene as a bottleneck. This is visible by the results of configuration , which on average used no more than twice the time of . For the remainder of these analyses, we selected this configuration, which establishes a good balance between the exploitation of the capabilities of B&S and the computational effort.
Impact of the tunneling strategy. Finally, Table 1 compares the performance of UHGSBS without and with the tunneling strategy. The range parameter has been set to , and . The leftmost group of columns reports the instance name, number of customers average number of visits per route in the BKS for each instance. Then, the next columns present, for each approach, the average solution value over ten runs, best solution value, average CPU time, percentage of routes which have been successfully queried from the global memory without a reevaluation, number of executions of the B&S optimization procedure, total number of iterations.
#  Instance  Without tunneling strategy  With tunneling strategy  

Time  Avg.  Best  Cache  B&S  Iters  Time  Avg.  Best  Cache  B&S  Iters  
21  Xn195k51  3.8  44289.5  44225  60.6%  1.4E+08  34577  5.2  44272.6  44225  60.2%  1.8E+08  41713 
22  Xn200k36  5.3  58618.5  58578  97.8%  1.1E+07  43664  6.2  58625.7  58578  97.8%  1.0E+07  41788 
23  Xn204k19  5.8  19569.0  19565  91.0%  2.7E+07  31315  6.6  19568.5  19565  92.6%  2.2E+07  31118 
24  Xn209k16  15.5  30673.1  30656  77.6%  1.1E+08  35319  17.8  30676.8  30656  81.0%  9.1E+07  36174 
25  Xn214k11  34.8  10873.2  10856  70.5%  2.1E+08  48054  32.3  10872.0  10856  76.2%  1.5E+08  46133 
26  Xn219k73  12.9  117602.7  117595  3.7%  5.0E+08  30949  11.0  117603.8  117595  3.6%  5.0E+08  30773 
27  Xn223k34  6.2  40490.3  40437  95.4%  2.3E+07  38698  7.7  40489.0  40437  95.7%  2.1E+07  39268 
28  Xn228k23  14.5  25784.1  25743  87.9%  8.5E+07  48303  14.6  25779.5  25743  88.6%  6.9E+07  42435 
29  Xn233k16  10.9  19305.3  19230  87.1%  4.7E+07  31274  9.5  19309.0  19230  90.4%  3.2E+07  29652 
30  Xn237k14  16.6  27053.8  27050  75.9%  9.7E+07  26927  17.5  27047.2  27042  81.1%  8.1E+07  30294 
31  Xn242k48  8.7  82922.0  82792  94.3%  4.4E+07  49734  12.9  82919.6  82792  93.6%  6.6E+07  62591 
32  Xn247k47  14.6  37380.7  37278  49.2%  6.1E+08  63084  13.7  37393.1  37281  49.1%  5.1E+08  52872 
33  Xn251k28  12.7  38767.4  38699  90.8%  6.4E+07  48107  11.4  38793.5  38699  91.2%  4.5E+07  35265 
34  Xn256k16  8.8  18880.0  18880  85.9%  4.3E+07  22991  9.4  18875.9  18839  89.6%  3.3E+07  24930 
35  Xn261k13  43.0  26628.7  26579  71.9%  2.3E+08  46952  38.9  26620.5  26586  75.1%  1.7E+08  39770 
36  Xn266k58  14.3  75703.6  75558  98.9%  1.5E+07  69013  17.8  75737.9  75646  98.9%  1.5E+07  72435 
37  Xn270k35  6.6  35310.5  35303  96.8%  1.5E+07  34796  8.5  35318.5  35303  97.1%  1.4E+07  35654 
38  Xn275k28  11.7  21256.9  21245  89.0%  5.9E+07  36470  12.8  21255.0  21245  91.0%  4.4E+07  34488 
39  Xn280k17  69.1  33614.0  33593  66.7%  4.1E+08  53528  86.3  33600.8  33584  70.1%  4.3E+08  63405 
40  Xn284k15  82.2  20306.0  20255  70.0%  4.3E+08  73605  67.4  20283.3  20238  73.9%  3.0E+08  61772 
41  Xn289k60  18.3  95501.4  95395  69.6%  4.5E+08  75465  20.7  95475.6  95245  69.5%  4.4E+08  72940 
42  Xn294k50  9.3  47272.9  47240  96.8%  2.6E+07  47162  11.6  47289.7  47239  96.3%  3.3E+07  53125 
43  Xn298k31  7.1  34282.8  34231  94.6%  2.3E+07  27625  8.1  34280.8  34231  95.2%  2.0E+07  26203 
44  Xn303k21  22.3  21839.1  21744  86.7%  1.0E+08  43974  21.1  21833.3  21744  88.9%  7.9E+07  39415 
45  Xn308k13  73.6  25893.2  25864  68.5%  3.3E+08  48194  63.6  25911.8  25866  71.8%  2.6E+08  42804 
46  Xn313k71  15.8  94270.5  94169  57.8%  4.7E+08  54350  16.8  94289.7  94192  57.6%  4.9E+08  58766 
47  Xn317k53  29.7  78389.3  78372  97.5%  3.3E+07  60402  31.6  78405.0  78380  97.7%  3.2E+07  65092 
48  Xn322k28  10.9  29918.8  29880  93.3%  3.7E+07  36519  17.4  29894.6  29834  93.9%  4.3E+07  48316 
49  Xn327k20  31.2  27618.2  27560  81.3%  1.5E+08  44348  35.5  27591.7  27560  83.4%  1.5E+08  47634 
50  Xn331k15  71.6  31138.6  31103  65.2%  3.7E+08  47298  63.2  31128.3  31103  67.3%  3.0E+08  41354 
51  Xn336k84  34.3  139522.2  139303  39.0%  1.4E+09  86713  40.9  139572.7  139339  39.7%  1.4E+09  89964 
52  Xn344k43  11.7  42161.9  42086  97.1%  2.3E+07  42768  22.7  42136.7  42066  97.4%  3.0E+07  65989 
53  Xn351k40  26.1  26010.1  25958  94.6%  7.9E+07  72726  32.8  25991.2  25947  94.7%  8.3E+07  76997 
54  Xn359k29  65.3  51676.6  51608  84.9%  3.1E+08  84963  56.1  51684.5  51569  85.9%  2.4E+08  67011 
55  Xn367k17  92.0  22870.1  22814  72.2%  4.0E+08  60550  66.8  22894.6  22814  71.1%  2.8E+08  39572 
56  Xn376k94  61.3  147736.8  147718  99.2%  1.6E+07  64604  53.1  147738.1  147718  99.2%  1.4E+07  61378 
57  Xn384k52  26.0  66254.7  66095  97.7%  4.2E+07  77380  29.1  66184.1  66133  97.6%  3.8E+07  68679 
58  Xn393k38  22.7  38361.0  38269  92.7%  7.2E+07  47398  26.9  38305.5  38260  93.7%  6.5E+07  48965 
59  Xn401k29  79.6  66357.8  66269  83.9%  3.4E+08  86787  69.7  66347.6  66225  86.1%  2.7E+08  78713 
60  Xn411k19  110.9  19753.4  19725  72.8%  4.9E+08  65729  74.8  19748.3  19719  77.9%  2.9E+08  51034 
Average values:  30.4  0.17%  0.02%  80.2%  2.1E+08  51058  29.3  0.15%  0.01%  81.5%  1.8E+08  49912 
The tunneling strategy leads to an average Gap(%) of , compared to without tunneling. Although this is only a difference of , reducing the gap indeed becomes harder as the solution quality approaches known BKS and optimal solutions. This improvement in solution quality also comes with a reduction of the overall CPU time, as the tunneling stimulates a faster convergence towards local minima and allows a better management of the global memory, by storing at most one route per customer set. As a consequence, the chances of successful queries in the memory is sensibly higher (81.5% compared to 80.2% on average), thereby reducing on average the number of B&S executions as well as the total number of iterations. Based on these observations, tunneling is beneficial without any other visible counterpart. We will use this mechanism for the final tests on the complete set of benchmark instances, in the next section.
4.3 Comparison with recent stateoftheart algorithms
Finally, this section reports detailed results of UHGSBS, using the baseline configuration and the tunneling strategy, on the complete set of 100 instances proposed by Uchoa et al. (2017). The results of UHGSBS are compared to that of the current stateoftheart algorithms: the Hybrid Iterated Local Search (HILS) proposed by Subramanian et al. (2013), and the original UHGS of Vidal et al. (2012, 2014a), which were executed 50 times for each instance. To keep the total computational effort within reasonable limits, UHGSBS was executed 10 times for each instance. The maximum number of consecutive iterations without improvements was set to to evaluate UHGSBS in the same conditions as UHGS on this set of instances (Uchoa et al. 2017). A hard runtime limit of 24 hours was imposed for each run.
Tables 2 and 3 report the results obtained with UHGSBS, in comparison with UHGS and HILS. The columns present the average CPU time in minutes, the average solution value and the best solution value for all approaches. The best results are highlighted in boldface in the table, with indicating an improvement over the best known solution from http://vrp.atdlab.inf.pucrio.br/, and the average gap to the best solution produced by the considered methods (including those generated during this research) is presented for each algorithm in each table’s final row.
#  Instance  ILS  UHGS  UHGSBS  
Time  Average  Best  Time  Average  Best  Time  Average  Best  
1  Xn101k25  0.1  27591.0  27591  1.4  27591.0  27591  2.4  27591.0  27591  
2  Xn106k14  2.0  26375.9  26362  4.0  26381.8  26378  17.6  26374.3  26362  
3  Xn110k13  0.2  14971.0  14971  1.6  14971.0  14971  3.9  14971.0  14971  
4  Xn115k10  0.2  12747.0  12747  1.8  12747.0  12747  6.4  12747.0  12747  
5  Xn120k6  1.7  13337.6  13332  2.3  13332.0  13332  38.3  13332.0  13332  
6  Xn125k30  1.4  55673.8  55539  2.7  55542.1  55539  6.1  55540.0  55539  
7  Xn129k18  1.9  28998.0  28948  2.7  28948.5  28940  8.4  28940.0  28940  
8  Xn134k13  2.1  10947.4  10916  3.3  10934.9  10916  20.2  10916.0  10916  
9  Xn139k10  1.6  13603.1  13590  2.3  13590.0  13590  8.9  13590.0  13590  
10  Xn143k7  1.6  15745.2  15726  3.1  15700.2  15700  33.1  15700.0  15700  
11  Xn148k46  0.8  43452.1  43448  3.2  43448.0  43448  5.1  43448.0  43448  
12  Xn153k22  0.5  21400.0  21340  5.5  21226.3  21220  15.2  21225.6  21225  
13  Xn157k13  0.8  16876.0  16876  3.2  16876.0  16876  28.4  16876.0  16876  
14  Xn162k11  0.5  14160.1  14138  3.3  14141.3  14138  12.2  14138.0  14138  
15  Xn167k10  0.9  20608.7  20562  3.7  20563.2  20557  36.1  20557.0  20557  
16  Xn172k51  0.6  45616.1  45607  3.8  45607.0  45607  5.7  45607.0  45607  
17  Xn176k26  1.1  48249.8  48140  7.6  47957.2  47812  16.2  47830.7  47812  
18  Xn181k23  1.6  25571.5  25569  6.3  25591.1  25569  13.9  25569.4  25569  
19  Xn186k15  1.7  24186.0  24145  5.9  24147.2  24145  20.2  24145.0  24145  
20  Xn190k8  2.1  17143.1  17085  12.1  16987.9  16980  161.9  16985.3  16980  
21  Xn195k51  0.9  44234.3  44225  6.1  44244.1  44225  9.3  44283.8  44225  
22  Xn200k36  7.5  58697.2  58626  8.0  58626.4  58578  12.0  58615.1  58578  
23  Xn204k19  1.1  19625.2  19570  5.4  19571.5  19565  15.7  19567.0  19565  
24  Xn209k16  3.8  30765.4  30667  8.6  30680.4  30656  35.7  30671.3  30656  
25  Xn214k11  2.3  11126.9  10985  10.2  10877.4  10856  52.3  10872.1  10856  
26  Xn219k73  0.9  117595.0  117595  7.7  117604.9  117595  18.2  117600.5  117595  
27  Xn223k34  8.5  40533.5  40471  8.3  40499.0  40437  18.6  40478.4  40437  
28  Xn228k23  2.4  25795.8  25743  9.8  25779.3  25742  29.0  25768.0  25743  
29  Xn233k16  3.0  19336.7  19266  6.8  19288.4  19230  34.0  19276.5  19230  
30  Xn237k14  3.5  27078.8  27042  8.9  27067.3  27042  32.2  27048.8  27042  
31  Xn242k48  17.8  82874.2  82774  12.4  82948.7  82804  18.1  82920.9  82751  
32  Xn247k47  2.1  37507.2  37289  20.4  37284.4  37274  27.7  37388.9  37274  
33  Xn251k28  10.8  38840.0  38727  11.7  38796.4  38699  20.2  38778.7  38684  
34  Xn256k16  2.0  18883.9  18880  6.5  18880.0  18880  23.0  18867.7  18839  
35  Xn261k13  6.7  26869.0  26706  12.7  26629.6  26558  48.6  26618.1  26558  
36  Xn266k58  10.0  75563.3  75478  21.4  75759.3  75517  29.9  75710.7  75478  
37  Xn270k35  9.1  35363.4  35324  11.3  35367.2  35303  18.9  35314.6  35303  
38  Xn275k28  3.6  21256.0  21245  12.0  21280.6  21245  22.7  21255.0  21245  
39  Xn280k17  9.6  33769.4  33624  19.1  33605.8  33505  136.2  33587.9  33503  
40  Xn284k15  8.6  20448.5  20295  19.9  20286.4  20227  97.7  20282.1  20228  
41  Xn289k60  16.1  95450.6  95315  21.3  95469.5  95244  41.7  95447.2  95211  
42  Xn294k50  12.4  47254.7  47190  14.7  47259.0  47171  27.0  47272.7  47161  
43  Xn298k31  6.9  34356.0  34239  10.9  34292.1  34231  20.7  34276.3  34231  
44  Xn303k21  14.2  21895.8  21812  17.3  21850.9  21748  48.4  21811.2  21744  
45  Xn308k13  9.5  26101.1  25901  15.3  25895.4  25859  112.8  25897.3  25861  
46  Xn313k71  17.5  94297.3  94192  22.4  94265.2  94093  30.6  94280.4  94045  
47  Xn317k53  8.6  78356.0  78355  22.4  78387.8  78355  50.3  78385.3  78355  
48  Xn322k28  14.7  29991.3  29877  15.2  29956.1  29870  27.7  29892.5  29834  
49  Xn327k20  19.1  27812.4  27599  18.2  27628.2  27564  68.7  27590.8  27532  
50  Xn331k15  15.7  31235.5  31105  24.4  31159.6  31103  102.1  31126.7  31103  
Average gap:  0.37%  0.13%  0.14%  0.02%  0.10%  0.00% 
#  Instance  ILS  UHGS  UHGSBS  
Time  Average  Best  Time  Average  Best  Time  Average  Best  
51  Xn336k84  21.4  139461.0  139197  38.0  139534.9  139210  66.0  139460.1  139303  
52  Xn344k43  22.6  42284.0  42146  21.7  42208.8  42099  39.7  42156.1  42056  
53  Xn351k40  25.2  26150.3  26021  33.7  26014.0  25946  51.5  25981.8  25938  
54  Xn359k29  48.9  52076.5  51706  34.9  51721.7  51509  112.0  51640.7  51555  
55  Xn367k17  13.1  23003.2  22902  22.0  22838.4  22814  117.3  22876.2  22814  
56  Xn376k94  7.1  147713.0  147713  28.3  147750.2  147717  70.3  147740.5  147714  
57  Xn384k52  34.5  66372.5  66116  40.2  66270.2  66081  56.8  66170.3  65997  
58  Xn393k38  20.8  38457.4  38298  28.7  38374.9  38269  49.3  38309.3  38260  
59  Xn401k29  60.4  66715.1  66453  49.5  66365.4  66243  110.2  66359.0  66212  
60  Xn411k19  23.8  19954.9  19792  34.7  19743.8  19718  126.0  19736.7  19721  
61  Xn420k130  22.2  107838.0  107798  53.2  107924.1  107798  87.7  107913.7  107798  
62  Xn429k61  38.2  65746.6  65563  41.5  65648.5  65501  65.6  65661.6  65470  
63  Xn439k37  39.6  36441.6  36395  34.6  36451.1  36395  57.1  36410.1  36395  
64  Xn449k29  59.9  56204.9  55761  64.9  55553.1  55378  132.6  55432.7  55330  
65  Xn459k26  60.6  24462.4  24209  42.8  24272.6  24181  92.9  24226.0  24145  
66  Xn469k138  36.3  222182.0  221909  86.7  222617.1  222070  142.3  222427.5  222235  
67  Xn480k70  50.4  89871.2  89694  67.0  89760.1  89535  73.1  89744.7  89513  
68  Xn491k59  52.2  67226.7  66965  71.9  66898.0  66633  81.9  66794.1  66607  
69  Xn502k39  80.8  69346.8  69284  63.6  69328.8  69253  177.7  69277.1  69247  
70  Xn513k21  35.0  24434.0  24332  33.1  24296.6  24201  99.4  24256.2  24201  
71  Xn524k137  27.3  155005.0  154709  80.7  154979.5  154774  207.3  155038.1  154787  
72  Xn536k96  62.1  95700.7  95524  107.5  95330.6  95122  144.5  95335.4  95112  
73  Xn548k50  64.0  86874.1  86710  84.2  86998.5  86822  136.6  86881.0  86778  
74  Xn561k42  68.9  43131.3  42952  60.6  42866.4  42756  77.2  42860.0  42733  
75  Xn573k30  112.0  51173.0  51092  188.2  50915.1  50780  782.4  50876.9  50801  
76  Xn586k159  78.5  190919.0  190612  175.3  190838.0  190543  234.3  190752.4  190442  
77  Xn599k92  73.0  109384.0  109056  125.9  109064.2  108813  166.9  108993.3  108576  
78  Xn613k62  74.8  60444.2  60229  117.3  59960.0  59778  103.6  59859.7  59654  
79  Xn627k43  162.7  62905.6  62783  239.7  62524.1  62366  543.1  62442.9  62254  
80  Xn641k35  140.4  64606.1  64462  158.8  64192.0  63839  304.4  64105.6  63859  
81  Xn655k131  47.2  106782.0  106780  150.5  106899.1  106829  253.2  106855.6  106804  
82  Xn670k126  61.2  147676.0  147045  264.1  147222.7  146705  267.7  147663.9  147163  
83  Xn685k75  73.9  68988.2  68646  156.7  68654.1  68425  177.0  68596.0  68496  
84  Xn701k44  210.1  83042.2  82888  253.2  82487.4  82293  368.0  82409.2  82174  
85  Xn716k35  225.8  44171.6  44021  264.3  43641.4  43525  437.2  43599.9  43498  
86  Xn733k159  111.6  137045.0  136832  244.5  136587.6  136366  334.2  136607.4  136424  
87  Xn749k98  127.2  78275.9  77952  313.9  77864.9  77715  308.3  77862.8  77605  
88  Xn766k71  242.1  115738.0  115443  383.0  115147.9  114683  330.5  115115.9  114812  
89  Xn783k48  235.5  73722.9  73447  269.7  73009.6  72781  351.2  72892.4  72738  
90  Xn801k40  432.6  74005.7  73830  289.2  73731.0  73587  424.0  73651.6  73466  
91  Xn819k171  148.9  159425.0  159164  374.3  158899.3  158611  675.6  158849.0  158592  
92  Xn837k142  173.2  195027.0  194804  463.4  194476.5  194266  634.9  194504.0  194356  
93  Xn856k95  153.7  89277.6  89060  288.4  89238.7  89118  314.6  89220.0  89020  
94  Xn876k59  409.3  100417.0  100177  495.4  99884.1  99715  543.1  99780.3  99610  
95  Xn895k37  410.2  54958.5  54713  321.9  54439.8  54172  500.2  54407.4  54254  
96  Xn916k207  226.1  330948.0  330639  560.8  330198.3  329836  1082.5  330153.2  329866  
97  Xn936k151  202.5  134530.0  133592  531.5  133512.9  133140  1022.2  133729.3  133376  
98  Xn957k87  311.2  85936.6  85697  432.9  85822.6  85672  307.9  85681.5  85555  
99  Xn979k58  687.2  120253.0  119994  554.0  119502.1  119194  928.4  119527.7  119188  
100  Xn1001k43  792.8  73985.4  73776  549.0  72956.0  72742  952.8  72903.3  72629  
Average gap:  0.74%  0.42%  0.30%  0.06%  0.24%  0.03% 
These tables highlight how the local search applied on resulted in several improvements over the best solutions obtained by the stateoftheart methods considered, with an average gap of 0.10% and 0.24% on the medium and large instances, respectively, in comparison with 0.14% and 0.30% for a classical search in . This improvement in solution quality comes at the price of an overall twofold increase of CPU time.
Another notable observation of these experiments is that the search in finds solutions which are structurally different. Indeed, we obtained some new best known solutions for unexpectedly small instances, with 256, 294, 322 and 327 customers (marked with a ). This is not a coincidence, given that the BKS listed at http://vrp.atdlab.inf.pucrio.br/ already originate from various previous articles and methods, over a large cumulated amount of test runs and parameter settings. For the particular case of the instance Xn256k16, one interesting characteristic was observed: only 16 vehicles are used, with a total capacity usage of 99.6%.
5 Conclusions and future work
In this article, we investigated decisionset decompositions for the classic CVRP. Our experiments show that decomposing the problem into Assignment and Sequencing decisions, and conducting local search in the Assignment space () generates consistently better results than heuristically searching in the complete search space . When doing so, each solution is systematically decoded by finding optimal routes (Sequencing decisions) with the Concorde TSP solver. However, the extra CPU dedicated to solution decoding is prohibitively high to employ this technique within stateoftheart metaheuristics. To circumvent this issue, the B&S neighborhood was employed to define an intermediate search space () which can be more efficiently searched while retaining a high solution quality.
Different techniques were proposed and evaluated for efficiently searching in space : neighborhood reduction, dynamic move filters, concatenation techniques and efficient memory structures. Moreover, tunneling techniques were employed to reshape into a search space more similar to as the search progresses. The combination of these techniques within the UHGS solver resulted in an significant improvement of solution accuracy. Multiple instances from the literature had their best known solution improved, and new stateoftheart results were obtained for the CVRP.
This improvement of search space, however, still results in some extra computational effort. Therefore, many research perspectives are open about how to fully exploit a larger search space such as without any time compromise. One possibility could be to employ the search in only in exceptional circumstances, for a selected subset of promising solutions (e.g., each new best incumbent solutions during the search). Other open research possibilities relate to the search space . Indeed, even with efficient memory structures and neighborhood reduction techniques, using Concorde for each solution evaluation remains impracticable. This effort could be mitigated if good and fast lower bounds were proposed for the cost of the routes, therefore permitting to filter a large proportion of moves as in Vidal (2017). Concorde is also not optimized to handle millions of small cardinality routes issued from a local search, such that a dedicated TSP solution procedure exploiting information from the current incumbent solution represents another promising research line. Finally, the proposed approach can be naturally evaluated for other variants of the CVRP, in the presence of different types of attributes that have an impact on the Assignment and Sequencing decision classes. These are all promising perspectives for future research at the crossroads of dynamic programming, integer programming, and metaheuristics.
Acknowledgements
This research has been funded by the Belgian Science Policy Office (BELSPO) in the Interuniversity Attraction Pole COMEX (http://comex.ulb.ac.be), the National Council for Scientific and Technological Development (CNPQ – grant number 308498/20151) and FAPERJ in Brazil (grant number E26/203.310/2016). The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Hercules Foundation and the Flemish Government – department EWI. This support is gratefully acknowledged. In addition, we would like to thank Luke Connolly for providing editorial consultation.