# Multiobjective Quantum Evolutionary Algorithm for the Vehicle Routing Problem with Customer Satisfaction

Authors: Jingling Zhang; Wanliang Wang; Yanwei Zhao; Carlo Cattani
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The multiobjective vehicle routing problem considering customer satisfaction (MVRPCS) involves the distribution of orders from several depots to a set of customers over a time window. This paper presents a self-adaptive grid multi-objective quantum evolutionary algorithm (MOQEA) for the MVRPCS, which takes into account customer satisfaction as well as travel costs. The degree of customer satisfaction is represented by proposing an improved fuzzy due-time window, and the optimization problem is modeled as a mixed integer linear program. In the MOQEA, nondominated solution set is constructed by the Challenge Cup rules. Moreover, an adaptive grid is designed to achieve the diversity of solution sets; that is, the number of grids in each generation is not fixed but is automatically adjusted based on the distribution of the current generation of nondominated solution set. In the study, the MOQEA is evaluated by applying it to classical benchmark problems. Results of numerical simulation and comparison show that the established model is valid and the MOQEA is effective for MVRPCS.


---

## Body

## 1. Introduction

The vehicle routing problem (VRP) is one of the most important and widely studied combinatorial optimization problems, with many real-world applications in logistic distribution and transportation [1]. Since the VRP was firstly proposed by Dantzig and Ramser in 1959 [2], it has been focused in the field of operational research and combinatorial optimization [3–5].The aim of VRP is to find optimal routes for a fleet of vehicles serving a set of customers with known demands. Each customer is serviced exactly once and must be assigned a satisfactory vehicle without exceeding vehicle capacities. A solution for this problem is to find out a set of minimum cost routes that are used to represent vehicles distribution and clients’ permutation. However, current studies on VRP [6, 7] are mainly focused on the single objective problem and the objective is to optimize the number of vehicles dispatched and the travel distance, that is, reducing the service costs of the provider.Actually, to achieve competitive advantage, a service provider needs to consider not only service costs but also service quality that can determine customers’ satisfaction. Most of the research on multiobjective VRP (MOVRP) does not take into account this objective, only focusing on the traditional objectives of minimum costs and the length of the longest route [8, 9]. Hong and Park [10] constructed a linear goal programming (GP) model for the biobjective vehicle routing with time window constraints (BVRPTW) and proposed a heuristic algorithm to relieve the computational burden inherent to the application of the GP model. Zitzler and Thiele [11] proposed a multiobjective evolutionary algorithm based on the Pareto approach for VRP. Lim and wang [12] proposed a method to deal with MOVRP by assigning different weights of the objectives. Tan et al. [13] proposed a hybrid multiobjective evolutionary algorithm (HMOEA) that incorporates various heuristics for local exploitation in the evolutionary search and the concept of Pareto’s optimality for solving multiobjective optimization in vehicle routing problem with time window constraints (VRPTW). Garcia-Najera and Bullinaria [14] studied an improved multiobjective evolutionary algorithm for VRP with time windows.The VRPTW is developed from VRP and has been widely studied in the last decade [15–19]. In the VRPTW, each customer has a time window with values about the deadline and the earliest time constraints for the service he/she requires. Thus this problem involves a routing combination and scheduling component. Routes must be designed to minimize the total cost, but, at the same time, scheduling must be performed to ensure time feasibility.In practice, this time window actually does not well describe customers’ satisfaction. A major reason is that customers are asked to provide a fixed time window for service, but in reality they really hope to be served at a desired time. Cheng and Gen [20, 21] called such a desired time the due-time and proposed to use the concept of fuzzy due-time to replace this time window because, as they claimed, it can describe customers’ satisfaction better. Thus customers’ satisfaction can be also described as a convex fuzzy number [22–24].Cheng and Gen [20, 21] introduced the fuzzy due-time, used triangle fuzzy numbers to describe customers’ satisfaction, and solved the VRP by genetic algorithms (GAs). Zhang et al. [25] proposed a multiobjective fuzzy VRP and used trapezoid fuzzy numbers to describe customers’ satisfaction.  Jia [26] used multiobjective hybrid GA for this problem.  Wu [27] studied the open VRP based on customers’ satisfaction. Lin [28] proposed a GA-based multiobjective decision-making method for optimal vehicle transportation, which is focused on a fuzzy vehicle routing and scheduling problem (FVRSP) based on five attributes, namely, space utility, service satisfaction, waiting time, delay time, and transportation distance. Wang and Li [29] proposed a hybrid algorithm based on GA and incorporated some methods based on greedy algorithms to solve the MOP model for which depot desires and clients expectations are considered simultaneously.The above studies use the weighted sums of objectives to solve the multiobjective problem; the higher an objective’s importance, the larger its corresponding weight coefficient. In general, no single solution can attain the optimum of all objectives at the same time. Therefore, it is desirable to obtain a set of Pareto optimal solution, that is, the Pareto set. The points in the objective space that correspond to the results in the set are usually called the Pareto front.In this paper, a self-adaptive grid multiobjective quantum evolutionary algorithm (MOQEA) is proposed to solve the MVRPCS problem. In particular, the quantum evolutionary algorithm (QEA) is used in the MOQEA due to its high efficiency, convergence speed, strong full-searching optimization ability [30]. With the MOQEA, an optimal or a nearly optimal set of vehicle routes solution with the minimal total travel cost and maximal customers’ satisfaction can be obtained by decoding the chromosome and simultaneously obtains several solution sets. This method can support the dispatcher to more efficiently determine how to distribute the shipment to serve customers by available vehicles.The remainder of the paper is organized as follows. Section2 presents the mathematical model for the MVRPCS with consideration of customer satisfaction. Section 3 describes the proposed MOQEA in detail. In Section 4, the application of the proposed algorithm to a classic problem is introduced and the simulation results are discussed and compared with an early algorithm. Finally, conclusions are given in Section 5.


## 2.  Model for the MVRPCS

### 2.1. Representation of Customer Satisfaction

