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

Authors: Xuhao Zhang; Shiquan Zhong; Yiliu Liu; Xuelian Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

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


---

## Body

## 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 k
=
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 thatG
=
(
V
,
A
) is a complete graph, V
=
(
1
,
…
,
n
) is the set of vertexes in the chart and A is the set of arcs. V
c
=
(
1,2
,
…
,
n
) represents vertexes, V
d
=
(
n
+
1
,
…
,
n
+
m
) represents depots, and q
i is the delivery amount to vertex i. The capacity of vehicles from depot j is W
j, j
=
n
+
1
,
…
,
n
+
m, and c
i
j is the delivery cost from vertex i to j. In this paper, MDVRP is described as how to establish multiple paths in the chart G, 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 depotj is no more than W
j, j
=
n
+
1
,
…
,
n
+
m.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 Za
¨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 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 Section2, 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 Figures1 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.Figure 1
The route with a solution value 3731.38.Figure 2
The route with a solution value 3647.22.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 Figures1 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 Figure3.Figure 3
The basic structure of FL based tabu algorithm for LSMDVRP.


### 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 iteratedN
d times;
(b)
the local optimal solution keeps unchanged forN
i generations;
(c)
all local and global optimums are less thanB
q after N
q 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Δ
s
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 A, B, and C are s
A, s
B, and s
C, respectively. B is the optimal vertex in this period, and (
s
A
-
s
B
)
/
s
B
<
Δ
s
max
⁡ and (
s
C
-
s
B
)
/
s
B
<
Δ
s
max
⁡.
(b)
Recognize all local optimal vertexes and the local worst solution before the optimum occurs. As shown in Figure4, the local worst solution vertexes before A, B, and C are A
′, B
′, and C
′, 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Ω
k
,
i, i
=
1,2
,
…
,
n
L
k, as the set of different links between the ith local optimal route and its corresponding worst route in the optimization period k and n
L
k as the number of local optimal solutions (including the periodical optimum) in the period k. For example, in Figure 4, comparing routes of A and A
′, different links between them can be found. Similarly, different links between B and B
′ and between C and C
′ also can be obtained. Consider Ω
k
=
{
Ω
k
,
1
,
Ω
k
,
2
,
…
,
Ω
k
,
n
L
k
}.
(d)
Set the corresponding link set of the optimal solution in the optimization periodk as Ω
B
k and the set of candidate links to the FL base as Ω
k
′, Ω
k
′
=
{
l
∣
l
∈
Ω
k
,
l
∈
Ω
B
k
}. Different parts between A and A
′, B and B
′, and C and C
′ are extracted and compared, and links in the route of B are the candidate links to be added in the FL base.Figure 4
Generation of framing links.


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

