# A New Improved Quantum Evolution Algorithm with Local Search Procedure for Capacitated Vehicle Routing Problem

Authors: Ligang Cui; Lin Wang; Jie Deng; Jinlong Zhang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The capacitated vehicle routing problem (CVRP) is the most classical vehicle routing problem (VRP); many solution techniques are proposed to find its better answer. In this paper, a new improved quantum evolution algorithm (IQEA) with a mixed local search procedure is proposed for solving CVRPs. First, an IQEA with a double chain quantum chromosome, new quantum rotation schemes, and self-adaptive quantum Not gate is constructed to initialize and generate feasible solutions. Then, to further strengthen IQEA's searching ability, three local search procedures 1-1 exchange, 1-0 exchange, and 2-OPT, are adopted. Experiments on a small case have been conducted to analyze the sensitivity of main parameters and compare the performances of the IQEA with different local search strategies. Together with results from the testing of CVRP benchmarks, the superiorities of the proposed algorithm over the PSO, SR-1, and SR-2 have been demonstrated. At last, a profound analysis of the experimental results is presented and some suggestions on future researches are given.


---

## Body

## 1. Introduction

Vehicle routing problem (VRP) has been presented more than fifty years by Dantzig and Ramser [1], and it is still drawing many researchers’ attention [2, 3]. In a VRP, a vertex usually denotes the distribution center and other vertexes are used to denote the clients on demand; the distance between the vertexes is usually considered as the traveling distance or traveling cost; a classical VRP [4, 5] may include two principles: each client needs to be visited only once  and the starting point and terminal point are both at the distribution center. So the objective of a VRP is to minimize the total distances.It is well known that VRP is an NP-hard problem [6]. While CVRP is the most basic problem of VRPs and studying CVRP has a fundamental meaning to learn other advanced problems in this field and develop its solution methodologies [7]. Laporte [8] reviewed the VRP algorithms development from 1959 to 2009; these algorithms were classified as exact algorithm, classical heuristic algorithm, and metaheuristic algorithm. As it was stated in his paper, solving an NP-hard CVRP will require more computational time but the outcome may not be improved a lot.In recent years, the solving methods of VRPs are mainly on finding a proper intelligent or heuristic algorithm with problem-oriented solving mechanisms to find near optimal solutions [9]. Vera et al. [10] used the variable neighborhood search (VNS) heuristic algorithm solving the periodic VRP (PVRP) without time window. Barker and Ayechew [11] proposed a genetic algorithm for VRP and Rego and Alidaee [12] reviewed some of the most important tabu search heuristics for the VRPs. Marinakis et al. [13] proposed a hybrid particle swarm optimization algorithm for the VRP, while Szeto et al. [14] used a bee colony algorithm for the CVRP. These algorithms have been proven successfully in solving VRPs. They are widely accepted as  being more efficient than some exact algorithms [4] to some extent.In 1996, new evolutionary algorithm called quantum-inspired algorithm was proposed [15], which combined the quantum mechanism and the evolution computing characters together, so Quantum-inspired algorithm is also called the quantum evolution algorithm (QEA). QEA has been widely used in a variety of fields, such as SVM parameter optimization [16] and nuclear reactor reload optimization problem [17].QEA was born to solve the TSP; there have been only a few satisfactory results in dealing  with the VRP, such as Zhang et al. [18] and Zhang et al. [19]. Current QEA-based methods for solving VRPs have three main shortages. (1) The first one is the complex encoding process. Classical QEA encoding procession need to be converted among quantum encoding-binary encoding-real number encoding-integer encoding, which increases the computational complexity. Nevertheless, this conversion has lost half of quantum encoding carrying information. (2) Another disadvantage of classical QEA is its initialization largely based on the encoding method, which has separated the relations between visiting customer and candidate vehicles. (3) The updating operation is stylized as it relies on the lookup table. The non problem-oriented process also impacts QEA searching efficiency.Therefore, an improved quantum evolution algorithm (IQEA) with a mixed local search procedure called IQEA+ is proposed for solving CVRPs, in which a new double chain based encoding method to deal with the relations between vehicles and clients, a revised quantum rotation mechanism to simplify the updating processions, and an adaptive mutation operator are adopted. To further strengthen the searching abilities of IQEA, three local search procedures, 1-1 exchange, 1-0 exchange, and 2-OPT are also adopted.The remaining contents are organized as follows. In Section2, we review the fundamental knowledge and the improvements of QEA. Section 3 mainly discusses the encoding mechanism, the revised updating method, and mutation operations. Moreover, the framework of the IQEA is shown in this section. While the numerical experiments and results are presented in Section 4. In Section 5, conclusions and future research directions are provided.


## 2. Literature Review of QEA

### 2.1. Basic Knowledge of QEA

Quantum-inspired algorithm was born to solve the TSP [15]. QEA is also a probability algorithm, which includes some similar features of the evolutionary algorithm (EA), such as the representation of the individuals, the fitness functions, the updating processes, and the mutations [27]. However, it is not quite the same as a classical EA for it is mainly constructed with the concept of quantum computation [28] that reflects microparticle behaviors based on quantum mechanics. In QEA, an individual of a chromosome is usually expressed by a certain number of quantum bits (Q-bits), which is represented as a linear superposition mode. The state of a Q-bit can be expressed as follow:
(1)|Ψ〉=α|0〉+β|1〉,
where α and β are complex numbers. In QEA, they are called probability amplitudes of |0〉 and |1〉, respectively. The smallest unit of information stored in a Q-bit may be in |0〉, |1〉, or any superposition of the two states, while |α|2 and |β|2 denote the probabilities of the Q-bit in “0” state and in “1” state and they satisfy the normalization equation:
(2)|α|2+|β|2=1.Then, we can construct a population ofQ-bit individuals in the tth generation like that Q(t)={q1t,q2t,…,qnt}, where n is the size of the population and qtj is the jth chromosome of a certain population in the tth generation. Therefore, a string of n  Q-bits can be expressed as follows:
(3)qjt=[α1tα2t⋯αntβ1tβ2t⋯βnt],
where t is the tth generation, m is the length of qtj, j is the jth chromosome, and j=1,2,…,n. The expression can also be described as another form:
(4)qjt=[sin(θ)1tsin(θ)2t⋯sin(θ)ntcos(θ)1tcos(θ)2t⋯cos(θ)nt],
where θ∈[0,2π].A basic QEA procedure is as follows.(
1
) Set the maximum generation and the population size. Initialize the population Q(t) where t=0. A general method is as follows: let all 2×m×n probability amplitudes be equal to 1/2, which means that all the possible linear superposition states are endowed with the same probability 1/2m; that is,
(5)|Ψt=0qj〉=∑k=12m‍12m|sk〉,
where Sk is a binary string (x1,x2,…,xm) denotes the kth state of the quantum chromosome, and xi=0,1; i=1,2,…,m.(
2
) Observe the Q(t) to obtain its binary solution set P(t)=(x1t,x2t,…,xmt). Each xtj(j=1,2,…,n) is a binary string (X1t,X2t,…,Xmt) with a length of m and Xit=0 or 1(i=1,2,…,m) according to certain conditions: randomly generate a number p  which lies in [0,1] if p is larger than the square of probability amplitude of the ith Q-bit, then we take “1” in the binary state; otherwise, “0” is adopted in the binary state.(
3
) Convert the binary solutions to real number solutions corresponding to specific questions. Use a fitness function to evaluate every individual and store the best solutions; if the object reached, stop the algorithm; otherwise, turn to step 4.(
4
) Update and execute mutation operations. Different from recombination and mutation operations of conventional EA, quantum-gate (Q-gate) is adopted in QEA to update the current solutions. There are several kinds of Q-gates, such as the NOT gate, controlled-Not gate, Hadamard gate, R & Nϵ gate and Phase gate [29], but all of them must meet requirement of the Schrödinger equation; that is,
(6)UUT=I,
where U is a Q-gate and I is a unit matrix. Because QEA’s updating and mutation operations are mainly constructed based on a single Q-bit, the single quantum bit gate—quantum rotation gate—has been widely used. A rotation gate U(θi) as an updating operator is employed in updating a Q-bit [αi,βi]T:
(7)[αi′βi′]=U(θi)[αiβi]=[cos(θ)-sin(θ)sin(θ)cos(θ)][αiβi],
where θi is the rotation angle and its magnitude is obtained by the equation θi=s(αi,βi)Δθi. s(αi,βi) is the sign of the rotation angle whose lookup table is shown in Table 1. In Table 1, bi is the best solution B of ith Q-bit and xi is the binary solution X of the ith bit, respectively.Table 1
The lookup table of rotation angle.
x
i
b
i
f
(
x
i
)
≥
f
(
b
i
)
Δ
θ
i
s
(
α
i
,
β
i
)
α
i
β
i
>
0
α
i
β
i
<
0
α
i
=
0
β
i
=
0
0
0
False
0
0
0
0
0
0
0
True
0
0
0
0
0
0
1
False
0
0
0
0
0
0
1
True
0.05
π
−1
+1
±1
0
1
0
False
0.01
π
−1
+1
±1
0
1
0
True
0.025
π
+1
−1
0
±1
1
1
False
0.005
π
+1
−1
0
±1
1
1
True
0.025
π
+1
−1
0
±1While in the mutation operation, the quantum NOT gate is commonly used. The essence of the NOT gate is to swap the positions of aQ-bit. A NOT gate U is expressed as follows:
(8)U=[1001].
The mutation procedure is as follows: firstly, set a mutation probability p1 and generate a p2 in [0,1]; secondly, compare p1 with p2; if p1 is larger than p2, select a certain chromosome, from which a number of Q-bits are chosen to carry out the mutation operation; otherwise, do nothing.(
5
) Check whether the maximum generation is reached. If not, t=t+1, go back to step 2; otherwise, stop the iteration.The flow chart of the original QEA is shown in Figure1.Figure 1
Flow chart of original QEA.


### 2.2. The Development of QEA

Compared with other evolutionary algorithms, better performance of QEA has been showed in balancing the exploration and exploitation [22]. But, for all that, it shows great potential in solving some specific problems, especially for some intractable problems [30], such as in the quantum turing machine and quantum automata, quantum circuits, adiabatic quantum computation and topological quantum computation.


#### 2.2.1. Hybrid Quantum Evolution Algorithms (HQEA)

Combining quantum mechanism with other EA algorithms is widely accepted by researchers; here we can just call them the hybrid quantum evolution algorithms (HQEA). Zheng and Yamashiro [31] proposed a novel quantum differential evolutionary algorithm (QDEA) with a strategy of largest rotating angle value rule to determine job sequence based on job’s quantum information. Su et al. [32] used a quantum-inspired differential evolution (DE) algorithm to discover the classification rules. Similarly, Su and Yang [33] proposed DE/QDE to learn the Takagi-Sugeno (T-S) fuzzy model.Typically speaking, the quantum mechanism is the most favored to be integrated with the swarm algorithm. For instance, Sun et al. [34] proposed quantum delta potential well based particle swarm optimization algorithm on the hypothesis that the PSO system is a quantum system, each particle of which is of quantum state formulated by the wave function. To deal with the problem of prematurity of particle swarm optimization in global searching, a global convergence guaranteed method called quantum behaved PSO (QPSO) was developed and a vector evaluated QPSO was also presented [35].


#### 2.2.2. Improved QEA (IQEA)

Other researchers have paid attention to the improvement of the QEA itself performances, such as redesigning the encoding mechanism, optimizing the rotation angle, and setting properly rotate directions; we can just call these algorithms the improved QEA (IQEA).P. L. Li and S. Y. Li [21] presented a double chain quantum evolution algorithm (DCQEA), in which each quantum chromosome contained two gene chains and each gene chain represented a feasible solution. While in their another paper, a QEA algorithm under the Bloch coordinates of the Q-bits (BQEA) encoding is proposed [36]. Gao and Wang [23] adopted a Q-bit representation to solve the multiobjective optimization problems (MOOPs) and a chaos-based approach is suggested to initialize the population. Based on paralleled quantum chromosomes, Gu et al. [37] used a novel parallel quantum genetic algorithm for stochastic job shop scheduling. In order to improve the efficiency of genetic algorithm in solving the above question, a novel parallel quantum genetic algorithm was presented by Gu et al. [38].QEA is also a probabilistic algorithm; its solution spaces update is both guided by the best individual and the quantum rotation gate mechanism. Researchers have found that quantum rotation gate has significant impact on its convergence speed, fitness value, and robustness [22]. Zhang and Gao [20] have introduced two types of rotation angles, which can provide a good balance between “fine search” and “coarse search.” By referring to the gradient of fitness function, P. L. Li and S. Y. Li [21] presented a new quantum rotation gate. While Gao and Wang [23] proposed a chaos-based rotation gate which depended on the rank of the parent cell in the population and a chaos equation. Lee et al. [25] constructed a self-adaptive operator which would adjust the size of the grid and was very simple to implement. Some of the above improvement measures are listed in Table 2.Table 2
Different quantum rotation angles of QEA.
No.
Different quantum rotation angles
Authors
1
Type I:Δθj,it=sgn(αj,itβj,it)θ0exp(|βj,it|/|αj,it|)+γ), where θ0 is an initial rotation angle and γ is a scale parameter; type II: αj,it+1=αj,it+1/(fix(k/5)+1);βj,it+1=1-(αj,it+1)2, where fix(·) is a round function and k denotes the number of continuously invalid evolution alleles
Zhang and Gao [20]
2
θ
=
-
sgn
(
A
)
×
θ
0
×
(
∥
∇
fixmax
-
∇
fixmin
∥
/
∥
∇
fix
(
X
)
-
∇
fixmin
∥
), where A is sin(θ1-θ0), θ0 and θ1 are the initial and the current best rotation angle, ∥·∥ is a certain norm, and ∇fix(·) is the gradient function of fitness at X.
P. L. Li and S. Y. Li [21]
3
θ
i
=
x
Δ
θ
i
(
1
+
x
λ
i
); Δθi=∑n=1K‍(ωbin-ωi)/K; x=sgn(αiβi(bi-0.5)); λ0=f(Xi); λt+1=μλt(1-λt), where ωi represents the ith rotation angle and ωbin the nth best individual's ith rotation angle; the last two equations are the logistic equations
Hu and Wu [22]
4
θ
i
=
λ
·
Φ
k
j
·
Φ
k
j; Φkj=8Φk-1j(1-Φk-1j)-1; λ=γλ0((1-K)/K), where the second equation is a modified chaos logistic equation, λ0 is a control parameter, γ is a temporary parameter, and generally γ=1.
Gao and Wang [23]
5
Δ
θ
=
min
Δ
θ
+
(
max
Δ
θ
-
min
Δ
θ
)
(
ω
1
·
(
∥
x
-
b
∥
/
n
)
+
ω
2
·
(
|
f
(
x
)
-
f
(
b
)
|
/
max
{
f
(
x
)
,
f
(
b
)
}
)
), where maxΔθ and minΔθ represent the maximum and minimum rotation angle separately, ω1 and ω2 are two predetermined constant numbers, and n is the length of its binary chromosome.
Chen et al. [24]
6
θ
=
k
·
f
(
α
i
,
β
i
); k=π·exp(-t/itermax), where t is the current iteration of the algorithm and itermax is the max iteration of the algorithm.
Lee et al. [25]In the encoding stage, the quantum encoding is often transferred into a binary encoding of a certain evolutionary algorithm, for example, the encoding method of the genetic algorithm (GA). The fact is that there are a few papers trying to solve the VRPs or VRPs variants with QEA. However, their methods of the evolutionary process are mainly on the basic QEA. So, in Section3, a new improved QEA is introduced to deal with CVRPs with new rotation mechanisms.