In traditional VRP, customers’ time constraints are represented by time windows as shown in Figure1. For this method, if the customer is serviced at a time within the window, then the satisfaction degree is 100%, otherwise the satisfaction degree is 0. It is really unrealistic to measure the degree of satisfaction accordingly based on time window because customer satisfaction is not necessarily the same if they are serviced at different times within the window. In fact, service time can be divided into two categories; namely, the service time can be tolerated and the desirable service time.Figure 1
Traditional time windows.Fuzzy due-time windows have been introduced to describe different degrees of satisfaction. Generally, the tolerable service time for customeri can be described as [Ei,Li], where Ei is the earliest time and Li is the latest time. For example, in [20, 25] customers satisfaction is represented by the triangular fuzzy number, as shown in Figure 2. When the desirable service time is taken into account, customer satisfaction can be described using trapezoidal fuzzy number. In [26], this method is used and the desirable service time is described by the due-time [ai,bi], as shown in Figure 3. In this case, customer satisfaction is zero no matter vehicles arrive early or late if the expected service time slot is not achieved.Figure 2
Triangular fuzzy number.Figure 3
Trapezoidal fuzzy number.In this paper, an improved fuzzy due-time window is proposed, as shown in Figure4. In this method, if a customer is served at the desired time, the degree of satisfaction is 1; otherwise, the degree of satisfaction decreases as the difference between the actual service time and the desired time increases. The degree of satisfaction for the cases in which vehicles arrive early than the earliest expected time is not zero, but equals to that of the case when the earliest expected time is met. The degree of customer satisfaction is represented by the membership function of the improved fuzzy due-time window, that is, an improved trapezoidal fuzzy number. For customer i, mark his/her satisfaction as μi(ti) for a given service time ti. Then the degree of satisfaction can be calculated using the membership function as follows:
(2.1)μi(ti)={expk(Ei-ai),ti<Eiexpk(ti-ai),Ei<ti<ai1,ai<ti<biexpk(bi-ti),bi<ti<Li0,ti>Li.Figure 4
Improved trapezoidal fuzzy number.


### 2.2. Mathematics Model

The MVRPCS can be described as follows: there areM depots each of which has Km(m=1,2,…,M) vehicles with a capacity of bk. These vehicles will be dispatched to L customers to meet their demands.Mark the demand of customeri as di(i=1,2,…,L), and assume di<bk. Each customer can be served by any vehicles from a depot, but the service will only take place once. In addition, each vehicle can complete the shipping task without having to return to the original depot. Thus an appropriate vehicle scheduling program is required to meet the needs of all customers. The meanings of variables used in this research are described as follows.Customer number is1,2,…,L. Depot number is L+1,L+2,…,L+M.Fixed cost of sending a vehicle isFk(m∈{L+1,L+2,…,L+M};k∈{1,2,…,Km}).Distribution cost from customeri to customer j is cij, and i,j∈(1,2,…,L,L+1,…,L+M).Time window of customeri is [Ei,Li]; the arrival time of customer i is ti; the travel time from i to j is tij,siis the service time of customer i and μi(ti) is the degree of satisfaction for customer i.If vehiclek travels directly from customer i to customer j and arrives too early at j, it will wait; wj(tj) is the waiting time of customer j for a vehicle. Thus, the MVRPCS model can be established, as discussed in detail in the following paragraphs. Two decision variables are defined as follows:
(2.2)xijkm={1vehiclekfromitojofdepotm,0other,ykm={1vehiclekdispatchedfromdepotm,0other.
This problem has two objectives that are determined by both the service quality and the service costs. Considering the two objectives can help the dispatcher make better decision without compromising any of the two objectives. The dispatcher can still have preferences on different objectives by selecting different parameters.The service quality objective is to maximize average customers satisfaction:(2.3)max1L∑i=1Lμi(ti).This objective is equivalent to minimizing the average customer dissatisfaction:(2.4)min(1-1L∑i=1Lμi(ti)).The other objective of the service costs is to minimize travel cost, fixed cost and waiting cost. For this objective, the fixed cost of sending a vehicle is considered because vehicles in operation have depreciation and fuel consumption. Also the fixed cost is related to the number of vehicle, that is, the more vehicles, the higher fixed cost. To the best of our knowledge, no previous work has been done to take into account this fixed cost when solving multiobjective VRP with fuzzy due-time. Based on the above discussion, the mathematical model for the MVRPCS can be established as follows:(2.5)minZ=∑i=1L+M∑j=1L+M∑k=1Km∑m=1Mxijkmcij+∑k=1Km∑m=1MykmFk+∑i=1Lwi(ti)s.t.,(2.6)∑i=1L∑j=1L+Mxijkmdi≤bk,m∈{L+1,L+2,…,L+M};k∈{1,2,…Km},(2.7)∑j=1L∑k=1Kmxijkm≤Km,i=m∈{L+1,L+2,…,L+M},(2.8)∑j=1L+M∑k=1Km∑m=1Mxijkm=1,i∈{1,2,…,L},(2.9)∑i=1L+M∑k=1Km∑m=1Mxijkm=1,j∈{1,2,…,L},(2.10)tj={ti+si+tij,ti≥Ei,Ei+si+tij,ti<Ei,(2.11)tj<Lj,(2.12)wj(tj)=Max{0,Ej-tj}.
In the model, formula (2.6) assures the carrying capacity of each vehicle; formula (2.7) assures the number of vehicles that are dispatched from each depot does not exceed the capacity of the depot; formulas (2.8) and (2.9) assure each client is served only by one vehicle; formulas (2.10), (2.11), and (2.12) assure customers can be served within time windows; and from (2.10) and (2.12), the waiting cost for the vehicle can be obtained.


## 2.1. Representation of Customer Satisfaction

In traditional VRP, customers’ time constraints are represented by time windows as shown in Figure1. For this method, if the customer is serviced at a time within the window, then the satisfaction degree is 100%, otherwise the satisfaction degree is 0. It is really unrealistic to measure the degree of satisfaction accordingly based on time window because customer satisfaction is not necessarily the same if they are serviced at different times within the window. In fact, service time can be divided into two categories; namely, the service time can be tolerated and the desirable service time.Figure 1
Traditional time windows.Fuzzy due-time windows have been introduced to describe different degrees of satisfaction. Generally, the tolerable service time for customeri can be described as [Ei,Li], where Ei is the earliest time and Li is the latest time. For example, in [20, 25] customers satisfaction is represented by the triangular fuzzy number, as shown in Figure 2. When the desirable service time is taken into account, customer satisfaction can be described using trapezoidal fuzzy number. In [26], this method is used and the desirable service time is described by the due-time [ai,bi], as shown in Figure 3. In this case, customer satisfaction is zero no matter vehicles arrive early or late if the expected service time slot is not achieved.Figure 2
Triangular fuzzy number.Figure 3
Trapezoidal fuzzy number.In this paper, an improved fuzzy due-time window is proposed, as shown in Figure4. In this method, if a customer is served at the desired time, the degree of satisfaction is 1; otherwise, the degree of satisfaction decreases as the difference between the actual service time and the desired time increases. The degree of satisfaction for the cases in which vehicles arrive early than the earliest expected time is not zero, but equals to that of the case when the earliest expected time is met. The degree of customer satisfaction is represented by the membership function of the improved fuzzy due-time window, that is, an improved trapezoidal fuzzy number. For customer i, mark his/her satisfaction as μi(ti) for a given service time ti. Then the degree of satisfaction can be calculated using the membership function as follows:
(2.1)μi(ti)={expk(Ei-ai),ti<Eiexpk(ti-ai),Ei<ti<ai1,ai<ti<biexpk(bi-ti),bi<ti<Li0,ti>Li.Figure 4
Improved trapezoidal fuzzy number.


## 2.2. Mathematics Model

The MVRPCS can be described as follows: there areM depots each of which has Km(m=1,2,…,M) vehicles with a capacity of bk. These vehicles will be dispatched to L customers to meet their demands.Mark the demand of customeri as di(i=1,2,…,L), and assume di<bk. Each customer can be served by any vehicles from a depot, but the service will only take place once. In addition, each vehicle can complete the shipping task without having to return to the original depot. Thus an appropriate vehicle scheduling program is required to meet the needs of all customers. The meanings of variables used in this research are described as follows.Customer number is1,2,…,L. Depot number is L+1,L+2,…,L+M.Fixed cost of sending a vehicle isFk(m∈{L+1,L+2,…,L+M};k∈{1,2,…,Km}).Distribution cost from customeri to customer j is cij, and i,j∈(1,2,…,L,L+1,…,L+M).Time window of customeri is [Ei,Li]; the arrival time of customer i is ti; the travel time from i to j is tij,siis the service time of customer i and μi(ti) is the degree of satisfaction for customer i.If vehiclek travels directly from customer i to customer j and arrives too early at j, it will wait; wj(tj) is the waiting time of customer j for a vehicle. Thus, the MVRPCS model can be established, as discussed in detail in the following paragraphs. Two decision variables are defined as follows:
(2.2)xijkm={1vehiclekfromitojofdepotm,0other,ykm={1vehiclekdispatchedfromdepotm,0other.
This problem has two objectives that are determined by both the service quality and the service costs. Considering the two objectives can help the dispatcher make better decision without compromising any of the two objectives. The dispatcher can still have preferences on different objectives by selecting different parameters.The service quality objective is to maximize average customers satisfaction:(2.3)max1L∑i=1Lμi(ti).This objective is equivalent to minimizing the average customer dissatisfaction:(2.4)min(1-1L∑i=1Lμi(ti)).The other objective of the service costs is to minimize travel cost, fixed cost and waiting cost. For this objective, the fixed cost of sending a vehicle is considered because vehicles in operation have depreciation and fuel consumption. Also the fixed cost is related to the number of vehicle, that is, the more vehicles, the higher fixed cost. To the best of our knowledge, no previous work has been done to take into account this fixed cost when solving multiobjective VRP with fuzzy due-time. Based on the above discussion, the mathematical model for the MVRPCS can be established as follows:(2.5)minZ=∑i=1L+M∑j=1L+M∑k=1Km∑m=1Mxijkmcij+∑k=1Km∑m=1MykmFk+∑i=1Lwi(ti)s.t.,(2.6)∑i=1L∑j=1L+Mxijkmdi≤bk,m∈{L+1,L+2,…,L+M};k∈{1,2,…Km},(2.7)∑j=1L∑k=1Kmxijkm≤Km,i=m∈{L+1,L+2,…,L+M},(2.8)∑j=1L+M∑k=1Km∑m=1Mxijkm=1,i∈{1,2,…,L},(2.9)∑i=1L+M∑k=1Km∑m=1Mxijkm=1,j∈{1,2,…,L},(2.10)tj={ti+si+tij,ti≥Ei,Ei+si+tij,ti<Ei,(2.11)tj<Lj,(2.12)wj(tj)=Max{0,Ej-tj}.
In the model, formula (2.6) assures the carrying capacity of each vehicle; formula (2.7) assures the number of vehicles that are dispatched from each depot does not exceed the capacity of the depot; formulas (2.8) and (2.9) assure each client is served only by one vehicle; formulas (2.10), (2.11), and (2.12) assure customers can be served within time windows; and from (2.10) and (2.12), the waiting cost for the vehicle can be obtained.


## 3. MOQEA for the MVRPCS

In this research, the multiobjective optimization method of Pareto optimal solution [31] is used. Its main advantage is to approximate the Pareto front in order to provide a set of equivalent solutions to the decision maker [32, 33]. The algorithm to solve multiobjective optimization problem of Pareto optimality involves two main questions:(i)
how to construct a Pareto optimal solution set, namely non-dominated solutions set, and make it close to the Pareto optimal front as much as possible?(ii)
how to attain the diversity and variety of solutions?To address these two issues, a self-adaptive grid multiobjective quantum evolutionary algorithm (MOQEA) to solve the MVRPCS problem is proposed. The method of constructing non-dominated solution set and attaining the diversity and variety of solutions is described in the following sections.


### 3.1. Constructing Nondominated Solution Set

In this paper, the Challenge Cup rule [34] is used to construct non-dominated solution set. Assuming P is an evolution group and Q is a constructed set, let Q=P initially. Assume Ndset is a non-dominated set which is empty initially. The basic idea of this rule is to take any individual x from Q, followed by comparison with all other individuals y in Q: if x≻y, then clear y from Q; if y≻x, then use y to replace x, and then y is the new champion and continues to be compared with other individuals. After a comparison, cluster(x)={y∣x≻y,x,y∈P} is formed, where x is the smallest element. Add x to Ndset and continue the comparison until Q is empty.


### 3.2. Method of Attaining the Diversity and Variety of Set Based on Self-Adaptive Grid

To attain the variety of the set, the individual space is divided into several small areas each of which is a called a grid, as shown in Figure5. Thus each individual is associated with a grid in the figure, and the number of individuals in each grid can be defined as extrusion coefficient.Figure 5
Individual space divided by grid.A grid is used in many different ways to maintain the diversity and variety of solutions. Knowles and Crone [35] proposed a pareto archived evolution strategy for pareto multiobjective optimization. Corne et al. [36] proposed a pareto envelope-based selection algorithm for multiobjective optimization.When a grid contains more than one individual, these individuals are treated as the same solution. As such, the size of grid is very important. When the grid is too large, multiple individuals will exist in the same grid, and the resultant solution distribution is not accurate. When the grid is too small, it is likely that there are no individuals in some grids, and so it takes longer computation time though the resultant solution distribution is accurate. Therefore, computation time and accuracy must be traded off when determining the grid size.There are two objective functions in this optimization problem. The range of the customers’ satisfaction is[0,1], and the range of travel costs and waiting costs is changed with the experimental data. The boundaries of the grid in the target space can be determined by the range of the above two objective functions.In this paper, the number of grids is not fixed in each generation but automatically adjusted based on the distribution of the current generation of non-dominated solution set. The grid boundary is a fixed value. In the process of each evolutionary generation, the number of grid is adjusted by theD-value between the maximum and minimum values of each dimension in the non-dominated solution set. The method of self-adaptive grid is designed as follows.The two objections can be described byk={1,2}. Mark the number of each dimension grid in generation 1 as N1,k, Ndset1 as the non-dominated solutions set in generation 1, and Ndsett as the non-dominated solutions set in generation t. If ft,k is the objection of each dimension, then the D-value of each dimension can be described as follows:generation 1:∀k={1,2}(3.1)Diffence1,k=max{f1,k∣f1∈Ndset1}-min{f1,k∣f1∈Ndset1},generationt:∀k={1,2}(3.2)Diffencet,k=max{ft,k∣ft∈Ndsett}-min{ft,k∣ft∈Ndsett}.The number of each dimension grid in generationt is
(3.3)Nt,k=[N1,k*Diffence1,kDiffencet,k],where[]meansbeingrounded.                     To keep the diversity and variety of the non-dominated solution set, choose the individual with the maximum extrusion coefficient and delete it form the non-dominated solution set.


### 3.3. Representation

Quantum evolutionary algorithm (QEA) [37] is based on the concept and principles of quantum computing such as a quantum bit and superposition of states. Instead of binary and numeric representation, QEA uses a Q-bit chromosome as a probabilistic representation and a Q-bits individual is modeled by a string of Q-bits.The smallest unit of information stored in two-state quantum computer is called a Q-bit, which may be in the “1” state, or in the “0” state, or in any superposition of the two. The state of a Q-bit can be represented as follows(3.4)|ψ〉=α|0〉+β|1〉,
where α and β are complex numbers. |α|2 and |β|2 donate the probabilities that the Q-bit will be found in the “0” state and in the “1” state, respectively. Normalization of the state to unity is used to meet |α|2+|β|2=1.So a Q-bit individual with a string of m Q-bits can be expressed as follows(3.5)[α1β1|α2β2|⋯αmβm],
where |αi|2+|βi|2=1,i=1,2,…,m.The main advantage of the representation is that any linear superposition of solutions can be represented. For example, a three-Q-bit system can contain the information of eight states. QEA with Q-bit representation has a better characteristic of population diversity than other representations, as it can represent linear superposition of state’s probabilities.In this paper, a method of converting integer representation to Q-bit representation is designed. For the MVRPCS withL customers, the representation of customers route is described as the permutation of 1~L. Note the permutation of 1~L is (g1,g2,…,gL), and represent each gene gj(j=1,2,…,L) as a string ofr-Q-bit, then L  groups n-Q-bit are obtained. So, a quantum individual is described as a rL×2 Q-bit matrix. Here r={m}, where {m} is the smallest integer not less than m and 2m≥L, thus we can get m≥log2L.


### 3.4. Decoding Method

The “Customers permutation Route First, Vehicles distribution Cluster Second” rule is adopted for decoding.(1)
Firstly, get the customers permutation route. The solution of MVRPCS is a permutation of all customers and Q-bit representation cannot be evaluated directly. So it should be converted to permutation for evaluation.The Q-bit string is firstly converted to binary stringγ. Specifically, a random number η between [0,1] is generated, if the ith bit αi of Q-bit string satisfied |αi|2>η, then let the corresponding bit γi of the binary string γ be 1, otherwise let it be 0. Then the binary representation is converted to integer representation, which is viewed as random key representation [38], and customer permutation is constructed based on the generated random key. If two random key values are different, the smaller random key denotes the customer with smaller number; otherwise, let the one that first appears denote the customer with smaller number.(2)
Secondly, distribute the vehicles and get the subroute. A vehicle is dispatched to service customers according to the customers’ permutation route, if the vehicle cannot serve the next customer when it cannot meet the time window or loading capacity constraints, a new vehicle will be dispatched. For example, the customers’ permutation route is [859341267], the 3 depots notes [101112], use this decoding method, the subroute is: Route 1: 11-8-5-9; Route 2: 10-3-4-1; Route 3: 12-2-6-7.


### 3.5. Strategy of Updating by Q-Gate

In the MOQEA, a Q-gate is an evolution operator which is the same as the QEA in [39]. A rotation gate is often used to update a Q-bit individual, as shown in (3.6)
(3.6)[αi′βi′]=U(θi)[αiβi]=[cos(θi)-sin(θi)sin(θi)cos(θi)][αiβi],
where [αi,βi]T is the ith Q-bit and θi is the rotation angle of each Q-bit. θi=s(αi,βi)Δθi.The lookup table ofθi is shown in Table 1. In this paper, a non-dominated solution is randomly selected as the current objective solution from the non-dominated solution set. In the multiobjective optimization, it is unable to find the optimal value with all objectives met. So we need to choose an objective solution y for each individual r only to find the non-dominated solution set.Table 1
Lookup table of rotation angle.
r
i
y
i
f
(
y
)
≻
f
(
r
)
Δ
θ
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
θ
+1
−1
0
±
1
0
1
False
0
0
0
0
0
0
1
True
θ
−1
+1
±
1
0
1
0
False
θ
−1
+1
±
1
0
1
0
True
θ
−1
+1
0
±
1
1
1
False
0
0
0
0
0
1
1
True
θ
−1
+1
0
±
1In the above table,Δθi is the magnitude of rotation angle. s(αi,βi) is the sign of θi that determines the direction. ri and yi are the ith Q-bit of the binary solution in individual r and the objective solution y respectively. f(∙) is the objective. θ is the rotation angle of size, which affects the convergence speed and search capability. In this paper, a method of dynamically adjusting the rotation angle is proposed, that is; θ will be changed with the extrusion coefficient as discussed in Section 3:
(3.7)θ=0.5snπ,
where n is the extrusion coefficient. s is a constant in the range [0,1]. From (3.7), we can see when the extrusion coefficient n is small, the rotation angle θ is big to accelerate the convergence speed; when the Extrusion Coefficient n is big, and the search step size will be reduced to enhance the solution diversity.


### 3.6. Procedure of MOQEA

The flow chart of the MOQEA for this problem is illustrated in Figure6.Figure 6
Flow chart of MOQEA.The detailed procedure of the MOQEA is as follows.Step 1.
Lett=0 and randomly generate an initial population Q(t)={q1t,q2t⋯qnt}={α1tβ1t|α2tβ2t|⋯αntβnt}, that is, randomly generate any value in [0,1] for αi and βi, where qjt denotes the jth individual in the tth generation. At the same time, construct initially empty external set O(t) with size of  m.Step 2.
ConvertQ(t) to binary population R(t), then convert it to integer population P(t).Step 3.
According to the decoding method to get the subroute, evaluate the objectives to get the MVRPCS solution setM(t)={(f11t,f21t)⋯(f1nt,f2nt)}, where (f1jt,f2jt) is the jth value of the two objectives in the tth generation. And let M(t) be the construction set.Step 4.
Use the formulas (2.4) and (2.5) to evaluate the domination of each (f1jt,f2jt) in construction set M(t) and use the method of Challenge Cup rules to construct non-dominated solution set NDset(t).Step 5.
Whent=0, reproduce NDset(0) to the external set O(t). When t>0, if a certain individual in NDset(t) dominates one solution in O(t) delete the solution and join the individual into the O(t), and if a solution in O(t) dominates a certain individual in NDset(t), the solutions in O(t) does not change; otherwise, join the individual of NDset(t) into the O(t).Step 6.
Adjust the size ofO(t) to the number of m and satisfy the distribution of the non-dominated solution set. The method is discussed in detail as follows. If the size of O(t) is less than m, randomly select individual of NDset(t) into the O(t) until the size of O(t) is m; otherwise, use the self-adaptive grid method to keep the diversity and variety of O(t), divide the individual space of O(t) into several small grids, choose the grid with the maximum extrusion coefficient, and randomly delete a solution set from it.Step 7.
If the stopping condition is satisfied, then output the Pareto set; otherwise, go on to the following steps.Step 8.
Randomly select some individuals from the Q-bitQ(t), which is instated of the individuals from Q(t) corresponding to O(t).Step 9.
Use (3.6) to perform rotation operation for Q(t) to generate Q(t+1).Step 10.
Lett=t+1 and go back to Step 2.


## 3.1. Constructing Nondominated Solution Set

In this paper, the Challenge Cup rule [34] is used to construct non-dominated solution set. Assuming P is an evolution group and Q is a constructed set, let Q=P initially. Assume Ndset is a non-dominated set which is empty initially. The basic idea of this rule is to take any individual x from Q, followed by comparison with all other individuals y in Q: if x≻y, then clear y from Q; if y≻x, then use y to replace x, and then y is the new champion and continues to be compared with other individuals. After a comparison, cluster(x)={y∣x≻y,x,y∈P} is formed, where x is the smallest element. Add x to Ndset and continue the comparison until Q is empty.


## 3.2. Method of Attaining the Diversity and Variety of Set Based on Self-Adaptive Grid

To attain the variety of the set, the individual space is divided into several small areas each of which is a called a grid, as shown in Figure5. Thus each individual is associated with a grid in the figure, and the number of individuals in each grid can be defined as extrusion coefficient.Figure 5
Individual space divided by grid.A grid is used in many different ways to maintain the diversity and variety of solutions. Knowles and Crone [35] proposed a pareto archived evolution strategy for pareto multiobjective optimization. Corne et al. [36] proposed a pareto envelope-based selection algorithm for multiobjective optimization.When a grid contains more than one individual, these individuals are treated as the same solution. As such, the size of grid is very important. When the grid is too large, multiple individuals will exist in the same grid, and the resultant solution distribution is not accurate. When the grid is too small, it is likely that there are no individuals in some grids, and so it takes longer computation time though the resultant solution distribution is accurate. Therefore, computation time and accuracy must be traded off when determining the grid size.There are two objective functions in this optimization problem. The range of the customers’ satisfaction is[0,1], and the range of travel costs and waiting costs is changed with the experimental data. The boundaries of the grid in the target space can be determined by the range of the above two objective functions.In this paper, the number of grids is not fixed in each generation but automatically adjusted based on the distribution of the current generation of non-dominated solution set. The grid boundary is a fixed value. In the process of each evolutionary generation, the number of grid is adjusted by theD-value between the maximum and minimum values of each dimension in the non-dominated solution set. The method of self-adaptive grid is designed as follows.The two objections can be described byk={1,2}. Mark the number of each dimension grid in generation 1 as N1,k, Ndset1 as the non-dominated solutions set in generation 1, and Ndsett as the non-dominated solutions set in generation t. If ft,k is the objection of each dimension, then the D-value of each dimension can be described as follows:generation 1:∀k={1,2}(3.1)Diffence1,k=max{f1,k∣f1∈Ndset1}-min{f1,k∣f1∈Ndset1},generationt:∀k={1,2}(3.2)Diffencet,k=max{ft,k∣ft∈Ndsett}-min{ft,k∣ft∈Ndsett}.The number of each dimension grid in generationt is
(3.3)Nt,k=[N1,k*Diffence1,kDiffencet,k],where[]meansbeingrounded.                     To keep the diversity and variety of the non-dominated solution set, choose the individual with the maximum extrusion coefficient and delete it form the non-dominated solution set.


## 3.3. Representation

Quantum evolutionary algorithm (QEA) [37] is based on the concept and principles of quantum computing such as a quantum bit and superposition of states. Instead of binary and numeric representation, QEA uses a Q-bit chromosome as a probabilistic representation and a Q-bits individual is modeled by a string of Q-bits.The smallest unit of information stored in two-state quantum computer is called a Q-bit, which may be in the “1” state, or in the “0” state, or in any superposition of the two. The state of a Q-bit can be represented as follows(3.4)|ψ〉=α|0〉+β|1〉,
where α and β are complex numbers. |α|2 and |β|2 donate the probabilities that the Q-bit will be found in the “0” state and in the “1” state, respectively. Normalization of the state to unity is used to meet |α|2+|β|2=1.So a Q-bit individual with a string of m Q-bits can be expressed as follows(3.5)[α1β1|α2β2|⋯αmβm],
where |αi|2+|βi|2=1,i=1,2,…,m.The main advantage of the representation is that any linear superposition of solutions can be represented. For example, a three-Q-bit system can contain the information of eight states. QEA with Q-bit representation has a better characteristic of population diversity than other representations, as it can represent linear superposition of state’s probabilities.In this paper, a method of converting integer representation to Q-bit representation is designed. For the MVRPCS withL customers, the representation of customers route is described as the permutation of 1~L. Note the permutation of 1~L is (g1,g2,…,gL), and represent each gene gj(j=1,2,…,L) as a string ofr-Q-bit, then L  groups n-Q-bit are obtained. So, a quantum individual is described as a rL×2 Q-bit matrix. Here r={m}, where {m} is the smallest integer not less than m and 2m≥L, thus we can get m≥log2L.


## 3.4. Decoding Method

The “Customers permutation Route First, Vehicles distribution Cluster Second” rule is adopted for decoding.(1)
Firstly, get the customers permutation route. The solution of MVRPCS is a permutation of all customers and Q-bit representation cannot be evaluated directly. So it should be converted to permutation for evaluation.The Q-bit string is firstly converted to binary stringγ. Specifically, a random number η between [0,1] is generated, if the ith bit αi of Q-bit string satisfied |αi|2>η, then let the corresponding bit γi of the binary string γ be 1, otherwise let it be 0. Then the binary representation is converted to integer representation, which is viewed as random key representation [38], and customer permutation is constructed based on the generated random key. If two random key values are different, the smaller random key denotes the customer with smaller number; otherwise, let the one that first appears denote the customer with smaller number.(2)
Secondly, distribute the vehicles and get the subroute. A vehicle is dispatched to service customers according to the customers’ permutation route, if the vehicle cannot serve the next customer when it cannot meet the time window or loading capacity constraints, a new vehicle will be dispatched. For example, the customers’ permutation route is [859341267], the 3 depots notes [101112], use this decoding method, the subroute is: Route 1: 11-8-5-9; Route 2: 10-3-4-1; Route 3: 12-2-6-7.


## 3.5. Strategy of Updating by Q-Gate

In the MOQEA, a Q-gate is an evolution operator which is the same as the QEA in [39]. A rotation gate is often used to update a Q-bit individual, as shown in (3.6)
(3.6)[αi′βi′]=U(θi)[αiβi]=[cos(θi)-sin(θi)sin(θi)cos(θi)][αiβi],
where [αi,βi]T is the ith Q-bit and θi is the rotation angle of each Q-bit. θi=s(αi,βi)Δθi.The lookup table ofθi is shown in Table 1. In this paper, a non-dominated solution is randomly selected as the current objective solution from the non-dominated solution set. In the multiobjective optimization, it is unable to find the optimal value with all objectives met. So we need to choose an objective solution y for each individual r only to find the non-dominated solution set.Table 1
Lookup table of rotation angle.
r
i
y
i
f
(
y
)
≻
f
(
r
)
Δ
θ
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
θ
+1
−1
0
±
1
0
1
False
0
0
0
0
0
0
1
True
θ
−1
+1
±
1
0
1
0
False
θ
−1
+1
±
1
0
1
0
True
θ
−1
+1
0
±
1
1
1
False
0
0
0
0
0
1
1
True
θ
−1
+1
0
±
1In the above table,Δθi is the magnitude of rotation angle. s(αi,βi) is the sign of θi that determines the direction. ri and yi are the ith Q-bit of the binary solution in individual r and the objective solution y respectively. f(∙) is the objective. θ is the rotation angle of size, which affects the convergence speed and search capability. In this paper, a method of dynamically adjusting the rotation angle is proposed, that is; θ will be changed with the extrusion coefficient as discussed in Section 3:
(3.7)θ=0.5snπ,
where n is the extrusion coefficient. s is a constant in the range [0,1]. From (3.7), we can see when the extrusion coefficient n is small, the rotation angle θ is big to accelerate the convergence speed; when the Extrusion Coefficient n is big, and the search step size will be reduced to enhance the solution diversity.


## 3.6. Procedure of MOQEA

The flow chart of the MOQEA for this problem is illustrated in Figure6.Figure 6
Flow chart of MOQEA.The detailed procedure of the MOQEA is as follows.Step 1.
Lett=0 and randomly generate an initial population Q(t)={q1t,q2t⋯qnt}={α1tβ1t|α2tβ2t|⋯αntβnt}, that is, randomly generate any value in [0,1] for αi and βi, where qjt denotes the jth individual in the tth generation. At the same time, construct initially empty external set O(t) with size of  m.Step 2.
ConvertQ(t) to binary population R(t), then convert it to integer population P(t).Step 3.
According to the decoding method to get the subroute, evaluate the objectives to get the MVRPCS solution setM(t)={(f11t,f21t)⋯(f1nt,f2nt)}, where (f1jt,f2jt) is the jth value of the two objectives in the tth generation. And let M(t) be the construction set.Step 4.
Use the formulas (2.4) and (2.5) to evaluate the domination of each (f1jt,f2jt) in construction set M(t) and use the method of Challenge Cup rules to construct non-dominated solution set NDset(t).Step 5.
Whent=0, reproduce NDset(0) to the external set O(t). When t>0, if a certain individual in NDset(t) dominates one solution in O(t) delete the solution and join the individual into the O(t), and if a solution in O(t) dominates a certain individual in NDset(t), the solutions in O(t) does not change; otherwise, join the individual of NDset(t) into the O(t).Step 6.
Adjust the size ofO(t) to the number of m and satisfy the distribution of the non-dominated solution set. The method is discussed in detail as follows. If the size of O(t) is less than m, randomly select individual of NDset(t) into the O(t) until the size of O(t) is m; otherwise, use the self-adaptive grid method to keep the diversity and variety of O(t), divide the individual space of O(t) into several small grids, choose the grid with the maximum extrusion coefficient, and randomly delete a solution set from it.Step 7.
If the stopping condition is satisfied, then output the Pareto set; otherwise, go on to the following steps.Step 8.
Randomly select some individuals from the Q-bitQ(t), which is instated of the individuals from Q(t) corresponding to O(t).Step 9.
Use (3.6) to perform rotation operation for Q(t) to generate Q(t+1).Step 10.
Lett=t+1 and go back to Step 2.


## 4. Experiment Results and Comparisons

### 4.1. Experimental Data

There are few studies on multi-VRP taking into account customers’ satisfaction. Among the studies that have taken into account customer satisfaction, most of them are evaluated using randomly generated test cases. Therefore, there is no standard test cases library. The tests data used in this research is from the benchmark problems in the standard example library of MDVRPTW (multiple depot vehicle routing problem with time windows), and all examples can be downloaded fromhttp://neo.lcc.uma.es/radi-aeb/WebVRP/. [Ei,Li] is the time windows in initial data, as the tolerable time in this paper. [ai,bi] is the desirable service time, which can be computed using the following formula [27]:
(4.1)ai=Ei+rand()*0.5*(Li-Ei),bi=Li-rand()*0.5*(Li-Ei).


### 4.2.                       Parameters Discussion of MOQEA

The parameters involved in the MOQEA include coefficientk in formula (2.1) and constant s in formula (3.7). The proposed MOQEA for different parameters was discussed and analyzed results are shown in Tables 2 and 3. VN is the vehicle numbers. CS is the customer nonsatisfaction. These tables only list the solution that the total cost is the least. From the tables we can see when s=0.2 and k=0.05, the resultant vehicle number is the smallest.Table 2
Optimization results of coefficientk (s=0.2).
Problem
pr02 (96, 4)
pr07 (72, 6)
k
VN
CS
Cost
VN
CS
Cost
0.08
16
0.402
4539
15
0.492
4116
0.06
15
0.378
4045
14
0.424
3904
0.05
14
0.355
3906
13
0.399
3331
0.02
15
0.364
3998
14
0.431
3969Table 3
Optimization results of constants (k=0.06).
Problem
pr02 (96, 4)
pr07 (72, 6)
s
VN
CS
Cost
VN
CS
Cost
0.1
15
0.381
4013
14
0.435
3391
0.2
14
0.357
3902
13
0.389
3367
0.3
15
0.361
3926
15
0.486
4124
0.4
16
0.411
4625
15
0.498
4102


### 4.3. Simulation Experiments

All the programs in this research are developed using the JAVA language and run on a PC with Dual 2.8 GHz CPU and 1.0 GB of memory. A manufacturing company has 4 warehouses and provides goods to 48 vendors. The actual distribution process can be attributed to the open, capacity constraints, and multidepot VRP. The capacity is 20 t. The proposed MOQEA is used to solve this problem, and the distance and demand of each client and depot are shown in Tables4 and 5. CN is the customer NO. ST is the service time. De is the demand. TW is the time windows.Table 4
The distance and demand of each client.
CN
Coordinate/km
ST
De/t
TW
1
−29.730
64.136
2
12
[399525]
2
−30.664
5.463
7
8
[121299]
3
51.642
5.469
21
16
[389483]
4
−13.171
69.336
24
5
[204304]
5
−67.413
68.323
1
12
[317458]
6
48.907
6.274
17
5
[160257]
7
5.243
22.260
6
13
[170287]
8
−65.002
77.234
5
20
[215321]
9
−4.175
−1.569
7
13
[80233]
10
23.029
11.639
1
18
[90206]
11
25.482
6.287
4
7
[397525]
12
−42.615
−26.392
10
6
[271 420]
13
−76.672
99.341
2
9
[108 266]
14
−20.673
57.892
16
9
[340 462]
15
−52.039
6.567
23
4
[226 377]
16
−41.376
50.824
18
25
[446 604]
17
−91.943
27.588
3
5
[444 566]
18
−65.118
30.212
15
17
[434 557]
19
18.597
96.716
13
3
[319 460]
20
−40.942
83.209
10
16
[192 312]
21
−37.756
−33.325
4
25
[414 572]
22
23.767
29.083
23
21
[371 462]
23
−43.030
20.453
20
14
[378 472]
24
−35.297
−24.896
10
19
[308 477]
25
−54.755
14.368
4
14
[329 444]
26
−49.329
33.374
2
6
[269 377]
27
57.404
23.822
23
16
[398 494]
28
−22.754
55.408
6
9
[257 416]
29
−56.622
73.340
8
20
[198 294]
30
−38.562
−3.705
10
13
[375 467]
31
−16.779
19.537
7
10
[200 338]
32
−11.560
11.615
1
16
[456 632]
33
−46.545
97.974
21
19
[72  179]
34
16.229
9.320
6
22
[182 282]
35
1.294
7.349
4
14
[159 306]
36
−26.404
29.529
13
10
[321 500]
37
4.352
14.685
9
11
[322 430]
38
−50.665
−23.126
22
15
[443 564]
39
−22.833
−9.814
22
13
[207 348]
40
−71.100
−18.616
18
15
[457 588]
41
−7.849
32.074
10
8
[203 382]
42
11.877
−24.933
25
22
[75 167]
43
−18.927
−23.730
23
24
[459 598]
44
−11.920
11.755
4
3
[174 332]
45
29.840
11.633
9
25
[130 225]
46
12.268
−55.811
17
19
[169 283]
47
−37.933
−21.613
10
21
[115 232]
48
42.883
−2.966
17
10
[414 531]Table 5
The distance of each depot.
Depot no.
Coordinate/km
Service time
Demand/t
Time windows
49
4.163
13.559
0
0
[0 1000]
50
21.387
17.105
0
0
[0 1000]
51
−36.118
49.097
0
0
[0 1000]
52
−31.201
0.235
0
0
[0 1000]The results obtained are shown in Table6 and Figure 7. Specifically, the number of iterations is 2000, and the population size is 30. The coefficient k is set as 0.06. The constant s in formula (4.1) is 0.2. The Pareto optimal solution set is (630, 3523), (593, 3879), (556, 4164), (543, 4203), (538, 4261), (525, 4302), (516, 4365), (495, 4498), (472, 4532). The value of non-satisfaction magnified 1000 times. It can be seen from Figure 7 that the Pareto front that is close to the axis forms a more satisfactory solution set. Comparing the leftmost and the rightmost points, we can see that in this instance, it would be possible to decrease the total cost by 20% at the expense of an increase in the non-satisfaction which is about 25%.Table 6
Optimization results of pr01.
Pareto set
The solution of maximum satisfaction
The solution of minimizing travel cost, fixed cost, and waiting cost
(472, 4532)
(630, 3523)
VN
9
11
Route
52 21 36 31 27 6 24
52 43 21 14 13
51 28 13 44 25
52 30 45 44 38
51 5 39 17
50 6 40 48 47 17
49 32 19 20 18 42
49 7 37 11 32 26 16 36 42
51 4 11 41 26 1 29 22
52 12 41 15 2 4
51 14 12 43 16 33
51 33 18 3
50 3 45 8 40 10 37
51 1 35 20 39
52 23 34 46 38 30 15 47
51 8 29 9 22
49 7 9 35 48 2
51 5 34 24 23
51 19 27 28 10 31 25
49 46Figure 7
Pareto optimal solution set.


### 4.4. Comparison and Discussion

In order to evaluate the performance of the algorithm, the proposed MOQEA is compared with the hybrid multiobjective evolutionary algorithm (HMOEA) developed in [34]. In the HMOEA, feasible individuals are constructed as the initial population by using the push-forward insertion heuristic (PFIH), and the GA is used to update these populations to obtain the new subpopulation and to improve the individuals of the subpopulation by the local search method of λ-interchange with variable probability, then non-dominated solution set is constructed by using the Challenge Cup rule.Table7 shows the comparison of the results obtained from the two algorithms. Because the calculation result is a solution set, this table only lists the solutions in which the total cost is the least. VN is the vehicle numbers. CS is the customer non-satisfaction. From the table we can see for these two multidepot VRPs obtained from the proposed MOQEA the total cost is smaller than that from the HMOEA, and the customers’ satisfaction obtained from the MOQEA is greater than that from the HMOEA.Table 7
Comparisons of the MOQEA to the HMOEA in [34].
Algorithm
Project 1
Project 2
MOQEA
HMOEA in [34]
Problems
VN
CS
Cost
VN
CS
Cost
pr02 (96, 4)
14
0.355
3906
15
0.391
4005
pr07 (72, 6)
13
0.399
3331
13
0.402
3397


## 4.1. Experimental Data

There are few studies on multi-VRP taking into account customers’ satisfaction. Among the studies that have taken into account customer satisfaction, most of them are evaluated using randomly generated test cases. Therefore, there is no standard test cases library. The tests data used in this research is from the benchmark problems in the standard example library of MDVRPTW (multiple depot vehicle routing problem with time windows), and all examples can be downloaded fromhttp://neo.lcc.uma.es/radi-aeb/WebVRP/. [Ei,Li] is the time windows in initial data, as the tolerable time in this paper. [ai,bi] is the desirable service time, which can be computed using the following formula [27]:
(4.1)ai=Ei+rand()*0.5*(Li-Ei),bi=Li-rand()*0.5*(Li-Ei).


## 4.2.                       Parameters Discussion of MOQEA

The parameters involved in the MOQEA include coefficientk in formula (2.1) and constant s in formula (3.7). The proposed MOQEA for different parameters was discussed and analyzed results are shown in Tables 2 and 3. VN is the vehicle numbers. CS is the customer nonsatisfaction. These tables only list the solution that the total cost is the least. From the tables we can see when s=0.2 and k=0.05, the resultant vehicle number is the smallest.Table 2
Optimization results of coefficientk (s=0.2).
Problem
pr02 (96, 4)
pr07 (72, 6)
k
VN
CS
Cost
VN
CS
Cost
0.08
16
0.402
4539
15
0.492
4116
0.06
15
0.378
4045
14
0.424
3904
0.05
14
0.355
3906
13
0.399
3331
0.02
15
0.364
3998
14
0.431
3969Table 3
Optimization results of constants (k=0.06).
Problem
pr02 (96, 4)
pr07 (72, 6)
s
VN
CS
Cost
VN
CS
Cost
0.1
15
0.381
4013
14
0.435
3391
0.2
14
0.357
3902
13
0.389
3367
0.3
15
0.361
3926
15
0.486
4124
0.4
16
0.411
4625
15
0.498
4102


## 4.3. Simulation Experiments

All the programs in this research are developed using the JAVA language and run on a PC with Dual 2.8 GHz CPU and 1.0 GB of memory. A manufacturing company has 4 warehouses and provides goods to 48 vendors. The actual distribution process can be attributed to the open, capacity constraints, and multidepot VRP. The capacity is 20 t. The proposed MOQEA is used to solve this problem, and the distance and demand of each client and depot are shown in Tables4 and 5. CN is the customer NO. ST is the service time. De is the demand. TW is the time windows.Table 4
The distance and demand of each client.
CN
Coordinate/km
ST
De/t
TW
1
−29.730
64.136
2
12
[399525]
2
−30.664
5.463
7
8
[121299]
3
51.642
5.469
21
16
[389483]
4
−13.171
69.336
24
5
[204304]
5
−67.413
68.323
1
12
[317458]
6
48.907
6.274
17
5
[160257]
7
5.243
22.260
6
13
[170287]
8
−65.002
77.234
5
20
[215321]
9
−4.175
−1.569
7
13
[80233]
10
23.029
11.639
1
18
[90206]
11
25.482
6.287
4
7
[397525]
12
−42.615
−26.392
10
6
[271 420]
13
−76.672
99.341
2
9
[108 266]
14
−20.673
57.892
16
9
[340 462]
15
−52.039
6.567
23
4
[226 377]
16
−41.376
50.824
18
25
[446 604]
17
−91.943
27.588
3
5
[444 566]
18
−65.118
30.212
15
17
[434 557]
19
18.597
96.716
13
3
[319 460]
20
−40.942
83.209
10
16
[192 312]
21
−37.756
−33.325
4
25
[414 572]
22
23.767
29.083
23
21
[371 462]
23
−43.030
20.453
20
14
[378 472]
24
−35.297
−24.896
10
19
[308 477]
25
−54.755
14.368
4
14
[329 444]
26
−49.329
33.374
2
6
[269 377]
27
57.404
23.822
23
16
[398 494]
28
−22.754
55.408
6
9
[257 416]
29
−56.622
73.340
8
20
[198 294]
30
−38.562
−3.705
10
13
[375 467]
31
−16.779
19.537
7
10
[200 338]
32
−11.560
11.615
1
16
[456 632]
33
−46.545
97.974
21
19
[72  179]
34
16.229
9.320
6
22
[182 282]
35
1.294
7.349
4
14
[159 306]
36
−26.404
29.529
13
10
[321 500]
37
4.352
14.685
9
11
[322 430]
38
−50.665
−23.126
22
15
[443 564]
39
−22.833
−9.814
22
13
[207 348]
40
−71.100
−18.616
18
15
[457 588]
41
−7.849
32.074
10
8
[203 382]
42
11.877
−24.933
25
22
[75 167]
43
−18.927
−23.730
23
24
[459 598]
44
−11.920
11.755
4
3
[174 332]
45
29.840
11.633
9
25
[130 225]
46
12.268
−55.811
17
19
[169 283]
47
−37.933
−21.613
10
21
[115 232]
48
42.883
−2.966
17
10
[414 531]Table 5
The distance of each depot.
Depot no.
Coordinate/km
Service time
Demand/t
Time windows
49
4.163
13.559
0
0
[0 1000]
50
21.387
17.105
0
0
[0 1000]
51
−36.118
49.097
0
0
[0 1000]
52
−31.201
0.235
0
0
[0 1000]The results obtained are shown in Table6 and Figure 7. Specifically, the number of iterations is 2000, and the population size is 30. The coefficient k is set as 0.06. The constant s in formula (4.1) is 0.2. The Pareto optimal solution set is (630, 3523), (593, 3879), (556, 4164), (543, 4203), (538, 4261), (525, 4302), (516, 4365), (495, 4498), (472, 4532). The value of non-satisfaction magnified 1000 times. It can be seen from Figure 7 that the Pareto front that is close to the axis forms a more satisfactory solution set. Comparing the leftmost and the rightmost points, we can see that in this instance, it would be possible to decrease the total cost by 20% at the expense of an increase in the non-satisfaction which is about 25%.Table 6
Optimization results of pr01.
Pareto set
The solution of maximum satisfaction
The solution of minimizing travel cost, fixed cost, and waiting cost
(472, 4532)
(630, 3523)
VN
9
11
Route
52 21 36 31 27 6 24
52 43 21 14 13
51 28 13 44 25
52 30 45 44 38
51 5 39 17
50 6 40 48 47 17
49 32 19 20 18 42
49 7 37 11 32 26 16 36 42
51 4 11 41 26 1 29 22
52 12 41 15 2 4
51 14 12 43 16 33
51 33 18 3
50 3 45 8 40 10 37
51 1 35 20 39
52 23 34 46 38 30 15 47
51 8 29 9 22
49 7 9 35 48 2
51 5 34 24 23
51 19 27 28 10 31 25
49 46Figure 7
Pareto optimal solution set.


## 4.4. Comparison and Discussion

In order to evaluate the performance of the algorithm, the proposed MOQEA is compared with the hybrid multiobjective evolutionary algorithm (HMOEA) developed in [34]. In the HMOEA, feasible individuals are constructed as the initial population by using the push-forward insertion heuristic (PFIH), and the GA is used to update these populations to obtain the new subpopulation and to improve the individuals of the subpopulation by the local search method of λ-interchange with variable probability, then non-dominated solution set is constructed by using the Challenge Cup rule.Table7 shows the comparison of the results obtained from the two algorithms. Because the calculation result is a solution set, this table only lists the solutions in which the total cost is the least. VN is the vehicle numbers. CS is the customer non-satisfaction. From the table we can see for these two multidepot VRPs obtained from the proposed MOQEA the total cost is smaller than that from the HMOEA, and the customers’ satisfaction obtained from the MOQEA is greater than that from the HMOEA.Table 7
Comparisons of the MOQEA to the HMOEA in [34].
Algorithm
Project 1
Project 2
MOQEA
HMOEA in [34]
Problems
VN
CS
Cost
VN
CS
Cost
pr02 (96, 4)
14
0.355
3906
15
0.391
4005
pr07 (72, 6)
13
0.399
3331
13
0.402
3397


## 5. Conclusions

This paper presents the modeling of vehicle scheduling problem that takes into account customer satisfaction and the development of the MVRPCS. Specifically, an improved trapezoidal fuzzy number is proposed to represent customers’ satisfaction and the MOQEA for this problem is developed. The MOQEA can get multiple solutions, namely, the Pareto optimal solution set, according to his own expectations. These solutions will be used by the decision maker to choose the best one on the basis of different preferences on satisfaction maximization and travel costs minimization. In the MOQEA, the Challenge Cup rule is constructed for non-dominated solution set and a method for attaining keeping the variety of the solution set, is designed, based on self-adaptive grid. Simulation results and comparisons show that the MOQEA is an effective method. In our future work, we will focus on improving the algorithm and test it on other datasets.



---

*Source: 879614-2012-12-10.xml*