A Hybrid GRASP+VND Heuristic for the Two-Echelon Vehicle Routing Problem Arising in City Logistics

The two-echelon vehicle routing problem (2E-VRP) is a variant of the classical vehicle routing problem (VRP) arising in twolevel transportation systems such as those encountered in the context of city logistics. In the 2E-VRP, freight from a depot is compulsorily delivered through intermediate depots, named satellites. The first echelons are routes that distribute freight from depot to satellites, and the second are those from satellites to customers. This problem is solved by a hybrid heuristic which is composed of a greedy randomized adaptive search procedure (GRASP) with a route-first cluster-second procedure embedded and a variable neighborhood descent (VND), called GRASP+VND hereafter. Firstly, an extended split algorithm in the GRASP continuously splits randomly generated permutations of all customers and assigns customers to satellites reasonably until a feasible assignment appears, and a complete 2E-VRP feasible solution is obtained by solving the first echelon problem subsequently and, secondly, a VND phase attempts to improve this solution until no more improvements can be found. The process above is iterated until the maximum number of iterations is reached. Computational tests conducted on three sets of benchmark instances from the literature show that our algorithm is both effective and efficient and outperforms the best existing heuristics for the 2E-VRP.


Introduction
The transportation of freight constitutes an extremely important activity taking place in urban areas, but it is also very disturbing.The increase in the number of freight transportation trucks using urban roads makes a more and more significant contribution to traffic congestion and many associated negative environmental impacts, such as air pollution and noise.For the purpose of preventing the urban environment from getting worse, many municipalities place restrictions on these big trucks to keep them out of their city centers by creating peripheral intermediate facilities, called satellites.External carriers need to supply these satellites from central depots and then smaller and environmentally friendly vehicles would distribute the freight downtown from these satellites [1][2][3].Therefore, two distribution echelons are involved in city logistics.With the customer demands, satellite capacities, and vehicle capacities from the two levels known in advance, the two-echelon vehicle routing problem (2E-VRP) [4,5] consists in building a set of the least-cost trips (the fewest number of vehicles used and least vehicle traveling cost) for the two echelons.
Perboli et al. [4] introduced a family of two-echelon vehicle routing problems and proposed a three-index flow-based formulation for the 2E-VRP.The authors also introduced some valid inequalities and two math-heuristics based on the 2E-VRP model, which were used within a branch-and-cut framework.They were able to solve to optimality instances containing up to 21 customers.
Perboli and Tadei [5] proposed several new classes of valid inequalities based on the traveling salesman problem (TSP) and the VRP and strengthened the previous 2E-VRP formulation with new cuts (including capacity cuts), which allowed their algorithm to solve seven new instances to optimality and reduce the optimality gap on several other instances.
Jepsen et al. [6] presented an edge flow based model for the 2E-VRP and employed a specialized branching scheme to branch on infeasible integer solutions in their branch-and-cut algorithm to obtain feasible solutions.Their algorithm was able to solve 47 instances to optimality, surpassing previous exact algorithms.They found that the coupling between the two echelons in the 2E-VRP would pose a challenge to incorporate.
Baldacci et al. [7] proposed a new mathematical formulation of the 2E-VRP (used to derive valid lower bounds) and a new exact method.They decomposed the 2E-VRP into a limited set of multidepot vehicle routing problems (MDVRPs) with side constraints.Computational results on extensive benchmark instances showed that their exact algorithm outperformed the state-of-the-art exact methods in terms of size, number of problems solved to optimality, and computing time.
Since exact algorithms are usually computationally expensive for large-scale combinatorial optimization problems, approximate solutions with sufficient accuracy that can be obtained fast are often desired in practice.
Crainic et al. [8] developed a family of multistart heuristics, based on separating the depot-to-satellite transfer and the satellite-to-customer delivery by iteratively solving the two resulting routing subproblems, while adjusting the satellite workloads that linked them.Besides, an intensification phase aiming at improving feasible solutions by local search was followed by a diversification phase to avoid local optimum in their algorithms.Their ideas are very useful to handle the 2E-VRP; however, they obtained relatively poor results even in small-scale instances.
Meihua et al. [9] proposed a hybrid ant colony optimization algorithm which combined three heuristics for the 2E-VRP.They firstly divided the problem into several VRPs by a separation strategy and then applied improved ant colony optimization with multiple neighborhood descent to build better feasible solutions.Computational tests on 22 benchmark instances from the literature showed that their algorithm was better than previous published algorithms.
Hemmelmayr et al. [10] developed an adaptive large neighborhood search (ALNS) heuristic for the 2E-VRP, based on the destroy-and-repair principle in which two different sets of operators (destroy operators and repair operators) are alternated.They used existing operators and new operators designed specifically for the 2E-VRP.Their algorithm was shown to provide better solutions and outperform existing heuristic methods for the 2E-VRP but is much complicated to implement in practice because of adopting too many kinds of operators (8 kinds of destroy operators, 5 kinds of repair operators, and 5 kinds of local search operators) and many parameters.
Crainic et al. [11] proposed a heuristic algorithm based on greedy randomized adaptive search procedure (GRASP) combined with path relinking to address the 2E-VRP.The problem was treated by separating the depot-to-satellite transfer and the satellite-to-customer delivery and iteratively solving the two resulting routing subproblems, while adjusting the satellite workloads that link them; this idea is the same as Crainic et al. [8].The path relinking procedure with feasibility search was applied to link current solution to elite one.The computational results on instances with up to 50 customers and 5 satellites showed that the proposed heuristic can improve literature results, in both efficiency and accuracy, but their solution quality cannot outperform the ALNS of Hemmelmayr et al. [10].
According to the published computational results of these heuristic methods mentioned above, we can retrieve that the best existing heuristic for the 2E-VRP is the ALNS algorithm of Hemmelmayr et al. [10].
Greedy randomized adaptive search procedure (GRASP) is one of the most well-known multistart heuristics for combinatorial optimization problems.It was introduced by Feo and Resende [12].Each GRASP iteration consists basically of constructing a feasible solution and then applying a local search procedure to improve it until a local optimum is found, and the best overall solution is kept as the final result [12,13].
Variable neighborhood descent (VND) is a deterministic version of variable neighborhood search (VNS), originally proposed by Mladenović and Hansen [14].VNS is a metaheuristic for solving combinatorial optimization problems, whose basic idea is a systematic change of neighborhoods, both within a descent phase to find a local optimum and in a perturbation phase to get out of the corresponding valley [14,15].VNS explores an ordered list of neighborhoods.It starts with a given neighborhood and switches to the next one in the list when it finds a local minimum.The search is reinitialized from the first neighborhood whenever a new better solution is found or when all neighborhoods have been checked.VND also changes the neighborhood operators once the search is stuck in a local optimum, but it differs from VNS that no random perturbation is applied.
The hybrid of GRASP and VND has successfully solved several kinds of routing problems, such as pickup-anddelivery traveling salesman problem [16], truck and trailer routing problem [17], traveling repairman problem [18], three-dimensional bin packing problem [19], and school bus routing problem [20].Because the combinations of GRASP and VND are not only effective and efficient in solving combinatorial optimization problems but also very simple to implement, we design a hybrid GRASP+VND heuristic for the 2E-VRP based on the characteristics of this problem to meet the practical requirements arising in city logistics.
The remainder of this paper is organized as follows.Section 2 first describes the 2E-VRP and its mathematical model.Section 3 gives the details of our hybrid heuristic for the 2E-VRP.Section 4 presents computational results on three sets of instances.Section 5 contains the discussion and Section 6 gives the conclusion.

