A hybrid metaheuristic for solving asymmetric distance-constrained vehicle routing problem

The Asymmetric Distance-Constrained Vehicle Routing Problem (ADVRP) is NP-hard as it is a natural extension of the NP-hard Vehicle Routing Problem. In ADVRP problem, each customer is visited exactly once by a vehicle; every tour starts and ends at a depot; and the traveled distance by each vehicle is not allowed to exceed a predetermined limit. We propose a hybrid metaheuristic algorithm combining the Randomized Variable Neighborhood Search (RVNS) and the Tabu Search (TS) to solve the problem. The combination of multiple neighborhoods and tabu mechanism is used for their capacity to escape local optima while exploring the solution space. Furthermore, the intensification and diversification phases are also included to deliver optimized and diversified solutions. Extensive numerical experiments and comparisons with all the state-of-the-art algorithms show that the proposed method is highly competitive in terms of solution quality and computation time, providing new best solutions for a number of instances.

special case of Asymmetric VRP.In Asymmetric VRPs, the real distance will depend upon the specific location of the nodes in the territory and also on the structure of the road network that communicates them.Moreover, when considering oriented networks, real distances might not have to be symmetric.
To the best of our knowledge, the ADVRP has not been addressed well in the literature before.There are a few papers which proposed exact algorithms [8,7] that could solve large instances.However, it could not always find a feasible solution due to the lack of memory.Moreover, when distance constraint is tight, solving the problem becomes harder, and these exact methods were stopped before it could find any feasible solution.For these reasons, developing metaheuristics to find a good feasible solution in a short computation time is a suitable approach.However, there exists only a metaheuristic [8] reported in the literature for ADVRP.The algorithm is based on the principles of Randomized Variable Neighborhood Search (RVNS).Unfortunately, it can be trapped into cycles as it returns to the points previously explored in the solution space.In this case, the algorithm gets stuck into local optima.Therefore, in this paper, we propose a hybrid metaheuristic combining RVNS and tabu search to overcome all these issues.
Besides ADVRP, we use algorithm to address also the Asymmetric Capacitated VRP (ACVRP) and Multiple Traveling Repairman Problem With Distance Constraints (MTRPD) to demonstrate its performance.Extensive numerical experiments on benchmark instances of ADVRP show that our algorithm could be comparable with the previous state-of-the-art algorithms in terms of solution quality and running time.Interestingly, in many cases, our proposed method is able to improve the best-known solutions available from the literature.
The remainder of the paper is structured as follows.We define the ADVRP in next section.Our proposed metaheuristic and its components are presented next, followed by the section dedicated to the experimental results, and the conclusion.

Problem description
The ADVRP is a generalization of the VRP, it is, thus, also a NP-hard problem.The ADVRP is defined on a complete graph K n with the vertex set V = {v 1 , v 2 , ..., v n } , an asymmetric distance matrix C = {c(v i , v j ) | i, j = 1, 2, ..., n} , where c(v i , v j ) is the distance between two vertices v i and v j , (c(v i , v j ) = c(v j , v i )) .A fleet of k identical vehicles is based at the depot v 1 .Suppose that the tour T = {R 1 , ..., R l , ..., R k } is a set of k routes obtained from these k vehicles, respectively.The route performed by the vehicle lth made up of a sequence of vertices, starting and ending at the depot v 1 .The length of the route R l is defined as the total traveled distances of the vehicle lth and its length cannot exceed the predetermined limit D max : The length of the tour T is the sum of its all k vehicles' total traveled distances: The ADVRP aims to determine the tour T with minimal length.

