#### Abstract

A framing link (FL) based tabu search algorithm is proposed in this paper for a large-scale multidepot vehicle routing problem (LSMDVRP). Framing links are generated during continuous great optimization of current solutions and then taken as skeletons so as to improve optimal seeking ability, speed up the process of optimization, and obtain better results. Based on the comparison between pre- and postmutation routes in the current solution, different parts are extracted. In the current optimization period, links involved in the optimal solution are regarded as candidates to the FL base. Multiple optimization periods exist in the whole algorithm, and there are several potential FLs in each period. If the update condition is satisfied, the FL base is updated, new FLs are added into the current route, and the next period starts. Through adjusting the borderline of multidepot sharing area with dynamic parameters, the authors define candidate selection principles for three kinds of customer connections, respectively. Link split and the roulette approach are employed to choose FLs. 18 LSMDVRP instances in three groups are studied and new optimal solution values for nine of them are obtained, with higher computation speed and reliability.

#### 1. Introduction

Nowadays, logistic cost is considered to be influential in fierce business competitions. The logistical delivery is a process from pickup to drop-off of goods, connecting vendors, shippers, and customers. In related studies, Vehicle Routing Problem (VRP) is always regarded as a significant issue, due to its impact on optimization of delivery vehicles scheduling and then on profit of logistic providers.

In common sense, Vehicle Routing Problem is defined as follows: how to determine appropriate delivery routes between series of collection and reception terminals and guarantee delivery vehicles in a proper order, so as to satisfy some requirements (e.g., shortest distance, minimum cost, delivery time, and needed vehicles) within kinds of constraints, such as amount of goods, sending time, vehicle capacity, mileage restriction, and time limitation. Currently, VRPs have been identified in many applications, such as products outbound distribution scheduling [1], home care crew scheduling [2], newspaper delivery [3], school bus routing [4], cargo routing [5], airline crew scheduling [6], waste collection scheduling [7], service system design [8], and computer system integration [9].

It is difficult to solve large-scale VRPs due to their complexities. Lenstra and Rinnooy Kan [10] proved that capacity constraint VRP (CVRP) was a NP-hard problem and recently Hassin and Rubinstein [11] verified the availability of polynomial-time algorithm in the cases where or 4 for the k-VRP issue. Imai and his partners [12] highlighted that the Vehicle Routing Problem with full container load was also NP-hard. Furthermore, Solomon [13] realized that VRP with time window constraints was more complicated, and Hashimoto et al. [14] confirmed the NP-hard characteristics of VRP with soft time windows. In the paper of Savelsbergh [15], the author proved that it is a NP-complete problem to decide whether a feasible solution existed or not for the TSP with time windows. Lenstra and Rinnooy Kan [10] proved that all types of VRPs are NP-hard problems.

From basic capacity constraint VRP, researchers have recognized varied VRP problems, for example, VRP with time window (VRPTW), periodical VRP (PVRP), VRP with pickups and deliveries (VRPPD), and multidepot VRP (MDVRP). Compared with VRP with a single depot, MDVRP with more than one depot is more complicated. In MDVRP, not only the delivery sequence but also the vehicle type and amount for customers at each distribution center need to be determined.

Given that is a complete graph, is the set of vertexes in the chart and is the set of arcs. represents vertexes, represents depots, and is the delivery amount to vertex . The capacity of vehicles from depot is , , and is the delivery cost from vertex to . In this paper, MDVRP is described as how to establish multiple paths in the chart , so as to satisfy the following:(a)each path starts and ends at the same depot;(b)all vertexes are allocated to paths, and each vertex is allocated once;(c)the total delivery amount of depot is no more than , .A heuristic approach is adapted by most researchers in dealing with MDVRP. It can be divided into the classical heuristic and the metaheuristic. Initially, the classical heuristic was more acceptable; for example, Gillett and Johnson [16] applied sweep heuristic in MDVRP and Golden et al. [17] used borderline search strategy and improved saving algorithm. A three-level heuristic algorithm was proposed by Salhi and Sari [18], where feasible solutions were generated on level 1 and delivery routes were optimized on levels 2 and 3. Sumichras and Markham [19] developed a C-W saving method. Wasner and Zpfel [20] proposed a local search strategy with a series of feedback loops, and Nagy and Sahli [21] put forward multiple enhanced optimization strategies for MDVRP with pickups and deliveries problem (MDVRPPD). Lim and Wang [22] offered two solution methodologies—one-stage and two-stage approaches—to solve MDVRP with fixed distribution of vehicles. Recently, three hybrid heuristics were proposed by Mirabi et al. [23] for MDVRP, which were based on deterministic, stochastic techniques, and simulated annealing (SA) methods. Contardo and Martinelli [24] designed a new exact method to solve the MDVRP based upon the vehicle-flow formulation and the set-partitioning formulation.

