# A Hybrid Bat Algorithm with Path Relinking for Capacitated Vehicle Routing Problem

Authors: Yongquan Zhou; Jian Xie; Hongqing Zheng
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The capacitated vehicle routing problem (CVRP) is an NP-hard problem with wide engineering and theoretical background. In this paper, a hybrid bat algorithm with path relinking (HBA-PR) is proposed to solve CVRP. The HBA-PR is constructed based on the framework of continuous bat algorithm; the greedy randomized adaptive search procedure (GRASP) and path relinking are effectively integrated into bat algorithm. Moreover, in order to further improve the performance, the random subsequences and single-point local search are operated with certain loudness (probability). In order to verify the validity of the method in this paper, and it's efficiency and with other existing methods, several classical CVRP instances from three classes of CVRP benchmarks are selected to tested. Experimental results and comparisons show that the HBA-PR is effective for CVRP.


---

## Body

## 1. Introduction

The vehicle routing problem (VRP) is a classical combinatorial optimization problem that was proposed in the late 1950s, and it is still a research hotspot in the field of Operations Research. The capacitated vehicle routing problem (CVRP) is introduced by Dantzig and Ramser in 1959 [1], which designs a set of customer demands that have to be served with a fleet of vehicles from a depot or central node, each vehicle has the uniform capacity, and each customer has a certain demand that must be satisfied at minimal cost. These costs usually represent distances, traveling times, number of vehicles employed, or a combination of these factors.It is known that the CVRP is an NP-hard problem [2]. Various approaches have been presented to solve the CVRP during the last decades, such as linear programming [3], several metaheuristics [4–6], and many hybrid heuristics with variable neighborhood search or constructive heuristic methods [7–10]. The overview of methods presented in [6] shows that at least 29 different methods for the CVRP exist; all achieve more or less comparable performance. Although several methods can produce good solutions, the computational time is long when the scale of instances is large. Meanwhile, an abundance of methods for the CVRP is a population-based algorithm and the parameter setting of the algorithm, is pretty important; however, the parameter setting of many metaheuristics is not considered in the literature.Bat algorithm (BA) is a fairly new metaheuristics proposed by Yang [11] in 2010, which inspired by the intelligent echolocation behavior of microbats when they forage. As we know, many new metaheuristics have been widely used and successfully applied to solve the CVRP notwithstanding BA has not yet been applied to solve the CVRP, and BA has been applied to solve other problems with great success. For example, Gandom et al. focus on solving constrained optimization tasks [12]. Yang and Gandom apply bat algorithm to solve many global engineering optimizations [13]. A classification mode is proposed by Mishra et al. using bat algorithm to update the weights of a functional link artificial neural network (FLANN) classifier [14]. Meanwhile, some researchers have improved bat algorithm and applied it to various optimization problems. Xie et al. proposed a DLBA bat algorithm based on differential operator and Lévy flights trajectory to solve function optimization and nonlinear equations [15]. Wang et al. proposed a new bat algorithm with mutation (BAM) to solve the uninhabited combat air vehicle (UCAV) path planning problem [16]. In this paper, we propose a hybrid bat algorithm (HBA-PR) to solve capacitated vehicle routing problem.The rest of this paper is organized as follows. The problem description of CVRP and original bat algorithm are described in Section2. The hybrid bat algorithm (HBA-PR) for CVRP is described detailedly in Section 3. The experimental results of the HBA-PR and comparisons with other previous algorithms are shown in Section 4. In the last section, we conclude this paper and point out some future work in Section 5.


## 2. Problem Descriptions and Bat Algorithm

### 2.1. Capacitated Vehicle Routing Problem

The CVRP is considered to be the classical version of the VRP, which designs a set of customer demands that have to be served with a fleet of vehicles from a depot or central node, each vehicle has the uniform capacity, and each customer has a certain demand. The objective makes the expended cost as minimal as possible. LetG=(V,E) be a complete graph with a set of vertices V={0,1,…,k}, where the vertex {0} represents the depot and the remaining ones represent the customers. Each edge eij={i,j}∈E has a nonnegative cost cij, and each customer i∈V′=V∖{0} has a demand di. Let S=(1,2,…,m) be the set of homogeneous vehicles with capacity Q. The CVRP consists in constructing a set of up to k routes in such a way that (1) every route starts and ends at the depot; (2) all demands are accomplished; (3) the vehicle’s capacity is not exceeded; (4) a customer is visited by only a single vehicle; (5) the sum of costs is minimized. The mathematical formulas are defined as follows [10]:
(1)minZ=∑i=0k∑j=0k∑s=0mcijeijs,(2)s.t.∑i=0kdiyis≤Q,s=1,2,…,m,∑s=1myis=1,i=1,2,…,k,∑i=0keijs=yis,j=1,2,…,k;s=1,2,…,m,∑j=0keijs=yis,i=1,2,…,k;s=1,2,…,m,
where s denotes the number of vehicles; eijs=1 if vehicle s is from i to j, otherwise eijs=0; the yis=1 if vehicle s is loading (active), otherwise yis=0.


### 2.2. Basic Bat Algorithm

The basic bat algorithm (BA) is a metaheuristic first introduced by Yang in 2010. In simulations of BA, under several under idealized rules, the updated rules of bat’s positionsxi and velocities vi in a D-dimensional search space are defined. The new solutions xit and velocities vit at generation t are given by
(3)fri=frmin+(frmax-frmin)β,vit=vit-1+(xit-x*)fri,xit=xit-1+vit,
where β∈[0,1] is a random vector drawn from a uniform distribution and fri denotes frequency of each bat. Generally, the frequency fri∈[frmin,frmax]. Here x* is the current global best location (solution) which is located after comparing all the solutions among all the n bats.After the position updating of bat, a random number is generated; if the random number is greater than the pulse emission rateri, a new position will be generated around the current best solutions; it can be represented by
(4)x=x*+ε×Ldt,
where ε∈[-1,1] is a random number, while Ldt=〈Ldit〉 is the average loudness of all the bats at current generation t.Furthermore, the loudnessLdi and the pulse emission rate ri will be updated, and a solution will be accepted if a random number is less than loudness Ldi and f(xi)<f(x*). Ldi and ri are updated by
(5)Ldit+1=α×Ldit,rit+1=ri0×[1-exp(-γ×t)],
where α, γ are constants and f(·) is fitness function. The algorithm repeat until the termination criterion is reached.