The proposed hybrid metaheuristic
The efficient hybrid metaheuristic we proposed brings together the components of Greedy Randomized Adaptive Search Procedures (GRASP) [9], Tabu search (TS) [10] and Randomized VNS [11].It, thus, consists of two phases performed iteratively: (1) the first phase uses the GRASP combining with k-means [12] to generate initial solutions; (2) the second phase then enhances solutions consecutively by means of an RVNS with multiple neighborhoods combined with a short-term memory mechanism in the spirit of TS.
An outline of our proposed algorithm (HVT algorithm) is shown in Algorithm 1.In Step 1, the initial solution is generated using the cluster first-route second scheme.This initial solution is then enhanced in Step 2. We use the idea of RVNS based on the principle of systematically exploring several different neighborhoods.Moreover, to avoid cycling, we use tabu lists of the recent types of moves in the solution space to prohibit reversing these moves.Thus, at each iteration of the algorithm, one neighborhood is selected randomly, then the selected neighborhood is explored, and the best move is chosen.This move must not be tabu unless it improves the current best solution T * .In Step 3, the new solution that has just been created from Step 2 is then added to a promising solution set P if its fitness value is not worse more than 10% to that of the T * .Step 2 is executed for a number of times to find good promising solutions to insert to the set P.
Step 4 uses elite solutions in P to exploiting the current solution space, thus enhances the solution quality.To explore more solution space, a diversification phase is added in Step 5.The search then restarts from the perturbed solution produced in Step 5.The search is stopped after maxITER number of iterations without improvement on the current best solution T * . (1) (3) In the following, the main components of our proposed algorithm (Algorithm 1) will be detailed.The search space is first presented in "Search space" section.The way to generate initial solution in Step 1 is then displayed in "The initial solution", "Neighborhoods of RVNS", "Tabu lists and tabu durations" section dedicate to explain the implementation of Step 2. They, thus, describe the set of eight neighborhood structures used for RVNS and the tabu mechanism for each neighborhood, respectively.We detail in "Intensification and diversification" section the intensification (Step 4) and diversification (Step 5) phases.

Search space
For a given solution T, let L(T) denote the total length of its routes and let V(T) denote the total violation of vehicle length.The total vehicle-length violation V(T) is computed on a route basis with respect to the value D max , thus is equal to R l ∈T max{D max − L(R l ), 0}.
Solutions are then evaluated according to the weighted fitness function L ′ (T ) = L(T ) + ρ * V (T ) , where ρ is the penalty parameter.

The initial solution
At Step 1 of Algorithm 1, we generate intial solution using the GRASP with clustering.The detail is described in Algorithm 2. The k-means [12] is first used to cluster the ADVRP problem to k smaller ADVRP problems.Specifically, the K n is converted into k smaller complete graphs K 1 n , K 2 n , ..., K k n .Routing is then performed iteratively on each cluster K l n (1 ≤ l ≤ k) to build a correspond route R l .All routes are initialized with the main depot v 1 .Each vertex of the K n is then added to the routes by using its Restricted Candidate List (RCL).The RCL of a vertex includes a number of vertices which are closest to it.At an iteration, assuming vertex v e be the current last vertex of the route R l , an unvisited vertex v is picked randomly from its RCL to add to R l such that this addition does not make the constraint violation.If there does not exist any such vertex v, we then accept the infeasibility.Thus, a not-yet-routed vertex v ′ is picked up randomly to add to the last position of the route R l .A solution is generated when all vertices of K n are routed.
The procedure then returns the feasible solution if any.Otherwise, for added randomness in routing, it tries to generate n solutions, then the one with the minimum fitness value will be returned.

Neighborhoods of RVNS
The RVNS in Step 2 of Algorithm 1 exploits eight neighborhoods including different intra-and inter-route moves as following: Five intra-route neighborhoods: -Remove-insert move each vertex is shifted from its current location to the last position of the same route.-Swap-adjacent move two adjacent vertices in the same route are exchanged.
-Swap move two vertices in the same route are exchanged.
-2-opt move for each pair of vertices v i , v j in the same route, the edges emanating from them

Three inter-route neighborhoods:
-Exchange-route move the two vertices of two different routes are exchanged.
-Insert-route move a vertex is shifted from its current position to another position of the other route.-Cross-route move two segments {v i , v k } and {v j , v l } on two different routes are exchanged.

Tabu lists and tabu durations
The tabu lists for the eight moves described above are included in our algorithm.The solution elements receiving a tabu status following an intra-route move are: -Remove-insert move the position of vertex v i just moved at the end of the route can- not be changed by the same type of move while it is tabu; -Swap-adjacent, swap move vertices v i and v j just swapped cannot be swapped again while they are tabu; -2-opt move a 2-opt move applied to vertices v i and v j cannot be applied again to the same vertices.-3-opt a 3-opt move applied to vertices v i , v j and v k cannot be applied again to the same vertices.
while for inter-route moves: -Exchange-route move vertices v i in route R l and v j in route R h cannot be swapped again by the same type of move while they are tabu.-Insert-route move the position of vertex v i in the route R l just inserted to the position jth of the route R h cannot be changed by the same type of move while it is tabu.-Cross-route move two segments {v i , v k } on R l and {v j , v l } on R h cannot be swapped again by the same type of move while they are tabu.
A tabu status is assigned to an element for θ iterations, where θ is randomly selected from an uniform interval [13,14].At each iteration, a neighborhood is selected to explore.Its best neighbor T ′ is then accepted as the next starting solution if the corresponding move is non-tabu and T ′ is better than the current starting solution T. Otherwise, T ′ must be better than the current best solution T * .