Since the 1980s, some innovative optimization methods, such as genetic algorithm, simulated annealing, tabu search, and ant colony algorithm, have been developed greatly and acted as creative roles in tackling VRP. Cordeau et al. [25] proposed a tabu search heuristic effective for three well-known routing problems: PVRP, the periodic traveling salesman problem (PTSP), and MDVRP. Such a method generalized a tabu search in solutions to VRP. This author [26] then improved the above algorithm, to solve the periodic and the multidepot vehicle routing problems with time windows (PVRPTW and MDVRPTW). Similar methods occurred in the studies of Renaud et al. [27] and Crevier et al. [28]. Renaud et al. [27] divided the algorithm into two stages and employed a tabu search to optimize the feasible solutions generated with heuristic. Crevier et al. [28] combined adaptive memory principle, a tabu search method, and the integer programming. Belhaiza et al. [29] presented a new hybrid variable neighborhood-tabu search heuristic for VRP with multiple time windows. Genetic algorithm is also acceptable by researchers. Bae et al. [30] developed an integrated VRP solver based on an improved genetic algorithm. Ho et al. [31] designed two hybrid genetic algorithms: one generated initial routes randomly and another created initial routes with Wright saving method and the nearest neighbor heuristic. In addition, Wang et al. [32] used genetic algorithm to study more complicated MDVRPs with constraints of time windows, limited numbers of vehicles, and multitype vehicles. Besides, the ant colony algorithm [33] and the simulated annealing approach [34, 35] were also applied in solving MDVRP.

However, few existing literatures pay attention to large-scale MDVRP (LSMDVRP, with customers more than 150), and the effectiveness of its solving methods should be discussed further. Framing link (FL) introduced in the following parts is helpful to reduce the searching space effectively and is a new direction in optimizing LSMDVRP. In fact, some scholars have accepted similar concepts in studying VRP; for example, Tarantilis and Kiranoudis [36] presented an adaptive memory based method for solving the Capacitated Vehicle Routing Problem (CVRP), called bone route. Zhong [37] proposed the concept and principium of kernel route, and a tabu search algorithm was designed to solve open Vehicle Routing Problem (OVRP) with capacity and distance limits. In this paper, the authors will combine the FL with a tabu search in solving LSMDVRP.

This paper is organized as follows. In Section 2, the authors propose the principle of FL for LSMDVRP and and the structure of tabu algorithm, and procedure of the method is described in Section 3. Then, the authors measure algorithm sensitivity and the relationship between FL and optimization results and then compare the result of the proposed method with those in references, which shows that the former is performed. Conclusions occur at the end.

#### 2. LSMDVRP Framing Link

##### 2.1. Principle of FL and Structure of Optimization

When an intelligent algorithm is applied in the optimization of vehicle routes, iteration is conducted based on the previous generation solutions. As a result, there are numerous links between solutions of neighbor generations, and these links are updated continuously with iterations; good links (in optimized solutions) are kept and bad ones are decomposed or combined with others. In the procedure of a so-called “good” algorithm to VRP, more good links are generated. As shown in Figures 1 and 2, the routes are achieved through iterations of a tabu search for MDVRP case p10 in Gillett and Johnson [16], and their solution values are 3714.28 and 3647.22, respectively.

The two solutions shown above occur at the generation 1233 and generation 1671, respectively, in the tabu search. Although there are over 400 generations between the two solutions, their structures have many similar parts as illustrated with red lines in Figures 1 and 2. In those routes with less generation gap, more similar links can be found. It is possible to optimize the VRP solutions based on these kinds of links as the skeleton, so as to obtain better results. In this paper, links that occur in both the optimal solution and suboptimal solutions are named as FLs.