Problem Description and Mathematical Formulation
The definition of 2E-VRP provided here follows those in [4,7,10].2E-VRP can be defined on a weighted, complete, and undirected graph  = (, ) with node set  and arc set . Node set  is composed of three kinds of nodes: a depot  0 , a subset   of   satellites, and a subset   of   customers.
In the arc set , each arc (, ) represents the shortest path in the actual road network, with a known traveling cost   .Costs in both levels are assumed to satisfy the triangle inequalities.Each customer  ∈   has an associated demand   , known in advance, and cannot be split.The demand should not be delivered by direct shipping from the depot but must be consolidated in a satellite.The deliveries to the satellites on the first echelon can be split.Each vehicle has a capacity constraint that has to be respected.This capacity is the same for all vehicles belonging to the same level but may differ for each level.The capacities of the first-and the second-level vehicles are denoted by  1 and  2 , respectively.The total number of vehicles available is given by  1 for the first level and  2 for the second.Each satellite  ∈   has a capacity   that limits the total amount of customers' demands that can be delivered to it by first-level routes.Moreover, there are   first-level vehicles available at each satellite  ∈   .No additional limitation on the route size, neither in length nor in number of visited customers, is introduced.Figure 1 illustrates a 2E-VRP solution with a depot, three satellites, and nine customers.The satellites are represented by five-pointed stars and customers by circles and the depot by a square.The routes that serve the satellites from the depot are called the first-level routes.The second-level routes are those that start from a satellite, visit the customers, and return to the same satellite.Vehicle routes of the first and the second echelon are represented by continuous thick lines and dashed thin lines, respectively.