## 2.1. Basic Knowledge of QEA

Quantum-inspired algorithm was born to solve the TSP [15]. QEA is also a probability algorithm, which includes some similar features of the evolutionary algorithm (EA), such as the representation of the individuals, the fitness functions, the updating processes, and the mutations [27]. However, it is not quite the same as a classical EA for it is mainly constructed with the concept of quantum computation [28] that reflects microparticle behaviors based on quantum mechanics. In QEA, an individual of a chromosome is usually expressed by a certain number of quantum bits (Q-bits), which is represented as a linear superposition mode. The state of a Q-bit can be expressed as follow:
(1)|Ψ〉=α|0〉+β|1〉,
where α and β are complex numbers. In QEA, they are called probability amplitudes of |0〉 and |1〉, respectively. The smallest unit of information stored in a Q-bit may be in |0〉, |1〉, or any superposition of the two states, while |α|2 and |β|2 denote the probabilities of the Q-bit in “0” state and in “1” state and they satisfy the normalization equation:
(2)|α|2+|β|2=1.Then, we can construct a population ofQ-bit individuals in the tth generation like that Q(t)={q1t,q2t,…,qnt}, where n is the size of the population and qtj is the jth chromosome of a certain population in the tth generation. Therefore, a string of n  Q-bits can be expressed as follows:
(3)qjt=[α1tα2t⋯αntβ1tβ2t⋯βnt],
where t is the tth generation, m is the length of qtj, j is the jth chromosome, and j=1,2,…,n. The expression can also be described as another form:
(4)qjt=[sin(θ)1tsin(θ)2t⋯sin(θ)ntcos(θ)1tcos(θ)2t⋯cos(θ)nt],
where θ∈[0,2π].A basic QEA procedure is as follows.(
1
) Set the maximum generation and the population size. Initialize the population Q(t) where t=0. A general method is as follows: let all 2×m×n probability amplitudes be equal to 1/2, which means that all the possible linear superposition states are endowed with the same probability 1/2m; that is,
(5)|Ψt=0qj〉=∑k=12m‍12m|sk〉,
where Sk is a binary string (x1,x2,…,xm) denotes the kth state of the quantum chromosome, and xi=0,1; i=1,2,…,m.(
2
) Observe the Q(t) to obtain its binary solution set P(t)=(x1t,x2t,…,xmt). Each xtj(j=1,2,…,n) is a binary string (X1t,X2t,…,Xmt) with a length of m and Xit=0 or 1(i=1,2,…,m) according to certain conditions: randomly generate a number p  which lies in [0,1] if p is larger than the square of probability amplitude of the ith Q-bit, then we take “1” in the binary state; otherwise, “0” is adopted in the binary state.(
3
) Convert the binary solutions to real number solutions corresponding to specific questions. Use a fitness function to evaluate every individual and store the best solutions; if the object reached, stop the algorithm; otherwise, turn to step 4.(
4
) Update and execute mutation operations. Different from recombination and mutation operations of conventional EA, quantum-gate (Q-gate) is adopted in QEA to update the current solutions. There are several kinds of Q-gates, such as the NOT gate, controlled-Not gate, Hadamard gate, R & Nϵ gate and Phase gate [29], but all of them must meet requirement of the Schrödinger equation; that is,
(6)UUT=I,
where U is a Q-gate and I is a unit matrix. Because QEA’s updating and mutation operations are mainly constructed based on a single Q-bit, the single quantum bit gate—quantum rotation gate—has been widely used. A rotation gate U(θi) as an updating operator is employed in updating a Q-bit [αi,βi]T:
(7)[αi′βi′]=U(θi)[αiβi]=[cos(θ)-sin(θ)sin(θ)cos(θ)][αiβi],
where θi is the rotation angle and its magnitude is obtained by the equation θi=s(αi,βi)Δθi. s(αi,βi) is the sign of the rotation angle whose lookup table is shown in Table 1. In Table 1, bi is the best solution B of ith Q-bit and xi is the binary solution X of the ith bit, respectively.Table 1
The lookup table of rotation angle.
x
i
b
i
f
(
x
i
)
≥
f
(
b
i
)
Δ
θ
i
s
(
α
i
,
β
i
)
α
i
β
i
>
0
α
i
β
i
<
0
α
i
=
0
β
i
=
0
0
0
False
0
0
0
0
0
0
0
True
0
0
0
0
0
0
1
False
0
0
0
0
0
0
1
True
0.05
π
−1
+1
±1
0
1
0
False
0.01
π
−1
+1
±1
0
1
0
True
0.025
π
+1
−1
0
±1
1
1
False
0.005
π
+1
−1
0
±1
1
1
True
0.025
π
+1
−1
0
±1While in the mutation operation, the quantum NOT gate is commonly used. The essence of the NOT gate is to swap the positions of aQ-bit. A NOT gate U is expressed as follows:
(8)U=[1001].
The mutation procedure is as follows: firstly, set a mutation probability p1 and generate a p2 in [0,1]; secondly, compare p1 with p2; if p1 is larger than p2, select a certain chromosome, from which a number of Q-bits are chosen to carry out the mutation operation; otherwise, do nothing.(
5
) Check whether the maximum generation is reached. If not, t=t+1, go back to step 2; otherwise, stop the iteration.The flow chart of the original QEA is shown in Figure1.Figure 1
Flow chart of original QEA.


## 2.2. The Development of QEA

Compared with other evolutionary algorithms, better performance of QEA has been showed in balancing the exploration and exploitation [22]. But, for all that, it shows great potential in solving some specific problems, especially for some intractable problems [30], such as in the quantum turing machine and quantum automata, quantum circuits, adiabatic quantum computation and topological quantum computation.


### 2.2.1. Hybrid Quantum Evolution Algorithms (HQEA)

Combining quantum mechanism with other EA algorithms is widely accepted by researchers; here we can just call them the hybrid quantum evolution algorithms (HQEA). Zheng and Yamashiro [31] proposed a novel quantum differential evolutionary algorithm (QDEA) with a strategy of largest rotating angle value rule to determine job sequence based on job’s quantum information. Su et al. [32] used a quantum-inspired differential evolution (DE) algorithm to discover the classification rules. Similarly, Su and Yang [33] proposed DE/QDE to learn the Takagi-Sugeno (T-S) fuzzy model.Typically speaking, the quantum mechanism is the most favored to be integrated with the swarm algorithm. For instance, Sun et al. [34] proposed quantum delta potential well based particle swarm optimization algorithm on the hypothesis that the PSO system is a quantum system, each particle of which is of quantum state formulated by the wave function. To deal with the problem of prematurity of particle swarm optimization in global searching, a global convergence guaranteed method called quantum behaved PSO (QPSO) was developed and a vector evaluated QPSO was also presented [35].


### 2.2.2. Improved QEA (IQEA)

Other researchers have paid attention to the improvement of the QEA itself performances, such as redesigning the encoding mechanism, optimizing the rotation angle, and setting properly rotate directions; we can just call these algorithms the improved QEA (IQEA).P. L. Li and S. Y. Li [21] presented a double chain quantum evolution algorithm (DCQEA), in which each quantum chromosome contained two gene chains and each gene chain represented a feasible solution. While in their another paper, a QEA algorithm under the Bloch coordinates of the Q-bits (BQEA) encoding is proposed [36]. Gao and Wang [23] adopted a Q-bit representation to solve the multiobjective optimization problems (MOOPs) and a chaos-based approach is suggested to initialize the population. Based on paralleled quantum chromosomes, Gu et al. [37] used a novel parallel quantum genetic algorithm for stochastic job shop scheduling. In order to improve the efficiency of genetic algorithm in solving the above question, a novel parallel quantum genetic algorithm was presented by Gu et al. [38].QEA is also a probabilistic algorithm; its solution spaces update is both guided by the best individual and the quantum rotation gate mechanism. Researchers have found that quantum rotation gate has significant impact on its convergence speed, fitness value, and robustness [22]. Zhang and Gao [20] have introduced two types of rotation angles, which can provide a good balance between “fine search” and “coarse search.” By referring to the gradient of fitness function, P. L. Li and S. Y. Li [21] presented a new quantum rotation gate. While Gao and Wang [23] proposed a chaos-based rotation gate which depended on the rank of the parent cell in the population and a chaos equation. Lee et al. [25] constructed a self-adaptive operator which would adjust the size of the grid and was very simple to implement. Some of the above improvement measures are listed in Table 2.Table 2
Different quantum rotation angles of QEA.
No.
Different quantum rotation angles
Authors
1
Type I:Δθj,it=sgn(αj,itβj,it)θ0exp(|βj,it|/|αj,it|)+γ), where θ0 is an initial rotation angle and γ is a scale parameter; type II: αj,it+1=αj,it+1/(fix(k/5)+1);βj,it+1=1-(αj,it+1)2, where fix(·) is a round function and k denotes the number of continuously invalid evolution alleles
Zhang and Gao [20]
2
θ
=
-
sgn
(
A
)
×
θ
0
×
(
∥
∇
fixmax
-
∇
fixmin
∥
/
∥
∇
fix
(
X
)
-
∇
fixmin
∥
), where A is sin(θ1-θ0), θ0 and θ1 are the initial and the current best rotation angle, ∥·∥ is a certain norm, and ∇fix(·) is the gradient function of fitness at X.
P. L. Li and S. Y. Li [21]
3
θ
i
=
x
Δ
θ
i
(
1
+
x
λ
i
); Δθi=∑n=1K‍(ωbin-ωi)/K; x=sgn(αiβi(bi-0.5)); λ0=f(Xi); λt+1=μλt(1-λt), where ωi represents the ith rotation angle and ωbin the nth best individual's ith rotation angle; the last two equations are the logistic equations
Hu and Wu [22]
4
θ
i
=
λ
·
Φ
k
j
·
Φ
k
j; Φkj=8Φk-1j(1-Φk-1j)-1; λ=γλ0((1-K)/K), where the second equation is a modified chaos logistic equation, λ0 is a control parameter, γ is a temporary parameter, and generally γ=1.
Gao and Wang [23]
5
Δ
θ
=
min
Δ
θ
+
(
max
Δ
θ
-
min
Δ
θ
)
(
ω
1
·
(
∥
x
-
b
∥
/
n
)
+
ω
2
·
(
|
f
(
x
)
-
f
(
b
)
|
/
max
{
f
(
x
)
,
f
(
b
)
}
)
), where maxΔθ and minΔθ represent the maximum and minimum rotation angle separately, ω1 and ω2 are two predetermined constant numbers, and n is the length of its binary chromosome.
Chen et al. [24]
6
θ
=
k
·
f
(
α
i
,
β
i
); k=π·exp(-t/itermax), where t is the current iteration of the algorithm and itermax is the max iteration of the algorithm.
Lee et al. [25]In the encoding stage, the quantum encoding is often transferred into a binary encoding of a certain evolutionary algorithm, for example, the encoding method of the genetic algorithm (GA). The fact is that there are a few papers trying to solve the VRPs or VRPs variants with QEA. However, their methods of the evolutionary process are mainly on the basic QEA. So, in Section3, a new improved QEA is introduced to deal with CVRPs with new rotation mechanisms.


## 2.2.1. Hybrid Quantum Evolution Algorithms (HQEA)

Combining quantum mechanism with other EA algorithms is widely accepted by researchers; here we can just call them the hybrid quantum evolution algorithms (HQEA). Zheng and Yamashiro [31] proposed a novel quantum differential evolutionary algorithm (QDEA) with a strategy of largest rotating angle value rule to determine job sequence based on job’s quantum information. Su et al. [32] used a quantum-inspired differential evolution (DE) algorithm to discover the classification rules. Similarly, Su and Yang [33] proposed DE/QDE to learn the Takagi-Sugeno (T-S) fuzzy model.Typically speaking, the quantum mechanism is the most favored to be integrated with the swarm algorithm. For instance, Sun et al. [34] proposed quantum delta potential well based particle swarm optimization algorithm on the hypothesis that the PSO system is a quantum system, each particle of which is of quantum state formulated by the wave function. To deal with the problem of prematurity of particle swarm optimization in global searching, a global convergence guaranteed method called quantum behaved PSO (QPSO) was developed and a vector evaluated QPSO was also presented [35].