Generation of framing links and update are keys of the FL based algorithm for LSMDVRP. The authors construct a FL base, which consists of links with higher frequencies in good routes in iterations. Addition and deletion of links are determined with the update condition. It is necessary to allocate vertexes to depots so as to generate LSMDVRP FLs, but, for those collection nodes close to several depots, possible allocations of them to different depots in optimization can result in generating unstable FLs and hence prevent them from entering the FL base. Consequently, the authors specify a FL tabu area among depots, where the generating principles of FLs are more rigorous. The basic structure of FL based tabu algorithm for LSMDVRP is described in Figure 3.

##### 2.2. Update Condition of Framing Links

If the requirement is satisfied, an update occurs in the FL base and the generated links are added into the base; meanwhile, those links without high usage frequency are deleted. In this study, the update condition of FLs includes:(a)the update of FLs has iterated times;(b)the local optimal solution keeps unchanged for generations;(c)all local and global optimums are less than after generations since the last update.

##### 2.3. Generation of New-Entry Links

###### 2.3.1. Generation of Candidate Links

FLs are generated with the continuous and large-scale optimization of current solutions in the procedure of the tabu algorithm. Compared with the premutation route, the postmutation one has different potential FLs. If these links are in the optimal solution during the current optimization period (one optimization period means the duration from updating the parameters of FLs to satisfying the next FL update condition), they are candidate links to enter the FL base. The entry procedure includes the following steps.(a)Identify the optimal vertex and those local optimal vertexes which have difference of less than with the current optimal vertex in solution value. Figure 4 shows the profile of parts of current solutions in an optimization period, where the optimal values of vertexes , , and are , , and , respectively. is the optimal vertex in this period, and and .(b)Recognize all local optimal vertexes and the local worst solution before the optimum occurs. As shown in Figure 4, the local worst solution vertexes before , , and are , , and , respectively.(c)Compare routes of local optimums and the periodical optimum with routes of their corresponding worst vertexes, and then choose different links among them. Define , , as the set of different links between the th local optimal route and its corresponding worst route in the optimization period and as the number of local optimal solutions (including the periodical optimum) in the period . For example, in Figure 4, comparing routes of and , different links between them can be found. Similarly, different links between and and between and also can be obtained. Consider .(d)Set the corresponding link set of the optimal solution in the optimization period as and the set of candidate links to the FL base as , . Different parts between and , and , and and are extracted and compared, and links in the route of are the candidate links to be added in the FL base.

###### 2.3.2. Selection of Links in the Sharing Area for the Base

In Figure 5, the area with depot as the center and (, ) as the radius is defined, and is the total collection of all depots. The sharing area of two depots is , , , as shown in blue, red, and green shadow area. For vertexes outside the sharing area of two depots, the rule of generating candidate links in Section 2.3.1 is inappropriate.

Principles to generate candidate links for the FL base in the sharing area include the following.(a)If all vertexes in a link are located in the sharing area of , , (Link Type 1), this link cannot be taken as a candidate link, as shown in Figure 6.(b)If vertexes of a link are spread in all following areas including , , , , (Link Type 2), this link cannot be taken as a candidate link, as shown in Figure 7.(c)If vertexes of a link are located in and , respectively, , (Link Type 3), this link can be a candidate link, as shown in Figure 8. If a link of Type 3 is added into the FL base, it can survive until the next th generation and then will be inspected whether to be kept in the FL base or not. Meanwhile, the borderline of the sharing area will be adjusted with the rules in Section 2.3.3.

###### 2.3.3. Adjusting Rules of Sharing Area

After the sharing area is initialized, its size and location will change with iterations, following the adhering rules.(a)Sharing area initialization: , , , is the borderline parameter of the sharing area, is the distance from depot to depot , and is the collection of all depots.(b)Adjusting rules of sharing area: since and remain constant, the area of is actually determined by and together. If there is a link of Type 3, the borderline of the sharing area will be adjusted. That means that if vertexes are located in the two areas and , , , and this link is added to the current route after generations and the current solution value is better than that before the th generation, and is the updating step of the sharing area borderline coefficient; if the link is added to the current route after generations but the current solution value is no better than that before the th generation, , as shown in Figure 9.

###### 2.3.4. Selection of Qualified Links into the FL Base

