Abstract
The twoechelon vehicle routing problem (2EVRP) 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 2EVRP, 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 routefirst clustersecond 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 2EVRP 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 2EVRP.
1. 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–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 twoechelon vehicle routing problem (2EVRP) [4, 5] consists in building a set of the leastcost trips (the fewest number of vehicles used and least vehicle traveling cost) for the two echelons.
The 2EVRP is a new twoechelon variant of the wellknown vehicle routing problem (VRP). Several models have been proposed and different kinds of algorithms have been designed, including both exact algorithms [4–7] and heuristic algorithms [8–11].
Perboli et al. [4] introduced a family of twoechelon vehicle routing problems and proposed a threeindex flowbased formulation for the 2EVRP. The authors also introduced some valid inequalities and two mathheuristics based on the 2EVRP model, which were used within a branchandcut 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 2EVRP 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 2EVRP and employed a specialized branching scheme to branch on infeasible integer solutions in their branchandcut 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 2EVRP would pose a challenge to incorporate.
Baldacci et al. [7] proposed a new mathematical formulation of the 2EVRP (used to derive valid lower bounds) and a new exact method. They decomposed the 2EVRP 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 stateoftheart exact methods in terms of size, number of problems solved to optimality, and computing time.
Since exact algorithms are usually computationally expensive for largescale 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 depottosatellite transfer and the satellitetocustomer 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 2EVRP; however, they obtained relatively poor results even in smallscale instances.
Meihua et al. [9] proposed a hybrid ant colony optimization algorithm which combined three heuristics for the 2EVRP. 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 2EVRP, based on the destroyandrepair 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 2EVRP. Their algorithm was shown to provide better solutions and outperform existing heuristic methods for the 2EVRP 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 2EVRP. The problem was treated by separating the depottosatellite transfer and the satellitetocustomer 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 2EVRP is the ALNS algorithm of Hemmelmayr et al. [10].
Greedy randomized adaptive search procedure (GRASP) is one of the most wellknown 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 pickupanddelivery traveling salesman problem [16], truck and trailer routing problem [17], traveling repairman problem [18], threedimensional 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 2EVRP 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 2EVRP and its mathematical model. Section 3 gives the details of our hybrid heuristic for the 2EVRP. Section 4 presents computational results on three sets of instances. Section 5 contains the discussion and Section 6 gives the conclusion.
2. Problem Description and Mathematical Formulation
The definition of 2EVRP provided here follows those in [4, 7, 10]. 2EVRP 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 , 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 secondlevel vehicles are denoted by and , respectively. The total number of vehicles available is given by for the first level and for the second. Each satellite has a capacity that limits the total amount of customers’ demands that can be delivered to it by firstlevel routes. Moreover, there are firstlevel 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 2EVRP solution with a depot, three satellites, and nine customers. The satellites are represented by fivepointed stars and customers by circles and the depot by a square. The routes that serve the satellites from the depot are called the firstlevel routes. The secondlevel 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.
2.1. Mathematical Formulation
Based on the above description, the 2EVRP can be formulated as follows [7]. Let be the index set of all firstlevel routes, and let be the subset of firstlevel routes serving satellite . Let and be the subset of satellites visited and subset of arcs traversed by a firstlevel route , respectively. Let be the index of the secondlevel routes passing through satellite , and let be the subset of routes passing through satellite and customer . The subset of customers visited and the subset of arcs traversed by a secondlevel route are represented by and , respectively. Let be a binary variable equal to 1 if and only if route is in 2EVRP solution, a binary variable equal to 1 if and only if route of satellite is in 2EVRP solution, and a nonnegative integer variable representing the quantity delivered by firstlevel route to satellite . And the mathematical model of the 2EVRP can be stated as follows [7]:
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 secondlevel route. Constraints (3), (4), and (6) impose the upper bounds on the number of first and secondlevel routes in the solution. Constraints (5) specify the capacities of each satellite. The balance between the quantity delivered by firstlevel 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 firstlevel vehicles is not exceeded.
Each secondlevel route must begin and end at the same satellite, and each customer must be served by exactly one secondlevel 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 secondlevel routes. Besides, any change to the customertosatellite assignment affects the firstlevel routing and therefore has an impact on the firstlevel transportation costs. The objective function to be minimized is the global transportation costs in both levels. 2EVRP can be easily seen to be a reduction to the vehicle routing problem (VRP), which is a special case of 2EVRP arising when just one satellite is considered, so it is also NPhard [10].
3. 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.
3.1. 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 secondlevel 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 2EVRP 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 fullload direct trips from the depot to each satellite as possible until the remaining demand of each satellite is less than the capacity of the firstlevel 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].
3.2. 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 routefirst clustersecond heuristic for the VRP. A random permutation of all customers is first cut into secondlevel 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 secondlevel routes, that is, the number of secondlevel 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 secondlevel vehicles used up to ), sets at the beginning, and tests before (the cost of the shortest path up to ) to minimize first the number of secondlevel 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 secondlevel 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 th customer in and records the predecessor of this customer on the shortest path. represents the additional penalty for each secondlevel route according to the distance between its assigned satellite and the depot. is a nonnegative value, it can be adjusted to change the impact of the firstlevel routes on the complete initial solution, and, as it increases, the firstlevel routes are given more importance. It is set to 1 as a default value to obtain a compromised consideration on routes of both echelons.