## 2.2.2. Improved QEA (IQEA)

Other researchers have paid attention to the improvement of the QEA itself performances, such as redesigning the encoding mechanism, optimizing the rotation angle, and setting properly rotate directions; we can just call these algorithms the improved QEA (IQEA).P. L. Li and S. Y. Li [21] presented a double chain quantum evolution algorithm (DCQEA), in which each quantum chromosome contained two gene chains and each gene chain represented a feasible solution. While in their another paper, a QEA algorithm under the Bloch coordinates of the Q-bits (BQEA) encoding is proposed [36]. Gao and Wang [23] adopted a Q-bit representation to solve the multiobjective optimization problems (MOOPs) and a chaos-based approach is suggested to initialize the population. Based on paralleled quantum chromosomes, Gu et al. [37] used a novel parallel quantum genetic algorithm for stochastic job shop scheduling. In order to improve the efficiency of genetic algorithm in solving the above question, a novel parallel quantum genetic algorithm was presented by Gu et al. [38].QEA is also a probabilistic algorithm; its solution spaces update is both guided by the best individual and the quantum rotation gate mechanism. Researchers have found that quantum rotation gate has significant impact on its convergence speed, fitness value, and robustness [22]. Zhang and Gao [20] have introduced two types of rotation angles, which can provide a good balance between “fine search” and “coarse search.” By referring to the gradient of fitness function, P. L. Li and S. Y. Li [21] presented a new quantum rotation gate. While Gao and Wang [23] proposed a chaos-based rotation gate which depended on the rank of the parent cell in the population and a chaos equation. Lee et al. [25] constructed a self-adaptive operator which would adjust the size of the grid and was very simple to implement. Some of the above improvement measures are listed in Table 2.Table 2
Different quantum rotation angles of QEA.
No.
Different quantum rotation angles
Authors
1
Type I:Δθj,it=sgn(αj,itβj,it)θ0exp(|βj,it|/|αj,it|)+γ), where θ0 is an initial rotation angle and γ is a scale parameter; type II: αj,it+1=αj,it+1/(fix(k/5)+1);βj,it+1=1-(αj,it+1)2, where fix(·) is a round function and k denotes the number of continuously invalid evolution alleles
Zhang and Gao [20]
2
θ
=
-
sgn
(
A
)
×
θ
0
×
(
∥
∇
fixmax
-
∇
fixmin
∥
/
∥
∇
fix
(
X
)
-
∇
fixmin
∥
), where A is sin(θ1-θ0), θ0 and θ1 are the initial and the current best rotation angle, ∥·∥ is a certain norm, and ∇fix(·) is the gradient function of fitness at X.
P. L. Li and S. Y. Li [21]
3
θ
i
=
x
Δ
θ
i
(
1
+
x
λ
i
); Δθi=∑n=1K‍(ωbin-ωi)/K; x=sgn(αiβi(bi-0.5)); λ0=f(Xi); λt+1=μλt(1-λt), where ωi represents the ith rotation angle and ωbin the nth best individual's ith rotation angle; the last two equations are the logistic equations
Hu and Wu [22]
4
θ
i
=
λ
·
Φ
k
j
·
Φ
k
j; Φkj=8Φk-1j(1-Φk-1j)-1; λ=γλ0((1-K)/K), where the second equation is a modified chaos logistic equation, λ0 is a control parameter, γ is a temporary parameter, and generally γ=1.
Gao and Wang [23]
5
Δ
θ
=
min
Δ
θ
+
(
max
Δ
θ
-
min
Δ
θ
)
(
ω
1
·
(
∥
x
-
b
∥
/
n
)
+
ω
2
·
(
|
f
(
x
)
-
f
(
b
)
|
/
max
{
f
(
x
)
,
f
(
b
)
}
)
), where maxΔθ and minΔθ represent the maximum and minimum rotation angle separately, ω1 and ω2 are two predetermined constant numbers, and n is the length of its binary chromosome.
Chen et al. [24]
6
θ
=
k
·
f
(
α
i
,
β
i
); k=π·exp(-t/itermax), where t is the current iteration of the algorithm and itermax is the max iteration of the algorithm.
Lee et al. [25]In the encoding stage, the quantum encoding is often transferred into a binary encoding of a certain evolutionary algorithm, for example, the encoding method of the genetic algorithm (GA). The fact is that there are a few papers trying to solve the VRPs or VRPs variants with QEA. However, their methods of the evolutionary process are mainly on the basic QEA. So, in Section3, a new improved QEA is introduced to deal with CVRPs with new rotation mechanisms.


## 3. The Proposed IQEA with Local Search Procedure for CVRP