Selection of qualified links into the FL base follows several principles.(a)Amount restriction: a number of vertexes to enter the FL base are restricted into the range of ; ones out of this range are excluded.(b)Minimized split: for those routes to be added in the base, except that the complete routes are kept, they also should be decomposed into links connecting vertexes. Considering the case 3-7-4-1, , this route needs to be decomposed into six FLs: 3-7-4-1, 3-7-4, 7-4-1, 3-7, 7-4, and 4-1.(c)Backward generation: links in the FL base occur in pairs with opposite directions, so, once a link is selected into the FL base, its opposite link is generated synchronously. For example, if the link 3-7-4-1 is added, then a link 1-4-7-3 is generated.(d)Entry: based on the minimized split results of an optimized route, If this route does not exist in the FL base, it will be added in.(e)Update: once an achieved FL has more than vertexes, it should be updated, as well as its subroutes. In above case, the occurrence frequencies and corresponding optimal solution values of 3-7-4-1, 3-7-4, 7-4-1, 3-7, 7-4, and 4-1 should be updated.

##### 2.4. Deletion Principle in the Framing Link Base

(a)Comparison between links split from the current route is conducted every iterations. After generations, links before generation have been compared for times. Set the usage frequency of a link as . If , this link is deleted from the FL base. Here, is the lower limit of FL usage, as a function of comparison times .(b)If the link is included in the current optimal solution, even if its usage is less than , it cannot be deleted.

##### 2.5. Adding Framing Links into the Current Route

In the FL base, some link parameters are set: the optimal solution value of the route generated based on current links, the value of the worst solution, the average value of solutions, the usage times, the usage frequency, and the depot. Steps of adding FLs into the current route include the following.(a)Computation of the link adaptability , : is the amount of links in the base, and is the optimal solution value of the route where the th link lies in. Set as the collection of links in the FL base, and , .(b)Judgment on : if it is equal to , the procedure ends.(c)Selection of link with a roulette according to the value of , : set .(d)Calculation with , where is a vertex included in the link: if , go back to step 2.

#### 3. Framing Link Based LSMDVRP Tabu Search Algorithm

In this study, the algorithm for MDVRP consists of two parts: initial optimization and follow-up optimization. Initial optimization is based on links extracted from the FL base, and extracted link is regarded as a vertex, not to be decomposed, so as to maximize FLs’ advantages in generating optimized routes. On the other hand, follow-up optimization splits links individually, excluding non-FLs so as to avoid inferior solutions.

##### 3.1. Initial Optimization

###### 3.1.1. -Neighborhood

The -neighborhood of a vertex: The vertexes or links in the nearest collection of vertex are the -neighborhood of , denoted by . If the element in is a vertex, it is illustrated with , and if it is a link, it is described with the link in the neighborhood, where represents the start vertex of link and represents the ending vertex. Furthermore, that means that is not in the neighborhood of vertex ; otherwise, . As showed in Figure 10, if , and .

The -neighborhood of a link: vertexes or links in the nearest collection of the start vertex of link are the -neighborhood of link , represented with . At the same time, vertexes or links in the nearest collection of link are the -neighborhood of the ending vertex of link , represented with .

###### 3.1.2. Insertion Method

For MDVRP with the coexistence of vertexes and links, there are three kinds of insertions: (a) insertion between two vertexes, (b) insertion between two links, and (c) insertion between a vertex and a link. The insertion method is similar to traditional insertion method and the only difference is that this method treats the link as a node. The specific content can refer to Solomon [38].

###### 3.1.3. Generation of the Initial Solution

The generation of the initial solution includes the following steps:

*Step 1. *Allocate vertexes and links to initial depots. For a point, the nearest depot is regarded as the initial depot; for a link, the initial depot is the one where its initial route is included. represents the collection of all unallocated vertexes of the depot , represents the total delivery amount of the th route of the th depot , represents the delivery limit of vehicles in the th route of the depot , represents the delivery amount of vertex or link , and is the amount of depots.

*Step 2. *Set .

*Step 3. *Consider .

*Step 4. *If , go to Step 8.

*Step 5. *Set as the collection of all vertexes or links belonging to the depot , , and . For the th route of the depot , ---. , .

*Step 6. *For all two neighboring vertexes and links - in the route , . If , update and turn to Step 5. Otherwise, insert to form a route of --.

*Step 7. *Update , . If , turn to Step 3; otherwise, turn to Step 6.

*Step 8. *End.

###### 3.1.4. Construction of Neighborhood