Mathematical Formulation.
Based on the above description, the 2E-VRP can be formulated as follows [7].Let M be the index set of all first-level routes, and let M  ⊆ M be the subset of first-level routes serving satellite  =∈   .Let   and () be the subset of satellites visited and subset of arcs traversed by a first-level route  ∈ M, respectively.Let R  be the index of the second-level routes passing through satellite  ∈   , and let R  ⊆ R  be the subset of routes passing through satellite  ∈   and customer i ∈   .The subset of customers visited and the subset of arcs traversed by a second-level route  ∈ R  are represented by   and (), respectively.Let   be a binary variable equal to 1 if and only if route  ∈ M is in 2E-VRP solution,   a binary variable equal to 1 if and only if route  ∈ R  of satellite  ∈   is in 2E-VRP solution, and   a nonnegative integer variable representing the quantity delivered by first-level route  ∈ M to satellite  ∈   .And the mathematical model of the 2E-VRP can be stated as follows [7]: subject to: The objective function (1) aims to minimize the total cost of two echelons.Constraints (2) ensure that each customer  ∈   must be visited by exactly one second-level route.Constraints (3), ( 4), and ( 6) impose the upper bounds on the number of first-and second-level routes in the solution.Constraints (5) specify the capacities of each satellite.The balance between the quantity delivered by first-level routes to every satellite and the customers demands supplied from this satellite is ensured by constraints (7).Constraints (8) impose that the vehicle capacity of the first-level vehicles is not exceeded.
Each second-level route must begin and end at the same satellite, and each customer must be served by exactly one second-level vehicle.The demand of each satellite is the total demand of its assigned customers, so each satellite must receive enough freight from the depot to satisfy the customers of its second-level routes.Besides, any change to the customer-to-satellite assignment affects the first-level routing and therefore has an impact on the first-level transportation costs.The objective function to be minimized is the global transportation costs in both levels.2E-VRP can be easily seen to be a reduction to the vehicle routing problem (VRP), which is a special case of 2E-VRP arising when just one satellite is considered, so it is also NP-hard [10].

The Hybrid GRASP+VND Heuristic
The proposed algorithm is a memoryless multistart heuristic method in which each iteration consists of two phases: a greedy randomized adaptive search procedure (GRASP) construction phase and a variable neighborhood descent (VND) improvement phase.Since the solutions generated by the GRASP construction phase are not guaranteed to be locally optimal, it may be very beneficial to apply a local search to further improve each constructed solution [12,17].We use the term GRASP+VND to symbolize this hybrid algorithm.GRASP+VND independently creates and improves a number of initial solutions and in the end returns the best solution obtained during the entire search.