Intensification and diversification
We use a promising solution set P as a pool of high-quality solutions found during the algorithm.This set is initialized empty and limited in size.Each solution T produced by the RVNS+TS in Step 2 of Algorithm 1 is inserted into the set P if its fitness is not worse more than 10% to that of the current best solution T * .Step 2 is thus executed for a number of times to find such solutions T for inserting them into the set P. Once the set P is full, the intensification phase (Step 4 of the Algorithm 1) will be triggered to enhance the quality of the solutions in P. Specifically, only the RVNS of Step 2 is applied to each solution in P. Without the tabu list restriction, the RVNS starts by applying a neighborhood selected randomly in the set NL of eight neighborhoods described in the "Neighborhoods of RVNS" section.The selected neighborhood is then searched on all possible moves, and the best neighbor is returned.The current solution will be modified if it is worse than its best neighbor.Otherwise, this selected neighborhood is removed from the set NL as it does not produce any improvement.The process is repeated until the set NL is empty.The best solution T produced by this intensification process will be served as the input for the diversification phase (Step 5 of the Algorithm 1).Using only intensification, the search could get stuck in the local optima.Thus, diversification then proceeds to perturb the search that starts from the solution T by per- forming a combined shakings of the solution (see Algorithm 3).Creating a new working solution from the T helps to conserve the best solution characteristics encountered so far.On the other hand, implementing shaking process to perturb a new working solution provides a certain level of diversity to the search.In this work, there are two shaking mechanisms used to diversify the solution: -Shaking intra-route: We use the shaking mechanism, called double-bridge, was originally developed by [15].The structure of the double-bridge move derives from a special 4-opt neighborhood where edges added and dropped need not be successively adjacent.This mechanism can also be seen as a permutation of two disjoint segments of a route.The double-bridge shaking is described in the Algorithm 4.
-Shaking inter-route: We randomly choose two routes R x and R y in the solution, then exchange a number of vertices between them.The detailed description of the implementation is given in the Algorithm 5.

Evaluation
Our proposed algorithm is implemented in C++.Experiments are conducted on an Intel Pentium core i7 duo 2.10 Ghz CPU, 8 GB RAM.The performance of the proposed algorithm is evaluated through comparison with published results on the instances provided in the literature of ADVRP, ACVRP, MTRPD.Through preliminary experiments, we observed that the values α = 5 , num = 5 , ρ = 100 , |E| = 5 , and maxITER = 100 resulted in a good trade-off between solution quality and running time.

Comparison with ADVRP's algorithms
The performance of our proposed algorithm (HVT) is evaluated through comparison with published results on the Group 1 instances of ADVRP provided in [16].These instances are divided into three sets, denoted as D max(1) , D max(2) , D max(3) , in descending order of the upper bound value D max .The number of customers is in the range of [40,1000].All pub- lished results used for competing are taken from: -M5SBB : The exact method (Multi Start Branch and Bound) of [16].
We report our best results over 10 runs.For concision sake, only aggregated results are provided in this section.A detailed comparison, instance by instance, may be found in Tables 4, 5, 6 of the Annex.Table 1 sums up comparison over all three sets.As the results of VNS for the D max(3) set have not been published, we use the hyphen "-" to indicate its unavailable in Table 1.The first two rows respectively show for each competing algorithm, the percentage of times that algorithm could obtain the optimal solutions (%opt) and could not find feasible solutions (% of no feasible).Note that all current published optimal solutions and best known solutions (BKS) are produced by M5SBB [16].The row Gap 1 [%] indicates the average gaps of best solutions obtained by each algorithm with repect to the optimal solutions (OPT) for instances whose optimal solutions have been published.Similarly, the row Gap 2 [%] displays the average gaps relative to BKS for instances whose the optimal solutions have not been found yet.More specifically, Gap 1 [%] and Gap 2 [%] in percentage on each instance are calculated as follows: One observes that only our algorithm could find feasible solutions in all cases.Indeed, for the easiest D max(1) set with the largest value of upper bound D max , all algorithms could find feasible solutions.However, for D max (2) set with smaller value of D max , the problem becomes harder to solve.As the result, the number of instances that the M5SBB could not find feasible solutions is three, and this value is even up to 19 (14.29%) for the VNS metaheuristic (see Table 5).On the hardest D max (3) instances, the M5SBB could not find feasible solutions for three cases, while the results for the VNS has not been published due to its worse results on the D max(2) .Especially, there are 11 intances, though the M5SBB cannot reach the optimal solutions due to lack of memory, the HVT obtains better solutions (see Table 6).
Experimental results illustrated clearly the superior performance of our HVT algorithm compared to the VNS metaheuristic of [16].Table 1 shows that the HVT provides high-quality solutions, with an average gap of 4.24% to the optimal solutions, compared to 10.97% of the VNS on the D max(1) , and D max(2) dataset.Overall 131 instances, the HVT produces 124 better solutions than those of VNS, and 41 optimal solutions while the VNS could not find any.The VNS was executed on intel Core(TM) 2 CPU 6600 2.4GHz with 3.24 GB of RAM which has different characteristics than our machine; therefore, comparing the running time between two algorithms is evaluated relatively.As reported in [16], our running time is comparable with that of the VNS.