CVRP is usually formulated as follows: given there areK(k=1,2,…,K) candidate vehicles providing the services for L(i=1,2,…,L) customers and i = 0  represented the distribution center. The capacity of each vehicle is bk(k=1,2,…,K), while the demands of each customer are di(i=1,2,…,L). Transport distance/cost from customer i to customer j is ci,j. The CVRP model is
(9)minz=∑k=1K‍∑i=0L‍∑j=0L‍cijxi,j,k,(10)xi,j,k={1vehiclekvisitingfromitoj,0other,∀i,j,k,(11)yi,k={1customeriisservedbyk,0other,∀i,k,(12)∑i=1L‍diyi,k≤bk,∀k,(13)∑k=1K‍yi,k=1,∀i,(14)∑i=1L‍xi,j,k=yj,k,∀j,k,(15)∑j=1L‍xi,j,k=yi,k,∀i,k,(16)∑‍i,j∈S×Sxi,j,k≤|S|-1,S⊂{1,2,…,L},S≠Φ,∀k,
where (9) is the objective function; here we want to obtain the minimized total routing distances; xi,j,k in (10) and yi,k in (11) are both the {0,1} decision variables and xi,j,k means that, when the vehicle k visits the customer j from i, xi,j,k is equal to 1; otherwise, xi,j,k is 0 and yi,j,k means that the vehicle k will visit customer i; (12) is capacity constraints, which means that the maximum loading of vehicle k can not bypass bk; (13) is to guarantee that there is only has one vehicle visiting customer i; (14) and (15) are used to guarantee that one vehicle will serve only one customer and (16) is used to eliminate subcircle or subroute. The essence of the CVRP is to minimize the total distribution distances under the constraints of capacity.


### 3.1. The Encoding Mechanism and Solution Presentation

In the encoding stage of classical QEA, quantum codes need to be translated into binary codes, and then the binary codes are continued to be translated into real numbers, which result in encoding procedure that is too complex and trivial to operate. Therefore, a new encoding mechanism is proposed to skip the intermediate encoding steps and take full advantage of quantum information rather than only care about theQ-bit with a bigger probability amplitude. Further, a double chain quantum chromosome is also adopted to improve the searching efficiency. Suppose there are l clients to be served by one distribution center and k vehicles are asked to serve them. Obviously, there are k routes but we do not know which client will be served by which vehicle.


#### 3.1.1. Solution Representations

In this part, in order to construct routes,k groups were set up, and each of them was assigned a string of l numbers which lie in (0,1). Definitely, the vehicle with bigger number will get more chances to serve the client. The encoding steps are listed below.(1)
According to the number of vehicles and clients, a quantum chromosome (Q-chrome) length L is randomly generated, where L=l×k, l is the number of clients, and k is the number of vehicles.(2)
Divide theQ-chrome into k groups; each group has l members, which just correspond to the clients that need to be served.(3)
Match the given number (quantum amplitude) in each group to the potential client. Make sure the client number is in ascending order (or in descending order).(4)
Select a number (quantum amplitude) everyl length from one chain of Q-chrome, and store the smallest number and its position. The procedure will iterate l times till all l biggest numbers are selected out and all the clients are allocated to the corresponding vehicles. The whole solution representation procedures are shown in Figure 3.The essence of the four-step assignment is that, if the position of the biggest number is fixed, the represented vehicle will render the service and the client will win the chance to receive the service accordingly. Still, we should remember that there are 2 chains in a Q-chrome, so all the groups and all processing procedures in step 4 will be doubled.


#### 3.1.2. Construct the Vehicle Routes

As the visiting clients in each route have been determined, visiting orders of the clients in each route will be arranged. By far, the smallest numbers we have got in the former part have been all divided intok groups. In one group, the numbers are sorted in descending order. Because the clients’ sequence is related to the value of Q-bits in each group, the service (visiting) orders are sorted according to its corresponding value.To illustrate the encoding and solution presentation clearly, an example with 6 clients 2 vehicles is given below, so a chromosome length 12 (6×2) is randomly generated. The whole process that generates the initial routes and visiting sequences is shown in Figure 4. There are 6 steps to finish the whole routing construction process.From Figure4, we could get the sequence of clients assigned to two vehicles. If we define the distribution center as 0, then the complete routes are [0-6-2-0] and [0-1-5-3-4-0]. And, if we consider those operations in the lower chain at the same time, obviously, the searching efficiency of the IQEA will be doubled.


### 3.2. The Proposed Fitness Function

The minimized objective function is adopted as the fitness function, and a big punishment number is used to deal with unfeasible routes; the fitness function equation can be described as follows:(17)f={zfeasiblesolutions,z+Motherwise,
where f is the fitness function, z is the objective function of the CVRP model, see (9), and M is a big positive number (103–106). Equation (17) means that, if the solutions we obtained are all feasible, the fitness function is equal to the objective function; if there are unfeasible solutions, a big number will be added to the objective function to eliminate unfeasible routes.


### 3.3. The Updating Mechanism of IQEA

There are many QEA updating methods; see Table2. Zhang and Gao’s [20] rotation angle is not simple enough, and Lee et al. [25] presented a method that is easy to be understood and implemented, but it is not a problem oriented method. Therefore, a new quantum rotation gate mechanism is used here. As we all know, there are two parts of a whole updating procedure: the magnitudes of quantum rotation angle and the rotation directions which need to be considered in chromosome updating. So these two part procedures are provided below.


#### 3.3.1. Rotation Angle ofQ-Bit

In order to let currentQ-bits search towards the best Q-bits, the magnitude of current quantum amplitudes and the best quantum amplitudes should be considered together. Because the biggest number is acquired from the comparisons between different groups in a single Q-chrome, the magnitude of rotation angle is derived from its own chromosome. So the rotation angle calculation method is different from the method in the work of Hu and Wu [22], whose rotation angle mainly depends upon the current and the current best rotation angle. The proposed rotation angle is formulated as follows:(18a)θ-i=1k∑j=0k-1‍θi+j×L,(18b)Δθi=1k∑j=0k-1‍(θi-θ-i)2,wherei=1,2,…,L and θi+j×L is the (i+j×L)th rotation angle, θ-i is the mean value of each group, and Δθi is the variance of numbers in each group. It is quite different from traditional strategies with fixed value in Table 1 that each group has a different deviate value Δθi.


#### 3.3.2. Rotation Direction ofQ-Bit

The rotate direction is as important as the magnitude of rotation angle to QEA’s efficiency, which is used to control the current rotation angle to rotate towards the current best rotation angle. Ifα1 and β1 are the probability amplitudes of current optimal solutions, α2 and β2 are the probability amplitudes of current solutions, E1=tan(β1/α1), ϵ1=|E1|, θ1=tan-1(β1/α1), E2=tan(β2/θ2), ϵ2=|E2|, and θ2=tan-1(β2/α2), then one has the following.Theorem 1.
Rotation directions selection strategies,(1)
whenϵ1≥ϵ2 and if E2>0, the symbol (of θi) is positive, or else, if E2<0, the symbol is negative;(2)
whenϵ1<ϵ2 and if E2>0, the symbol (of θi) is negative, or else, if E2<0, the symbol is positive.For the proof of Theorem1, please refer to Appendix A.


### 3.4. Mutation Operation

Even though mutation process can increase the diversity of the searching population, effective strategies of mutation operation for QEA have not drawn too much attention. So, to some extent, the quantum Not gate may have played a role in local searches, but it still depends on the mutation rate. If the mutation probability is too small, only a small part ofQ-bits is exchanged and the local search ability is not fully fulfilled. However, too large mutation rate will result in unstable searches. So an adaptive mutation strategy is adopted: in the initial stages, the mutation rate is set large to guide the algorithm searching in the whole space and keep solutions diversity; in the last stages, a small mutation rate is adopted to keep the searching process stable and efficient.The pseudo codes of Algorithm1 of the adaptive mutation operation are given.Algorithm 1:Pseudo codes of adaptive mutation operation.
(1)Settings: Pm=1-gen/maxgen, N=ceil(Pm×popsize);
(2)iF  N> popsize then
(3)N = popsize;
(4)end if
(5)Randomly generate the positions Npos of N Q-chroms, Npos ∈[1, popsize];
(6)Randomly generate K different intergers in each group, K∈Mutpo;
(7)if  N≠0  then
(8)for  i=1 to N  do
(9)for  k=1 to K  do
(10)     Use the quantum Not-gate mutating the specificQ-bit;
(11)end for
(12)end for
(13)end ifIn the pseudo codes, themaxgen means the maximum iterations, gen is the current iterated generation, N is the number mutated, and Npos is the position of N  Q-chrom. K is the number of vehicles, which means each group will select one Q-bit mutated. So, as iteration goes on, the Pm is getting smaller and smaller, but the mutation position will not be fixed, which can help keep the searching towards the best solution in one aspect and maintain population diversities.


### 3.5. Local Search Algorithm

Mutation operation can provide diversities to the quantum population, but we are quite not sure that the solutions, have been found are the global optimized solutions. That is to say, the proposed IQEA still has a risk of being trapped in local minimum solutions. Therefore, some local search operators should be employed to increase algorithm convergence and qualities of solutions. To improve the local search abilities and jump out the local minimum points, many local search procedures and operators are proposed in kinds of evolution algorithms, such ask-exchange [39] and λ-exchange neighborhood [40] and ICROSS, IRP, IOPT, and O-OPT [41]. In order to improve the local optimal both for a single route and multiple routes, Ai and Kachitvichyanukul [7] used three local search procedures: 2-OPT, exchange 1-1, and 1-0 exchange. But they were just operated separately because the combination mechanism of the three procedures was not set up. So a new local search mechanism that combines procedures of 2-OPT, exchange 1-1, and 1-0 exchange together is proposed.


#### 3.5.1. Three Local Search Procedures

As depicted in the work of Ai and Kachitvichyanukul [7], 2-OPT is used to improve a single route that is determined by IQEA procedures by exchanging two different clients’ visiting sequence in that route: if distance of the new route is smaller than the former route, the new one is accepted. Continue the procedure until all the clients are all checked and a local optimal is obtained (see Figure 5: 2-OPT).The next two procedures 1-1 exchange and 1-0 exchange are adopted to deal with two adjacent candidate routes that have not reached their local optimum. 1-1 exchange is tried to interchange one client from the first route with another client in the other route (see Figure5: 1-0 exchange) when the distance of the two clients is lower than a predetermined number according to a certain condition. If the total distance of the modified routes is better than that of the former ones, the bad routes will be replaced by new ones.While in 1-0 exchange, a client is removed from the first route and inserted into the second route. If the total distance of new routes is lower than that of the old routes, the routes will be updated till all the clients in the first route are checked. The basic principal of 1-0 exchange is shown in Figure5: 1-0 exchange.More detailed information on the three local search procedures is depicted in AppendixB.


#### 3.5.2. Design Local Search Algorithm

Before the design of the scheme of local searches, we should learn that there is a basic principle to follow, which is also our aim to do local searches: to convert an unfeasible route into a feasible route, and improve the feasible route to a local optimal route. Hence, we presented three local search schemes using the three above operators (there still are other schemes, the explanations on why we did not adopt them will be given later), and they are listed in Table3.Table 3
Schemes of the local search operators.
Scheme
Description
(1)
Positive sequence: 1-1 exchange→ 1-0 exchange→ 2-OPT
(2)
Positive sequence: 1-0 exchange→ 1-1 exchange → 2-OPT
(3)
Mixed sequence of 1-0 exchange, 1-1 exchange, and 2-OPTThere are some characters about three local search operators:(1) 1-1 exchange procedure deals with two client vertexes corresponding to two routes, without impacting the sequence of the two routes; (2) 1-0 exchange procedure is to extract one vertex from one route and insert it into another route, while the relative positions of remaining vertexes in two candidate routes are without changes (which still need to be visited by corresponding vehicles); (3) 2-OPT procedure only arranges the vertexes sequence of one candidate route and has nothing to do with other routes. Moreover, if we execute the 2-OPT firstly and then do the other two procedures, we will face the problem that, once the local search operation finished, the improved routes we got may still not be a local optimum; it still needs another 2-OPT operation to find its local optimal. Therefore, the 2-OPT procedure must be put in the end of the execution sequence.How to arrange 1-0 exchange and 1-1 exchange procedures is another problem. From a problem-oriented point of view, the two procedures are all used to deal with two routes. On their processing objects, there exist 3 cases: both of the two routes are unfeasible, one is feasible and the other one is unfeasible, and both of the two routes are feasible. According to the local search principle stated above, an unfeasible route is improved to a feasible route. In dealing the routes are all unfeasible ones, the 1-0 exchange operator which is not completed for a vertex that moved from one route to another may not get two feasible routes, while the 1-1 exchange operator has more chances to improve the situation. In the case with of feasible route and one unfeasible route, 1-0 exchange can generate two feasible routes by removing a vertex from the unfeasible route and inserting it into the other route. In the third situation, both of the operators can be used to get better routes.The pseudo codes of local search Algorithm2 are as follows.Algorithm 2:Pseudo codes of local search algorithm.
(1)for  i=1 to K  do
(2)  Randomly select 2 candidate routes;
(3)  Calculate the actual carrying capacity (AC) of each vehicle.
(4)if AC < Burden Rating (BR) then
(5)  Run 1-1 exchange and 1-0 exchange 5 times.
(6)if solutions improved and constraints meet then
(7)  Run 2-OPT 5 times
(8)else
(9)  Roll back to other routes
(10)end if
(11)else
(12)Run 1-1 exchange and 1-0 exchange 5 times.
(13)while solutions improved do
(14)Run 2-OPT 5 times
(15)end while
(16)  end if
(17)end forThe local search procedure should be performed considering the coordination of the searching efficiency and searching effectiveness, so Algorithm2 is conducted after the QEA searching has been finished. In order to separate the former IQEA, we named the IQEA with local search scheme 3 the IQEA+ (Figure 2).Figure 2
Flow chart of IQEA+.Figure 3
Solution representations.Figure 4
Vehicle routes construct process.Figure 5
Three local search procedures.


## 3.1. The Encoding Mechanism and Solution Presentation

In the encoding stage of classical QEA, quantum codes need to be translated into binary codes, and then the binary codes are continued to be translated into real numbers, which result in encoding procedure that is too complex and trivial to operate. Therefore, a new encoding mechanism is proposed to skip the intermediate encoding steps and take full advantage of quantum information rather than only care about theQ-bit with a bigger probability amplitude. Further, a double chain quantum chromosome is also adopted to improve the searching efficiency. Suppose there are l clients to be served by one distribution center and k vehicles are asked to serve them. Obviously, there are k routes but we do not know which client will be served by which vehicle.


### 3.1.1. Solution Representations

In this part, in order to construct routes,k groups were set up, and each of them was assigned a string of l numbers which lie in (0,1). Definitely, the vehicle with bigger number will get more chances to serve the client. The encoding steps are listed below.(1)
According to the number of vehicles and clients, a quantum chromosome (Q-chrome) length L is randomly generated, where L=l×k, l is the number of clients, and k is the number of vehicles.(2)
Divide theQ-chrome into k groups; each group has l members, which just correspond to the clients that need to be served.(3)
Match the given number (quantum amplitude) in each group to the potential client. Make sure the client number is in ascending order (or in descending order).(4)
Select a number (quantum amplitude) everyl length from one chain of Q-chrome, and store the smallest number and its position. The procedure will iterate l times till all l biggest numbers are selected out and all the clients are allocated to the corresponding vehicles. The whole solution representation procedures are shown in Figure 3.The essence of the four-step assignment is that, if the position of the biggest number is fixed, the represented vehicle will render the service and the client will win the chance to receive the service accordingly. Still, we should remember that there are 2 chains in a Q-chrome, so all the groups and all processing procedures in step 4 will be doubled.


### 3.1.2. Construct the Vehicle Routes

As the visiting clients in each route have been determined, visiting orders of the clients in each route will be arranged. By far, the smallest numbers we have got in the former part have been all divided intok groups. In one group, the numbers are sorted in descending order. Because the clients’ sequence is related to the value of Q-bits in each group, the service (visiting) orders are sorted according to its corresponding value.To illustrate the encoding and solution presentation clearly, an example with 6 clients 2 vehicles is given below, so a chromosome length 12 (6×2) is randomly generated. The whole process that generates the initial routes and visiting sequences is shown in Figure 4. There are 6 steps to finish the whole routing construction process.From Figure4, we could get the sequence of clients assigned to two vehicles. If we define the distribution center as 0, then the complete routes are [0-6-2-0] and [0-1-5-3-4-0]. And, if we consider those operations in the lower chain at the same time, obviously, the searching efficiency of the IQEA will be doubled.


## 3.1.1. Solution Representations

In this part, in order to construct routes,k groups were set up, and each of them was assigned a string of l numbers which lie in (0,1). Definitely, the vehicle with bigger number will get more chances to serve the client. The encoding steps are listed below.(1)
According to the number of vehicles and clients, a quantum chromosome (Q-chrome) length L is randomly generated, where L=l×k, l is the number of clients, and k is the number of vehicles.(2)
Divide theQ-chrome into k groups; each group has l members, which just correspond to the clients that need to be served.(3)
Match the given number (quantum amplitude) in each group to the potential client. Make sure the client number is in ascending order (or in descending order).(4)
Select a number (quantum amplitude) everyl length from one chain of Q-chrome, and store the smallest number and its position. The procedure will iterate l times till all l biggest numbers are selected out and all the clients are allocated to the corresponding vehicles. The whole solution representation procedures are shown in Figure 3.The essence of the four-step assignment is that, if the position of the biggest number is fixed, the represented vehicle will render the service and the client will win the chance to receive the service accordingly. Still, we should remember that there are 2 chains in a Q-chrome, so all the groups and all processing procedures in step 4 will be doubled.


## 3.1.2. Construct the Vehicle Routes

As the visiting clients in each route have been determined, visiting orders of the clients in each route will be arranged. By far, the smallest numbers we have got in the former part have been all divided intok groups. In one group, the numbers are sorted in descending order. Because the clients’ sequence is related to the value of Q-bits in each group, the service (visiting) orders are sorted according to its corresponding value.To illustrate the encoding and solution presentation clearly, an example with 6 clients 2 vehicles is given below, so a chromosome length 12 (6×2) is randomly generated. The whole process that generates the initial routes and visiting sequences is shown in Figure 4. There are 6 steps to finish the whole routing construction process.From Figure4, we could get the sequence of clients assigned to two vehicles. If we define the distribution center as 0, then the complete routes are [0-6-2-0] and [0-1-5-3-4-0]. And, if we consider those operations in the lower chain at the same time, obviously, the searching efficiency of the IQEA will be doubled.


## 3.2. The Proposed Fitness Function

The minimized objective function is adopted as the fitness function, and a big punishment number is used to deal with unfeasible routes; the fitness function equation can be described as follows:(17)f={zfeasiblesolutions,z+Motherwise,
where f is the fitness function, z is the objective function of the CVRP model, see (9), and M is a big positive number (103–106). Equation (17) means that, if the solutions we obtained are all feasible, the fitness function is equal to the objective function; if there are unfeasible solutions, a big number will be added to the objective function to eliminate unfeasible routes.


## 3.3. The Updating Mechanism of IQEA

There are many QEA updating methods; see Table2. Zhang and Gao’s [20] rotation angle is not simple enough, and Lee et al. [25] presented a method that is easy to be understood and implemented, but it is not a problem oriented method. Therefore, a new quantum rotation gate mechanism is used here. As we all know, there are two parts of a whole updating procedure: the magnitudes of quantum rotation angle and the rotation directions which need to be considered in chromosome updating. So these two part procedures are provided below.


### 3.3.1. Rotation Angle ofQ-Bit

In order to let currentQ-bits search towards the best Q-bits, the magnitude of current quantum amplitudes and the best quantum amplitudes should be considered together. Because the biggest number is acquired from the comparisons between different groups in a single Q-chrome, the magnitude of rotation angle is derived from its own chromosome. So the rotation angle calculation method is different from the method in the work of Hu and Wu [22], whose rotation angle mainly depends upon the current and the current best rotation angle. The proposed rotation angle is formulated as follows:(18a)θ-i=1k∑j=0k-1‍θi+j×L,(18b)Δθi=1k∑j=0k-1‍(θi-θ-i)2,wherei=1,2,…,L and θi+j×L is the (i+j×L)th rotation angle, θ-i is the mean value of each group, and Δθi is the variance of numbers in each group. It is quite different from traditional strategies with fixed value in Table 1 that each group has a different deviate value Δθi.


### 3.3.2. Rotation Direction ofQ-Bit

The rotate direction is as important as the magnitude of rotation angle to QEA’s efficiency, which is used to control the current rotation angle to rotate towards the current best rotation angle. Ifα1 and β1 are the probability amplitudes of current optimal solutions, α2 and β2 are the probability amplitudes of current solutions, E1=tan(β1/α1), ϵ1=|E1|, θ1=tan-1(β1/α1), E2=tan(β2/θ2), ϵ2=|E2|, and θ2=tan-1(β2/α2), then one has the following.Theorem 1.
Rotation directions selection strategies,(1)
whenϵ1≥ϵ2 and if E2>0, the symbol (of θi) is positive, or else, if E2<0, the symbol is negative;(2)
whenϵ1<ϵ2 and if E2>0, the symbol (of θi) is negative, or else, if E2<0, the symbol is positive.For the proof of Theorem1, please refer to Appendix A.


## 3.3.1. Rotation Angle ofQ-Bit

In order to let currentQ-bits search towards the best Q-bits, the magnitude of current quantum amplitudes and the best quantum amplitudes should be considered together. Because the biggest number is acquired from the comparisons between different groups in a single Q-chrome, the magnitude of rotation angle is derived from its own chromosome. So the rotation angle calculation method is different from the method in the work of Hu and Wu [22], whose rotation angle mainly depends upon the current and the current best rotation angle. The proposed rotation angle is formulated as follows:(18a)θ-i=1k∑j=0k-1‍θi+j×L,(18b)Δθi=1k∑j=0k-1‍(θi-θ-i)2,wherei=1,2,…,L and θi+j×L is the (i+j×L)th rotation angle, θ-i is the mean value of each group, and Δθi is the variance of numbers in each group. It is quite different from traditional strategies with fixed value in Table 1 that each group has a different deviate value Δθi.


## 3.3.2. Rotation Direction ofQ-Bit

The rotate direction is as important as the magnitude of rotation angle to QEA’s efficiency, which is used to control the current rotation angle to rotate towards the current best rotation angle. Ifα1 and β1 are the probability amplitudes of current optimal solutions, α2 and β2 are the probability amplitudes of current solutions, E1=tan(β1/α1), ϵ1=|E1|, θ1=tan-1(β1/α1), E2=tan(β2/θ2), ϵ2=|E2|, and θ2=tan-1(β2/α2), then one has the following.Theorem 1.
Rotation directions selection strategies,(1)
whenϵ1≥ϵ2 and if E2>0, the symbol (of θi) is positive, or else, if E2<0, the symbol is negative;(2)
whenϵ1<ϵ2 and if E2>0, the symbol (of θi) is negative, or else, if E2<0, the symbol is positive.For the proof of Theorem1, please refer to Appendix A.


## 3.4. Mutation Operation

Even though mutation process can increase the diversity of the searching population, effective strategies of mutation operation for QEA have not drawn too much attention. So, to some extent, the quantum Not gate may have played a role in local searches, but it still depends on the mutation rate. If the mutation probability is too small, only a small part ofQ-bits is exchanged and the local search ability is not fully fulfilled. However, too large mutation rate will result in unstable searches. So an adaptive mutation strategy is adopted: in the initial stages, the mutation rate is set large to guide the algorithm searching in the whole space and keep solutions diversity; in the last stages, a small mutation rate is adopted to keep the searching process stable and efficient.The pseudo codes of Algorithm1 of the adaptive mutation operation are given.Algorithm 1:Pseudo codes of adaptive mutation operation.
(1)Settings: Pm=1-gen/maxgen, N=ceil(Pm×popsize);
(2)iF  N> popsize then
(3)N = popsize;
(4)end if
(5)Randomly generate the positions Npos of N Q-chroms, Npos ∈[1, popsize];
(6)Randomly generate K different intergers in each group, K∈Mutpo;
(7)if  N≠0  then
(8)for  i=1 to N  do
(9)for  k=1 to K  do
(10)     Use the quantum Not-gate mutating the specificQ-bit;
(11)end for
(12)end for
(13)end ifIn the pseudo codes, themaxgen means the maximum iterations, gen is the current iterated generation, N is the number mutated, and Npos is the position of N  Q-chrom. K is the number of vehicles, which means each group will select one Q-bit mutated. So, as iteration goes on, the Pm is getting smaller and smaller, but the mutation position will not be fixed, which can help keep the searching towards the best solution in one aspect and maintain population diversities.


## 3.5. Local Search Algorithm

Mutation operation can provide diversities to the quantum population, but we are quite not sure that the solutions, have been found are the global optimized solutions. That is to say, the proposed IQEA still has a risk of being trapped in local minimum solutions. Therefore, some local search operators should be employed to increase algorithm convergence and qualities of solutions. To improve the local search abilities and jump out the local minimum points, many local search procedures and operators are proposed in kinds of evolution algorithms, such ask-exchange [39] and λ-exchange neighborhood [40] and ICROSS, IRP, IOPT, and O-OPT [41]. In order to improve the local optimal both for a single route and multiple routes, Ai and Kachitvichyanukul [7] used three local search procedures: 2-OPT, exchange 1-1, and 1-0 exchange. But they were just operated separately because the combination mechanism of the three procedures was not set up. So a new local search mechanism that combines procedures of 2-OPT, exchange 1-1, and 1-0 exchange together is proposed.


### 3.5.1. Three Local Search Procedures

As depicted in the work of Ai and Kachitvichyanukul [7], 2-OPT is used to improve a single route that is determined by IQEA procedures by exchanging two different clients’ visiting sequence in that route: if distance of the new route is smaller than the former route, the new one is accepted. Continue the procedure until all the clients are all checked and a local optimal is obtained (see Figure 5: 2-OPT).The next two procedures 1-1 exchange and 1-0 exchange are adopted to deal with two adjacent candidate routes that have not reached their local optimum. 1-1 exchange is tried to interchange one client from the first route with another client in the other route (see Figure5: 1-0 exchange) when the distance of the two clients is lower than a predetermined number according to a certain condition. If the total distance of the modified routes is better than that of the former ones, the bad routes will be replaced by new ones.While in 1-0 exchange, a client is removed from the first route and inserted into the second route. If the total distance of new routes is lower than that of the old routes, the routes will be updated till all the clients in the first route are checked. The basic principal of 1-0 exchange is shown in Figure5: 1-0 exchange.More detailed information on the three local search procedures is depicted in AppendixB.


### 3.5.2. Design Local Search Algorithm

Before the design of the scheme of local searches, we should learn that there is a basic principle to follow, which is also our aim to do local searches: to convert an unfeasible route into a feasible route, and improve the feasible route to a local optimal route. Hence, we presented three local search schemes using the three above operators (there still are other schemes, the explanations on why we did not adopt them will be given later), and they are listed in Table3.Table 3
Schemes of the local search operators.
Scheme
Description
(1)
Positive sequence: 1-1 exchange→ 1-0 exchange→ 2-OPT
(2)
Positive sequence: 1-0 exchange→ 1-1 exchange → 2-OPT
(3)
Mixed sequence of 1-0 exchange, 1-1 exchange, and 2-OPTThere are some characters about three local search operators:(1) 1-1 exchange procedure deals with two client vertexes corresponding to two routes, without impacting the sequence of the two routes; (2) 1-0 exchange procedure is to extract one vertex from one route and insert it into another route, while the relative positions of remaining vertexes in two candidate routes are without changes (which still need to be visited by corresponding vehicles); (3) 2-OPT procedure only arranges the vertexes sequence of one candidate route and has nothing to do with other routes. Moreover, if we execute the 2-OPT firstly and then do the other two procedures, we will face the problem that, once the local search operation finished, the improved routes we got may still not be a local optimum; it still needs another 2-OPT operation to find its local optimal. Therefore, the 2-OPT procedure must be put in the end of the execution sequence.How to arrange 1-0 exchange and 1-1 exchange procedures is another problem. From a problem-oriented point of view, the two procedures are all used to deal with two routes. On their processing objects, there exist 3 cases: both of the two routes are unfeasible, one is feasible and the other one is unfeasible, and both of the two routes are feasible. According to the local search principle stated above, an unfeasible route is improved to a feasible route. In dealing the routes are all unfeasible ones, the 1-0 exchange operator which is not completed for a vertex that moved from one route to another may not get two feasible routes, while the 1-1 exchange operator has more chances to improve the situation. In the case with of feasible route and one unfeasible route, 1-0 exchange can generate two feasible routes by removing a vertex from the unfeasible route and inserting it into the other route. In the third situation, both of the operators can be used to get better routes.The pseudo codes of local search Algorithm2 are as follows.Algorithm 2:Pseudo codes of local search algorithm.
(1)for  i=1 to K  do
(2)  Randomly select 2 candidate routes;
(3)  Calculate the actual carrying capacity (AC) of each vehicle.
(4)if AC < Burden Rating (BR) then
(5)  Run 1-1 exchange and 1-0 exchange 5 times.
(6)if solutions improved and constraints meet then
(7)  Run 2-OPT 5 times
(8)else
(9)  Roll back to other routes
(10)end if
(11)else
(12)Run 1-1 exchange and 1-0 exchange 5 times.
(13)while solutions improved do
(14)Run 2-OPT 5 times
(15)end while
(16)  end if
(17)end forThe local search procedure should be performed considering the coordination of the searching efficiency and searching effectiveness, so Algorithm2 is conducted after the QEA searching has been finished. In order to separate the former IQEA, we named the IQEA with local search scheme 3 the IQEA+ (Figure 2).Figure 2
Flow chart of IQEA+.Figure 3
Solution representations.Figure 4
Vehicle routes construct process.Figure 5
Three local search procedures.


## 3.5.1. Three Local Search Procedures

As depicted in the work of Ai and Kachitvichyanukul [7], 2-OPT is used to improve a single route that is determined by IQEA procedures by exchanging two different clients’ visiting sequence in that route: if distance of the new route is smaller than the former route, the new one is accepted. Continue the procedure until all the clients are all checked and a local optimal is obtained (see Figure 5: 2-OPT).The next two procedures 1-1 exchange and 1-0 exchange are adopted to deal with two adjacent candidate routes that have not reached their local optimum. 1-1 exchange is tried to interchange one client from the first route with another client in the other route (see Figure5: 1-0 exchange) when the distance of the two clients is lower than a predetermined number according to a certain condition. If the total distance of the modified routes is better than that of the former ones, the bad routes will be replaced by new ones.While in 1-0 exchange, a client is removed from the first route and inserted into the second route. If the total distance of new routes is lower than that of the old routes, the routes will be updated till all the clients in the first route are checked. The basic principal of 1-0 exchange is shown in Figure5: 1-0 exchange.More detailed information on the three local search procedures is depicted in AppendixB.


## 3.5.2. Design Local Search Algorithm

Before the design of the scheme of local searches, we should learn that there is a basic principle to follow, which is also our aim to do local searches: to convert an unfeasible route into a feasible route, and improve the feasible route to a local optimal route. Hence, we presented three local search schemes using the three above operators (there still are other schemes, the explanations on why we did not adopt them will be given later), and they are listed in Table3.Table 3
Schemes of the local search operators.
Scheme
Description
(1)
Positive sequence: 1-1 exchange→ 1-0 exchange→ 2-OPT
(2)
Positive sequence: 1-0 exchange→ 1-1 exchange → 2-OPT
(3)
Mixed sequence of 1-0 exchange, 1-1 exchange, and 2-OPTThere are some characters about three local search operators:(1) 1-1 exchange procedure deals with two client vertexes corresponding to two routes, without impacting the sequence of the two routes; (2) 1-0 exchange procedure is to extract one vertex from one route and insert it into another route, while the relative positions of remaining vertexes in two candidate routes are without changes (which still need to be visited by corresponding vehicles); (3) 2-OPT procedure only arranges the vertexes sequence of one candidate route and has nothing to do with other routes. Moreover, if we execute the 2-OPT firstly and then do the other two procedures, we will face the problem that, once the local search operation finished, the improved routes we got may still not be a local optimum; it still needs another 2-OPT operation to find its local optimal. Therefore, the 2-OPT procedure must be put in the end of the execution sequence.How to arrange 1-0 exchange and 1-1 exchange procedures is another problem. From a problem-oriented point of view, the two procedures are all used to deal with two routes. On their processing objects, there exist 3 cases: both of the two routes are unfeasible, one is feasible and the other one is unfeasible, and both of the two routes are feasible. According to the local search principle stated above, an unfeasible route is improved to a feasible route. In dealing the routes are all unfeasible ones, the 1-0 exchange operator which is not completed for a vertex that moved from one route to another may not get two feasible routes, while the 1-1 exchange operator has more chances to improve the situation. In the case with of feasible route and one unfeasible route, 1-0 exchange can generate two feasible routes by removing a vertex from the unfeasible route and inserting it into the other route. In the third situation, both of the operators can be used to get better routes.The pseudo codes of local search Algorithm2 are as follows.Algorithm 2:Pseudo codes of local search algorithm.
(1)for  i=1 to K  do
(2)  Randomly select 2 candidate routes;
(3)  Calculate the actual carrying capacity (AC) of each vehicle.
(4)if AC < Burden Rating (BR) then
(5)  Run 1-1 exchange and 1-0 exchange 5 times.
(6)if solutions improved and constraints meet then
(7)  Run 2-OPT 5 times
(8)else
(9)  Roll back to other routes
(10)end if
(11)else
(12)Run 1-1 exchange and 1-0 exchange 5 times.
(13)while solutions improved do
(14)Run 2-OPT 5 times
(15)end while
(16)  end if
(17)end forThe local search procedure should be performed considering the coordination of the searching efficiency and searching effectiveness, so Algorithm2 is conducted after the QEA searching has been finished. In order to separate the former IQEA, we named the IQEA with local search scheme 3 the IQEA+ (Figure 2).Figure 2
Flow chart of IQEA+.Figure 3
Solution representations.Figure 4
Vehicle routes construct process.Figure 5
Three local search procedures.


## 4. Experiments and Discussions

The experiment is firstly conducted with a simple small case. For large-scale VRPs, parts of the standard benchmark are based on the work of Christofides et al. [42] and part of data sets are obtained from the website: http://branchandcut.org/VRP/data/. All the experiments were carried out in Windows system using Matlab 7.1 and performed on an Intel Pentium Dual-Core at 3.2 GHz and 2 GB RAM computer.


### 4.1. A Small CVRP Case and Parameters Analysis

#### 4.1.1.  A Contrastive Experiment: Small Case of CVRP

To illustrate the algorithm clearly, a small case of CVRP was introduced from Zhang et al. [18] and Chen et al. [43]. The content of the small case is as follows: 8 clients are served by 2 vehicles, the load capacity of which is 8 ton and the depot is node 0, and the optimal value of objective function is 67.5. The distance of every pair of clients and the demand of each client  are shown in Table 4. All algorithms were run 10 times, and the max iteration number is 100. For IQEA and IQEA+, the population size was set to 30 and the mutation probability was set to 0.025π. From the results in Table 5, we can see that the IQEA+ has better performance than IQEA.Table 4
Distance matrix and the demand quantities of each client.
Client
0
1
2
3
4
5
6
7
8
Demand
0
0
4
6
7.5
9
20
10
16
8
—
1
4
0
6.5
4
10
5
7.5
11
10
1
2
6
6.5
0
7.5
10
10
7.5
7.5
7.5
2
3
7.5
4
7.5
0
10
5
9
9
15
1
4
9
10
10
10
0
10
7.5
7.5
10
2
5
20
5
10
5
10
0
7
9
7.5
1
6
10
7.5
7.5
9
7.5
7
0
7
10
4
7
16
11
7.5
9
7.5
9
7
0
10
2
8
8
10
7.5
15
10
7.5
10
10
0
2Table 5
Comparison of different algorithms.
Orders
DGA
PSO-SA
IQEA
IQEA + 2-OPT
IQEA + 1-1 exchange
IQEA + 1-0 exchange
IQEA+
1
69.0
67.5
67.5
69
67.5
69.5
67.5
2
70.5
67.5
69.0
69
67.5
69
67.5
3
67.5
67.5
67.5
67.5
69
69
67.5
4
70.5
67.5
67.5
67.5
67.5
67.5
67.5
5
69.0
67.5
70.5
69
69.5
67.5
67.5
6
71.0
67.5
69.0
67.5
69
67.5
67.5
7
69.5
67.5
69.5
69
67.5
69.5
67.5
8
70.0
67.5
67.5
67.5
69.5
67.5
67.5
9
67.5
67.5
67.5
67.5
69
67.5
67.5
10
70.0
67.5
67.5
69
67.5
67.5
67.5
AVG
69.45
67.5
68.3
68.25
68.35
68.20
67.5


#### 4.1.2. Parameters Tuning Test for IQEA and IQEA+

Before further test of the algorithm on large scale CVRPs, we should first investigate the parameters of the IQEA+ such as population size, rotation angles, and mutation probabilities to find their impacts on the performances (Table11).(1)  Population Size Test. Experiments were performed on the above problem in Section 4.1.1 with the population size of 20, 30, and 50, respectively. Each case was run 30 times, the iteration step was set to 500, and the average of steps and the average of running time of IQEA+ that reached the best answer were given in Table 6. We can see that, as the population size increases, the running time cost increases linearly, especially when the population size goes to 30, and both of the running steps and running time reach a stable state.Table 6
Outputs in finding the best solutions under different population size.
Pop size
IQEA
IQEA + 2-OPT
IQEA + 1-1 exchange
IQEA + 1-0 exchange
IQEA+
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
20
0.72
113
1.69
271
167
263
1.42
220
1.39
209
30
0.95
102
2.09
213
1.89
194
1.52
156
2.26
234
50
2.12
140
2.61
170
3.42
226
2.20
138
3.59
232(2)  The Impacts of Rotation Angles and Mutation Probabilities. Each case was run 30 times, while the population size was set to 30. The average running time in finding the best solutions are listed in Tables 7 and 8 under different rotation angels and mutation probabilities, respectively. Table 7 shows that the IQEA and IQEA+ can find optimal solution faster when the rotation angle equals 0.01 pi and 0.025 pi., while our proposed rotation angle scheme can find the best solutions faster than any other approaches.Table 7
CPU times (s) under different rotation angles.
Rotation angle
IQEA
IQEA+
0.005 pi [26]
2.06
2.18
0.01 pi [26]
1.88
1.99
0.025 pi [26]
1.82
1.92
0.05 pi [26]
1.95
2.07
Std/5
2.04
2.16
Std/10
1.92
2.03
Std/(10 * selfsgn)
1.38
1.49Table 8
Optimal value under different mutation probabilities.
Probabilities
IQEA
IQEA+
0.01 [26]
71.8
70.0
0.025 [26]
70.7
68.9
0.05 [26]
72.3
70.2
0.1 [26]
69.5
68.4
0.2
69.1
68.4
0.5
72.3
70.9
Self-adapted P
67.8
67.5Table8 shows the average computation times under different mutation probabilities. When the mutation probability is smaller than the 0.01, IQEA+ has spent too much computation time  and is very hard to find the optimal solutions. We can find that the mutation probability is very important for IQEA+. It is partly because the client vertex in the candidate route is selected out by comparing between different Q-bits, and it will very easily fall into local minimum points with a small mutation probability. However, if the mutation probability was set too large (see the mutation probability = 0.5), IQEA+ cannot find the optimal solution for the unstable searching which makes the algorithm jump out of the spaces where the best solution may lie.Still, there is a need to ascertain how to arrange three local search procedures to test their best performance, so we select theA-n32-k5.vrp, the optimal value of which is 784 to test their performance. Each case was run 5 times, and the comparison results were shown in Table 9, where A is 2-OPT, B is 1-1 exchange and C is 1-0 exchange. We can easily find that the C-B-C-A scheme has the best result, so this scheme was used to test the CVRP data sets.Table 9
Different local search schemes for IQEA+.
No.
Scheme
Best objective
AVE
AVE.CPU time (s)
1
A
→
B
→
C
798
821.6
23.56
2
B
→
A
→
C
806
842
20.59
3
C
→
A
→
B
→
C
784
828.8
21.61
4
C
→
B
→
A
→
C
789
826
18.51
5
C
→
B
→
C
→
A
784
791.6
15.74


### 4.2. Large Scale CVRP Benchmarks Tests and Analysis

#### 4.2.1. Parameters Setting Analysis

(1)  Experiments of Medium Scale CVRPs. In order to make IQEA+ suitable for solving large scale CVRP benchmarks, we firstly test it by solving some medium scale cases. The same parameters settings as mentioned in Section 4.1 are still adopted. Each case is run 20 times to obtain the average value. The results are shown in Table 10.Table 10
Results of some medium-scale CVRP cases by IQEA+.
No.
Problem instance
Object function
CPU time (s)
Type
Clients
Vehicles
Optimal
Best
AVE
1
A-n32-k5
31
5
784
784
791.6
15.7
2
A-n34-k5
33
5
778
780
794.4
17.4
3
B-n31-k5
30
5
672
672
679.8
11.9
4
B-n34-k5
33
5
788
788
798.8
16.4
5
P-n19-k2
18
2
212
212
213.6
4.8
6
P-n20-k2
19
2
216
216
217.4
5.2
7
P-n40-k5
39
5
458
458
469.6
13.1
8
E-n23-k3
22
3
569
574
575.4
8.0
9
E-n22-k4
21
4
375
375
387.6
9.5
10
F-n45-k4
44
4
724
730
763.8
20.2Table 11
The impacts of the parameters of IQEA+ under two large cases.
Pn78k10
Pn101k4
Best value
AVE
AVE of CPU time
Best value
AVE
AVE of CPU time
Pop size
20
1507
1651
78.5
852
1303
85
30
1256
1295
72.2
788
1322
136
50
1358
1704
192
984
1306
219
Rotation angle
0.005 pi
1574
1628
233
1252
1337
133
0.01 pi
1562
1618
244
1145
1221
138
0.025 pi
1607
1582
250
1076
1123
136
0.05 pi
1469
1610
234
1021
1068
140
Std/5
1602
1557
181
904
937
137
Std/10
1487
1615
136
897
955
144
Δ1
1321
1465
132
885
992
136
Mutation probability
0.01
2016
2048
130
1138
1230
143
0.025
1970
2066
133
1062
1217
144
0.05
1885
1918
144
1007
1193
137
0.1
1775
1848
137
913
1097
148
0.2
1539
1679
146
1282
1326
140
0.5
1493
1523
138
1223
1302
141
Δ2
1385
1505
134
895
1202
141
Δ1 = Std/(10 * selfsgn), Δ2 = Self-adapted P.Results of Table10 show that the parameters setting of IQEA+ mentioned in Section 4.1 is relatively reasonable.(2)  Parameter Settings for IQEA+ for Large Scale CVRPs. In fact, related parameters of the most intelligent algorithms are confirmed according to the specific problems. Finding a universal parameter setting for intelligence algorithms is rather difficult. As discussed in Section 4.1, the parameters of IQEA+ are also confirmed by parameters setting test and/or experiences of similar problems.Related parameters of IQEA+ for solving larger-scale CVRP should be also discussed. To simplify our analysis, we only retest two large scale cases “Bn78k10” and “Pn101k”. Also, these two cases with different number of vehicles can also show some performance information on local searching procedures. Those experiments are tested with ten runs and fixed 500 iterations each without any break intervention. The computation results show us the following.(a)
Three important parameters of IQEA+, population size, rotation angle, and mutation probability, can cause great impacts of IQEA+’s performance.(b)
The main influences of population size lies in the searching time but not in the quality of the searching results. Different rotation angles and mutation probabilities have great influences of IQEA+ on finding the optimal values, but the searching times in different settings are quite approximative to each other.(c)
A large mutation probability seems to have better performance than that with a smaller mutation probability.(d)
The longer the length of a single patrolling route, the worse the searching ability of the local search procedure. In “Bn78k10” and “Pn101k” analyses, the average lengths of the patrolling routes of the former case and the latter case are 7.8 and 25.3, respectively. So, the former case has better optimal results.Finding a universal parameter setting of IQEA+ is rather challengeable and costly under several types of CVRP benchmarks. Therefore, based on results of IQEA+ for solving the medium scale CVRPs and two large scale CVRP cases and the experiences of [19, 22], the similar parameters’ setting of IQEA+ in Section 4.1 is also be applied to solve large scale CVRPs.


#### 4.2.2. Experiments for Large Scale CVRP Benchmarks

Large scale VRP benchmarks are also obtained from the website (http://branchandcut.org/VRP/data/), which were also introduced by Chen et al. [43] and Ai and Kachitvichyanukul [7]. There are sixteen CVRP cases, and the total number of clients is varying from 29 to 134, and the vehicle number is varying from 3 to 10. The settings of the IQEA+ for large scale CVRPs are as follows: the iteration is 500 times and the population size is 30. The setting of the rotation angle and mutation probability are the same as in Section 4.1. The searching performance of IQEA+ is compared with that from Chen’s et al. [43] DPSO algorithm and Ai and Kachitvichyanukul’s [7] SR-2 algorithm; all the comparison results are shown in Table 12.Table 12
Results of large-scale CVRPs with different algorithms.
Problem
Clients
Vehicles
Optimal value
Best objective function
AVE. IQEA+
CPU time (s)
Chen
SR-1
SR-2
IQEA+
Chen
SR-1
SR-2
IQEA+
An33k5
32
5
661
661
661
661
661
669.2
32
11
13
15.7
An46k7
45
7
914
914
914
914
914
970.6
129
19
23
72.0
An60k9
59
9
1354
1354
1366
1355
1363
1410.4
309
28
40
148.1
Bn35k5
34
5
955
955
955
955
955
967.6
38
12
14
21.2
Bn45k5
44
5
751
751
751
751
751
789.4
134
17
20
50.3
Bn68k9
67
9
1272
1272
1278
1274
1281
1349.8
344
33
50
65.4
Bn78k10
77
10
1221
1239
1239
1223
1256
1295.4
429
41
64
72.2
En30k3
29
3
534
534
541
534
534
551.2
28
11
16
13.4
En51k5
50
5
521
528
521
521
524
566
301
32
22
20.5
En76k7
75
7
682
688
691
682
703
739.4
527
38
60
150.0
Fn72k4
71
4
237
244
237
237
237
293.4
398
58
53
64.5
Fn135k7
134
7
1162
1215
1184
1162
1202
1414.2
1526
178
258
159.2
Mn101k10
100
10
820
824
821
820
956
1001.6
874
60
114
210.4
Mn121k7
120
7
1034
1038
1041
1036
1082
1207.4
1734
88
89
128.5
Pn76k4
75
4
593
602
599
594
594
612
496
51
48
64.9
Pn101k4
100
4
681
694
686
683
681
689.4
978
99
86
135.0


### 4.3. Discussions

In order to test the proposed IQEA+, firstly, it has been used to solve a simple CVRP to find different performances of the IQEA and IQEA+; then we test different body parts of IQEA to best play out its advantages. Three local search procedures combination schemes were also examined, and all of them were used to test the CVRP dataset.Experiments results show that the proposed IQEA+ has potential in solving the CVRP to some extent with relatively high precision and speed; see Table10 and several cases in Table 12. However, experiments also show some useful findings.(1)
IQEA+ outperforms DPSO, SR-1, and SR-2 in searching precision but not in the searching efficiency; see casePn101k4 in Table 12;(2)
Sometimes, IQEA+ does not show many advantages on searching performance; see caseFn135k7 in Table 12.(3)
IQEA can search quickly, while the local search procedures take more than a half of the total computing time, especially when problem scale becomes larger.Three main reasons can explain our findings. (a) There is a weakness of the IQEA+ in the length of quantum chromosome which originated from the encoding mechanism. When the client vertexes and vehicle number increase, especially as the length of quantum chromosome reached to thousand level, algorithm searching efficiency becomes very low. (b) The tightness (demand/capacity) of demand and capacity is another important factor that impacts the searching performance; see caseFn72k4 in Table 12. Different tightness may result in very different results: the more intensive are the tightness and the worse is the results, and the more intensive are the tightness, the greater chances of unfeasible routes are generated out. (c) Local search procedures are used to deal with the point-to-point operation and evaluation, which make it cost too much time.


## 4.1. A Small CVRP Case and Parameters Analysis

### 4.1.1.  A Contrastive Experiment: Small Case of CVRP

To illustrate the algorithm clearly, a small case of CVRP was introduced from Zhang et al. [18] and Chen et al. [43]. The content of the small case is as follows: 8 clients are served by 2 vehicles, the load capacity of which is 8 ton and the depot is node 0, and the optimal value of objective function is 67.5. The distance of every pair of clients and the demand of each client  are shown in Table 4. All algorithms were run 10 times, and the max iteration number is 100. For IQEA and IQEA+, the population size was set to 30 and the mutation probability was set to 0.025π. From the results in Table 5, we can see that the IQEA+ has better performance than IQEA.Table 4
Distance matrix and the demand quantities of each client.
Client
0
1
2
3
4
5
6
7
8
Demand
0
0
4
6
7.5
9
20
10
16
8
—
1
4
0
6.5
4
10
5
7.5
11
10
1
2
6
6.5
0
7.5
10
10
7.5
7.5
7.5
2
3
7.5
4
7.5
0
10
5
9
9
15
1
4
9
10
10
10
0
10
7.5
7.5
10
2
5
20
5
10
5
10
0
7
9
7.5
1
6
10
7.5
7.5
9
7.5
7
0
7
10
4
7
16
11
7.5
9
7.5
9
7
0
10
2
8
8
10
7.5
15
10
7.5
10
10
0
2Table 5
Comparison of different algorithms.
Orders
DGA
PSO-SA
IQEA
IQEA + 2-OPT
IQEA + 1-1 exchange
IQEA + 1-0 exchange
IQEA+
1
69.0
67.5
67.5
69
67.5
69.5
67.5
2
70.5
67.5
69.0
69
67.5
69
67.5
3
67.5
67.5
67.5
67.5
69
69
67.5
4
70.5
67.5
67.5
67.5
67.5
67.5
67.5
5
69.0
67.5
70.5
69
69.5
67.5
67.5
6
71.0
67.5
69.0
67.5
69
67.5
67.5
7
69.5
67.5
69.5
69
67.5
69.5
67.5
8
70.0
67.5
67.5
67.5
69.5
67.5
67.5
9
67.5
67.5
67.5
67.5
69
67.5
67.5
10
70.0
67.5
67.5
69
67.5
67.5
67.5
AVG
69.45
67.5
68.3
68.25
68.35
68.20
67.5


### 4.1.2. Parameters Tuning Test for IQEA and IQEA+

Before further test of the algorithm on large scale CVRPs, we should first investigate the parameters of the IQEA+ such as population size, rotation angles, and mutation probabilities to find their impacts on the performances (Table11).(1)  Population Size Test. Experiments were performed on the above problem in Section 4.1.1 with the population size of 20, 30, and 50, respectively. Each case was run 30 times, the iteration step was set to 500, and the average of steps and the average of running time of IQEA+ that reached the best answer were given in Table 6. We can see that, as the population size increases, the running time cost increases linearly, especially when the population size goes to 30, and both of the running steps and running time reach a stable state.Table 6
Outputs in finding the best solutions under different population size.
Pop size
IQEA
IQEA + 2-OPT
IQEA + 1-1 exchange
IQEA + 1-0 exchange
IQEA+
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
20
0.72
113
1.69
271
167
263
1.42
220
1.39
209
30
0.95
102
2.09
213
1.89
194
1.52
156
2.26
234
50
2.12
140
2.61
170
3.42
226
2.20
138
3.59
232(2)  The Impacts of Rotation Angles and Mutation Probabilities. Each case was run 30 times, while the population size was set to 30. The average running time in finding the best solutions are listed in Tables 7 and 8 under different rotation angels and mutation probabilities, respectively. Table 7 shows that the IQEA and IQEA+ can find optimal solution faster when the rotation angle equals 0.01 pi and 0.025 pi., while our proposed rotation angle scheme can find the best solutions faster than any other approaches.Table 7
CPU times (s) under different rotation angles.
Rotation angle
IQEA
IQEA+
0.005 pi [26]
2.06
2.18
0.01 pi [26]
1.88
1.99
0.025 pi [26]
1.82
1.92
0.05 pi [26]
1.95
2.07
Std/5
2.04
2.16
Std/10
1.92
2.03
Std/(10 * selfsgn)
1.38
1.49Table 8
Optimal value under different mutation probabilities.
Probabilities
IQEA
IQEA+
0.01 [26]
71.8
70.0
0.025 [26]
70.7
68.9
0.05 [26]
72.3
70.2
0.1 [26]
69.5
68.4
0.2
69.1
68.4
0.5
72.3
70.9
Self-adapted P
67.8
67.5Table8 shows the average computation times under different mutation probabilities. When the mutation probability is smaller than the 0.01, IQEA+ has spent too much computation time  and is very hard to find the optimal solutions. We can find that the mutation probability is very important for IQEA+. It is partly because the client vertex in the candidate route is selected out by comparing between different Q-bits, and it will very easily fall into local minimum points with a small mutation probability. However, if the mutation probability was set too large (see the mutation probability = 0.5), IQEA+ cannot find the optimal solution for the unstable searching which makes the algorithm jump out of the spaces where the best solution may lie.Still, there is a need to ascertain how to arrange three local search procedures to test their best performance, so we select theA-n32-k5.vrp, the optimal value of which is 784 to test their performance. Each case was run 5 times, and the comparison results were shown in Table 9, where A is 2-OPT, B is 1-1 exchange and C is 1-0 exchange. We can easily find that the C-B-C-A scheme has the best result, so this scheme was used to test the CVRP data sets.Table 9
Different local search schemes for IQEA+.
No.
Scheme
Best objective
AVE
AVE.CPU time (s)
1
A
→
B
→
C
798
821.6
23.56
2
B
→
A
→
C
806
842
20.59
3
C
→
A
→
B
→
C
784
828.8
21.61
4
C
→
B
→
A
→
C
789
826
18.51
5
C
→
B
→
C
→
A
784
791.6
15.74


## 4.1.1.  A Contrastive Experiment: Small Case of CVRP

To illustrate the algorithm clearly, a small case of CVRP was introduced from Zhang et al. [18] and Chen et al. [43]. The content of the small case is as follows: 8 clients are served by 2 vehicles, the load capacity of which is 8 ton and the depot is node 0, and the optimal value of objective function is 67.5. The distance of every pair of clients and the demand of each client  are shown in Table 4. All algorithms were run 10 times, and the max iteration number is 100. For IQEA and IQEA+, the population size was set to 30 and the mutation probability was set to 0.025π. From the results in Table 5, we can see that the IQEA+ has better performance than IQEA.Table 4
Distance matrix and the demand quantities of each client.
Client
0
1
2
3
4
5
6
7
8
Demand
0
0
4
6
7.5
9
20
10
16
8
—
1
4
0
6.5
4
10
5
7.5
11
10
1
2
6
6.5
0
7.5
10
10
7.5
7.5
7.5
2
3
7.5
4
7.5
0
10
5
9
9
15
1
4
9
10
10
10
0
10
7.5
7.5
10
2
5
20
5
10
5
10
0
7
9
7.5
1
6
10
7.5
7.5
9
7.5
7
0
7
10
4
7
16
11
7.5
9
7.5
9
7
0
10
2
8
8
10
7.5
15
10
7.5
10
10
0
2Table 5
Comparison of different algorithms.
Orders
DGA
PSO-SA
IQEA
IQEA + 2-OPT
IQEA + 1-1 exchange
IQEA + 1-0 exchange
IQEA+
1
69.0
67.5
67.5
69
67.5
69.5
67.5
2
70.5
67.5
69.0
69
67.5
69
67.5
3
67.5
67.5
67.5
67.5
69
69
67.5
4
70.5
67.5
67.5
67.5
67.5
67.5
67.5
5
69.0
67.5
70.5
69
69.5
67.5
67.5
6
71.0
67.5
69.0
67.5
69
67.5
67.5
7
69.5
67.5
69.5
69
67.5
69.5
67.5
8
70.0
67.5
67.5
67.5
69.5
67.5
67.5
9
67.5
67.5
67.5
67.5
69
67.5
67.5
10
70.0
67.5
67.5
69
67.5
67.5
67.5
AVG
69.45
67.5
68.3
68.25
68.35
68.20
67.5


## 4.1.2. Parameters Tuning Test for IQEA and IQEA+

Before further test of the algorithm on large scale CVRPs, we should first investigate the parameters of the IQEA+ such as population size, rotation angles, and mutation probabilities to find their impacts on the performances (Table11).(1)  Population Size Test. Experiments were performed on the above problem in Section 4.1.1 with the population size of 20, 30, and 50, respectively. Each case was run 30 times, the iteration step was set to 500, and the average of steps and the average of running time of IQEA+ that reached the best answer were given in Table 6. We can see that, as the population size increases, the running time cost increases linearly, especially when the population size goes to 30, and both of the running steps and running time reach a stable state.Table 6
Outputs in finding the best solutions under different population size.
Pop size
IQEA
IQEA + 2-OPT
IQEA + 1-1 exchange
IQEA + 1-0 exchange
IQEA+
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
Time (s)
Steps
20
0.72
113
1.69
271
167
263
1.42
220
1.39
209
30
0.95
102
2.09
213
1.89
194
1.52
156
2.26
234
50
2.12
140
2.61
170
3.42
226
2.20
138
3.59
232(2)  The Impacts of Rotation Angles and Mutation Probabilities. Each case was run 30 times, while the population size was set to 30. The average running time in finding the best solutions are listed in Tables 7 and 8 under different rotation angels and mutation probabilities, respectively. Table 7 shows that the IQEA and IQEA+ can find optimal solution faster when the rotation angle equals 0.01 pi and 0.025 pi., while our proposed rotation angle scheme can find the best solutions faster than any other approaches.Table 7
CPU times (s) under different rotation angles.
Rotation angle
IQEA
IQEA+
0.005 pi [26]
2.06
2.18
0.01 pi [26]
1.88
1.99
0.025 pi [26]
1.82
1.92
0.05 pi [26]
1.95
2.07
Std/5
2.04
2.16
Std/10
1.92
2.03
Std/(10 * selfsgn)
1.38
1.49Table 8
Optimal value under different mutation probabilities.
Probabilities
IQEA
IQEA+
0.01 [26]
71.8
70.0
0.025 [26]
70.7
68.9
0.05 [26]
72.3
70.2
0.1 [26]
69.5
68.4
0.2
69.1
68.4
0.5
72.3
70.9
Self-adapted P
67.8
67.5Table8 shows the average computation times under different mutation probabilities. When the mutation probability is smaller than the 0.01, IQEA+ has spent too much computation time  and is very hard to find the optimal solutions. We can find that the mutation probability is very important for IQEA+. It is partly because the client vertex in the candidate route is selected out by comparing between different Q-bits, and it will very easily fall into local minimum points with a small mutation probability. However, if the mutation probability was set too large (see the mutation probability = 0.5), IQEA+ cannot find the optimal solution for the unstable searching which makes the algorithm jump out of the spaces where the best solution may lie.Still, there is a need to ascertain how to arrange three local search procedures to test their best performance, so we select theA-n32-k5.vrp, the optimal value of which is 784 to test their performance. Each case was run 5 times, and the comparison results were shown in Table 9, where A is 2-OPT, B is 1-1 exchange and C is 1-0 exchange. We can easily find that the C-B-C-A scheme has the best result, so this scheme was used to test the CVRP data sets.Table 9
Different local search schemes for IQEA+.
No.
Scheme
Best objective
AVE
AVE.CPU time (s)
1
A
→
B
→
C
798
821.6
23.56
2
B
→
A
→
C
806
842
20.59
3
C
→
A
→
B
→
C
784
828.8
21.61
4
C
→
B
→
A
→
C
789
826
18.51
5
C
→
B
→
C
→
A
784
791.6
15.74


## 4.2. Large Scale CVRP Benchmarks Tests and Analysis

### 4.2.1. Parameters Setting Analysis

(1)  Experiments of Medium Scale CVRPs. In order to make IQEA+ suitable for solving large scale CVRP benchmarks, we firstly test it by solving some medium scale cases. The same parameters settings as mentioned in Section 4.1 are still adopted. Each case is run 20 times to obtain the average value. The results are shown in Table 10.Table 10
Results of some medium-scale CVRP cases by IQEA+.
No.
Problem instance
Object function
CPU time (s)
Type
Clients
Vehicles
Optimal
Best
AVE
1
A-n32-k5
31
5
784
784
791.6
15.7
2
A-n34-k5
33
5
778
780
794.4
17.4
3
B-n31-k5
30
5
672
672
679.8
11.9
4
B-n34-k5
33
5
788
788
798.8
16.4
5
P-n19-k2
18
2
212
212
213.6
4.8
6
P-n20-k2
19
2
216
216
217.4
5.2
7
P-n40-k5
39
5
458
458
469.6
13.1
8
E-n23-k3
22
3
569
574
575.4
8.0
9
E-n22-k4
21
4
375
375
387.6
9.5
10
F-n45-k4
44
4
724
730
763.8
20.2Table 11
The impacts of the parameters of IQEA+ under two large cases.
Pn78k10
Pn101k4
Best value
AVE
AVE of CPU time
Best value
AVE
AVE of CPU time
Pop size
20
1507
1651
78.5
852
1303
85
30
1256
1295
72.2
788
1322
136
50
1358
1704
192
984
1306
219
Rotation angle
0.005 pi
1574
1628
233
1252
1337
133
0.01 pi
1562
1618
244
1145
1221
138
0.025 pi
1607
1582
250
1076
1123
136
0.05 pi
1469
1610
234
1021
1068
140
Std/5
1602
1557
181
904
937
137
Std/10
1487
1615
136
897
955
144
Δ1
1321
1465
132
885
992
136
Mutation probability
0.01
2016
2048
130
1138
1230
143
0.025
1970
2066
133
1062
1217
144
0.05
1885
1918
144
1007
1193
137
0.1
1775
1848
137
913
1097
148
0.2
1539
1679
146
1282
1326
140
0.5
1493
1523
138
1223
1302
141
Δ2
1385
1505
134
895
1202
141
Δ1 = Std/(10 * selfsgn), Δ2 = Self-adapted P.Results of Table10 show that the parameters setting of IQEA+ mentioned in Section 4.1 is relatively reasonable.(2)  Parameter Settings for IQEA+ for Large Scale CVRPs. In fact, related parameters of the most intelligent algorithms are confirmed according to the specific problems. Finding a universal parameter setting for intelligence algorithms is rather difficult. As discussed in Section 4.1, the parameters of IQEA+ are also confirmed by parameters setting test and/or experiences of similar problems.Related parameters of IQEA+ for solving larger-scale CVRP should be also discussed. To simplify our analysis, we only retest two large scale cases “Bn78k10” and “Pn101k”. Also, these two cases with different number of vehicles can also show some performance information on local searching procedures. Those experiments are tested with ten runs and fixed 500 iterations each without any break intervention. The computation results show us the following.(a)
Three important parameters of IQEA+, population size, rotation angle, and mutation probability, can cause great impacts of IQEA+’s performance.(b)
The main influences of population size lies in the searching time but not in the quality of the searching results. Different rotation angles and mutation probabilities have great influences of IQEA+ on finding the optimal values, but the searching times in different settings are quite approximative to each other.(c)
A large mutation probability seems to have better performance than that with a smaller mutation probability.(d)
The longer the length of a single patrolling route, the worse the searching ability of the local search procedure. In “Bn78k10” and “Pn101k” analyses, the average lengths of the patrolling routes of the former case and the latter case are 7.8 and 25.3, respectively. So, the former case has better optimal results.Finding a universal parameter setting of IQEA+ is rather challengeable and costly under several types of CVRP benchmarks. Therefore, based on results of IQEA+ for solving the medium scale CVRPs and two large scale CVRP cases and the experiences of [19, 22], the similar parameters’ setting of IQEA+ in Section 4.1 is also be applied to solve large scale CVRPs.


### 4.2.2. Experiments for Large Scale CVRP Benchmarks

Large scale VRP benchmarks are also obtained from the website (http://branchandcut.org/VRP/data/), which were also introduced by Chen et al. [43] and Ai and Kachitvichyanukul [7]. There are sixteen CVRP cases, and the total number of clients is varying from 29 to 134, and the vehicle number is varying from 3 to 10. The settings of the IQEA+ for large scale CVRPs are as follows: the iteration is 500 times and the population size is 30. The setting of the rotation angle and mutation probability are the same as in Section 4.1. The searching performance of IQEA+ is compared with that from Chen’s et al. [43] DPSO algorithm and Ai and Kachitvichyanukul’s [7] SR-2 algorithm; all the comparison results are shown in Table 12.Table 12
Results of large-scale CVRPs with different algorithms.
Problem
Clients
Vehicles
Optimal value
Best objective function
AVE. IQEA+
CPU time (s)
Chen
SR-1
SR-2
IQEA+
Chen
SR-1
SR-2
IQEA+
An33k5
32
5
661
661
661
661
661
669.2
32
11
13
15.7
An46k7
45
7
914
914
914
914
914
970.6
129
19
23
72.0
An60k9
59
9
1354
1354
1366
1355
1363
1410.4
309
28
40
148.1
Bn35k5
34
5
955
955
955
955
955
967.6
38
12
14
21.2
Bn45k5
44
5
751
751
751
751
751
789.4
134
17
20
50.3
Bn68k9
67
9
1272
1272
1278
1274
1281
1349.8
344
33
50
65.4
Bn78k10
77
10
1221
1239
1239
1223
1256
1295.4
429
41
64
72.2
En30k3
29
3
534
534
541
534
534
551.2
28
11
16
13.4
En51k5
50
5
521
528
521
521
524
566
301
32
22
20.5
En76k7
75
7
682
688
691
682
703
739.4
527
38
60
150.0
Fn72k4
71
4
237
244
237
237
237
293.4
398
58
53
64.5
Fn135k7
134
7
1162
1215
1184
1162
1202
1414.2
1526
178
258
159.2
Mn101k10
100
10
820
824
821
820
956
1001.6
874
60
114
210.4
Mn121k7
120
7
1034
1038
1041
1036
1082
1207.4
1734
88
89
128.5
Pn76k4
75
4
593
602
599
594
594
612
496
51
48
64.9
Pn101k4
100
4
681
694
686
683
681
689.4
978
99
86
135.0


## 4.2.1. Parameters Setting Analysis

(1)  Experiments of Medium Scale CVRPs. In order to make IQEA+ suitable for solving large scale CVRP benchmarks, we firstly test it by solving some medium scale cases. The same parameters settings as mentioned in Section 4.1 are still adopted. Each case is run 20 times to obtain the average value. The results are shown in Table 10.Table 10
Results of some medium-scale CVRP cases by IQEA+.
No.
Problem instance
Object function
CPU time (s)
Type
Clients
Vehicles
Optimal
Best
AVE
1
A-n32-k5
31
5
784
784
791.6
15.7
2
A-n34-k5
33
5
778
780
794.4
17.4
3
B-n31-k5
30
5
672
672
679.8
11.9
4
B-n34-k5
33
5
788
788
798.8
16.4
5
P-n19-k2
18
2
212
212
213.6
4.8
6
P-n20-k2
19
2
216
216
217.4
5.2
7
P-n40-k5
39
5
458
458
469.6
13.1
8
E-n23-k3
22
3
569
574
575.4
8.0
9
E-n22-k4
21
4
375
375
387.6
9.5
10
F-n45-k4
44
4
724
730
763.8
20.2Table 11
The impacts of the parameters of IQEA+ under two large cases.
Pn78k10
Pn101k4
Best value
AVE
AVE of CPU time
Best value
AVE
AVE of CPU time
Pop size
20
1507
1651
78.5
852
1303
85
30
1256
1295
72.2
788
1322
136
50
1358
1704
192
984
1306
219
Rotation angle
0.005 pi
1574
1628
233
1252
1337
133
0.01 pi
1562
1618
244
1145
1221
138
0.025 pi
1607
1582
250
1076
1123
136
0.05 pi
1469
1610
234
1021
1068
140
Std/5
1602
1557
181
904
937
137
Std/10
1487
1615
136
897
955
144
Δ1
1321
1465
132
885
992
136
Mutation probability
0.01
2016
2048
130
1138
1230
143
0.025
1970
2066
133
1062
1217
144
0.05
1885
1918
144
1007
1193
137
0.1
1775
1848
137
913
1097
148
0.2
1539
1679
146
1282
1326
140
0.5
1493
1523
138
1223
1302
141
Δ2
1385
1505
134
895
1202
141
Δ1 = Std/(10 * selfsgn), Δ2 = Self-adapted P.Results of Table10 show that the parameters setting of IQEA+ mentioned in Section 4.1 is relatively reasonable.(2)  Parameter Settings for IQEA+ for Large Scale CVRPs. In fact, related parameters of the most intelligent algorithms are confirmed according to the specific problems. Finding a universal parameter setting for intelligence algorithms is rather difficult. As discussed in Section 4.1, the parameters of IQEA+ are also confirmed by parameters setting test and/or experiences of similar problems.Related parameters of IQEA+ for solving larger-scale CVRP should be also discussed. To simplify our analysis, we only retest two large scale cases “Bn78k10” and “Pn101k”. Also, these two cases with different number of vehicles can also show some performance information on local searching procedures. Those experiments are tested with ten runs and fixed 500 iterations each without any break intervention. The computation results show us the following.(a)
Three important parameters of IQEA+, population size, rotation angle, and mutation probability, can cause great impacts of IQEA+’s performance.(b)
The main influences of population size lies in the searching time but not in the quality of the searching results. Different rotation angles and mutation probabilities have great influences of IQEA+ on finding the optimal values, but the searching times in different settings are quite approximative to each other.(c)
A large mutation probability seems to have better performance than that with a smaller mutation probability.(d)
The longer the length of a single patrolling route, the worse the searching ability of the local search procedure. In “Bn78k10” and “Pn101k” analyses, the average lengths of the patrolling routes of the former case and the latter case are 7.8 and 25.3, respectively. So, the former case has better optimal results.Finding a universal parameter setting of IQEA+ is rather challengeable and costly under several types of CVRP benchmarks. Therefore, based on results of IQEA+ for solving the medium scale CVRPs and two large scale CVRP cases and the experiences of [19, 22], the similar parameters’ setting of IQEA+ in Section 4.1 is also be applied to solve large scale CVRPs.


## 4.2.2. Experiments for Large Scale CVRP Benchmarks

Large scale VRP benchmarks are also obtained from the website (http://branchandcut.org/VRP/data/), which were also introduced by Chen et al. [43] and Ai and Kachitvichyanukul [7]. There are sixteen CVRP cases, and the total number of clients is varying from 29 to 134, and the vehicle number is varying from 3 to 10. The settings of the IQEA+ for large scale CVRPs are as follows: the iteration is 500 times and the population size is 30. The setting of the rotation angle and mutation probability are the same as in Section 4.1. The searching performance of IQEA+ is compared with that from Chen’s et al. [43] DPSO algorithm and Ai and Kachitvichyanukul’s [7] SR-2 algorithm; all the comparison results are shown in Table 12.Table 12
Results of large-scale CVRPs with different algorithms.
Problem
Clients
Vehicles
Optimal value
Best objective function
AVE. IQEA+
CPU time (s)
Chen
SR-1
SR-2
IQEA+
Chen
SR-1
SR-2
IQEA+
An33k5
32
5
661
661
661
661
661
669.2
32
11
13
15.7
An46k7
45
7
914
914
914
914
914
970.6
129
19
23
72.0
An60k9
59
9
1354
1354
1366
1355
1363
1410.4
309
28
40
148.1
Bn35k5
34
5
955
955
955
955
955
967.6
38
12
14
21.2
Bn45k5
44
5
751
751
751
751
751
789.4
134
17
20
50.3
Bn68k9
67
9
1272
1272
1278
1274
1281
1349.8
344
33
50
65.4
Bn78k10
77
10
1221
1239
1239
1223
1256
1295.4
429
41
64
72.2
En30k3
29
3
534
534
541
534
534
551.2
28
11
16
13.4
En51k5
50
5
521
528
521
521
524
566
301
32
22
20.5
En76k7
75
7
682
688
691
682
703
739.4
527
38
60
150.0
Fn72k4
71
4
237
244
237
237
237
293.4
398
58
53
64.5
Fn135k7
134
7
1162
1215
1184
1162
1202
1414.2
1526
178
258
159.2
Mn101k10
100
10
820
824
821
820
956
1001.6
874
60
114
210.4
Mn121k7
120
7
1034
1038
1041
1036
1082
1207.4
1734
88
89
128.5
Pn76k4
75
4
593
602
599
594
594
612
496
51
48
64.9
Pn101k4
100
4
681
694
686
683
681
689.4
978
99
86
135.0


## 4.3. Discussions

In order to test the proposed IQEA+, firstly, it has been used to solve a simple CVRP to find different performances of the IQEA and IQEA+; then we test different body parts of IQEA to best play out its advantages. Three local search procedures combination schemes were also examined, and all of them were used to test the CVRP dataset.Experiments results show that the proposed IQEA+ has potential in solving the CVRP to some extent with relatively high precision and speed; see Table10 and several cases in Table 12. However, experiments also show some useful findings.(1)
IQEA+ outperforms DPSO, SR-1, and SR-2 in searching precision but not in the searching efficiency; see casePn101k4 in Table 12;(2)
Sometimes, IQEA+ does not show many advantages on searching performance; see caseFn135k7 in Table 12.(3)
IQEA can search quickly, while the local search procedures take more than a half of the total computing time, especially when problem scale becomes larger.Three main reasons can explain our findings. (a) There is a weakness of the IQEA+ in the length of quantum chromosome which originated from the encoding mechanism. When the client vertexes and vehicle number increase, especially as the length of quantum chromosome reached to thousand level, algorithm searching efficiency becomes very low. (b) The tightness (demand/capacity) of demand and capacity is another important factor that impacts the searching performance; see caseFn72k4 in Table 12. Different tightness may result in very different results: the more intensive are the tightness and the worse is the results, and the more intensive are the tightness, the greater chances of unfeasible routes are generated out. (c) Local search procedures are used to deal with the point-to-point operation and evaluation, which make it cost too much time.


## 5. Conclusions and Future Research

Although the CVRP has been presented for tens of years, still much attention has been put on it. Due to its NP nature, a lot of heuristic and metaheuristic algorithms have been proposed to provide better solutions for this problem. This paper proposes a new IQEA with the mixed local search procedure (named IQEA+) to solve the CVRPs: IQEA is the improvements of original QEA in the encoding mechanism for solving the CVRP. Moreover, in order to make full use of the information embedded in quantum chromosome, a double chain quantum chromosome is adopted. Local search procedures are to enhance the IQEA local search ability. The computational results show that the proposed IQEA+ is more effective in solving the CVRP than the DPSO and also has some advantages over other widely used SR-1 and SR-2 when the length of chromosome is within a certain range.There are some obvious advantages of the proposed IQEA+:(1)
less parameters to be set make the proposed IQEA+ very simple and easy to be implemented;(2)
unlike the classic genetic algorithm, the updating and mutation operations were replaced by easier quantum rotation gate and the quantum Not gate procedures;(3)
the comprehensive improvement of QEA with the double chain, new rotation angle, and rotation direction and self-adaptive Not gate can simplify its searching procedure and improve its efficiency. Moreover, three local search procedures help the IQEA+ jump out of local minimum points and continue its searching toward the right path.Future work can be focused on finding a new encoding mechanism to further increase the efficiency and keep the diversity of the quantum population. There still are improvement spaces in the IQEA+ framework, such as the significance of the rotation angle is designed as the problem-oriented and the mutation procedure is set as the fitness function adapted. Moreover, the proposed IQEA+ can be used to solve some VRP variants.



---

*Source: 159495-2013-12-09.xml*