## 2.1. Capacitated Vehicle Routing Problem

The CVRP is considered to be the classical version of the VRP, which designs a set of customer demands that have to be served with a fleet of vehicles from a depot or central node, each vehicle has the uniform capacity, and each customer has a certain demand. The objective makes the expended cost as minimal as possible. LetG=(V,E) be a complete graph with a set of vertices V={0,1,…,k}, where the vertex {0} represents the depot and the remaining ones represent the customers. Each edge eij={i,j}∈E has a nonnegative cost cij, and each customer i∈V′=V∖{0} has a demand di. Let S=(1,2,…,m) be the set of homogeneous vehicles with capacity Q. The CVRP consists in constructing a set of up to k routes in such a way that (1) every route starts and ends at the depot; (2) all demands are accomplished; (3) the vehicle’s capacity is not exceeded; (4) a customer is visited by only a single vehicle; (5) the sum of costs is minimized. The mathematical formulas are defined as follows [10]:
(1)minZ=∑i=0k∑j=0k∑s=0mcijeijs,(2)s.t.∑i=0kdiyis≤Q,s=1,2,…,m,∑s=1myis=1,i=1,2,…,k,∑i=0keijs=yis,j=1,2,…,k;s=1,2,…,m,∑j=0keijs=yis,i=1,2,…,k;s=1,2,…,m,
where s denotes the number of vehicles; eijs=1 if vehicle s is from i to j, otherwise eijs=0; the yis=1 if vehicle s is loading (active), otherwise yis=0.


## 2.2. Basic Bat Algorithm

The basic bat algorithm (BA) is a metaheuristic first introduced by Yang in 2010. In simulations of BA, under several under idealized rules, the updated rules of bat’s positionsxi and velocities vi in a D-dimensional search space are defined. The new solutions xit and velocities vit at generation t are given by
(3)fri=frmin+(frmax-frmin)β,vit=vit-1+(xit-x*)fri,xit=xit-1+vit,
where β∈[0,1] is a random vector drawn from a uniform distribution and fri denotes frequency of each bat. Generally, the frequency fri∈[frmin,frmax]. Here x* is the current global best location (solution) which is located after comparing all the solutions among all the n bats.After the position updating of bat, a random number is generated; if the random number is greater than the pulse emission rateri, a new position will be generated around the current best solutions; it can be represented by
(4)x=x*+ε×Ldt,
where ε∈[-1,1] is a random number, while Ldt=〈Ldit〉 is the average loudness of all the bats at current generation t.Furthermore, the loudnessLdi and the pulse emission rate ri will be updated, and a solution will be accepted if a random number is less than loudness Ldi and f(xi)<f(x*). Ldi and ri are updated by
(5)Ldit+1=α×Ldit,rit+1=ri0×[1-exp(-γ×t)],
where α, γ are constants and f(·) is fitness function. The algorithm repeat until the termination criterion is reached.


## 3. Hybrid Bat Algorithm with Path Relinking for CVRP

### 3.1. Solution Representation in HBA-PR

Since standard BA is a continuous optimization algorithm, the standard continuous encoding scheme of BA cannot be used to solve CVRP directly. In order to apply BA to solve CVRP, the first step is to devise a suitable representation for the candidate solutions in designing a hybrid bat algorithm for a particular problem. Each individual is a sequence with integer number which is the order of visiting these customers; the 0 represents the depot. For example, if we have the following three routes:(6) the coded individual with integer is 0→3→5→8→4→9→2→6→1→7→0, and the bat individual is represented as 3→5→8→4→9→2→6→1→7.


### 3.2. Hybrid Bat Algorithm