Comparison with ACVRP's algorithms
In this section, we present the performance comparison between our HVT algorithm and published results on the ACVRP instances provided by Pessoa et al. in [17].The ACVRP is the Capacitated VRP in which the matrix cost is asymmetric.Table 2 displays the comparison between the best solutions obtained by our HVT algorithm and those obtained by two other algorithms: -VL: a metaheuristic combines some heuristic concepts with compact mixed-integer linear programming (MILP) formulation, proposed by Valeria Leggieria and Mohamed Haouari in [18].-AS: a column generation approach uses metaheuristic to generate columns and a sequence of set partitioning models for the Mixed-Integer Programming solver, proposed by A. Subramanian et al. in [19].
for each instance of the ACVRP.
The first three columns give the instance name, the number of vertices and vehicles, respectively.The column OPT indicates the optimal solutions reported by [20,6].We report for each algorithm the best solutions (Best.Sol column) found.The last three columns give our average results over 10 runs, the gap between our best found solution and the optimal solution, and our running time in seconds, respectively.One could observe that our HVT algorithm is comparable with the others in terms of solution quality and running time.It produces very near-optimal solutions, with an average gap of 0.01% to the optimal solutions.

Comparison with MTRPD's algorithms
In the MTRPD, a fleet of homogeneous vehicles is dispatched to serve a set of customers.Each customer is serviced exactly once by a vehicle.Each vehicle which starts from and ends at the depot is not allowed to travel a distance longer than a predetermined limit.The objective is to minimize the total waiting time of all customers after the vehicles leave the depot.Our algorithm runs on the MTRPD dataset inherited several instances in [20].The optimal solutions for these instances are found using the exact algorithms in [20,21].The dataset includes six TSP instances from the TSPLIB such as brd14051, d15112, d18512, fnl4461, nrw1379, and pr1002.For each TSP instance, ten MTRPD instances are generated by randomly selecting ten subsets of n vertices, where n = 30, 40, 50, 60, 70, and 80.Let (d max ) be the distance to the farthest vertices from the depot.The distance con- straint gets the values 2 × d max , 2.5 × d max , and 3 × d max .The exact algorithms for these instances are extracted from [5,15].Therefore, in total, 900 MTRPD instances are used in our experiment.
For small instances of size 30-50 vertices, we compare our results to optimal solutions produced by [15].However, for large instances with more than 50 vertices, there does not exist published optimal solutions.We thus compare our results to the lower bounds of the MTRPD which is the optimal solutions of the MTRP obtained in [12].Table 3 displays our aggregated results over 900 instances in terms of CPU time (in seconds), gap to the optimal solution on small instances and gap to lower bound on large instances.Each row represents to a set of 10 instances.As shown in Table 3, our HVT algorithm is capable of finding the optimal solutions for all small instances in a reasonable amount of time (0.52 seconds on average).It is, thus, better than those obtained by GRASP+VNS in [13] which fails to find the optimal solutions for all instances with 50 vertices.Moreover, for most instances consisting of 60-80 vertices, our solutions fall into the range of 0.30-0.32% of the lower bound.Consequently, our HVT algorithm could produce much better results compared to that of [13] which ranges from 2.73 to 4.75%.

Conclusion
We have proposed a new effective metaheuristic algorithm for the ADVRP, which combines GRASP with clustering, tabu search and RVNS.Extensive computational experiments and comparisons with published algorithms on all benchmark instances show that the proposed algorithm is highly competitive.It could provide not only new best solutions but also first published feasible solutions for some instances with tight value of D max .Additionally, our algorithm provides better solutions than the state-of- the-art metaheuristic algorithm for 124 out of 131 cases.However, the efficiency and

Table 3 The average experimental results for MTRPD instances
running time for the large instances with up to 1000 vertices need to be improved.It is our aim for future research.