According to the characteristics of FL MDVRP, the authors introduce three neighborhood operators: insertion, interchange, and crossover.

*Insertion. * is a random number in the range of . In the route , if capacity constraint of is satisfied, vertexes or links are randomly selected in , as (or ) to (or ). is set as the collection of vertexes or links outside the route , and, according Section 3.1.2, is selected.

*Interchange. * is a random number in the range of , as well as in . Choose two routes: and , and set and as the collection of vertexes and links in and . vertexes and links in are selected randomly, and then is obtained based on the method in Section 3.1.2, . Similarly, vertexes and links in are selected. If neither capacity constraint of nor that of is unsatisfied after insertions, can be inserted into the corresponding place in , or in . Otherwise, and will be selected again.

*Crossover. * and are random numbers in the range of , where , and the random number is in the range of . In the route of , a segment from to is extracted. Meanwhile, segment with a distance of is selected in a neighboring route of , where represents the distance from the center of segment to the center of the whole route . If the exchange of and to the counterpart routes cannot lead to the capacity unconstraint of these two routes, this transform is allowed; otherwise, new and need to be selected.

##### 3.2. Follow-Up Optimization

Although the FL method prompts to generate more desired solutions, FLs in the initial optimization are not necessarily parts of routes in the optimized solution; consequently, follow-up operation is needed to iterate continuously after splitting links into individual vertexes. Once the optimization solution is not yet updated generations after the initial optimization, follow-up optimization should be introduced. The initial solution of the follow-up optimization is the optimal one of the initial optimization, and if this initial solution has occurred in previous periods, the second optimal solution of the initial optimization can be selected.

#### 4. Computational Experiments

The algorithm proposed above is programmed in MS Visual C++ 6.0 and tested with a PC with a AMD Athlon (tm) X2 2.0 GHz CPU and 2 GB RAM. In this paper, the authors solve problems in existing literature with this new algorithm and then solved different LSMDVRPs with different parameters, in order to identify the coefficient valuing discipline.

The authors test the algorithm based on 18 instances from literature. Instances of p8–p11 were provided in Gillett and Johnson [16], p15–p23 by Chao et al. [39], and pr04–pr6, pr9, and pr10 by Cordeau et al. [25]. Table 1 is the column headings for Tables 2, 3, 4, 5, 6, 7, 8, 9, and 10.

The characteristics of the instances are summarized in Table 2 and the complete data sets and best known results are given in the website http://neo.lcc.uma.es/radi-aeb/WebVRP/index.html. The costs of the best solutions found by each method are listed in Table 2.

##### 4.1. Sensitivity Analyses

###### 4.1.1.

Firstly, all other parameters are set as constants. is selected from (3, 4, 5, 6, 7) and is from (2, 3, 4), given . Each of the 18 instances is calculated for 14 times, and the results are listed in Table 3.

It is not hard to find that the value of has not much influence on solutions, and average gap between solutions is 5.52%. Inferior solutions always occur in the condition where is relatively large and ; for example, . However, the value of greatly affects the computation speed of the algorithm, with an average gap of 453.88%, since can control the upper limits of lengths of candidate links, and determines their lower limits as well as split numbers. With the decrease in , the number of splits increases with an exponential distribution. According to this case, the relation of and is revealed: a larger leads to larger and . When the value of is less than 40, the commendation value of is , when , the value is or , and when is more than 100, the value is .

###### 4.1.2. Borderline Parameter and Update Step

Keeping other parameters unchanged, the authors measure the effect of here. Firstly, , is in [0.5, 0.9], the interval between successive is 0.05, and then for 18 instances are obtained through 9 times of calculations for each instance. Next, is taken varied in [0.02, 0.12], are allocated to these instances, and the interval for is defined as 0.02. for 18 instances are achieved through 6 times of calculations for each one.

According to Table 4, given is a constant, difference between the best and the worst solution values of different is 11.48%. If is regarded as a constant, difference between the best and the worst solution values of different only is 1.84%. Such results illustrate that valuing of has a more obvious impact on the computation result of the algorithm and that of . Similarly, valuing of is much more influential on the computation efficiency of the algorithm (average difference is 39.25%) than that of (average difference is 13.83%). The values of and are proportional to the customer numbers allocated in each depot (). The approximate value range of and is shown in Table 5.

###### 4.1.3.