Search Space.
The search is restricted to feasible solutions.We do not allow any violations of the constraints on the vehicle capacity, number of vehicles available, and the capacity constraints of the satellites.The reason that we always maintain the feasibility of solutions is that it is quite difficult and time consuming to perform feasibility repairing (due to the constraints of the number and capacity of the second-level vehicles) after the solutions are infeasible in the search.
The initial solutions of GRASP+VND are generated from solutions encoded as random permutations (TSP tours covering the entire   customers).Any random permutation  can be converted into a 2E-VRP solution  (with respect to the orders of customers in ) using an extended version of the splitting procedure split of Nguyen et al. [21].Customers are assigned to satellites by this extended split according to their orders in  in a relatively reasonable manner and then the demand of each satellite is obtained by the summation of the demands of all its assigned customers.
With the obtained demands of all the satellites, the firstlevel routes are constructed as follows [10,22]: we first create as many full-load direct trips from the depot to each satellite as possible until the remaining demand of each satellite is less than the capacity  1 of the first-level vehicle and then solve a TSP tour using the saving algorithm of Clarke and Wright [23], as it was handled by Hemmelmayr et al. [10].

Initial Solution Construction
Phase.The goal of GRASP is to create initial feasible solutions for the second phase to further improve them.For the hybrid GRASP+VND algorithm to work well, it is essential that solutions of relatively high quality are constructed during the solution construction phase.
The idea of splitting a good TSP tour to a VRP solution was first introduced by Beasley [24], as a route-first clustersecond heuristic for the VRP.A random permutation  of all customers is first cut into second-level routes by an extended version of the split algorithm of Nguyen et al. [21].The first stage of our extended algorithm consists of minimizing the number of second-level routes, that is, the number of second-level vehicles used.This extension [25] of split defines a second label   for the number of arcs on the shortest path (i.e., the number of second-level vehicles used up to ), sets  0 = 0 at the beginning, and tests   before   (the cost of the shortest path up to ) to minimize first the number of second-level vehicles.
The adoption of this extension is because the number of routes is often considered as the primary objective of VRP and we also need to obtain feasible solutions concerning the upper bound on the number of the second-level vehicles.
Algorithm 1 gives the framework of the initial solution construction procedure in our GRASP+VND.
Algorithm 2 implements the extended split for a random permutation  of all the   customers.  records the assigned satellite of the jth customer in  and   records the predecessor of this customer on the shortest path.W⋅(,  0 ) represents the additional penalty for each second-level route according to the distance between its assigned satellite and the depot.W is a nonnegative value, it can be adjusted to change the impact of the first-level routes on the complete initial solution, and, as it increases, the first-level routes are given more importance.It is set to 1 as a default value to obtain a compromised consideration on routes of both echelons.