In Figure5, the area A with depot i as the center and R
i
,
j (i
,
j
∈
V
d, i
≠
j) as the radius is defined, and V
d is the total collection of all depots. The sharing area of two depots is A
i
∩
A
j, i
,
j
∈
V, i
≠
j, 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.Figure 5
Illustration of the sharing area.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 ofA
i
∩
A
j, i
,
j
∈
V
d, i
≠
j (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 includingA
i
∩
A
j, A
i, A
j, i
,
j
∈
V
d, i
≠
j (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 inA
i and A
i
∩
A
j, respectively, i
,
j
∈
V
d, i
≠
j (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 N
cth 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.Figure 6
Link Type 1.Figure 7
Link Type 2.Figure 8
Link Type 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:R
i
,
j
=
δ
i
,
j
d
i
,
j, i
,
j
∈
V
d, i
≠
j, δ
i
,
j is the borderline parameter of the sharing area, d
i
,
j is the distance from depot i to depot j, and V
d is the collection of all depots.
(b)
Adjusting rules of sharing area: sinced
i
,
j and d
j
,
i remain constant, the area of A
i
∩
A
j is actually determined by δ
i
,
j and δ
j
,
i 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 A
i and A
i
∩
A
j, i
,
j
∈
V
d, i
≠
j, and this link is added to the current route after N
c generations and the current solution value is better than that before the N
cth generation, δ
j
,
i
∶
=
δ
j
,
i
+
Δ
δ and Δ
δ is the updating step of the sharing area borderline coefficient; if the link is added to the current route after N
c generations but the current solution value is no better than that before the N
cth generation, δ
j
,
i
∶
=
δ
j
,
i
-
Δ
δ, as shown in Figure 9.Figure 9
Update of the public 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[
L
min
⁡
,
L
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 connectingL
min
⁡ vertexes. Considering the case 3-7-4-1, L
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 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 thanL
min
⁡ 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 k generations, links before generation d have been compared for c
=
[
k
/
η
] times. Set the usage frequency of a link as u
(
c
). If u
(
c
)
<
u
s
(
c
), this link is deleted from the FL base. Here, u
s
(
c
) is the lower limit of FL usage, as a function of comparison times c.
(b)
If the link is included in the current optimal solution, even if its usage is less thanu
s
(
c
), 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 adaptabilityf
i
=
1
/
s
i
*, i
=
1,2
,
…
,
N
l: N
i is the amount of links in the base, and s
i
* is the optimal solution value of the route where the ith link lies in. Set H as the collection of links in the FL base, and H
′
=
H, H
*
=
Φ.
(b)
Judgment onH
′: if it is equal to Φ, the procedure ends.
(c)
Selection of linkl
j with a roulette according to the value of f
i, l
j
∈
H
′: set H
*
=
H
*
+
{
l
j
}.
(d)
Calculation withL
=
{
l
h
∣
l
h
∈
H
′
,
there
is
v
∈
l
j
and
v
∈
l
h
}, where v is a vertex included in the link: if H
′
=
H
′
-
L, go back to step 2.


## 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 Figures1 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.Figure 1
The route with a solution value 3731.38.Figure 2
The route with a solution value 3647.22.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 Figures1 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 Figure3.Figure 3
The basic structure of FL based tabu algorithm for LSMDVRP.


## 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 iteratedN
d times;
(b)
the local optimal solution keeps unchanged forN
i generations;
(c)
all local and global optimums are less thanB
q after N
q 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Δ
s
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 A, B, and C are s
A, s
B, and s
C, respectively. B is the optimal vertex in this period, and (
s
A
-
s
B
)
/
s
B
<
Δ
s
max
⁡ and (
s
C
-
s
B
)
/
s
B
<
Δ
s
max
⁡.
(b)
Recognize all local optimal vertexes and the local worst solution before the optimum occurs. As shown in Figure4, the local worst solution vertexes before A, B, and C are A
′, B
′, and C
′, 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Ω
k
,
i, i
=
1,2
,
…
,
n
L
k, as the set of different links between the ith local optimal route and its corresponding worst route in the optimization period k and n
L
k as the number of local optimal solutions (including the periodical optimum) in the period k. For example, in Figure 4, comparing routes of A and A
′, different links between them can be found. Similarly, different links between B and B
′ and between C and C
′ also can be obtained. Consider Ω
k
=
{
Ω
k
,
1
,
Ω
k
,
2
,
…
,
Ω
k
,
n
L
k
}.
(d)
Set the corresponding link set of the optimal solution in the optimization periodk as Ω
B
k and the set of candidate links to the FL base as Ω
k
′, Ω
k
′
=
{
l
∣
l
∈
Ω
k
,
l
∈
Ω
B
k
}. Different parts between A and A
′, B and B
′, and C and C
′ are extracted and compared, and links in the route of B are the candidate links to be added in the FL base.Figure 4
Generation of framing links.


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

In Figure5, the area A with depot i as the center and R
i
,
j (i
,
j
∈
V
d, i
≠
j) as the radius is defined, and V
d is the total collection of all depots. The sharing area of two depots is A
i
∩
A
j, i
,
j
∈
V, i
≠
j, 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.Figure 5
Illustration of the sharing area.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 ofA
i
∩
A
j, i
,
j
∈
V
d, i
≠
j (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 includingA
i
∩
A
j, A
i, A
j, i
,
j
∈
V
d, i
≠
j (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 inA
i and A
i
∩
A
j, respectively, i
,
j
∈
V
d, i
≠
j (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 N
cth 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.Figure 6
Link Type 1.Figure 7
Link Type 2.Figure 8
Link Type 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:R
i
,
j
=
δ
i
,
j
d
i
,
j, i
,
j
∈
V
d, i
≠
j, δ
i
,
j is the borderline parameter of the sharing area, d
i
,
j is the distance from depot i to depot j, and V
d is the collection of all depots.
(b)
Adjusting rules of sharing area: sinced
i
,
j and d
j
,
i remain constant, the area of A
i
∩
A
j is actually determined by δ
i
,
j and δ
j
,
i 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 A
i and A
i
∩
A
j, i
,
j
∈
V
d, i
≠
j, and this link is added to the current route after N
c generations and the current solution value is better than that before the N
cth generation, δ
j
,
i
∶
=
δ
j
,
i
+
Δ
δ and Δ
δ is the updating step of the sharing area borderline coefficient; if the link is added to the current route after N
c generations but the current solution value is no better than that before the N
cth generation, δ
j
,
i
∶
=
δ
j
,
i
-
Δ
δ, as shown in Figure 9.Figure 9
Update of the public 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[
L
min
⁡
,
L
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 connectingL
min
⁡ vertexes. Considering the case 3-7-4-1, L
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 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 thanL
min
⁡ 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.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Δ
s
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 A, B, and C are s
A, s
B, and s
C, respectively. B is the optimal vertex in this period, and (
s
A
-
s
B
)
/
s
B
<
Δ
s
max
⁡ and (
s
C
-
s
B
)
/
s
B
<
Δ
s
max
⁡.
(b)
Recognize all local optimal vertexes and the local worst solution before the optimum occurs. As shown in Figure4, the local worst solution vertexes before A, B, and C are A
′, B
′, and C
′, 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Ω
k
,
i, i
=
1,2
,
…
,
n
L
k, as the set of different links between the ith local optimal route and its corresponding worst route in the optimization period k and n
L
k as the number of local optimal solutions (including the periodical optimum) in the period k. For example, in Figure 4, comparing routes of A and A
′, different links between them can be found. Similarly, different links between B and B
′ and between C and C
′ also can be obtained. Consider Ω
k
=
{
Ω
k
,
1
,
Ω
k
,
2
,
…
,
Ω
k
,
n
L
k
}.
(d)
Set the corresponding link set of the optimal solution in the optimization periodk as Ω
B
k and the set of candidate links to the FL base as Ω
k
′, Ω
k
′
=
{
l
∣
l
∈
Ω
k
,
l
∈
Ω
B
k
}. Different parts between A and A
′, B and B
′, and C and C
′ are extracted and compared, and links in the route of B are the candidate links to be added in the FL base.Figure 4
Generation of framing links.


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

In Figure5, the area A with depot i as the center and R
i
,
j (i
,
j
∈
V
d, i
≠
j) as the radius is defined, and V
d is the total collection of all depots. The sharing area of two depots is A
i
∩
A
j, i
,
j
∈
V, i
≠
j, 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.Figure 5
Illustration of the sharing area.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 ofA
i
∩
A
j, i
,
j
∈
V
d, i
≠
j (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 includingA
i
∩
A
j, A
i, A
j, i
,
j
∈
V
d, i
≠
j (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 inA
i and A
i
∩
A
j, respectively, i
,
j
∈
V
d, i
≠
j (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 N
cth 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.Figure 6
Link Type 1.Figure 7
Link Type 2.Figure 8
Link Type 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:R
i
,
j
=
δ
i
,
j
d
i
,
j, i
,
j
∈
V
d, i
≠
j, δ
i
,
j is the borderline parameter of the sharing area, d
i
,
j is the distance from depot i to depot j, and V
d is the collection of all depots.
(b)
Adjusting rules of sharing area: sinced
i
,
j and d
j
,
i remain constant, the area of A
i
∩
A
j is actually determined by δ
i
,
j and δ
j
,
i 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 A
i and A
i
∩
A
j, i
,
j
∈
V
d, i
≠
j, and this link is added to the current route after N
c generations and the current solution value is better than that before the N
cth generation, δ
j
,
i
∶
=
δ
j
,
i
+
Δ
δ and Δ
δ is the updating step of the sharing area borderline coefficient; if the link is added to the current route after N
c generations but the current solution value is no better than that before the N
cth generation, δ
j
,
i
∶
=
δ
j
,
i
-
Δ
δ, as shown in Figure 9.Figure 9
Update of the public 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[
L
min
⁡
,
L
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 connectingL
min
⁡ vertexes. Considering the case 3-7-4-1, L
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 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 thanL
min
⁡ 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 k generations, links before generation d have been compared for c
=
[
k
/
η
] times. Set the usage frequency of a link as u
(
c
). If u
(
c
)
<
u
s
(
c
), this link is deleted from the FL base. Here, u
s
(
c
) is the lower limit of FL usage, as a function of comparison times c.
(b)
If the link is included in the current optimal solution, even if its usage is less thanu
s
(
c
), 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 adaptabilityf
i
=
1
/
s
i
*, i
=
1,2
,
…
,
N
l: N
i is the amount of links in the base, and s
i
* is the optimal solution value of the route where the ith link lies in. Set H as the collection of links in the FL base, and H
′
=
H, H
*
=
Φ.
(b)
Judgment onH
′: if it is equal to Φ, the procedure ends.
(c)
Selection of linkl
j with a roulette according to the value of f
i, l
j
∈
H
′: set H
*
=
H
*
+
{
l
j
}.
(d)
Calculation withL
=
{
l
h
∣
l
h
∈
H
′
,
there
is
v
∈
l
j
and
v
∈
l
h
}, where v is a vertex included in the link: if H
′
=
H
′
-
L, 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.p-Neighborhood

Thep-neighborhood of a vertex: The p vertexes or links in the nearest collection A of vertex v are the p-neighborhood of v, denoted by N
p
(
v
,
A
). If the element in N
p
(
v
,
A
) is a vertex, it is illustrated with v
′, and if it is a link, it is described with the link v
l
s
(
z
l
,
v
s
)
-
v
l
2
(
z
l
,
v
2
)
-
v
l
3
(
z
l
,
v
3
)
-
⋯
-
v
l
e
(
z
l
,
v
e
) in the neighborhood, where v
l
s represents the start vertex of link l and v
l
e represents the ending vertex. Furthermore, that z
l
,
v
s
=
0 means that v
l
s is not in the p neighborhood of vertex v; otherwise, z
l
,
v
s
=
1. As showed in Figure 10, if p
=
2, N
2
(
4
)
=
{
8,3
(
1
)
-
2
(
1
)
-
1
(
0
)
} and N
2
(
7
)
=
{
9
(
1
)
-
10
(
1
)
-
11
(
0
)
,
12
(
1
)
-
13
(
0
)
}.Figure 10
p-neighborhoods of vertexes and links.Thep-neighborhood of a link: p vertexes or links in the nearest collection A of the start vertex v
l
s of link l are the p-neighborhood of link l, represented with N
p
s
(
l
,
A
). At the same time, p vertexes or links in the nearest collection A of link l are the p-neighborhood of the ending vertex v
l
e of link l, represented with N
p
e
(
l
,
A
).


#### 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.N
r
(
c
h
) represents the collection of all unallocated vertexes of the depot h, G
(
c
h
,
r
) represents the total delivery amount of the rth route of the hth depot c
h, G
max
⁡
(
c
h
,
r
) represents the delivery limit of vehicles in the rth route of the depot c
h, g
(
v
k
(
or
l
k
)
) represents the delivery amount of vertex v
k or link l
k, and N
C is the amount of depots.Step 2.
Seth
=
0.Step 3.
Considerh
:
=
h
+
1.Step 4.
Ifh
>
N
C, go to Step 8.Step 5.
SetA
c
h as the collection of all vertexes or links belonging to the depot c
h, ∀
v
i
(
or
l
i
)
∈
N
p
(
c
h
,
A
c
h
), and v
j
(
or
l
j
)
∈
N
p
(
c
h
,
A
c
h
). For the rth route R
(
c
h
,
r
) of the depot c
h, c
h-v
i
(
or
l
i
)-v
j
(
or
l
j
)-c
h. G
(
c
h
,
r
)
=
g
(
v
i
(
or
l
i
)
)
+
g
(
v
j
(
or
l
j
)
), N
r
(
c
h
)
:
=
N
r
(
c
h
)
-
v
i
(
or
l
i
)
-
v
j
(
or
l
j
).Step 6.
For all two neighboring vertexes and linksv
i
(
or
l
i
)-v
j
(
or
l
j
) in the route r, v
k
(
or
l
k
)
=
min
⁡
{
v
(
or
l
)
∣
Δ
d
(
v
i
(
or
l
i
)
,
v
(
or
l
)
,
v
j
(
or
l
j
)
). If G
(
c
h
,
r
)
+
g
(
v
k
(
or
l
k
)
)
>
G
max
⁡
(
h
,
r
), update r
:
=
r
+
1 and turn to Step 5. Otherwise, insert v
k
(
or
l
k
) to form a route of v
i
(
or
l
i
)-v
k
(
or
l
k
)-v
j
(
or
l
j
).Step 7.
UpdateG
(
c
h
,
r
)
:
=
G
(
c
h
,
r
)
+
g
(
v
k
(
or
l
k
)
), N
r
(
c
h
)
:
=
N
r
(
c
h
)
-
v
k
(
or
l
k
). If N
(
c
h
)
=
Φ, 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 [
α
min
⁡
,
α
max
⁡
]. In the route R
(
c
h
,
r
), if capacity constraint of R
(
c
h
,
r
) is satisfied, α vertexes or links are randomly selected in r, as v
i (or l
i) to v
j (or l
j). A
R
-
(
c
h
,
r
) is set as the collection of vertexes or links outside the route R
(
c
h
,
r
), and, according Section 3.1.2, v
k
(
or
l
k
)
∈
A
R
-
(
c
h
,
r
) is selected.Interchange.
θ is a random number in the range of [
θ
min
⁡
,
θ
max
⁡
], as well as φ in [
φ
min
⁡
,
φ
max
⁡
]. Choose two routes: R
(
c
h
,
r
1
) and R
(
c
h
,
r
2
), and set A
R
(
c
h
,
r
1
) and A
R
(
c
h
,
r
2
) as the collection of vertexes and links in R
(
c
h
,
r
1
) and R
(
c
h
,
r
2
). θ vertexes and links in R
(
c
h
,
r
1
) are selected randomly, and then v
m
(
or
l
m
)
∈
A
R
(
c
h
,
r
2
) is obtained based on the method in Section 3.1.2, m
=
1,2
,
…
,
θ. Similarly, φ vertexes and links in R
(
c
h
,
r
2
) are selected. If neither capacity constraint of R
(
c
h
,
r
1
) nor that of R
(
c
h
,
r
2
) is unsatisfied after insertions, v
m
(
or
l
m
) can be inserted into the corresponding place in A
R
(
c
h
,
r
1
), or v
n
(
or
l
n
) in A
R
(
c
h
,
r
2
). Otherwise, v
m
(
or
l
m
) and v
n
(
or
l
n
) will be selected again.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 R
(
c
h
,
r
1
), a segment s
1 from γ
1 to γ
2 is extracted. Meanwhile, segment s
*
=
{
s
∣
min
⁡
d
(
c
s
1
,
c
s
)
} with a distance of λ is selected in a neighboring route of R
(
c
h
,
r
1
), where d
(
c
s
1
,
c
s
) represents the distance from the center of segment s
1 to the center of the whole route s. If the exchange of s
1 and s
* 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.


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


## 3.1. Initial Optimization

### 3.1.1.p-Neighborhood

Thep-neighborhood of a vertex: The p vertexes or links in the nearest collection A of vertex v are the p-neighborhood of v, denoted by N
p
(
v
,
A
). If the element in N
p
(
v
,
A
) is a vertex, it is illustrated with v
′, and if it is a link, it is described with the link v
l
s
(
z
l
,
v
s
)
-
v
l
2
(
z
l
,
v
2
)
-
v
l
3
(
z
l
,
v
3
)
-
⋯
-
v
l
e
(
z
l
,
v
e
) in the neighborhood, where v
l
s represents the start vertex of link l and v
l
e represents the ending vertex. Furthermore, that z
l
,
v
s
=
0 means that v
l
s is not in the p neighborhood of vertex v; otherwise, z
l
,
v
s
=
1. As showed in Figure 10, if p
=
2, N
2
(
4
)
=
{
8,3
(
1
)
-
2
(
1
)
-
1
(
0
)
} and N
2
(
7
)
=
{
9
(
1
)
-
10
(
1
)
-
11
(
0
)
,
12
(
1
)
-
13
(
0
)
}.Figure 10
p-neighborhoods of vertexes and links.Thep-neighborhood of a link: p vertexes or links in the nearest collection A of the start vertex v
l
s of link l are the p-neighborhood of link l, represented with N
p
s
(
l
,
A
). At the same time, p vertexes or links in the nearest collection A of link l are the p-neighborhood of the ending vertex v
l
e of link l, represented with N
p
e
(
l
,
A
).


### 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.N
r
(
c
h
) represents the collection of all unallocated vertexes of the depot h, G
(
c
h
,
r
) represents the total delivery amount of the rth route of the hth depot c
h, G
max
⁡
(
c
h
,
r
) represents the delivery limit of vehicles in the rth route of the depot c
h, g
(
v
k
(
or
l
k
)
) represents the delivery amount of vertex v
k or link l
k, and N
C is the amount of depots.Step 2.
Seth
=
0.Step 3.
Considerh
:
=
h
+
1.Step 4.
Ifh
>
N
C, go to Step 8.Step 5.
SetA
c
h as the collection of all vertexes or links belonging to the depot c
h, ∀
v
i
(
or
l
i
)
∈
N
p
(
c
h
,
A
c
h
), and v
j
(
or
l
j
)
∈
N
p
(
c
h
,
A
c
h
). For the rth route R
(
c
h
,
r
) of the depot c
h, c
h-v
i
(
or
l
i
)-v
j
(
or
l
j
)-c
h. G
(
c
h
,
r
)
=
g
(
v
i
(
or
l
i
)
)
+
g
(
v
j
(
or
l
j
)
), N
r
(
c
h
)
:
=
N
r
(
c
h
)
-
v
i
(
or
l
i
)
-
v
j
(
or
l
j
).Step 6.
For all two neighboring vertexes and linksv
i
(
or
l
i
)-v
j
(
or
l
j
) in the route r, v
k
(
or
l
k
)
=
min
⁡
{
v
(
or
l
)
∣
Δ
d
(
v
i
(
or
l
i
)
,
v
(
or
l
)
,
v
j
(
or
l
j
)
). If G
(
c
h
,
r
)
+
g
(
v
k
(
or
l
k
)
)
>
G
max
⁡
(
h
,
r
), update r
:
=
r
+
1 and turn to Step 5. Otherwise, insert v
k
(
or
l
k
) to form a route of v
i
(
or
l
i
)-v
k
(
or
l
k
)-v
j
(
or
l
j
).Step 7.
UpdateG
(
c
h
,
r
)
:
=
G
(
c
h
,
r
)
+
g
(
v
k
(
or
l
k
)
), N
r
(
c
h
)
:
=
N
r
(
c
h
)
-
v
k
(
or
l
k
). If N
(
c
h
)
=
Φ, 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 [
α
min
⁡
,
α
max
⁡
]. In the route R
(
c
h
,
r
), if capacity constraint of R
(
c
h
,
r
) is satisfied, α vertexes or links are randomly selected in r, as v
i (or l
i) to v
j (or l
j). A
R
-
(
c
h
,
r
) is set as the collection of vertexes or links outside the route R
(
c
h
,
r
), and, according Section 3.1.2, v
k
(
or
l
k
)
∈
A
R
-
(
c
h
,
r
) is selected.Interchange.
θ is a random number in the range of [
θ
min
⁡
,
θ
max
⁡
], as well as φ in [
φ
min
⁡
,
φ
max
⁡
]. Choose two routes: R
(
c
h
,
r
1
) and R
(
c
h
,
r
2
), and set A
R
(
c
h
,
r
1
) and A
R
(
c
h
,
r
2
) as the collection of vertexes and links in R
(
c
h
,
r
1
) and R
(
c
h
,
r
2
). θ vertexes and links in R
(
c
h
,
r
1
) are selected randomly, and then v
m
(
or
l
m
)
∈
A
R
(
c
h
,
r
2
) is obtained based on the method in Section 3.1.2, m
=
1,2
,
…
,
θ. Similarly, φ vertexes and links in R
(
c
h
,
r
2
) are selected. If neither capacity constraint of R
(
c
h
,
r
1
) nor that of R
(
c
h
,
r
2
) is unsatisfied after insertions, v
m
(
or
l
m
) can be inserted into the corresponding place in A
R
(
c
h
,
r
1
), or v
n
(
or
l
n
) in A
R
(
c
h
,
r
2
). Otherwise, v
m
(
or
l
m
) and v
n
(
or
l
n
) will be selected again.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 R
(
c
h
,
r
1
), a segment s
1 from γ
1 to γ
2 is extracted. Meanwhile, segment s
*
=
{
s
∣
min
⁡
d
(
c
s
1
,
c
s
)
} with a distance of λ is selected in a neighboring route of R
(
c
h
,
r
1
), where d
(
c
s
1
,
c
s
) represents the distance from the center of segment s
1 to the center of the whole route s. If the exchange of s
1 and s
* 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.


## 3.1.1.p-Neighborhood

Thep-neighborhood of a vertex: The p vertexes or links in the nearest collection A of vertex v are the p-neighborhood of v, denoted by N
p
(
v
,
A
). If the element in N
p
(
v
,
A
) is a vertex, it is illustrated with v
′, and if it is a link, it is described with the link v
l
s
(
z
l
,
v
s
)
-
v
l
2
(
z
l
,
v
2
)
-
v
l
3
(
z
l
,
v
3
)
-
⋯
-
v
l
e
(
z
l
,
v
e
) in the neighborhood, where v
l
s represents the start vertex of link l and v
l
e represents the ending vertex. Furthermore, that z
l
,
v
s
=
0 means that v
l
s is not in the p neighborhood of vertex v; otherwise, z
l
,
v
s
=
1. As showed in Figure 10, if p
=
2, N
2
(
4
)
=
{
8,3
(
1
)
-
2
(
1
)
-
1
(
0
)
} and N
2
(
7
)
=
{
9
(
1
)
-
10
(
1
)
-
11
(
0
)
,
12
(
1
)
-
13
(
0
)
}.Figure 10
p-neighborhoods of vertexes and links.Thep-neighborhood of a link: p vertexes or links in the nearest collection A of the start vertex v
l
s of link l are the p-neighborhood of link l, represented with N
p
s
(
l
,
A
). At the same time, p vertexes or links in the nearest collection A of link l are the p-neighborhood of the ending vertex v
l
e of link l, represented with N
p
e
(
l
,
A
).


## 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.N
r
(
c
h
) represents the collection of all unallocated vertexes of the depot h, G
(
c
h
,
r
) represents the total delivery amount of the rth route of the hth depot c
h, G
max
⁡
(
c
h
,
r
) represents the delivery limit of vehicles in the rth route of the depot c
h, g
(
v
k
(
or
l
k
)
) represents the delivery amount of vertex v
k or link l
k, and N
C is the amount of depots.Step 2.
Seth
=
0.Step 3.
Considerh
:
=
h
+
1.Step 4.
Ifh
>
N
C, go to Step 8.Step 5.
SetA
c
h as the collection of all vertexes or links belonging to the depot c
h, ∀
v
i
(
or
l
i
)
∈
N
p
(
c
h
,
A
c
h
), and v
j
(
or
l
j
)
∈
N
p
(
c
h
,
A
c
h
). For the rth route R
(
c
h
,
r
) of the depot c
h, c
h-v
i
(
or
l
i
)-v
j
(
or
l
j
)-c
h. G
(
c
h
,
r
)
=
g
(
v
i
(
or
l
i
)
)
+
g
(
v
j
(
or
l
j
)
), N
r
(
c
h
)
:
=
N
r
(
c
h
)
-
v
i
(
or
l
i
)
-
v
j
(
or
l
j
).Step 6.
For all two neighboring vertexes and linksv
i
(
or
l
i
)-v
j
(
or
l
j
) in the route r, v
k
(
or
l
k
)
=
min
⁡
{
v
(
or
l
)
∣
Δ
d
(
v
i
(
or
l
i
)
,
v
(
or
l
)
,
v
j
(
or
l
j
)
). If G
(
c
h
,
r
)
+
g
(
v
k
(
or
l
k
)
)
>
G
max
⁡
(
h
,
r
), update r
:
=
r
+
1 and turn to Step 5. Otherwise, insert v
k
(
or
l
k
) to form a route of v
i
(
or
l
i
)-v
k
(
or
l
k
)-v
j
(
or
l
j
).Step 7.
UpdateG
(
c
h
,
r
)
:
=
G
(
c
h
,
r
)
+
g
(
v
k
(
or
l
k
)
), N
r
(
c
h
)
:
=
N
r
(
c
h
)
-
v
k
(
or
l
k
). If N
(
c
h
)
=
Φ, 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 [
α
min
⁡
,
α
max
⁡
]. In the route R
(
c
h
,
r
), if capacity constraint of R
(
c
h
,
r
) is satisfied, α vertexes or links are randomly selected in r, as v
i (or l
i) to v
j (or l
j). A
R
-
(
c
h
,
r
) is set as the collection of vertexes or links outside the route R
(
c
h
,
r
), and, according Section 3.1.2, v
k
(
or
l
k
)
∈
A
R
-
(
c
h
,
r
) is selected.Interchange.
θ is a random number in the range of [
θ
min
⁡
,
θ
max
⁡
], as well as φ in [
φ
min
⁡
,
φ
max
⁡
]. Choose two routes: R
(
c
h
,
r
1
) and R
(
c
h
,
r
2
), and set A
R
(
c
h
,
r
1
) and A
R
(
c
h
,
r
2
) as the collection of vertexes and links in R
(
c
h
,
r
1
) and R
(
c
h
,
r
2
). θ vertexes and links in R
(
c
h
,
r
1
) are selected randomly, and then v
m
(
or
l
m
)
∈
A
R
(
c
h
,
r
2
) is obtained based on the method in Section 3.1.2, m
=
1,2
,
…
,
θ. Similarly, φ vertexes and links in R
(
c
h
,
r
2
) are selected. If neither capacity constraint of R
(
c
h
,
r
1
) nor that of R
(
c
h
,
r
2
) is unsatisfied after insertions, v
m
(
or
l
m
) can be inserted into the corresponding place in A
R
(
c
h
,
r
1
), or v
n
(
or
l
n
) in A
R
(
c
h
,
r
2
). Otherwise, v
m
(
or
l
m
) and v
n
(
or
l
n
) will be selected again.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 R
(
c
h
,
r
1
), a segment s
1 from γ
1 to γ
2 is extracted. Meanwhile, segment s
*
=
{
s
∣
min
⁡
d
(
c
s
1
,
c
s
)
} with a distance of λ is selected in a neighboring route of R
(
c
h
,
r
1
), where d
(
c
s
1
,
c
s
) represents the distance from the center of segment s
1 to the center of the whole route s. If the exchange of s
1 and s
* 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.


## 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.Table 1
Column headings for Tables2–10.
I
Instance ID
n
Number of customers
m
Number of depots
Q
Vehicle capacity
D
Maximum duration of a route
c
(
s
*
)
Previous best known solution cost
c
(
s
*
*
)
Best known solution found by proposed tabu search algorithm on FL
L
min
⁡, L
max
⁡
Min and max amount limitation for vertexes to be added into the FL base
δ
i
,
j
Borderline parameter of the sharing area
Δ
δ
Updating step ofδ
i
,
j
Δ
s
max
⁡
Max difference between optimal vertex and local optimal vertexes in an optimization period
n
Y
Number of optimization periods
n
L
The total number of local optimal vertexes
Ω
′
The total number of candidate links
ρ
1
i
The proportion of customers in FLs to total customers in optimization periodi
ρ
2
j
The proportion of customers in FLs to total customers in local optimal vertexj
ω
1
i
The optimal solution value of optimization periodi
ω
2
j
The value local optimal solutionjTable 2
MDVRP instances and previous best known solution cost.
I
source
n
Q
D
m
GJ
CGW
CGL
c
(
s
*
)
p08
GJ
249
500
310
2
4832.0
4511.6
4437.68
4437.68
p09
GJ
249
500
310
3
4219.7
3950.9
3900.22
3900.22
p10
GJ
249
500
310
4
3822.0
3727.1
3663.02
3663.02
p11
GJ
249
500
310
5
3754.1
3670.2
3554.18
3554.18
p15
CGW
160
60
0
4
—
2610.3
2505.42
2505.42
p16
CGW
160
60
200
4
—
2605.3
2572.23
2572.23
p17
CGW
160
60
180
4
—
2816.6
2709.09
2709.09
p18
CGW
240
60
0
6
—
3877.4
3702.85
3702.85
p19
CGW
240
60
200
6
—
3863.9
3827.06
3827.06
p20
CGW
240
60
180
6
—
4272.0
4058.07
4058.07
p21
CGW
360
60
0
9
—
5791.5
5474.84
5474.84
p22
CGW
360
60
200
9
—
5857.4
5702.16
5702.16
p23
CGW
360
60
180
9
—
6494.6
6095.46
6095.46
pr04
CGL
192
185
440
4
—
—
2072.52
2072.52
pr05
CGL
240
180
420
4
—
—
2385.77
2385.77
pr06
CGL
288
175
400
4
—
—
2723.27
2723.27
pr09
CGL
216
180
450
6
—
—
2153.1
2153.1
pr10
CGL
288
170
425
6
—
—
2921.85
2921.85Table 3
Sensitivity data of parametersL
min
⁡ and L
max
⁡.
I
n
m
[
L
min
⁡
*
,
L
max
⁡
*
]
Solution value % gap
CPU time % gap
p08
249
2
[
3,7
]
3.37
403.23
p09
249
3
[
3,6
]
4.6
476.87
p10
249
4
[
3,6
]
7.62
322.33
p11
249
5
[
3,5
]
6.26
398.74
p15
160
4
[
2,4
]
5.06
366.86
p16
160
4
[
2,4
]
3.01
485.78
p17
160
4
[
2,4
]
3.1
573.27
p18
240
6
[
2,4
]
6.53
397.2
p19
240
6
[
2,4
]
4.61
676.61
p20
240
6
[
2,4
]
7.85
238.97
p21
360
9
[
2,4
]
8.18
291.3
p22
360
9
[
2,4
]
6.76
334.59
p23
360
9
[
2,4
]
6.14
618.93
pr04
192
4
[
2,5
]
7.57
474.76
pr05
240
4
[
3,6
]
3.17
356.43
pr06
288
4
[
3,6
]
4.31
649.79
pr09
216
6
[
2,4
]
7.55
414.62
pr10
288
6
[
2,4
]
3.62
689.48Table 4
Sensitivity data of parameterδ
i
,
j and its update step Δ
δ.
I
n
m
δ
i
,
j
*
Δ
δ
*
Solution value
CPU time
% gap 1
% gap 2
% gap 3
% gap 4
p08
249
2
0.85
0.10
11.54
1.51
11.81
29.2
p09
249
3
0.75
0.08
9.82
1.79
19.13
34.48
p10
249
4
0.65
0.04
4.35
1.57
10.83
30.79
p11
249
5
0.65
0.04
17.68
1.2
16.61
51.14
p15
160
4
0.65
0.04
3.24
1.69
7.99
26.85
p16
160
4
0.65
0.06
16.32
1.95
8.68
41.14
p17
160
4
0.65
0.04
9.49
1.43
14.09
45.51
p18
240
6
0.7
0.06
9.59
2.59
12.67
48.92
p19
240
6
0.65
0.04
18.23
1.01
19.57
37.92
p20
240
6
0.65
0.06
13.66
1.99
17.61
44.74
p21
360
9
0.65
0.04
6.95
2.84
4.9
22.24
p22
360
9
0.65
0.04
8.49
1.9
19.68
57.11
p23
360
9
0.65
0.04
10.8
2.55
15.73
21.33
pr04
192
4
0.65
0.04
18.3
2.52
8.75
50.53
pr05
240
4
0.7
0.06
11.59
1.59
14.52
25.42
pr06
288
4
0.75
0.08
12.44
1.13
13.44
35.39
pr09
216
6
0.6
0.04
6.96
1.73
19.96
46.54
pr10
288
6
0.65
0.04
17.1
2.07
12.92
57.18
% Gap 1: % Gap between the worst and the best solution value with differentδ
i
,
j and Δ
δ
=
0.04.
% Gap 2: % Gap between the worst and the best solution value with differentΔ
δ and δ
i
,
j
*.
% Gap 3: % Gap between the longest and the shortest CPU time with differentδ
i
,
j and Δ
δ
=
0.04.
% Gap 4: % Gap between the longest and the shortest CPU time with differentΔ
δ and δ
i
,
j
*.Table 5
The value range ofδ
i
,
j
* and Δ
δ
*.
n
/
m
δ
i
,
j
*
Δ
δ
*
≤50
[
0.6,0.7
]
[
0.04,0.06
]
≤100
(
0.7,0.8
]
(
0.06,0.08
]
>100
(
0.8,0.9
]
(
0.08,0.10
]Table 6
The calculation results of parameterΔ
s
max
⁡.
I
n
m
Δ
s
max
⁡
*
n
Y
n
L
Ω
′
p08
249
2
0.03
3
17
201
p09
249
3
0.04
5
18
193
p10
249
4
0.03
7
21
187
p11
249
5
0.03
7
42
422
p15
160
4
0.06
5
18
143
p16
160
4
0.07
5
18
159
p17
160
4
0.06
5
17
129
p18
240
6
0.05
8
31
255
p19
240
6
0.06
9
48
646
p20
240
6
0.05
10
47
531
p21
360
9
0.02
10
45
559
p22
360
9
0.02
11
46
315
p23
360
9
0.02
9
43
302
pr04
192
4
0.06
5
23
309
pr05
240
4
0.03
6
36
400
pr06
288
4
0.03
6
21
215
pr09
216
6
0.05
8
27
232
pr10
288
6
0.04
7
41
343Table 7
The value ofn
Y, n
L, and Ω
′ and the accuracies of Ω
′ with different Δ
s
max
⁡ in the instance p09.
Δ
s
max
⁡
*
n
Y
n
L
Ω
′
Ω
′ accuracy
0.01
3
5
58
37.77%
0.02
3
8
101
38.37%
0.03
4
11
131
40.71%
0.04
5
18
193
42.40%
0.05
9
23
287
32.05%
0.06
10
26
402
23.77%
0.07
12
27
606
21.96%
0.08
12
40
726
15.59%
0.09
13
55
958
10.02%
0.1
15
101
974
11.23%Table 8
Numbers of FLs in every optimization period.
I
Optimization period
1
2
3
4
5
6
7
8
9
10
11
p08
20
34
38
—
—
—
—
—
—
—
—
p09
13
16
24
28
30
—
—
—
—
—
—
p10
11
15
19
26
30
33
34
—
—
—
—
p11
15
17
19
23
28
30
33
—
—
—
—
p15
12
18
22
26
27
—
—
—
—
—
—
p16
15
17
25
31
31
—
—
—
—
—
—
p17
16
21
25
27
27
—
—
—
—
—
—
p18
15
18
19
23
28
32
32
33
—
—
—
p19
13
18
22
27
32
37
37
38
40
—
—
p20
10
13
12
17
23
27
32
34
36
36
—
p21
11
15
19
27
31
35
38
45
48
50
—
p22
12
17
21
27
33
39
46
51
52
52
53
p23
12
15
19
25
33
43
47
54
56
—
—
pr04
14
21
25
24
25
—
—
—
—
—
—
pr05
16
25
34
39
39
40
—
—
—
—
—
pr06
15
19
25
32
36
38
—
—
—
—
—
pr09
13
18
23
27
31
32
33
35
—
—
—
pr10
19
22
29
34
34
35
37
—
—
—
—Table 9
The correlation ofρ
1
i, ω
1
i, ρ
2
j, and ω
2
j.
I
n
Y
n
L
Correlation coefficient ofρ
1
i and ω
1
i
Correlation coefficient ofρ
2
j and ω
2
j
p08
3
17
-
0.92
-
0.91
p09
5
18
-
0.91
-
0.99
p10
7
21
-
0.97
-
0.96
p11
7
42
-
0.97
-
0.97
p15
5
18
-
0.89
-
0.89
p16
5
18
-
0.97
-
0.90
p17
5
17
-
0.88
-
0.93
p18
8
31
-
0.95
-
0.90
p19
9
48
-
0.92
-
0.95
p20
10
47
-
0.93
-
0.98
p21
10
45
-
0.88
-
0.95
p22
11
46
-
0.90
-
0.96
p23
9
43
-
0.91
-
0.93
pr04
5
23
-
0.96
-
0.91
pr05
6
36
-
0.90
-
0.89
pr06
6
21
-
0.96
-
0.95
pr09
8
27
-
0.90
-
0.92
pr10
7
41
-
0.92
-
0.90Table 10
The results found by proposed tabu search algorithm on FL.
I
n
c
(
s
*
)
c
(
s
*
*
)
Average CPU time (min)
% CPU time
Average solution value
% solution value
% above best known solution
p08
249
4437.68
4417.46
10.40
27.04
4465.69
1.08
0.46
p09
249
3900.22
3886.82
10.32
28.60
3937.76
1.29
0.34
p10
249
3663.02
3647.22
9.61
28.20
3684.78
1.02
0.43
p11
249
3554.18
3547.70
9.20
12.03
3611.16
1.76
0.18
p15
160
2505.42
2513.43
2.88
8.31
2528.42
0.59
-
0.32
p16
160
2572.23
2578.88
4.04
9.03
2591.07
0.47
-
0.26
p17
160
2709.09
2709.09
3.43
11.86
2710.19
0.04
0
p18
240
3702.85
3712.61
4.22
12.45
3719.78
0.19
-
0.26
p19
240
3827.06
3827.06
3.73
12.58
3853.29
0.68
0
p20
240
4058.07
4066.85
5.33
18.03
4078.59
0.29
-
0.22
p21
360
5474.84
5483.31
7.72
20.46
5518.1
0.63
-
0.15
p22
360
5702.16
5708.12
8.05
11.57
5743.57
0.62
-
0.10
p23
360
6095.46
6095.46
7.07
17.23
6136.23
0.66
0
pr04
192
2072.52
2058.80
9.63
32.54
2081.42
1.09
0.66
pr05
240
2385.77
2369.00
8.59
19.30
2405.35
1.51
0.70
pr06
288
2723.27
2718.11
13.25
19.24
2763.24
1.63
0.19
pr09
216
2153.1
2147.40
8.04
13.21
2182.06
1.59
0.26
pr10
288
2921.85
2865.31
10.88
22.79
2934.16
2.35
1.94
% 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.The characteristics of the instances are summarized in Table2 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.[
L
min
,
L
max
]

Firstly, all other parameters are set as constants.L
max
⁡ is selected from (3, 4, 5, 6, 7) and L
min
⁡ is from (2, 3, 4), given L
max
⁡
≥
L
min
⁡. 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[
L
min
⁡
,
L
max
⁡
] has not much influence on solutions, and average gap between solutions is 5.52%. Inferior solutions always occur in the condition where L
min
⁡ is relatively large and L
min
⁡
=
L
max
⁡; for example, [
L
min
⁡
,
L
max
⁡
]
=
[
5
,
5
]. However, the value of [
L
min
⁡
,
L
max
⁡
] greatly affects the computation speed of the algorithm, with an average gap of 453.88%, since L
max
⁡ can control the upper limits of lengths of candidate links, and L
min
⁡ determines their lower limits as well as split numbers. With the decrease in L
min
⁡, the number of splits increases with an exponential distribution. According to this case, the relation of [
L
min
⁡
*
,
L
max
⁡
*
] and n
/
m is revealed: a larger n
/
m leads to larger L
min
⁡
* and L
max
⁡
*. When the value of n
/
m is less than 40, the commendation value of [
L
min
⁡
*
,
L
max
⁡
*
] is [
2
,
4
], when 40
<
n
/
m
≤
100, the value is [
3
,
5
] or [
3
,
6
], and when n
/
m is more than 100, the value is [
3
,
7
].


#### 4.1.2. Borderline Parameterδ
i
,
j and Update Step Δ
δ

Keeping other parameters unchanged, the authors measure the effect ofδ
i
,
j here. Firstly, Δ
δ
=
0.04, δ
i
,
j is in [0.5, 0.9], the interval between successive δ
i
,
j is 0.05, and then δ
i
,
j
* for 18 instances are obtained through 9 times of calculations for each instance. Next, Δ
δ is taken varied in [0.02, 0.12], δ
i
,
j
* 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 Table4, given Δ
δ is a constant, difference between the best and the worst solution values of different δ
i
,
j is 11.48%. If δ
i
,
j 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 δ
i
,
j has a more obvious impact on the computation result of the algorithm and that of Δ
δ. Similarly, valuing of δ
i
,
j 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 δ
i
,
j
* and Δ
δare proportional to the customer numbers allocated in each depot (n
/
m). The approximate value range of δ
i
,
j
* and Δ
δ
* is shown in Table 5.


#### 4.1.3.Δ
s
max

Set[
L
min
⁡
,
L
max
⁡
]
=
[
L
min
⁡
*
,
L
max
⁡
*
], Δ
s
max
⁡ in [0.01, 0.10], and the interval between different Δ
s
max
⁡ as 0.01. The authors calculate 10 times for each instance, the number of optimization periods is n
Y, the number of local optimal vertexes in optimization period k is n
L
k, n
L
=
∑
k
=
1
n
Y
n
L
k, Ω
k
′ is the number of candidate links produced in optimization period k, and Ω
′
=
∑
k
=
1
n
Y
Ω
k
′. The calculation results are shown in Table 6.According to the analysis results,Δ
s
max
⁡ can control the number of local optimal vertexes n
L, and the value of Δ
s
max
⁡ is inversely proportional to n. This is because global optimization is difficult with the increase in n and more local optimal vertexes attained. As a result, a smaller Δ
s
max
⁡ is needed to restrict the number of n
L, 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 Δ
s
max
⁡ for computation scale and characteristics is needed. Table 7 lists the value of n
Y, n
L, and Ω
′ and the accuracies of Ω
′ with different Δ
s
max
⁡ in the instance p09.As shown in Table7, when Δ
s
max
⁡ is at its best value Δ
s
max
⁡
*, Ω
′ is the most accurate, so as to guarantee the most efficient and effective optimization. Based on this study, the recommendation value of Δ
s
max
⁡
* is [0.06, 0.07] when 150
<
n
≤
200, [0.03, 0.06] when 200
<
n
≤
300, and [0.02, 0.03] when 300
<
n
≤
400.


### 4.2. FLs and Optimization Results

With the best parameters in Section4.1, the generation process of FLs in the optimization is tracked. Set ρ
1
i as the proportion of customers in FLs to total customers in each optimization period, ω
1
i as the optimal solution value of optimization period i, i
=
1,2
,
…
,
n
Y, ρ
2
j as the proportion of customers in FLs to total customers in each local optimal vertex, and ω
2
j as the value local optimal solution j, j
=
1,2
,
…
,
n
L. Table 8 shows numbers of FLs in every optimization period.In Tables8 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 ρ
1
i and ω
1
i reaches −0.924, and the number between ρ
2
j and ω
2
j 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
i and ρ
2
j 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 Table10. Figure 11 describes the corresponding route of the optimal solution for p09.Figure 11
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
<
n
≤
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.


## 4.1. Sensitivity Analyses

### 4.1.1.[
L
min
,
L
max
]

Firstly, all other parameters are set as constants.L
max
⁡ is selected from (3, 4, 5, 6, 7) and L
min
⁡ is from (2, 3, 4), given L
max
⁡
≥
L
min
⁡. 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[
L
min
⁡
,
L
max
⁡
] has not much influence on solutions, and average gap between solutions is 5.52%. Inferior solutions always occur in the condition where L
min
⁡ is relatively large and L
min
⁡
=
L
max
⁡; for example, [
L
min
⁡
,
L
max
⁡
]
=
[
5
,
5
]. However, the value of [
L
min
⁡
,
L
max
⁡
] greatly affects the computation speed of the algorithm, with an average gap of 453.88%, since L
max
⁡ can control the upper limits of lengths of candidate links, and L
min
⁡ determines their lower limits as well as split numbers. With the decrease in L
min
⁡, the number of splits increases with an exponential distribution. According to this case, the relation of [
L
min
⁡
*
,
L
max
⁡
*
] and n
/
m is revealed: a larger n
/
m leads to larger L
min
⁡
* and L
max
⁡
*. When the value of n
/
m is less than 40, the commendation value of [
L
min
⁡
*
,
L
max
⁡
*
] is [
2
,
4
], when 40
<
n
/
m
≤
100, the value is [
3
,
5
] or [
3
,
6
], and when n
/
m is more than 100, the value is [
3
,
7
].


### 4.1.2. Borderline Parameterδ
i
,
j and Update Step Δ
δ

Keeping other parameters unchanged, the authors measure the effect ofδ
i
,
j here. Firstly, Δ
δ
=
0.04, δ
i
,
j is in [0.5, 0.9], the interval between successive δ
i
,
j is 0.05, and then δ
i
,
j
* for 18 instances are obtained through 9 times of calculations for each instance. Next, Δ
δ is taken varied in [0.02, 0.12], δ
i
,
j
* 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 Table4, given Δ
δ is a constant, difference between the best and the worst solution values of different δ
i
,
j is 11.48%. If δ
i
,
j 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 δ
i
,
j has a more obvious impact on the computation result of the algorithm and that of Δ
δ. Similarly, valuing of δ
i
,
j 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 δ
i
,
j
* and Δ
δare proportional to the customer numbers allocated in each depot (n
/
m). The approximate value range of δ
i
,
j
* and Δ
δ
* is shown in Table 5.


### 4.1.3.Δ
s
max

Set[
L
min
⁡
,
L
max
⁡
]
=
[
L
min
⁡
*
,
L
max
⁡
*
], Δ
s
max
⁡ in [0.01, 0.10], and the interval between different Δ
s
max
⁡ as 0.01. The authors calculate 10 times for each instance, the number of optimization periods is n
Y, the number of local optimal vertexes in optimization period k is n
L
k, n
L
=
∑
k
=
1
n
Y
n
L
k, Ω
k
′ is the number of candidate links produced in optimization period k, and Ω
′
=
∑
k
=
1
n
Y
Ω
k
′. The calculation results are shown in Table 6.According to the analysis results,Δ
s
max
⁡ can control the number of local optimal vertexes n
L, and the value of Δ
s
max
⁡ is inversely proportional to n. This is because global optimization is difficult with the increase in n and more local optimal vertexes attained. As a result, a smaller Δ
s
max
⁡ is needed to restrict the number of n
L, 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 Δ
s
max
⁡ for computation scale and characteristics is needed. Table 7 lists the value of n
Y, n
L, and Ω
′ and the accuracies of Ω
′ with different Δ
s
max
⁡ in the instance p09.As shown in Table7, when Δ
s
max
⁡ is at its best value Δ
s
max
⁡
*, Ω
′ is the most accurate, so as to guarantee the most efficient and effective optimization. Based on this study, the recommendation value of Δ
s
max
⁡
* is [0.06, 0.07] when 150
<
n
≤
200, [0.03, 0.06] when 200
<
n
≤
300, and [0.02, 0.03] when 300
<
n
≤
400.


## 4.1.1.[
L
min
,
L
max
]

Firstly, all other parameters are set as constants.L
max
⁡ is selected from (3, 4, 5, 6, 7) and L
min
⁡ is from (2, 3, 4), given L
max
⁡
≥
L
min
⁡. 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[
L
min
⁡
,
L
max
⁡
] has not much influence on solutions, and average gap between solutions is 5.52%. Inferior solutions always occur in the condition where L
min
⁡ is relatively large and L
min
⁡
=
L
max
⁡; for example, [
L
min
⁡
,
L
max
⁡
]
=
[
5
,
5
]. However, the value of [
L
min
⁡
,
L
max
⁡
] greatly affects the computation speed of the algorithm, with an average gap of 453.88%, since L
max
⁡ can control the upper limits of lengths of candidate links, and L
min
⁡ determines their lower limits as well as split numbers. With the decrease in L
min
⁡, the number of splits increases with an exponential distribution. According to this case, the relation of [
L
min
⁡
*
,
L
max
⁡
*
] and n
/
m is revealed: a larger n
/
m leads to larger L
min
⁡
* and L
max
⁡
*. When the value of n
/
m is less than 40, the commendation value of [
L
min
⁡
*
,
L
max
⁡
*
] is [
2
,
4
], when 40
<
n
/
m
≤
100, the value is [
3
,
5
] or [
3
,
6
], and when n
/
m is more than 100, the value is [
3
,
7
].


## 4.1.2. Borderline Parameterδ
i
,
j and Update Step Δ
δ

Keeping other parameters unchanged, the authors measure the effect ofδ
i
,
j here. Firstly, Δ
δ
=
0.04, δ
i
,
j is in [0.5, 0.9], the interval between successive δ
i
,
j is 0.05, and then δ
i
,
j
* for 18 instances are obtained through 9 times of calculations for each instance. Next, Δ
δ is taken varied in [0.02, 0.12], δ
i
,
j
* 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 Table4, given Δ
δ is a constant, difference between the best and the worst solution values of different δ
i
,
j is 11.48%. If δ
i
,
j 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 δ
i
,
j has a more obvious impact on the computation result of the algorithm and that of Δ
δ. Similarly, valuing of δ
i
,
j 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 δ
i
,
j
* and Δ
δare proportional to the customer numbers allocated in each depot (n
/
m). The approximate value range of δ
i
,
j
* and Δ
δ
* is shown in Table 5.


## 4.1.3.Δ
s
max

Set[
L
min
⁡
,
L
max
⁡
]
=
[
L
min
⁡
*
,
L
max
⁡
*
], Δ
s
max
⁡ in [0.01, 0.10], and the interval between different Δ
s
max
⁡ as 0.01. The authors calculate 10 times for each instance, the number of optimization periods is n
Y, the number of local optimal vertexes in optimization period k is n
L
k, n
L
=
∑
k
=
1
n
Y
n
L
k, Ω
k
′ is the number of candidate links produced in optimization period k, and Ω
′
=
∑
k
=
1
n
Y
Ω
k
′. The calculation results are shown in Table 6.According to the analysis results,Δ
s
max
⁡ can control the number of local optimal vertexes n
L, and the value of Δ
s
max
⁡ is inversely proportional to n. This is because global optimization is difficult with the increase in n and more local optimal vertexes attained. As a result, a smaller Δ
s
max
⁡ is needed to restrict the number of n
L, 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 Δ
s
max
⁡ for computation scale and characteristics is needed. Table 7 lists the value of n
Y, n
L, and Ω
′ and the accuracies of Ω
′ with different Δ
s
max
⁡ in the instance p09.As shown in Table7, when Δ
s
max
⁡ is at its best value Δ
s
max
⁡
*, Ω
′ is the most accurate, so as to guarantee the most efficient and effective optimization. Based on this study, the recommendation value of Δ
s
max
⁡
* is [0.06, 0.07] when 150
<
n
≤
200, [0.03, 0.06] when 200
<
n
≤
300, and [0.02, 0.03] when 300
<
n
≤
400.


## 4.2. FLs and Optimization Results

With the best parameters in Section4.1, the generation process of FLs in the optimization is tracked. Set ρ
1
i as the proportion of customers in FLs to total customers in each optimization period, ω
1
i as the optimal solution value of optimization period i, i
=
1,2
,
…
,
n
Y, ρ
2
j as the proportion of customers in FLs to total customers in each local optimal vertex, and ω
2
j as the value local optimal solution j, j
=
1,2
,
…
,
n
L. Table 8 shows numbers of FLs in every optimization period.In Tables8 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 ρ
1
i and ω
1
i reaches −0.924, and the number between ρ
2
j and ω
2
j 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
i and ρ
2
j 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 Table10. Figure 11 describes the corresponding route of the optimal solution for p09.Figure 11
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
<
n
≤
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.


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



---

*Source: 152494-2014-07-23.xml*