A Framing Link Based Tabu Search Algorithm for Large-Scale Multidepot Vehicle Routing Problems

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 preand 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.


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  = 3 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,  = (1, . . ., ) is the set of vertexes in the chart and  is the set of arcs.  = (1, 2, . . ., ) represents vertexes,   = ( + 1, . . .,  + ) represents depots, and   is the delivery amount to vertex .The capacity of vehicles from depot  is   ,  =  + 1, . . .,  + , 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   ,  =  + 1, . . .,  + .
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 Z ä pfel [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 approachesto 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 largescale 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.

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.28and 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. (a) Identify the optimal vertex and those local optimal vertexes which have difference of less than Δ max 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 (  −   )/  < Δ max and (  −   )/  < Δ max .

Update Condition of
(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 Ω , ,  = 1, 2, . . .,    , 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

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.    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.

Adjusting Rules of Sharing Area.
After the sharing area is initialized, its size and location will change with iterations, following the adhering rules.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. 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  min vertexes, it should be updated, as well as its subroutes.
Δd j,i  j,i d j,i i j

Postupdate Preupdate
Figure 9: Update of the public area.

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.

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   = 1/ *  ,  = 1, 2, . . .,   :   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  = { ℎ |  ℎ ∈   , there is V ∈   and V ∈  ℎ }, where V is a vertex included in the link: if   =   − , go back to step 2.

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.
The -neighborhood of a link:  vertexes or links in the nearest collection  of the start vertex V   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 V   of link , represented with    (, ).

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].

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  ℎ ,  max ( ℎ , ) represents the delivery limit of vehicles in the th route of the depot  ℎ , (V  (or   )) represents the delivery amount of vertex V  or link   , and   is the amount of depots.

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 [ min ,  max ].In the route ( ℎ , ), if capacity constraint of ( ℎ , ) is satisfied,  vertexes or links are randomly selected in , as V  (or   ) to V  (or   ). ( ℎ ,) is set as the collection of vertexes or links outside the route ( ℎ , ), and, according Section 3.1.2,V  (or   ) ∈  ( ℎ ,) is selected.
Crossover. 1 and  2 are random numbers in the range of [ min ,  max ], where  1 <  2 , and the random number  is in the range of [ min ,  max ].In the route of ( ℎ ,  1 ), a segment  1 from  1 to  2 is extracted.Meanwhile, segment  * = { | min (  1 ,   )} with a distance of  is selected in a neighboring route of ( ℎ ,  1 ), where (  1 ,   ) represents the distance from the center of segment  1 to the center of the whole route .If the exchange of  1 and  * to the counterpart routes cannot lead to the capacity unconstraint of these two routes, this transform is allowed; otherwise, new  1 and  2 need to be selected.

Follow-Up Optimization.
Although the FL method prompts to generate more desired solutions, FLs in the initial The proportion of customers in FLs to total customers in local optimal vertex      1   The optimal solution value of optimization period    2 The value local optimal solution  Once the optimization solution is not yet updated  generations after the initial optimization, follow-up optimization should be introduced.The initial solution of the followup 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.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.
It is not hard to find that the value of [ min ,  max ] has not much influence on solutions, and average gap between  [2,4], when 40 < / ≤ 100, the value is [3,5] or [3,6], and when / is more than 100, the value is [3,7].

Borderline Parameter 𝛿 𝑖,𝑗 and Update
Step Δ.Keeping other parameters unchanged, the authors measure the effect of  , here.Firstly, Δ = 0.04,  , 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 =1 Ω   .The calculation results are shown in Table 6.
According to the analysis results, Δ max can control the number of local optimal vertexes   , and the value of Δ max 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 Δ max 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 Δ max for computation scale and characteristics is needed.as the proportion of customers in FLs to total customers in each local optimal vertex, and   2 as the value local optimal solution ,  = 1, 2, . . .,   .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.% CPU time: the gap in percentage between the average CPU time and the fastest CPU time among the ten runs.% solution value: the gap in percentage between the average value of the solutions and the best solution value among the ten runs.
In the latter periods of optimization, low quality FLs are likely to be replaced by high quality ones.(d) the average correlation coefficient between   1 and   1 reaches −0.924, and the number between   2 and   2 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   1 and   2 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 (150 <  ≤ 360), 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.

Conclusions
The authors propose a new tabu algorithm to optimize largescale multidepot routes, using FLs as skeleton to improve search ability, speed up optimization, and then obtain better

Figure 1 :
Figure 1: The route with a solution value 3731.38.

Figure 2 :
Figure 2: The route with a solution value 3647.22.

Figure 3 :Figure 4 :
Figure 3: The basic structure of FL based tabu algorithm for LSMDVRP.
(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 Figure7.(c)If vertexes of a link are located in   and   ∩   , respectively, ,  ∈   ,  ̸ =  (Link Type 3), this link can be a candidate link, as shown in Figure8.If a link R

Figure 5 :
Figure 5: Illustration of the sharing area.

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 [ min ,  max ]; 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  min vertexes.Considering the case 3-7-4-1,  min = 2, 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

Figure 11 :
Figure 11: The corresponding route of the optimal solution for p09.

Table 1 :
Column headings for Tables2-10.The proportion of customers in FLs to total customers in optimization period

Table 2 :
MDVRP instances and previous best known solution cost.

Table 3 :
Sensitivity data of parameters  min and  max .

Table 4 :
Sensitivity data of parameter  , and its update step Δ.Gap between the worst and the best solution value with different  , and Δ = 0.04.%Gap 2: % Gap between the worst and the best solution value with different Δ and  , * .%Gap 3: % Gap between the longest and the shortest CPU time with different  , and Δ = 0.04.% Gap 4: % Gap between the longest and the shortest CPU time with different Δ and  , % Gap 1: % * .

Table 6 :
The calculation results of parameter Δ max .

Table 7 :
[5,5]alue of   ,   , and Ω  and the accuracies of Ω  with different Δ max in the instance p09.Inferior solutions always occur in the condition where  min is relatively large and  min =  max ; for example, [ min ,  max ] =[5,5].However, the value of [ min ,  max ] greatly affects the computation speed of the algorithm, with an average gap of 453.88%, since  max can control the upper limits of lengths of candidate links, and  min determines their lower limits as well as split numbers.With the decrease in  min , the number of splits increases with an exponential distribution.According to this case, the relation of [ * min ,  * max ] and / is revealed: a larger / leads to larger  * min and  * max .When the value of / is less than 40, the commendation value of [ * min ,  * max ] is

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 Table5.4.1.3.Δ max .Set [ min ,  max ] = [ * min ,  * max ],Δ max in [0.01, 0.10], and the interval between different Δ max 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    ,   = ∑ Table 7 lists the value of   ,   , and Ω  and the accuracies of Ω  with different Δ max in the instance p09.As shown in Table 7, when Δ max is at its best value Δ * max , Ω  is the most accurate, so as to guarantee the most efficient and effective optimization.Based on this study, the recommendation value of Δ * max is [0.06, 0.07] when 150 <  ≤ 200, [0.03, 0.06] when 200 <  ≤ 300, and [0.02, 0.03] when 300 <  ≤ 400.

Table 8 :
Numbers of FLs in every optimization period.

Table 9 :
The correlation of   1 , 1 ,  With the best parameters in Section 4.1, the generation process of FLs in the optimization is tracked.Set   1 as the proportion of customers in FLs to total customers in each optimization period,   1 as the optimal solution value of optimization period ,  = 1, 2, . . .,   ,

Table 10 :
The results found by proposed tabu search algorithm on FL.