Solution Improvement
Phase.After an initial feasible solution is generated, an attempt is made to improve it exhaustively.We apply a VND approach to perform the solution improvement procedure until no more improvements can be found.VND is a simplified variant of the VNS heuristic, in which the shaking phase is omitted.Therefore, VND is usually completely deterministic contrary to VNS.
The neighborhoods of our VND can be classified into two types: intersecond-level-route and intrasecond-level-route. Traditional VRP local search operators contain interroute and intraroute operators.Intraroute operators search and improve a single route at a time, while interroute operators deal with several routes simultaneously.When handling the 2E-VRP, we slightly modify the name of this classification for the second-level routes.The existing three interroute neighborhood structures Shift, Swap, and Cross and three intraroute neighborhood structures 2-Opt, Or-Opt, and Exchange of Subramanian et al. [26][27][28] are employed, and we also propose two new neighborhoods Satellites-Change and Satellites-Swap specially for the 2E-VRP.
The Satellites-Change neighborhood changes the assigned satellite of a second-level route at a time, while Satellites-Swap neighborhood swaps the assigned satellites of two distinct second-level routes if they belong to different satellites.Both of them may modify the demands of the satellites and hence may change the total cost of the 2E-VRP.The Satellites-Change can be seen as a special case of Satellites-Swap when a second-level route swaps its satellite with an empty secondlevel route.We define both of them as intersecond-levelroute neighborhoods.Neighborhoods of Satellites-Swap and Satellites-Change are shown in Figures 6 and 7, respectively.
Once a new initial feasible solution is generated, an intersecond-level-route move is performed to modify the second-level routes first and then an intrasecond-levelroute improvement procedure that uses 2-Opt, Or-Opt, and Exchange neighborhoods sequentially is triggered to improve the modified routes if they are feasible.Details of the three neighborhoods can be found in the survey of Bräysy and Gendreau [30].2-Opt moves delete two nonadjacent arcs and add two new arcs to generate a new route (see Figure 8); Or-Opt moves at most two adjacent customers back and forth in the current second-level route (see Figure 9), while Exchange exchanges the positions of two customers or the positions of a customer and two adjacent customers in the same second-level route (see Figure 10).Or-Opt used in this paper can be seen as the intraroute version of ℎ (1, 0) and ℎ (2, 0), while Exchange can be seen as the intraroute version of  (1, 1) and  (1,2).
All the neighborhoods are searched until no more improvements can be obtained, in a first-accept manner, and infeasible move against the capacity of each echelon vehicles is never allowed.When no more moves that improve the current solution can be found in a neighborhood, the search continues with the next neighborhood.VND ends when the current solution is a local optimum with respect to all the applied neighborhoods.The framework of VND [31,32] is briefly shown in Algorithm 3.Each  ℎ (ℎ < ℎ max ) is actually a combination of an intersecond-level-route neighborhood while ℎ ≤ ℎ max do (7)  ← LocalSearch(,  ℎ ); (8) if Cost() < Cost( * ) then (9) improve flag ← true; (10)  * ← ; (11) ℎ ← ℎ + 1; (12) end if (13) end while (14) until improve flag = false (15) return the local optimum  * found; Algorithm 3: VND for the improvement of a solution . (1) S * ← ; (7) end if (8) end for (9) return the best solution S * found; Algorithm 4: GRASP+VND for 2E-VRP outline.and the three intrasecond-level-route neighborhoods.But the neighborhood  ℎ max placed at the end is an exception and it is only composed of the above three intrasecond-levelroute neighborhoods, which are searched circularly until none of them can improve every second-level route of current solution.Every time the solution is modified due to a feasible intersecond-level-route move, the three intrasecond-levelroute operators are performed sequentially to improve the newly generated second-level routes.

Framework of the Resulting GRASP+VND.
Our hybrid GRASP+VND heuristic for the 2E-VRP combines a GRASP construction phase with a VND improvement phase.Both of them are iterated maxi times, and the best solution found of all iterations is kept as the final result.An outline of the proposed hybrid algorithm is shown in Algorithm 4. The term S * stands for the global best solution found.

Numeric Verification
We conducted computational study on three sets of benchmark instances in order to assess the proposed hybrid GRASP+VND algorithm with respect to solution quality and computing times.Our hybrid algorithm was coded in C++, compiled by Microsoft Visual C++ 6.0, and run on a single core of an Intel Pentium Dual-Core E5500 processor (2.8 GHz) and 2 GB of memory.

Instance Description.
The proposed hybrid algorithm was tested on the 2E-VRP benchmark instances Set 2, Set 3, and Set 4. Instances Set 2 and Set 3 were proposed by Perboli et al. [4] and they are based on the following VRP instances of Christofides and Eilon: E-n22-k4, E-n33-k4, and E-n51-k5.The cost matrix of each instance is given by the corresponding VRP instance.The capacity of the first-level vehicles is 2.5 times the capacity of the second.The capacity and the available number of the second-level vehicles are equal to the corresponding VRP instance.The satellites are located at several randomly chosen customers.Instances in this set range between 21 and 50 customers and consider 2 or 4 satellites.Set 4 was proposed by Crainic et al. [33] and contains 54 instances.Each instance has 50 customers and the number of available satellites is 2, 3, or 5.They were generated using three different customer distributions (Random, Centroids, and Quadrants) and three satellite location patterns (Random, Sliced, and Forbidden Zone).A summary of the characteristics of these instances can be found in Hemmelmayr et al. [10].

Parameter Setting.
The nonnegative weight W in Algorithm 2 is set to 1 as default, which has been determined  by experience.The maximum number of iterations maxi (i.e., the number of initial feasible solutions generated and improved by VND) is set to 500, and each instance is conducted on five independent runs, as it was done in [10].We record the best and average solutions, as well as the average time needed to find the best solution of five runs for each instance.