Set , in [0.01, 0.10], and the interval between different as 0.01. The authors calculate 10 times for each instance, the number of optimization periods is , the number of local optimal vertexes in optimization period is , , is the number of candidate links produced in optimization period , and . The calculation results are shown in Table 6.

According to the analysis results, can control the number of local optimal vertexes , and the value of is inversely proportional to . This is because global optimization is difficult with the increase in and more local optimal vertexes attained. As a result, a smaller is needed to restrict the number of , so as to control . Although a larger raises the possibilities to obtain FLs, it lowers the average quality of candidate links as well, and so some “bad” links corresponding to local optimal vertexes far from the periodical vertex are involved in the candidate link base. Once these unexpected links enter the FL base, the global optimization will be staggered and its results will be weakened. In this study, a proper for computation scale and characteristics is needed. Table 7 lists the value of , , and and the accuracies of with different in the instance p09.

As shown in Table 7, when is at its best value , is the most accurate, so as to guarantee the most efficient and effective optimization. Based on this study, the recommendation value of is [0.06, 0.07] when , [0.03, 0.06] when , and [0.02, 0.03] when .

##### 4.2. FLs and Optimization Results

With the best parameters in Section 4.1, the generation process of FLs in the optimization is tracked. Set as the proportion of customers in FLs to total customers in each optimization period, as the optimal solution value of optimization period , , as the proportion of customers in FLs to total customers in each local optimal vertex, and as the value local optimal solution , . Table 8 shows numbers of FLs in every optimization period.

In Tables 8 and 9, the following can be found. (a) Average numbers of FLs increase with the progress of optimization, to the most in the last optimization period, and the optimal solution is achieved at the same time. (b) Increasing speed of FL numbers falls down gradually in the optimization process. Due to the existence of sharing area, when the proportion of customers on FLs reaches a certain number, the increase in its absolute amount will be slowed. (c) The optimal solution keeps updating with optimization. Although FLs do not increase rapidly, the optimal search ability is not challenged. In the latter periods of optimization, low quality FLs are likely to be replaced by high quality ones. (d) the average correlation coefficient between and reaches −0.924, and the number between and is −0.933 for all 18 instances. Such a result apparently bridges the proportion of customers in FLs and their corresponding local optimal solutions, so as to prove that the improvement of and prompts to achieve better optimal solutions.

##### 4.3. Optimization Results

The 18 instances are calculated for 10 times based on parameters recommended above, respectively, and the optimization results are shown in Table 10. Figure 11 describes the corresponding route of the optimal solution for p09.

In 9 of all 18 instances, new optimal solutions are obtained. The average fluctuation value of optimal solution in 18 instances is 0.97%. This fluctuation value consolidates the independence of optimization results to initial solutions in the proposed algorithm, as well as its effectiveness and robustness.

Scales of all 18 instances are large enough (), but the average computation duration is only 7.58 min with a little fluctuation rate of 18.03% in 10 calculations. All these results show the improvement of this algorithm in optimization efficiency.

For p15–p23, no better route is attained, and values of the average computation time, fluctuation rate of computation time, and optimal solutions are lower than the average one in the total 18 instances. The reason is that regular distributions and similar requirements of customers in these instances reduce the optimization space and difficulty.

#### 5. Conclusions

The authors propose a new tabu algorithm to optimize large-scale multidepot routes, using FLs as skeleton to improve search ability, speed up optimization, and then obtain better solutions. The validity of the new algorithm has been verified by example of verification.

FL is a new concept proposed in view of the VRP. However, in terms of the present study level, it is hard to get an algorithm and recommended parameters applying to all VRP constraints because of the complexity of VRP. Studies of this paper only aimed at large-scale VRP and capacity constraints. For the VRP in other conditions, this paper has certain reference value on the way of thinking; however, the specific operation of operator, particularly the recommended value parameters, needs analysis case by case. In principle, the effectiveness of FL is the result of joint action of many factors and the designs and parameter values of different operators are not unique. As the research of FL under different constraint conditions of VRP and influence factors continues, more complete processing methods and recommended parameters will gradually form.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

Data of 18 instances is from the website http://neo.lcc.uma.es/radi-aeb/WebVRP/index.html. This research is supported by National Natural Science Foundation of China (no. 71371134). Meanwhile, the authors appreciate the help and comments from editors and reviewers.