Aimed at the capacitated vehicle routing problem, based on the idea of bat algorithm, a hybrid bat algorithm is proposed, which integrates greedy randomized adaptive search procedure (GRASP) heuristic and bat algorithm, and the path relinking as an intensification strategy to explore local trajectories connecting elite solutions obtained by the proposed algorithm. The hybrid bat algorithm with path relinking is named as HBA-PR.GRASP [17, 18] is a heuristic already applied to many optimization problems successfully [19–21]. GRASP consists of a two-phase iterative process: construction phase and local search phase. In the first phase, a greedy randomized solution is built. Due to this solution is not guaranteed to be locally optimal; a local search is performed in the second phase. The final result is simply the best solution found over all iterations.The construction phase can be described as a process which stepwise adds one element at a time to a partial (incomplete) solution. According to a greedy function, all elements are sorted, and the Restricted Candidate List (RCL) is constructed based on the order, and then from the list randomly select a element. In the second phase, a local search is initialized from these points, this iterative process is repeated until a termination criterion is met, and the best solution found over all iterations is taken as the result.RCL is created using a parameterα to restrict the size of this list. Candidate e∈C is sorted according to their greedy function value f(e). In a cardinality-based RCL, the latter is made up by the k top-ranked elements. In a value-based construction, the RCL consists of the elements in the set {e∈C:f*≤f(e)≤f*+α×(f*-f*)}, where f*=min(f(e):e∈C), f*=max{f(e):e∈C}, and α∈[0,1]. Since the best value for α is often difficult to determine, a random value is often assigned. The values for α adopted in the constructive heuristics are set using reactive strategies, which usually leads to better performance than using fixed values [22].In original bat algorithm, the bat individual randomly selects a certain range of frequency, its velocity is updated according to their selected frequency, and at last a new position is generated using its velocity and its own position. The idea is that the position of bat individual is updated by adjusting its frequency of sonic pulse. In this paper, the position updating of bat adopted GRASP to generate a new position, the frequency is used for restricting the size of CRL, frequency equivalent to parameterα in GRASP, and the frequency is variable value. In basic bat algorithm, the loudness Ld and the pulse emission rate r are updated according as the iterations proceed. As the loudness usually decreases, while the rate of pulse emission increases, it indicates that the bats approximate their prey (optimum solution). The pulse emission rate r is updated by
(7)r(t)=(1+exp(-5tmax×(t-tmax/2))-1),
where t denotes the tth generation and tmax is the maximal generation. The rate r is similar to Sigmoid function, and the frequency fr is determined according to the pulse emission rate r, which is represented by
(8)fr={1-max(0.2,min(0.8,r)),rand>r,max(0.2,min(0.8,r)),rand≤r,
where rand is a random number and 0.2 and 0.8 are experience parameters reference [18]. The frequency fr decreases gradually at firstly and then increases gradually, while the generation t increases. Figure 1 is the changing curve of rate r, Figure 2 is an example of frequency fr, and Algorithm 1 shows the pseudocode of greedy randomized construction with frequency fr.Algorithm 1:
Greedy_Randomized_Construction(fr).
S
←
R  rand_Chosen_Vertex(v);//    R is selected at random a vertex as initial solution
C
←
V
∖
v
; // The candidate set
C
is initialized
While
C
≠
⌀
do
IC←Evaluate_Incremental_Costs(S); //  the incremental costs are evaluated
icmin←min({ic∈IC});
icmax←max({ic∈IC});
RCL←{e∈C∣ic(e)≤icmin+fr×(icmax-icmin)};  //RCLis created
s←Select_Element (RCL); // a vertexsis randomly selected fromRCL
S←Obtain_Min_Solution(S,s); // the partial tour is updated by inserting the
vertexs
C←C∖{s}; // the candidate setCis updated
end
return
SFigure 1
Changing curve of rater.Figure 2
Changing curve of frequencyfr.The local search phase uses the 2-Opt heuristic for exchanging. Algorithm2 shows the local search procedure. The input parameter is an initial solution S obtained by the Greedy_Randomized_Construction procedure. The current route need is divided into m subroute according to the load, set the start, and end of suroute said to 0; for example, S={1,2,3,4,5,6,7}, subroute is 0→1→2→3→4→0, and 0→5→6→7→0.Algorithm 2:
Local_Search_Phase (S).
S
′
←
S
;
{
s
r
}
←
C
o
n
s
t
r
u
c
t
_Sub-route(S′);
for
each sub-route
s
r do
sr′←2-opt(sr); // carry out the 2-opt operation
end
S
′
←
C
o
n
s
t
r
u
c
t
_Indiviaual
(
{
s
r
′
}
);
if f(S′)<f(S)then
S←S′
end
return
S


### 3.3. Hybrid Bat Algorithm with Path Relinking

Path relinking was originally proposed by Glover [23]; Laguna and Martí [24] were the first to use path relinking within a GRASP strategy. Path relinking generates new solutions by exploring trajectories connecting an initial solution xs to an elite guiding solution xt. The path relinking procedure consists in selecting moves that introduce attributes contained in the guiding solution xt to the initial solution xs until the initial solution is completely transformed in the guiding solution xt. Path relinking may also be viewed as a constrained local search strategy applied to the initial solution xs. Furthermore, there are several alternatives that have been considered, which involve tradeoffs between computation time and solution quality. These alternatives include periodical relinking, forward relinking, backward relinking, back and forward relinking, mixed relinking, randomized relinking, and truncated relinking[18].One important issue in implementing a path relinking technique is the strategy to construct the elite set (ES). We adopted a fixed size elite set, and a solutionx is inserted into the ES as follows.A solutionx is always inserted into ES if it is not full. Otherwise, the generated solution x is inserted in ES only if its cost is better than the worst cost solution found in ES, and the worst cost solution is replaced by the solution x. Algorithm 3 shows the pseudocode for the algorithm to construct and maintain the elite set ES.Algorithm 3:
Elite_Set(x).
if cardinality(ES) ≠ps then  //psis the population size
E
S
←
E
S
∪
{
x
}
else
E
S
w
←
t
h
e
worst solution in
E
S
if f(x)<f(ESw)then
E
S
←
E
S
∖
{
E
S
w
}
ES←ES∪{x}
end
end
return
E
SAlgorithm4 shows the pseudocode of path relinking procedure; an instance is xs={1,2,3,4,5,6,7}, xt={1,3,4,2,5,6,7}, according to process of Algorithm 4, the first different element is position 2 in xs; after executing replace operation, the x1={1,3,4,2,5,6,7}; the second different element is position 3 in xs; after executing replace operation, the x2={1,2,4,3,5,6,7}; the third different element is position 4 in xs; after executing replace operation, the x3={1,4,3,2,5,6,7}; if the x2 is better than the xt, and the only, then return x2.Algorithm 4:
Path_Relinking (ES*, x).
x
s
←
x
; // initial solution
x
t
←
E
S
*
; //  guiding solution
f
(
x
t
)
←
f
(
E
S
*
);
Δ
←
D
i
f
f
e
r
e
n
c
e(xs,xt); // find out the difference position betweenxs xt
for i=1:  cardinality(Δ) do
j←Find_Position(xs,xt,Δi); // find the position of the element ofxtinΔi which inxs
xi←Replace(xs,Δi,j); // replace thejposition ofxswith the  element ofxsinΔi
xi←Replace(xs,xt,Δi);  //replace theΔiposition ofxswith the corresponding element ofxt
if
f
(
x
i
)
<
f
(
x
t
)
then
x
*
←
x
i
;
f(xt)←f(xi);
end
end
return
x
*


### 3.4. Subsequence and Single-Point Local Search

In this paper, the loudnessLdi of bat individual i is relative to its own fitness fiti, the better fitness, and the lower loudness. The loudness can be described by
(9)Ldi=(fiti-fit*+0.1)(fit*-fit*+0.1),
where constant 0.1 is used for avoiding the denominator is zero, and fiti is the fitness of individual i, fit* and fit* are the minimum and maximum fitness in current population, respectively. In HBA-PR, the loudness reflects the quality of individual. In this algorithm, there are two kinds of local search are embedded into HBA-PR to further improve the performance, random subsequence local search that and single-point local search. The random subsequence local search includes random subsequence inverse and random subsequence insert, and the single-point local search includes single-point insert and single-point swap.For random subsequence insert, firstly, randomly selected an origin of subsequence, and then randomly selected a length of subsequence which is less than length of individualS. Secondly, after determining the subsequence S1, randomly selected an insert point in remainder subsequence S2, S=S1∪S2; the S1 is inserted into S2 location in insert point. An example is shown in Figure 2. For random subsequence inverse, firstly, randomly selected a subsequence with random length, and then the inverse operation is performed. An example is shown in Figure 3.Figure 3
Random subsequence insert operation.For single-point swap, choose two different positions from a permutation randomly and swap them. For single-point insert, choose two different positions from a permutation randomly, and the element in first position is insert into the back of second element. Similarly, two instances are shown in Figures4 and 5.Figure 4
Random subsequence inverse operation.Figure 5
Single-point swap operation.In local search part, the random subsequence local search is performed before random single-point local search. The random subsequence insert and random subsequence inverse are preformed according to loudnessLd. In other words, if a random number is greater than the loudness Ldi, the random subsequence insert is performed; otherwise, the random subsequence inverse is performed. Similarly, the random single-point insert and random single-point swap are performed with the loudness Ld. If a random number is greater than the loudness Ldi, the insert operation is performed (Figure 6); otherwise, the swap operation is performed. Note that where local search is operated on the current optimal individual ES* in elite set ES, local search is performed for each individual. Algorithms 5 and 6 show the pseudocode of subsequence and single-point local search.Algorithm 5:
Sub-sequence_Local_Search (ES*).
for i=1:  ps do
if
rand
>
L
d
i then
x
←
S
u
b
-sequence_Insert(ES*); // perform random sub-sequence insert operation
else
x←Sub-sequence_Inverse(ES*); // perform random sub-sequence inverse operation
end
end
return
xAlgorithm 6:
Single-point_Local_Search (ES*).
for i=1:  psdo
if rand>Ldi then
x←Single-point_Insert(ES*);  //  perform random single-point insert operation
else
x←Single-point_Swap(ES*); // perform random ingle point swap operation
end
end
return
xFigure 6
Single-point insert operation.


### 3.5. HBA-PR Framework for CVRP

We propose in this work, to incorporate greedy randomized adaptive search procedure, path-relinking strategies, subsequence, and single-point local search to the bat algorithm by defining distinct ways to solve capacitated vehicle routing problem. This iterative process is repeated until the termination criterion is met; Algorithm7 shows the pseudocode of HBA-PR for CVRP.Algorithm 7:HBA-PR.
(1)Initialize theps,  bat population and other parameters;
(2)Evaluate fitness for each individual;
(3)ES←Elite_Set(x);
(4)while  t≤tmax  do
(5)Compute pulse emission rate by (7);
(6)Determine frequencyfrby (8);
(7)for i=1:  psdo
(8)x←Greedy_Randomized_Construction(fr);
(9)x←Local_Search_Phase(x);
(10)end
(11)Evaluate fitness for each individualx;
(12)ES←Elite_Set(x);
(13)ES*←Select_Best_Elite(ES); // select the best individual in elite setES
(14)for i=1:ps do
(15)x←Path_Relinking(ES*,  x);
(16)end
(17)Evaluate fitness for each individualx;
(18)ES←Elite_Set(x);
(19)Ldi←Compute_Loudness(f(x)); //Compute loudness of each individual by (9);
(20)ES*←Select_Best_Elite(ES);
(21)x←Sub-sequence_Local_Search(ES*)//carry out random sub-sequence local
search
(22)Evaluate fitness for each individualx;
(23)ES←Elite_Set(x);
(24)ES*←Select_Best_Elite(ES);
(25)x←Single-point_Local_Search(ES*)//carry out random single-point local search
(26)Evaluate fitness for each individualx;
(27)ES←Elite_Set(x);
(28)t=t+1
(29)end
(30)Output result and plot


## 3.1. Solution Representation in HBA-PR

Since standard BA is a continuous optimization algorithm, the standard continuous encoding scheme of BA cannot be used to solve CVRP directly. In order to apply BA to solve CVRP, the first step is to devise a suitable representation for the candidate solutions in designing a hybrid bat algorithm for a particular problem. Each individual is a sequence with integer number which is the order of visiting these customers; the 0 represents the depot. For example, if we have the following three routes:(6) the coded individual with integer is 0→3→5→8→4→9→2→6→1→7→0, and the bat individual is represented as 3→5→8→4→9→2→6→1→7.


## 3.2. Hybrid Bat Algorithm

Aimed at the capacitated vehicle routing problem, based on the idea of bat algorithm, a hybrid bat algorithm is proposed, which integrates greedy randomized adaptive search procedure (GRASP) heuristic and bat algorithm, and the path relinking as an intensification strategy to explore local trajectories connecting elite solutions obtained by the proposed algorithm. The hybrid bat algorithm with path relinking is named as HBA-PR.GRASP [17, 18] is a heuristic already applied to many optimization problems successfully [19–21]. GRASP consists of a two-phase iterative process: construction phase and local search phase. In the first phase, a greedy randomized solution is built. Due to this solution is not guaranteed to be locally optimal; a local search is performed in the second phase. The final result is simply the best solution found over all iterations.The construction phase can be described as a process which stepwise adds one element at a time to a partial (incomplete) solution. According to a greedy function, all elements are sorted, and the Restricted Candidate List (RCL) is constructed based on the order, and then from the list randomly select a element. In the second phase, a local search is initialized from these points, this iterative process is repeated until a termination criterion is met, and the best solution found over all iterations is taken as the result.RCL is created using a parameterα to restrict the size of this list. Candidate e∈C is sorted according to their greedy function value f(e). In a cardinality-based RCL, the latter is made up by the k top-ranked elements. In a value-based construction, the RCL consists of the elements in the set {e∈C:f*≤f(e)≤f*+α×(f*-f*)}, where f*=min(f(e):e∈C), f*=max{f(e):e∈C}, and α∈[0,1]. Since the best value for α is often difficult to determine, a random value is often assigned. The values for α adopted in the constructive heuristics are set using reactive strategies, which usually leads to better performance than using fixed values [22].In original bat algorithm, the bat individual randomly selects a certain range of frequency, its velocity is updated according to their selected frequency, and at last a new position is generated using its velocity and its own position. The idea is that the position of bat individual is updated by adjusting its frequency of sonic pulse. In this paper, the position updating of bat adopted GRASP to generate a new position, the frequency is used for restricting the size of CRL, frequency equivalent to parameterα in GRASP, and the frequency is variable value. In basic bat algorithm, the loudness Ld and the pulse emission rate r are updated according as the iterations proceed. As the loudness usually decreases, while the rate of pulse emission increases, it indicates that the bats approximate their prey (optimum solution). The pulse emission rate r is updated by
(7)r(t)=(1+exp(-5tmax×(t-tmax/2))-1),
where t denotes the tth generation and tmax is the maximal generation. The rate r is similar to Sigmoid function, and the frequency fr is determined according to the pulse emission rate r, which is represented by
(8)fr={1-max(0.2,min(0.8,r)),rand>r,max(0.2,min(0.8,r)),rand≤r,
where rand is a random number and 0.2 and 0.8 are experience parameters reference [18]. The frequency fr decreases gradually at firstly and then increases gradually, while the generation t increases. Figure 1 is the changing curve of rate r, Figure 2 is an example of frequency fr, and Algorithm 1 shows the pseudocode of greedy randomized construction with frequency fr.Algorithm 1:
Greedy_Randomized_Construction(fr).
S
←
R  rand_Chosen_Vertex(v);//    R is selected at random a vertex as initial solution
C
←
V
∖
v
; // The candidate set
C
is initialized
While
C
≠
⌀
do
IC←Evaluate_Incremental_Costs(S); //  the incremental costs are evaluated
icmin←min({ic∈IC});
icmax←max({ic∈IC});
RCL←{e∈C∣ic(e)≤icmin+fr×(icmax-icmin)};  //RCLis created
s←Select_Element (RCL); // a vertexsis randomly selected fromRCL
S←Obtain_Min_Solution(S,s); // the partial tour is updated by inserting the
vertexs
C←C∖{s}; // the candidate setCis updated
end
return
SFigure 1
Changing curve of rater.Figure 2
Changing curve of frequencyfr.The local search phase uses the 2-Opt heuristic for exchanging. Algorithm2 shows the local search procedure. The input parameter is an initial solution S obtained by the Greedy_Randomized_Construction procedure. The current route need is divided into m subroute according to the load, set the start, and end of suroute said to 0; for example, S={1,2,3,4,5,6,7}, subroute is 0→1→2→3→4→0, and 0→5→6→7→0.Algorithm 2:
Local_Search_Phase (S).
S
′
←
S
;
{
s
r
}
←
C
o
n
s
t
r
u
c
t
_Sub-route(S′);
for
each sub-route
s
r do
sr′←2-opt(sr); // carry out the 2-opt operation
end
S
′
←
C
o
n
s
t
r
u
c
t
_Indiviaual
(
{
s
r
′
}
);
if f(S′)<f(S)then
S←S′
end
return
S


## 3.3. Hybrid Bat Algorithm with Path Relinking

Path relinking was originally proposed by Glover [23]; Laguna and Martí [24] were the first to use path relinking within a GRASP strategy. Path relinking generates new solutions by exploring trajectories connecting an initial solution xs to an elite guiding solution xt. The path relinking procedure consists in selecting moves that introduce attributes contained in the guiding solution xt to the initial solution xs until the initial solution is completely transformed in the guiding solution xt. Path relinking may also be viewed as a constrained local search strategy applied to the initial solution xs. Furthermore, there are several alternatives that have been considered, which involve tradeoffs between computation time and solution quality. These alternatives include periodical relinking, forward relinking, backward relinking, back and forward relinking, mixed relinking, randomized relinking, and truncated relinking[18].One important issue in implementing a path relinking technique is the strategy to construct the elite set (ES). We adopted a fixed size elite set, and a solutionx is inserted into the ES as follows.A solutionx is always inserted into ES if it is not full. Otherwise, the generated solution x is inserted in ES only if its cost is better than the worst cost solution found in ES, and the worst cost solution is replaced by the solution x. Algorithm 3 shows the pseudocode for the algorithm to construct and maintain the elite set ES.Algorithm 3:
Elite_Set(x).
if cardinality(ES) ≠ps then  //psis the population size
E
S
←
E
S
∪
{
x
}
else
E
S
w
←
t
h
e
worst solution in
E
S
if f(x)<f(ESw)then
E
S
←
E
S
∖
{
E
S
w
}
ES←ES∪{x}
end
end
return
E
SAlgorithm4 shows the pseudocode of path relinking procedure; an instance is xs={1,2,3,4,5,6,7}, xt={1,3,4,2,5,6,7}, according to process of Algorithm 4, the first different element is position 2 in xs; after executing replace operation, the x1={1,3,4,2,5,6,7}; the second different element is position 3 in xs; after executing replace operation, the x2={1,2,4,3,5,6,7}; the third different element is position 4 in xs; after executing replace operation, the x3={1,4,3,2,5,6,7}; if the x2 is better than the xt, and the only, then return x2.Algorithm 4:
Path_Relinking (ES*, x).
x
s
←
x
; // initial solution
x
t
←
E
S
*
; //  guiding solution
f
(
x
t
)
←
f
(
E
S
*
);
Δ
←
D
i
f
f
e
r
e
n
c
e(xs,xt); // find out the difference position betweenxs xt
for i=1:  cardinality(Δ) do
j←Find_Position(xs,xt,Δi); // find the position of the element ofxtinΔi which inxs
xi←Replace(xs,Δi,j); // replace thejposition ofxswith the  element ofxsinΔi
xi←Replace(xs,xt,Δi);  //replace theΔiposition ofxswith the corresponding element ofxt
if
f
(
x
i
)
<
f
(
x
t
)
then
x
*
←
x
i
;
f(xt)←f(xi);
end
end
return
x
*


## 3.4. Subsequence and Single-Point Local Search

In this paper, the loudnessLdi of bat individual i is relative to its own fitness fiti, the better fitness, and the lower loudness. The loudness can be described by
(9)Ldi=(fiti-fit*+0.1)(fit*-fit*+0.1),
where constant 0.1 is used for avoiding the denominator is zero, and fiti is the fitness of individual i, fit* and fit* are the minimum and maximum fitness in current population, respectively. In HBA-PR, the loudness reflects the quality of individual. In this algorithm, there are two kinds of local search are embedded into HBA-PR to further improve the performance, random subsequence local search that and single-point local search. The random subsequence local search includes random subsequence inverse and random subsequence insert, and the single-point local search includes single-point insert and single-point swap.For random subsequence insert, firstly, randomly selected an origin of subsequence, and then randomly selected a length of subsequence which is less than length of individualS. Secondly, after determining the subsequence S1, randomly selected an insert point in remainder subsequence S2, S=S1∪S2; the S1 is inserted into S2 location in insert point. An example is shown in Figure 2. For random subsequence inverse, firstly, randomly selected a subsequence with random length, and then the inverse operation is performed. An example is shown in Figure 3.Figure 3
Random subsequence insert operation.For single-point swap, choose two different positions from a permutation randomly and swap them. For single-point insert, choose two different positions from a permutation randomly, and the element in first position is insert into the back of second element. Similarly, two instances are shown in Figures4 and 5.Figure 4
Random subsequence inverse operation.Figure 5
Single-point swap operation.In local search part, the random subsequence local search is performed before random single-point local search. The random subsequence insert and random subsequence inverse are preformed according to loudnessLd. In other words, if a random number is greater than the loudness Ldi, the random subsequence insert is performed; otherwise, the random subsequence inverse is performed. Similarly, the random single-point insert and random single-point swap are performed with the loudness Ld. If a random number is greater than the loudness Ldi, the insert operation is performed (Figure 6); otherwise, the swap operation is performed. Note that where local search is operated on the current optimal individual ES* in elite set ES, local search is performed for each individual. Algorithms 5 and 6 show the pseudocode of subsequence and single-point local search.Algorithm 5:
Sub-sequence_Local_Search (ES*).
for i=1:  ps do
if
rand
>
L
d
i then
x
←
S
u
b
-sequence_Insert(ES*); // perform random sub-sequence insert operation
else
x←Sub-sequence_Inverse(ES*); // perform random sub-sequence inverse operation
end
end
return
xAlgorithm 6:
Single-point_Local_Search (ES*).
for i=1:  psdo
if rand>Ldi then
x←Single-point_Insert(ES*);  //  perform random single-point insert operation
else
x←Single-point_Swap(ES*); // perform random ingle point swap operation
end
end
return
xFigure 6
Single-point insert operation.


## 3.5. HBA-PR Framework for CVRP

We propose in this work, to incorporate greedy randomized adaptive search procedure, path-relinking strategies, subsequence, and single-point local search to the bat algorithm by defining distinct ways to solve capacitated vehicle routing problem. This iterative process is repeated until the termination criterion is met; Algorithm7 shows the pseudocode of HBA-PR for CVRP.Algorithm 7:HBA-PR.
(1)Initialize theps,  bat population and other parameters;
(2)Evaluate fitness for each individual;
(3)ES←Elite_Set(x);
(4)while  t≤tmax  do
(5)Compute pulse emission rate by (7);
(6)Determine frequencyfrby (8);
(7)for i=1:  psdo
(8)x←Greedy_Randomized_Construction(fr);
(9)x←Local_Search_Phase(x);
(10)end
(11)Evaluate fitness for each individualx;
(12)ES←Elite_Set(x);
(13)ES*←Select_Best_Elite(ES); // select the best individual in elite setES
(14)for i=1:ps do
(15)x←Path_Relinking(ES*,  x);
(16)end
(17)Evaluate fitness for each individualx;
(18)ES←Elite_Set(x);
(19)Ldi←Compute_Loudness(f(x)); //Compute loudness of each individual by (9);
(20)ES*←Select_Best_Elite(ES);
(21)x←Sub-sequence_Local_Search(ES*)//carry out random sub-sequence local
search
(22)Evaluate fitness for each individualx;
(23)ES←Elite_Set(x);
(24)ES*←Select_Best_Elite(ES);
(25)x←Single-point_Local_Search(ES*)//carry out random single-point local search
(26)Evaluate fitness for each individualx;
(27)ES←Elite_Set(x);
(28)t=t+1
(29)end
(30)Output result and plot


## 4. Numerical Simulation Results and Comparisons

To test the performance of the proposed HBA-PR which is extensively investigated by a large number of experimental studies computational simulations are carried out with some well-studied problems taken from the webhttp://www.branchandcut.org/, a reference site which contains detailed information regarding a large number of benchmark instances. In this paper, 12 instances from three classes of benchmarks are selected. The first class is Augerat et al. Set A instances, the second class is Augerat et al. Set P instances, and the third class is Augerat et al. Set E instances. So far, these problems have been widely used as benchmarks to certify the performance of algorithms by many researchers.All computational experiments are conducted with MATLAB 2012a, and in our simulation, numerical experiments are run on a PC with AMD Athlon (tm) II X4 640 Processor 3.0 GHz and 2.0 GB memory. In the experiment, the termination criterion is set as maximum generationtmax=200. Each instance independently run15 times for comparison.


### 4.1. Parameter Analysis

In the subsection, parameters of HBA-PR are determined by experiments, and the impact of each parameter is analyzed. In particular, the HBA-PR has few parameters; we only need to test population size (ps) in HBA-PR. A small ps which may lead to insufficient population information is provided, and the diversity cannot be guaranteed. On the other side, a large one which indicates diversity is sufficient, but the computing time will increase and the precision of optimal solution may have lesser improvement. In order to evaluate the sensitivity of parameters ps, three benchmarks selected from different benchmark set are chosen to run 10 times. These benchmarks are A_n33_k5, E_n23_k3 and P_n19_k2, and the statistical result and convergence curves are shown in Figures 7, 8, 9, and 10.Figure 7
Box-and-whisker diagram of selected three instances.Figure 8
Convergence curve forA_n33_k5.Figure 9
Convergence curve forE_n23_k3.Figure 10
Convergence curve forP_n19_k2.The ordinate normalized fitness (log) in Figures8–10 is logarithm of normalized fitness; the aim is to show the convergence curves clearly. The normalization formula of fitness is (fit-fit*)/(fit*-fit*), where fit* is the best-known solution, and fit* is the initial fitness.Figure7 represents the relative error of test caseA_n33_k5, E_n23_k3, and P_n19_k2 after 10 times independent running, which shows the sensitivity of parameter ps. From the three test cases, the parameter ps can be determined. For A_n33_k5 and P_n19_k2, the performance of HBA-PR is better when ps=30. For E_n23_k3, that ps is equal to 50 is best, but the performance is good as well, while ps=30. From Figures 8–10, the information provided by population is sufficient when ps is greater than 20; however, the convergence rate is better when ps=30 among the three instances. Considering the tradeoff between the stability of algorithm and the rate of convergence, the parameter ps takes a compromise values, ps=30.


### 4.2. Comparisons of Simulation Results

In order to show the effective of HBA-PR, we carry out a simulation to compare HBA-PR with other state-of-art algorithms, that is, a parallel version of the classical Clarke and Wright Savings (CWS) heuristic and SR-GCWS (Simulation in Routing via the Generalized Clarke and Wright Savings heuristic) proposed by Juan et al. [9], CS-GRASP proposed by Zheng et al. [10]. Results of these simulations are summarized in Table 1, which contain the following information of each instance: vehicle capacity; tightness (demand/capacity); I_BKS is integral best-known solution (BKS) or “optimal” value according to the web http://www.branchandcut.org/; R_BSK is verifies real costs for the best-known solution according to [9]; CWS is the costs associated with the solution given by the parallel version of the CWS heuristic; SR-GCWS is the best solution obtained by SR-GCWS method; HBA-PR is our best solution, where “—” represents no records in the literature.Table 1
Comparison of results for Augerat et al. Set A, E, and P instances.
Instance
Capacity
Tightness
I_BKS
R_BSK
CWS
SR-GCWS
CS-GRASP
HBA-PR
A-n33-k5
100
0.89
661
662.76
712.05
662.11
662.1452
662.1101
A-n33-k6
100
0.9
742
742.83
776.26
742.69
743.5785
742.6933
A-n37-k5
100
0.81
669
672.59
707.26
672.47
672.4652
672.4652
A-n39-k6
100
0.88
831
833.20
863.08
833.20
—
835.2518
E-n23-k3
4500
0.75
569
—
—
—
569.7461
568.5625
E-n22-k4
6000
0.94
375
375.28
388.77
375.28
—
375.2798
E-n33-k4
8000
0.92
839
838.72
843.1
837.67
—
837.9253
E-n51-k5
160
0.97
521
524.94
584.64
524.61
—
524.6111
P-n19-k2
160
0.97
212
212.66
237.90
212.66
212.6569
212.6569
P-n20-k2
160
0.97
216
217.42
234.00
217.42
217.4156
217.4156
P-n22-k2
160
0.96
216
217.85
239.50
217.85
217.8522
217.8522
P-n51-k10
80
0.97
741
742.48
790.97
741.50
—
743.2648From the simulationresults obtained by testing HBA-PR, it demonstrates that using the proposed HBA-PR to solve the CVRP is effective, and the performance of HBA-PR is prominent. From Table1, all instances achieved a good quality solution, and solutions of 6 instances are better than the best-known solution (“optimal” value). HBA-PR has matched 10 of the 12 best-known solutions except for P_n51_k10 and A_n39_k6, and the average deviation from the real costs best-known solution is 0.057%. HBA-PR outperforms CWS for the 12 instances. Comparing with CS-GRASP and HBA-PR, HBA-PR outperforms CS-GRASP in terms of A_n33_k5, A_n33_k6, and E_n23_k3; the other instances have same results. The SR-GCWS and HBA-PR have similar results, and the gap is very small.Furthermore, Figure11 shows the best solution found so far by using our methodology for the A_n37_k5.vrp file, where the depot (using 1 instead of 0) is at the center. Analogously, Figures 12 and 13 show the best solution found by HBA-PR for the E-n51-k5.vrp file and P_n51_k10.vrp file. The convergence curve for some test problems has been shown in Figure 14. From Figure 14, for instance, A_n37_k5, it converges to an optimal solution after 43 generation, and it expends about 170 generations, for instance, E-n51-k5 andP_n51_k10 when algorithm converges to an optimal solution, which demonstrate that the HBA-PR has a faster convergence rate.Figure 11
Optimal routes ofA_n37_k5.Figure 12
Optimal routes ofE_n51_k5.Figure 13
Optimal routes ofP_n51_k10.Figure 14
Convergence curve of three instances.In general, the proposed HBA-PR can produce good solutions when compared with existing heuristics for solving the CVRP, and the convergence rate of HBA-PR is faster. These results seem to indicate that the hybrid bat algorithm with path relinking is an alternative to solve the capacitated vehicle routing problem.


## 4.1. Parameter Analysis

In the subsection, parameters of HBA-PR are determined by experiments, and the impact of each parameter is analyzed. In particular, the HBA-PR has few parameters; we only need to test population size (ps) in HBA-PR. A small ps which may lead to insufficient population information is provided, and the diversity cannot be guaranteed. On the other side, a large one which indicates diversity is sufficient, but the computing time will increase and the precision of optimal solution may have lesser improvement. In order to evaluate the sensitivity of parameters ps, three benchmarks selected from different benchmark set are chosen to run 10 times. These benchmarks are A_n33_k5, E_n23_k3 and P_n19_k2, and the statistical result and convergence curves are shown in Figures 7, 8, 9, and 10.Figure 7
Box-and-whisker diagram of selected three instances.Figure 8
Convergence curve forA_n33_k5.Figure 9
Convergence curve forE_n23_k3.Figure 10
Convergence curve forP_n19_k2.The ordinate normalized fitness (log) in Figures8–10 is logarithm of normalized fitness; the aim is to show the convergence curves clearly. The normalization formula of fitness is (fit-fit*)/(fit*-fit*), where fit* is the best-known solution, and fit* is the initial fitness.Figure7 represents the relative error of test caseA_n33_k5, E_n23_k3, and P_n19_k2 after 10 times independent running, which shows the sensitivity of parameter ps. From the three test cases, the parameter ps can be determined. For A_n33_k5 and P_n19_k2, the performance of HBA-PR is better when ps=30. For E_n23_k3, that ps is equal to 50 is best, but the performance is good as well, while ps=30. From Figures 8–10, the information provided by population is sufficient when ps is greater than 20; however, the convergence rate is better when ps=30 among the three instances. Considering the tradeoff between the stability of algorithm and the rate of convergence, the parameter ps takes a compromise values, ps=30.


## 4.2. Comparisons of Simulation Results

In order to show the effective of HBA-PR, we carry out a simulation to compare HBA-PR with other state-of-art algorithms, that is, a parallel version of the classical Clarke and Wright Savings (CWS) heuristic and SR-GCWS (Simulation in Routing via the Generalized Clarke and Wright Savings heuristic) proposed by Juan et al. [9], CS-GRASP proposed by Zheng et al. [10]. Results of these simulations are summarized in Table 1, which contain the following information of each instance: vehicle capacity; tightness (demand/capacity); I_BKS is integral best-known solution (BKS) or “optimal” value according to the web http://www.branchandcut.org/; R_BSK is verifies real costs for the best-known solution according to [9]; CWS is the costs associated with the solution given by the parallel version of the CWS heuristic; SR-GCWS is the best solution obtained by SR-GCWS method; HBA-PR is our best solution, where “—” represents no records in the literature.Table 1
Comparison of results for Augerat et al. Set A, E, and P instances.
Instance
Capacity
Tightness
I_BKS
R_BSK
CWS
SR-GCWS
CS-GRASP
HBA-PR
A-n33-k5
100
0.89
661
662.76
712.05
662.11
662.1452
662.1101
A-n33-k6
100
0.9
742
742.83
776.26
742.69
743.5785
742.6933
A-n37-k5
100
0.81
669
672.59
707.26
672.47
672.4652
672.4652
A-n39-k6
100
0.88
831
833.20
863.08
833.20
—
835.2518
E-n23-k3
4500
0.75
569
—
—
—
569.7461
568.5625
E-n22-k4
6000
0.94
375
375.28
388.77
375.28
—
375.2798
E-n33-k4
8000
0.92
839
838.72
843.1
837.67
—
837.9253
E-n51-k5
160
0.97
521
524.94
584.64
524.61
—
524.6111
P-n19-k2
160
0.97
212
212.66
237.90
212.66
212.6569
212.6569
P-n20-k2
160
0.97
216
217.42
234.00
217.42
217.4156
217.4156
P-n22-k2
160
0.96
216
217.85
239.50
217.85
217.8522
217.8522
P-n51-k10
80
0.97
741
742.48
790.97
741.50
—
743.2648From the simulationresults obtained by testing HBA-PR, it demonstrates that using the proposed HBA-PR to solve the CVRP is effective, and the performance of HBA-PR is prominent. From Table1, all instances achieved a good quality solution, and solutions of 6 instances are better than the best-known solution (“optimal” value). HBA-PR has matched 10 of the 12 best-known solutions except for P_n51_k10 and A_n39_k6, and the average deviation from the real costs best-known solution is 0.057%. HBA-PR outperforms CWS for the 12 instances. Comparing with CS-GRASP and HBA-PR, HBA-PR outperforms CS-GRASP in terms of A_n33_k5, A_n33_k6, and E_n23_k3; the other instances have same results. The SR-GCWS and HBA-PR have similar results, and the gap is very small.Furthermore, Figure11 shows the best solution found so far by using our methodology for the A_n37_k5.vrp file, where the depot (using 1 instead of 0) is at the center. Analogously, Figures 12 and 13 show the best solution found by HBA-PR for the E-n51-k5.vrp file and P_n51_k10.vrp file. The convergence curve for some test problems has been shown in Figure 14. From Figure 14, for instance, A_n37_k5, it converges to an optimal solution after 43 generation, and it expends about 170 generations, for instance, E-n51-k5 andP_n51_k10 when algorithm converges to an optimal solution, which demonstrate that the HBA-PR has a faster convergence rate.Figure 11
Optimal routes ofA_n37_k5.Figure 12
Optimal routes ofE_n51_k5.Figure 13
Optimal routes ofP_n51_k10.Figure 14
Convergence curve of three instances.In general, the proposed HBA-PR can produce good solutions when compared with existing heuristics for solving the CVRP, and the convergence rate of HBA-PR is faster. These results seem to indicate that the hybrid bat algorithm with path relinking is an alternative to solve the capacitated vehicle routing problem.


## 5. Conclusions

The capacitated vehicle routing problem is important in the fields of Operations Research, which is an NP-hard problem. Bat algorithm is a continuous metaheuristics; it cannot be used to solve the CVRP directly. In this paper, a hybrid bat algorithm with path relinking (HBA-PR) for solving the CVRP has been presented. This methodology, which does not require any particular fine tuning of parameters or configuration process, combines the classical greedy randomized adaptive search procedure (GRASP) with bat algorithm; the path relinking as an intensification strategy to explore local trajectories connecting elite solutions obtained by proposed algorithm; subsequence and single-point local search are effectively integrated into HBA-PR. Results show that our methodology is able to provide fine-quality solutions which can compete with the ones provided by some exact and heuristic approaches. Moreover, because of its simplicity and flexibility, we think that this methodology can easily be adapted to other variants of the vehicle routing problem and even to other combinatorial problems, for example, the vehicle routing problem with time window and traveling salesman problem, which is our further work.



---

*Source: 392789-2013-08-28.xml*