Computational Results.
In this section, we compare the computational results of our hybrid GRASP+VND heuristic with the best existing heuristic method for the 2E-VRP, that is, the ALNS algorithm of Hemmelmayr et al. [10].
The results of instances Set 2, Set 3, and Set 4 obtained by the two compared heuristics are shown in Tables 1, 2, and 3, respectively.The column instance name gives the names of these instances.The column BKS gives the best-known results of these 2E-VRP instances (published in Baldacci et al. [7]).The columns ALNS and GRASP+VND report the computational results of the two algorithms, respectively.The ALNS algorithm [10] and our GRASP+VND algorithm both report the best and average solutions (containing their percentage deviation from the best-know solutions), as well as the average time (in seconds) needed to find the best solutions of five runs.
The two heuristics were both implemented in C++ but used different workstations.The ALNS heuristic was tested on a single core of an AMD Opteron 275 processor (2.2 GHz).In order to make a fair comparison, we choose the running times of ALNS as a benchmark and scale the running times  of our GRASP+VND according to the CPU performances reported at http://www.cpubenchmark.net/cpulist.php.Our machine is 1.34 times faster than the AMD Opteron 275 processor (2.2 GHz) used by Hemmelmayr et al. [10]; therefore, our running times are multiplied by 1.34.The times shown in Tables 1, 2, and 3 are already the scaled value of actual running times.Values in bold fonts correspond to those that our GRASP+VND outperforms the ALNS, while those italic mean that our algorithm is worse than the ALNS.

Discussion
As seen from Tables 1, 2, and 3, our GRASP+VND and the ALNS can both find the best-known solutions for all the 39 instances in Set 2 and Set 3, and the ALNS can find the best-known 44 solutions for instances in Set 4, while the GRASP+VND can find 52; besides, there are ten better results in the best values found by our GRASP+VND.In terms of the average solutions, the two heuristics obtained the same results in Set 2 and Set 3, but we can see from Table 3 that our GRASP+VND is much better than the ALNS.In regard to    the deviations of the solutions from the best-known results, we can find that our GRASP+VND is as good as the ALNS in Set 2 and Set 3 but outperforms the ALNS in Set 4. With respect to the running times, our GRASP+VND is as fast as the ALNS in Set 2, a little faster than the ALNS in Set 3 and much faster in Set 4; besides, the running times of our GRASP+VND are limited and reasonable.Our GRASP+VND adopted 8 kinds of local search operators, while the ALNS used 8 kinds of destroy operators, 5 kinds of repair operators, and 5 kinds of local search operators; besides, our GRASP+VND used only 2 parameters, but the ALNS employed 9 parameters, so the GRASP+VND is much easier to implement and tune.From the comparisons above, we can retrieve that our GRASP+VND is both effective and efficient and outperforms the best existing heuristics for the 2E-VRP.

Conclusion
This paper has presented a very simple hybrid GRASP+VND heuristic to address the two-echelon vehicle routing problem (2E-VRP), a newly defined multiechelon variant of the classical vehicle routing problem (VRP).The heuristic is composed of a GRASP construction phase (embedding an extended split algorithm) to generate feasible and relatively good solutions and a VND phase to improve them.
Computational experiments on three sets of benchmark instances from the literature showed the effectiveness and efficiency of the proposed algorithm, and it outperformed the best existing heuristic methods for the 2E-VRP in both solutions quality and computing times.Moreover, the implementation of the hybrid heuristic is much easier than other heuristics.As a result, the proposed hybrid heuristic is more suitable for handling practical 2E-VRP and its relative variants arising in city logistics.

Figure 9 :
Figure 9: Two kinds of Or-Opt neighborhoods.

Figure 10 :
Figure 10: Two kinds of Exchange neighborhoods.

Table 1 :
Computational results comparison for Set 2 for the 2E-VRP.

Table 2 :
Computational results comparison for Set 3 for the 2E-VRP.

Table 3 :
Computational results comparison for Set 4 for the 2E-VRP.