Heuristics for vehicle routing problems: Sequence or set optimization?\SetKw
Continuecontinue \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 vehicle-to-customer “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 speed-up 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. Decision-set decompositions, Metaheuristics, Dynamic programming, Integer programming, Large neighborhood search, Vehicle routing problem
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: customer-to-vehicle Assignments, and Sequencing choices for each route (Vidal et al. 2013b), a duality that has left long-standing impressions in the literature, from the early developments of route-first cluster-second (Bodin and Berman 1979, Beasley 1983) and cluster-first route-second constructive methods (Fisher and Jaikumar 1981), all the way to the set-covering-based 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 intra-route neighborhood search: Relocate, Swap and 2-opt 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 stand-alone neighborhood, is unlikely to result in further improvements. For this reason, it is very uncommon to observe the use of larger intra-route neighborhoods (e.g., 3-Opt or beyond) in recent state-of-the-art 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, inter-route moves frequently lead to TSP-suboptimal 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 :
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 non-trivial 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 multi-start local search (MS-LS), 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 state-of-the-art 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 decision-set decompositions.
Sequencing and Assignment decisions were first optimized separately in early constructive heuristics, giving rise to different families of methods. Route-first cluster-second 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 cluster-first route-second 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 a-priori 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, 2-opt* and their close variants allow to optimize both decision subsets. The associated neighborhood search methods form the basis of the vast majority of state-of-the-art algorithms. Petal algorithms have withstood the test of time, and high-quality 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 state-of-the-art 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 service-mode choices (Vidal et al. 2015a, Vidal 2017).
Decision-set 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 decision-set decomposition does not result from supplementary problem attributes, but is instead used to define exponential-size polynomially-searchable neighborhoods and transform the search space. Exponential-size 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 NP-hard 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 polynomially-searchable 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 exponential-size 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 stand-alone 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 arc-routing problems with possible cluster constraints, and dial-a-ride problems. One common characteristic of these studies is that they employed B&S as a stand-alone 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 stand-alone 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 local-search 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 ‘’ 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 customer-route 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 ‘[1,3]’ or ‘’).
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 ‘’.
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 worst-case 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 post-optimization 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:
From an initial solution containing a B–optimal tour, a local search in the space explores only B–optimal tours.
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.
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 non-B-optimal 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 state-of-the-art 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 :
Second, consider the cost of the classical CVRP move to filter non-promising solutions in :
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 non-improving moves:
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 intra-route and inter-route Relocate and Swap moves, for single vertices or generalized to pairs of consecutive vertices, as well as the 2-Opt and 2-Opt* 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 large-scale 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.
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 non-trivial, given that it is an instance-dependent 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:
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 constant-time 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 Constant-Time Evaluations
We use the concatenation strategy of Vidal et al. (2014a, 2015b) to perform efficient cost- and load-feasibility 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 .
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:
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 non-filtered 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 high-quality tours without jeopardizing the discovery of better route configurations. It can be efficiently implemented with a refinement of the hashtable-based 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 re-introduces 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 customer-to-vehicle 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 UHGS-BS. Except this modification of the local search, all other parameters and procedures of UHGS-BS 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 UHGS-BS 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) E5-2680v3 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 UHGS-BS, we used the code base made available at https://github.com/vidalthi/HGS-CARP, 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 trade-off 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 state-of-the-art 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 UHGS-BS – 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 state-of-the-art metaheuristics translates back into significant quality improvements.
This section now analyses the performance of the UHGS-BS 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 state-of-the-art 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 UHGS-BS 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 (p-values ). 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 trade-off 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 non-improving moves (), which is equivalent to using B&S only as a post-optimization 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 non-improving 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 non-filtered 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 UHGS-BS 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 re-evaluation, number of executions of the B&S optimization procedure, total number of iterations.
|#||Instance||Without tunneling strategy||With tunneling strategy|
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 state-of-the-art algorithms
Finally, this section reports detailed results of UHGS-BS, using the baseline configuration and the tunneling strategy, on the complete set of 100 instances proposed by Uchoa et al. (2017). The results of UHGS-BS are compared to that of the current state-of-the-art 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, UHGS-BS was executed 10 times for each instance. The maximum number of consecutive iterations without improvements was set to to evaluate UHGS-BS 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 UHGS-BS, 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.atd-lab.inf.puc-rio.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.
These tables highlight how the local search applied on resulted in several improvements over the best solutions obtained by the state-of-the-art 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.atd-lab.inf.puc-rio.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 X-n256-k16, 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 decision-set 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 state-of-the-art 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 state-of-the-art 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.
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/2015-1) and FAPERJ in Brazil (grant number E-26/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.