3.3. 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: intersecondlevelroute and intrasecondlevelroute. 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 2EVRP, we slightly modify the name of this classification for the secondlevel routes. The existing three interroute neighborhood structures Shift, Swap, and Cross and three intraroute neighborhood structures 2Opt, OrOpt, and Exchange of Subramanian et al. [26–28] are employed, and we also propose two new neighborhoods SatellitesChange and SatellitesSwap specially for the 2EVRP.
Following the representation of Subramanian et al. [26–28], we consider , , and , . We restrict that Shift moves at most 3 adjacent customers from a secondlevel route to another; Swap swaps at most 3 adjacent customers from a secondlevel route with 3 adjacent customers from another. As a result, nine distinct neighborhood structures can be identified, that is, , , , , , , , , and . The neighborhoods of and are shown in Figures 2 and 3, respectively. In , consecutive customers in a secondlevel route are moved from a route to another route . In , consecutive customers from a secondlevel route are interchanged with consecutive customers from another secondlevel route . The Cross operator is also called 2 in the literature, firstly proposed by Potvin and Rousseau [29], consisting of replacing arc () from a route and arc () from a route by arcs () and () (see Figure 4) or by arcs () and () (see Figure 5).
The SatellitesChange neighborhood changes the assigned satellite of a secondlevel route at a time, while SatellitesSwap neighborhood swaps the assigned satellites of two distinct secondlevel 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 2EVRP. The SatellitesChange can be seen as a special case of SatellitesSwap when a secondlevel route swaps its satellite with an empty secondlevel route. We define both of them as intersecondlevelroute neighborhoods. Neighborhoods of SatellitesSwap and SatellitesChange are shown in Figures 6 and 7, respectively.
Once a new initial feasible solution is generated, an intersecondlevelroute move is performed to modify the secondlevel routes first and then an intrasecondlevelroute improvement procedure that uses 2Opt, OrOpt, 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]. 2Opt moves delete two nonadjacent arcs and add two new arcs to generate a new route (see Figure 8); OrOpt moves at most two adjacent customers back and forth in the current secondlevel 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 secondlevel route (see Figure 10). OrOpt used in this paper can be seen as the intraroute version of and , while Exchange can be seen as the intraroute version of and .
All the neighborhoods are searched until no more improvements can be obtained, in a firstaccept 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 is actually a combination of an intersecondlevelroute neighborhood and the three intrasecondlevelroute neighborhoods. But the neighborhood placed at the end is an exception and it is only composed of the above three intrasecondlevelroute neighborhoods, which are searched circularly until none of them can improve every secondlevel route of current solution. Every time the solution is modified due to a feasible intersecondlevelroute move, the three intrasecondlevelroute operators are performed sequentially to improve the newly generated secondlevel routes.

3.4. Framework of the Resulting GRASP+VND
Our hybrid GRASP+VND heuristic for the 2EVRP 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 stands for the global best solution found.

4. 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 DualCore E5500 processor (2.8 GHz) and 2 GB of memory.
4.1. Instance Description
The proposed hybrid algorithm was tested on the 2EVRP 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: En22k4, En33k4, and En51k5. The cost matrix of each instance is given by the corresponding VRP instance. The capacity of the firstlevel vehicles is 2.5 times the capacity of the second. The capacity and the available number of the secondlevel 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].
4.2. Parameter Setting
The nonnegative weight 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.
4.3. Computational Results
In this section, we compare the computational results of our hybrid GRASP+VND heuristic with the best existing heuristic method for the 2EVRP, 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 bestknown results of these 2EVRP 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 bestknow 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/cpu_list.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.
5. Discussion
As seen from Tables 1, 2, and 3, our GRASP+VND and the ALNS can both find the bestknown solutions for all the 39 instances in Set 2 and Set 3, and the ALNS can find the bestknown 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 bestknown 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 2EVRP.
6. Conclusion
This paper has presented a very simple hybrid GRASP+VND heuristic to address the twoechelon vehicle routing problem (2EVRP), 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 2EVRP 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 2EVRP and its relative variants arising in city logistics.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by the National Natural Science Foundation of China under Grant no. 71090404. The authors are extremely grateful to the anonymous referees for their helpful comments.