1. travel

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
IEEE TRANSACTIONS ON CYBERNETICS
1
Multiobjective Vehicle Routing Problems With
Simultaneous Delivery and Pickup and Time
Windows: Formulation, Instances,
and Algorithms
Jiahai Wang, Member, IEEE, Ying Zhou, Yong Wang, Member, IEEE, Jun Zhang, Senior Member, IEEE,
C. L. Philip Chen, Fellow, IEEE, and Zibin Zheng, Member, IEEE
Abstract—This paper investigates a practical variant of the
vehicle routing problem (VRP), called VRP with simultaneous
delivery and pickup and time windows (VRPSDPTW), in
the logistics industry. VRPSDPTW is an important logistics
problem in closed-loop supply chain network optimization.
VRPSDPTW exhibits multiobjective properties in realworld applications. In this paper, a general multiobjective
VRPSDPTW (MO-VRPSDPTW) with five objectives is first
defined, and then a set of MO-VRPSDPTW instances based
on data from the real-world are introduced. These instances
represent more realistic multiobjective nature and more
challenging MO-VRPSDPTW cases. Finally, two algorithms,
multiobjective local search (MOLS) and multiobjective
memetic algorithm (MOMA), are designed, implemented
and compared for solving MO-VRPSDPTW. The simulation
results on the proposed real-world instances and traditional
instances show that MOLS outperforms MOMA in most of
instances. However, the superiority of MOLS over MOMA
in real-world instances is not so obvious as in traditional
instances.
Manuscript received February 16, 2015; revised February 16, 2015;
accepted February 27, 2015. This work was supported in part by the National
High-Technology Research and Development Program (863 Program) of
China under Grant 2013AA01A212, in part by the National Natural
Science Foundation of China (NSFC) for Distinguished Young Scholars
under Grant 61125205, in part by the NSFC under Grant 61332002,
Grant 61300044, and Grant 61273314, and in part by the Program for
New Century Excellent Talents in University under Grant NCET-13-0596.
This paper was recommended by Associate Editor J. Wang.
J. Wang is with the Department of Computer Science, Sun Yat-sen
University, Guangzhou 510006, China (e-mail: wjiahai@hotmail.com).
Y. Zhou is with the Department of Computer Network Technology,
Shenzhen Institute of Information Technology, Shenzhen 518172,
China.
Y. Wang is with the School of Information Science and Engineering,
Central South University, Changsha 410083, China.
J. Zhang is with the Sun Yat-sen University, Guangzhou 510006, China, also
with the Key Laboratory of Machine Intelligence and Advanced Computing,
Ministry of Education, China, also with the Engineering Research Center
of Supercomputing Engineering Software, Ministry of Education, China, and
also with the Key Laboratory of Software Technology, Education Department
of Guangdong Province, Guangzhou 510006, China.
C. L. P. Chen is with the Faculty of Science and Technology, University
of Macau, Macau 99999, China.
Z. Zheng is with the Shenzhen Research Institute, The Chinese University
of Hong Kong, Shenzhen 518057, China.
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TCYB.2015.2409837
Index Terms—Bi-directional logistics, multiobjective optimization, simultaneous delivery and pickup, vehicle routing problem
with time windows (VRPTW).
I. I NTRODUCTION
ECENTLY, green manufacturing and logistics have
emerged as the new agenda item in supply chain
management [1]–[3]. They have received increasing attention from governments and business organization. One of
the actions taken by manufacturing companies toward green
manufacturing is to collect end-of-life products from customers for either reuse or proper disposal, which is known as
reverse logistics [3]–[5]. Economics, environmental laws, and
the environmental consciousness of consumers are the driving factors for adopting reverse logistics concepts in industry.
Depending on the nature of returned products, one option is
to design combined distribution-collection (delivery-pickup)
systems. For example, in the distribution system of grocery
store chains, each grocery store may have demand for both
delivery (fresh food or soft drinks) and pickup (outdated
items or empty bottles) and is served with a single stop
by the supplier. In the foundry industry, collection of used
sand and delivery of purified reusable sand at the same customer location are carried out with only a single stop. In the
printer manufacturing industry, full ink toners and cartridges
are delivered and empty ones are collected. In the photocopier manufacturing industry, manufacturers are required to
take back or properly dispose of end-of-life products. In these
cases, the utilization of vehicles increases significantly when
merging products brought to the customers (forward logistics)
with returning products brought back to the depot (reverse
logistics). Thus, the vehicle routing and the flows of freights
become more effective and balanced in the bi-directional
logistics [2].
In this paper, a more realistic and general variant of the
vehicle routing problem (VRP), called VRP with simultaneous delivery and pickup with time windows (VRPSDPTW), is
considered. VRPSPDTW considers simultaneous pickup and
delivery at each customer such that a customer is visited only
once within the specified time window and without violating
R
c 2015 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/
2168-2267 redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
2
the vehicle capacity constraints [4]–[6]. VRPSPDTW is a
challenging combinatorial optimization problem, containing
complex constraints not present in classic VRP [1], [7].
VRPSPDTW is NP-hard because it contains VRP as a
special case [7]. Thus, practical large-scale instances cannot
be solved by exact methodologies within acceptable computational time [5], and most researchers have focused on
metaheuristic approaches [8]–[10].
At present, studies on VRPSDPTW remain scarce because
the pickup and time window constraints make VRPSDPTW
more difficult. Angelelli and Mansini [11] made a first attempt
to solve VRPSPDTW via an exact algorithm which can solve
problems with up to 20 customers. Lai and Cao [6] proposed an improved differential evolution for VRPSDPTW
and carried out some experiments on their own small size
instances (with only 8 and 40 customers). Boubahri et al. [12]
proposed a multiagent colonies algorithm for VRPSDPTW,
but their method has not been tested on any instance.
Wang and Chen [5] proposed a co-evolution genetic algorithm
for VRPSPDTW, and developed 65 benchmark instances [13]
revised from the well-known Solomon benchmark for VRP
with time windows (VRPTW) [14]. Recently, Wang et al. [15]
and Kassem and Chen [4] adopted simulated annealing to deal
with VRPSPDTW. Deng et al. [16] proposed an improved
simulated annealing for VRPSDPTW with soft time windows.
In [17], a special VPRSDPTW with four types of demands
in home health care logistics was considered and solved by
genetic algorithm and tabu search.
Among the existing work mentioned, some researches study
on single-objective VRPSDPTW. For example, in [4] and [6],
the total travel distance is considered as a sole objective.
Considering VRPSDPTW with multiple objectives (often two
objectives including the number of vehicles and total distance),
most previous studies transform it to a single-objective optimization problem, and thus adopt single-objective approaches
to solve it and return a single solution. One of the most popular
techniques is to use scalar techniques. In [5], [15], and [16],
several algorithms are proposed to optimize a weighted linear
aggregation function of objectives. However, this kind of techniques needs to set the weights according to the importance
of the objectives, which is a difficult task [18].
Due to the constraints and problem structure of
VRPSDPTW, the optimization of one objective may lead to the
deterioration of other objectives, thus VRPSDPTW is essentially a multiobjective optimization problem (MOP) [18]–[20].
Since the decision maker’s preference is not known a priori,
multiobjective formulation is necessary for VRPSDPTW
to provide a set of solutions that represent the tradeoffs
among the objectives, rather than a single solution [18]–[20].
The feature of multiobjective formulation is to consider all
objectives with the same importance and obtain a set of
Pareto optimal solutions. To the authors’ best knowledge, no
previous work utilizes the multiobjective optimization method
for VRPSDPTW, which motivates this paper.
This paper first defines a general multiobjective
VRPSDPTW (MO-VRPSDPTW) with five objectives commonly used in the VRP literature. These objectives include:
number of vehicles, total travel distance, makespan
IEEE TRANSACTIONS ON CYBERNETICS
[i.e., travel time of the longest route (from/to depot)],
total waiting time, and total delay time [20]. Then this
paper proposes a set of MO-VRPSDPTW instances based
on data from a distribution company in Tenerife, Spain.
These MO-VRPSDPTW instances are quite different from
the revised Solomon instances in [5]. Finally, two algorithms, called multiobjective local search (MOLS) and
multiobjective memetic algorithm (MOMA), are designed,
implemented, and compared. The usefulness of the proposed
MO-VRPSDPTW formulation and algorithms are demonstrated by solving two sets of benchmark instances. The
proposed algorithms can be seen as benchmark algorithms
for the real-world MO-VRPSDPTW instances, which can be
used for comparison by future research.
The contribution of this paper is fourfold: 1) introducing a
five-objective variant of VRPSDPTW; 2) introducing a set of
realistic benchmark instances; 3) designing and testing multiobjective optimization algorithms to solve the five-objective
VRPSDPTW; and 4) performing extensive experiments to
evaluate the two proposed algorithms.
The remaining sections are organized as follows. Section II
presents problem formulation and benchmark instances of
MO-VRPSDPTW. Section III proposes two algorithms for
MO-VRPSDPTW. Section IV presents experimental results.
Finally, Section V presents the conclusion.
II. P ROBLEM F ORMULATION AND
B ENCHMARK I NSTANCES
A. MO-VRPSDPTW
Given a number of customers who require both forward
supply service and reverse recycling service within a given
time window, MO-VRPSDPTW considered in this paper aims
to send out a fleet of capacitated vehicles at a distribution
center to meet the requests with minimum costs or objectives.
Let us introduce the following nomenclature [21].
v = {0, . . . , N}
Vertex 0 is called depot and the
others are called customers.
C
Vehicle capacity.
Delivery demand of customer i.
gi
Pickup amount of recycling of
pi
customer i.
Service time of customer i.
si
[bi , ei ]
Time window of customer i,
where bi and ei denote the earliest and latest service time of
customer i, respectively.
Travel distance between cusdij
tomers i and j.
Travel time between customers
tij
i and j.
rj =< c(1, j), . . . , c(Nj , j)>Sequence of Nj customers
served in the jth route, where
c(i, j) denotes the ith customer
to be visited in the jth route.
For computational convenience,
let c(0, j) = c(Nj + 1, j) = 0
represents the depot.
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
WANG et al.: MO-VRPSDPTW: FORMULATION, INSTANCES, AND ALGORITHMS
Fig. 1.
Solution and its representation. (a) Solution. (b) Representation.
Note that the depot also has a time window [0, e0 ]. Delivery
demand and pickup amount of the depot are g0 = 0 and
p0 = 0, respectively. The aim of VRPSDTW is to design
a set of M routes (i.e., R = {r1 , . . . , rM }) with the lowest cost
such that each vehicle departs from and returns to the same
depot and each customer is served by exactly one vehicle.
As shown in Fig. 1(a), there are three routes (i.e., M = 3):
R = {r1 , r2 , r3 }. r1 = <c(1, 1), c(2, 1)> specifies the sequence
of two customers (customers 2 and 7) served in route 1, that
is, r1 = <2, 7>. In addition, the remaining two routes have
three and two customers served, respectively.
The total travel distance of the jth route is defined as
Distj =
Nj
dc(i,j)c(i+1,j) .
(1)
3
paper, soft time windows are considered as in [16] and [20].
VRPSDPTW with soft time windows can be viewed as a generalization of VRPSDPTW with hard time windows. It has
many practical applications. In many cases, relaxing the time
windows may be more appropriate [24], [25] because it may
result in lower cost solutions requiring fewer vehicles, shorter
travel distance, and less travel time. In VRPSDPTW with soft
time windows, a vehicle can arrive late within the maximum
allowed time. Arriving outwith the maximum allowed time is
not allowed. Let md denote the maximum allowed time that a
vehicle can arrive after the end of the time window.
The delay time of vehicle j at the ith vertex is
0,
if ac(i,j) ≤ ec(i,j)
(7)
dtc(i,j) =
ac(i,j) − ec(i,j) , otherwise
and the total delay time of this route is
Nj
DTj =
ac(i,j) = lc(i−1,j) + tc(i−1,j)c(i,j)
(2)
where lc(0,j) = 0, which indicates that vehicle j departs from
the depot at time 0. If a vehicle arrives at a customer before the
earliest service time, it will cause waiting time. The waiting
time of vehicle j at the ith customer can be described as
0,
if ac(i,j) ≥ bc(i,j)
(3)
wc(i,j) =
bc(i,j) − ac(i,j) , otherwise
f1 = |R| = M.
(4)
Hence, the total travel time of route rj is
Nj
Tj =
tc(i,j)c(i+1,j) + wc(i+1,j) + sc(i+1,j)
2) Total travel distance ( f2 )
(5)
and the total waiting time of this route is
Nj
wc(i,j) .
f2 =
M
Distj .
(10)
j=1
3) Makespan ( f3 ), i.e., travel time of the longest route
(from/to depot)
f3 = max{Tj | j = 1 . . . M}.
j
(11)
4) Total waiting time due to early arrivals ( f4 )
f4 =
M
Wj .
(12)
j=1
f5 =
M
DTj .
(13)
j=1
i=0
Wj =
(9)
5) Total delay time due to late arrivals ( f5 )
and the leaving time from the ith customer of vehicle j is
lc(i,j) = ac(i,j) + wc(i,j) + sc(i,j) .
(8)
Based on the introduction above, the objectives for
MO-VRPSDPTW can be defined as follows.
1) Number of vehicles ( f1 )
i=0
Let ac(i−1,j) be the arrival time of vehicle j at the (i − 1)th
vertex and lc(i−1,j) be the leaving time. Thus, the arrival time
of vehicle j at the ith vertex is
dtc(i,j) .
i=1
(6)
i=1
In VRPSDPTW with hard time windows [5], arriving after
the latest service time is not allowed. In practice, the hard
time windows are often relaxed because travel time is usually stochastic and cannot be precisely predictable [22], [23].
Additionally, a small time window violation is usually not
a critical breach of service requirements. Hence, in this
Minimization of f1 aims to reduce the fixed costs of buying
(or hiring) and maintaining vehicles. Variable costs considered in routing and distribution/collection problems are often
estimated by using a function of the total distance traveled
(denoted as f2 ). Hence, f1 and f2 can be considered as transportation costs (i.e., economic objectives). f3 can be considered
as a social objective. This objective has a twofold contribution toward social sustainability. On the one hand, it is used
to minimize the maximum working hours among all drivers,
thus enabling a balanced workload to promote equity among
the drivers. On the other hand, minimizing the maximum
working hours releases drivers to activities such as recycling
awareness campaigns or training, which is helpful for improving the career development and promoting versatility among
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
4
IEEE TRANSACTIONS ON CYBERNETICS
human resources [3]. Minimization of f4 improves work efficiency and avoids wasting working hours. Finally, f5 can be
considered as a service cost related to the satisfaction of
customers [22]. Needless to say, customers want to be served
just in time and cannot tolerate severe delays.
In addition, the constraints for MO-VRPSDPTW can be
defined as follows.
1) Vehicle capacity constraint: The total delivery demand of
each route rj should not exceed the vehicle capacity, that is
Nj
gc(i,j) ≤ C
∀j = 1, . . . , M.
(14)
i=1
When a vehicle arrives at customer c(i, j), its load en route
is denoted by Cc(i,j) , then the following constraint should be
satisfied:
Cc(i,j) − gc(i,j) + pc(i,j) ≤ C ∀i = 1, . . . , Nj
∀j = 1, . . . , M.
(15)
2) Travel time constraint: Delay time should not exceed the
maximum allowed delay time
dtc(i,j) ≤ md ∀i = 1, . . . , Nj ,
∀j = 1, . . . , M.
(16)
3) Return time constraint: Vehicles should return to the
depot before the closing time, that is
ac(Nj +1,j) ≤ ec(Nj +1,j)
∀j = 1, . . . , M.
(17)
Thus, MO-VRPSDPTW with the five objectives considered
in this paper can be summarized as min F, where
F = { f1 , f2 , f3 , f4 , f5 }
(18)
is subject to (14)–(17), and minimization of the vector function F is supposed here in the sense of Pareto optimization
(see Section III).
B. Real-World Instances
References [4] and [5] generated VRPSDPTW instances
revised from Solomon benchmark instances originally
designed for VRPTW [14]. But there are two disadvantages.
Firstly, the Solomon dataset relies on Euclidean distance for
both travel distance and travel time, and one unit of travel time
corresponds to one unit of travel distance. Thus, the distance
and time matrices are the same and symmetric. This is hardly
a realistic scenario because the travel time is often not directly
proportional to the travel distance. Triangle inequality holds
for the travel distance and the travel time in the Solomon
dataset. Secondly, Solomon dataset is not suitable to conduct a proper multiobjective study because weak dependency
relationships are presented among objectives in the dataset.
These disadvantages motivate us to propose new instances.
Our new instances are based on the real-world multiobjective VRPTW (MO-VRPTW) instances recently proposed
by Castro-Gutierrez et al. [20]. The MO-VRPTW instances
are based on data from a distribution company in Tenerife,
Spain, [26]. These MO-VRPTW instances are quite different
from the Solomon instances [14]. On the one hand, the distance and time matrices are distinct and nonsymmetric in these
real-world instance, hence, they represent a realistic trade-off
between travel distance and travel time. The effects of the
asymmetry on realistic VRP are studied in [27]. Moreover,
the triangle inequality violations for both travel distance
and travel time are prevalent in the real-world MO-VRPTW
instances. The effects of the triangle inequality violations on
VRP have also been studied in [28]. The triangle inequality violation for travel time has an effect on the algorithm
design in this paper because it directly affects the feasibility of time constraints of a modified route. On the other
hand, strong dependency relationships are presented among
objectives in the real-world MO-VRPTW instances. Thus the
real-world MO-VRPTW instances exhibit more realistic multiobjective nature and are suitable to assess the performance
of multiobjective optimization algorithms [20], [26].
To generate our real-world MO-VRPSDPTW instances, the
MO-VRPTW instances are modified by adding pickup quantity for the customers. The pickup quantity pi for customer i
is set to 10, 20, or 30, each with probability 1/3 in the
MO-VRPSDPTW instances. The proposed instances can be
downloaded from our website [29]. Details will be described
in Section IV. Although the modification seems to be minor,
a new real-world benchmark dataset for MO-VRPSDPTW is
proposed in this paper for the first time. The multiobjective
nature of MO-VRPSDPTW can be better investigated and multiobjective optimization algorithms could be tested in a more
reasonable manner using this realistic dataset.
III. P ROPOSED A LGORITHMS FOR MO-VRPSDPTW
A. Overview of Multiobjective Optimization
For better understanding our algorithms and empirical
studies, this section briefly reviews the basic concepts and
algorithms of MOPs. In general, a MOP can be defined as
follows:
min F(x) = { f1 (x), . . . , fm (x)}
(19)
subject to x ∈ , where is the solution space. F : → Rm
consists of m objective functions that often conflict with each
other. Let x, y ∈ , y is said to Pareto dominate x if and only if
fi (y) ≤ fi (x) for every i ∈ 1, . . . , m, and fj (y) < fj (x) for at least
one j ∈ 1, . . . , m. A solution x∗ is Pareto optimal if there is no
solution x ∈ such that x Pareto dominates x∗ . In this case, x∗
is also called a nondominated solution. F(x∗ ) is called a Pareto
optimal objective vector. The set of all Pareto optimal solutions
is called Pareto set, and the set of all Pareto optimal objective
vectors is called Pareto front. The goal of a multiobjective
optimization algorithm for a MOP is to seek a set of solutions
which perform well in terms of convergence and diversity.
That is, the solutions obtained by a multiobjective optimization
algorithm should be close to and well-distributed along the
Pareto front.
Over the years, a number of metaheuristics have
been extended to solve MOPs. Multiobjective evolutionary algorithms (MOEAs) strive to obtain an accurate and
well-distributed approximation of the true Pareto front.
Popular MOEAs include nondominated sorting genetic
algorithm II (NSGA-II) [30] and MOEA based on decomposition (MOEA/D) [31]. More studies about the development
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
WANG et al.: MO-VRPSDPTW: FORMULATION, INSTANCES, AND ALGORITHMS
and application of MOEAs can be found in [32]–[35]. Besides
MOEAs, local search-based algorithms, such as Pareto local
search [36] and multidirectional local search [37], are promising alternative approaches to solve MOPs. The merit of
local search-based algorithms is that problem-specific knowledge can be directly used to guide the search toward the
Pareto front. Thus, they are specially suitable for multiobjective combinatorial optimization problems. More details
about local search-based algorithms can be found in [38].
Moreover, problem-specific heuristics for local search and
the evolutionary algorithm framework for global search are
often combined to maintain a good balance between local
search (exploitation) and global search (exploration) in multiobjective optimization. This kind of algorithm is often called
memetic algorithm [38]–[40].
In this paper, MOLS is firstly proposed for
MO-VRPSDPTW, and then MOMA is proposed. Note
that, in response to the particularities of a MOP, different
multiobjective optimization algorithms may differ in the
encoding scheme (responsible for the characterization of the
search space), objective function, and operators that depend
on the kind of encoding scheme adopted. As a consequence,
the proposed algorithms in this paper are different from the
previous studies [31], [32], [37] since our algorithms are
composed of dedicated modules, for example, local search
operators, to solve MO-VRPSDPTW.
B. Solution Representation for MO-VRPSDPTW
In the proposed algorithms, the solution representation
for MO-VRPSDPTW is based on variable-length solution
representation [21]. This representation has been shown in
Fig. 1(b). A solution consists of several routes, and each route
has a sequence of customers to be served. Note that each route
should begin and end with a vertex 0 which denotes the depot.
C. MOLS for MO-VRPSDPTW
The general framework of MOLS is presented in
Algorithm 1. This algorithm framework is successfully used to
solve real-world MO-VRPTW and obtains better results than
NSGA-II in [21]. It is extended to solve MO-VRPSDPTW
in this paper, and expected to also obtain good results.
MO-VRPSDPTW generalizes MO-VRPTW and it is conceptually a harder problem. The critical feature of MO-VRPSDPTW
is that both pickup and delivery activities should be carried out
simultaneously by the same vehicle. Hence, a mechanism that
checks the fluctuating load on the vehicle at each customer
should be imposed to prevent vehicle overload.
The main idea of MOLS is to use different local search procedures, called objectivewise local searches, to optimize different objectives of a given solution in parallel [21], [37], [41].
More specifically, for each objective obj, an objectivewise
local search LSobj is defined. LSobj (x) is performed to improve
solution x with respect to objective obj. In MOLS, an archive
A is initialized by generating several nondominated solutions.
Then the main loop consists of: 1) selecting a solution from
archive A; 2) using objectivewise local searches to improve
each objective of the selected solution; and 3) updating archive
5
Algorithm 1 MOLS
1: archive A = ∅
2: generate several nondominated solutions to initialize A /*Initialization*/
3: while running time ≤ maximum computation time do
4:
x = randomly select a solution from archive A
5:
for obj = 1 to 5 do
6:
perform objectivewise local search LSobj (x) /*Objectivewise local
searches*/
7:
update archive A using neighbor solutions generated during the
objectivewise local search. /*Archive updating scheme*/
8:
end for
9: end while
10: return A
A to keep all nondominated solutions found during the
search [37].
The main components in MOLS are described in detail as
follows.
1) Initialization: A solution is randomly generated as
follows. Firstly, we randomly select a customer and create
a route. Then, another customer is randomly selected to be
inserted into this route. If a customer cannot be inserted into
any existing route, a new route is created. This procedure is
repeated until all customers are inserted properly.
Five (i.e., the number of objectives of MO-VRPSDPTW)
solutions are generated by the procedure described above.
Then the ith solution (i.e., xi ) is optimized by the objectivewise local search of the ith objective [i.e., LSi (xi )]. As a result,
different solutions are improved along different objectives or
directions. Finally, the nondominated solutions are stored into
archive A. Since different objectives (directions) are selected to
improve different solutions, the resultant nondominated solutions in the initial archive are expected to scatter over the
objective space to a certain extent.
2) Archive Updating Scheme: Archive A is updated using
the resultant solutions during the objectivewise local searches.
However, the archive may contain more and more nondominated solutions during the search process since the number of
nondominated solutions increases dramatically with objective
number m. It is important to control the size of the archive. In
MOLS, the concept of -dominance [42], [43] is adopted. The
size of the -dominance archive is controlled by a parameter .
Each solution in the archive is assigned an identification array
B = {B1 , . . . , Bm }, where Bi = log fi /log(1 + ) [42].
The identification array divides the whole objective space
into hyper-boxes. Each hyper-box can be occupied by only one
solution, thereby it provides two properties: 1) well-distributed
solutions can be maintained and 2) the final archive size is
bounded [43]. The original -dominance archive often loses
extreme solutions [44]. In MOLS, extreme solutions are stored
during the search.
In
3) Objectivewise
Local
Searches
LSobj (x):
MO-VRPSDPTW, different objectives have different physical
meanings and characteristics, hence, objectivewise local
searches are designed for these objectives as follows.
A specific local search, Algorithm 2, is designed to deal
with f1 . First, the route which has the fewest customers is
selected. Then, we enumerate all customers in the selected
route to try to insert them into other possible routes. The
customer should be inserted into the first position where the
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
6
Algorithm 2 Local Search LS1 (x) for f1
1: flag = false
2: while !flag do
3:
select the route with the fewest customers from solution x
4:
enumerate all customers in the selected route to try to insert them into
other routes
5:
if a customer cannot be inserted into other routes properly then
6:
flag = true
7:
update A with x
8:
break and return x
9:
end if
10:
if all customers in the selected route are inserted into other routes
successfully then
11:
one vehicle is reduced, i.e., f1 (x) = f1 (x) − 1
12:
continue
13:
end if
14: end while
insertion does not delay the start time of other customers. If all
customers in the selected route are inserted into other routes
successfully, one vehicle can be reduced and the local search
procedure proceeds to reduce one more vehicle. The local
search procedure stops when a customer cannot be inserted
into other routes properly. The resultant solution is used to
update archive.
Local search procedures for f2 , . . . , f5 are described in
Algorithm 3. Based on the solution representation, three
neighborhood operators, N1 , N2 , and N3 , are introduced for
f2 , . . . , f5 , which are presented below.
N1 removes a random customer from a selected route, and
then reinserts it into the best position among all possible positions of all possible routes. This operator involves a basic
function, selectRoute to select a route, and a definition of the
best position to insert a customer. Both function and definition
are also used in N2 and N3 . Since different objective-wise local
searches are used to optimize their corresponding objectives,
the function, selectRoute, and the definition of the best position
are both based on different objectives. Function selectRoute
selects a route with the longest travel time for f3 , and randomly selects a route for f2 , f4 , and f5 . The definitions of
the best position to be inserted are described as follows for
different objectives.
f2 : The position which makes the resultant solution after
insertion have the lowest total distance.
f3 : The position which makes the resultant solution after
insertion have the lowest travel time.
f4 : The position which makes the resultant solution after
insertion have the lowest total waiting time.
f5 : The position which makes the resultant solution after
insertion have the lowest total delay time.
Given the function and definition above, this neighborhood
operator can be described in detail as follows: first, a route
is selected using the selectRoute function. Then, a customer,
randomly chosen from the selected route, is removed from
this route. Finally, the removed customer is reinserted into the
best position. The criterion of selecting position to be inserted
is according to the corresponding objective. Since the best
position is selected for insertion, the best improvement strategy
is adopted in this neighborhood operator.
N2 firstly removes a random number of customers from a
selected route. For example, for a selected route with seven
IEEE TRANSACTIONS ON CYBERNETICS
Algorithm 3 Local Search LSobj (x) for the objth Objective
(obj = 2, ..., 5)/*Objectivewise Local Search*/
1: depth = 1
2: while depth < MaxDepth do
3:
randomly select a neighborhood operator from (N1 , N2 , N3 )
4:
x = generate a neighbor solution of x using the selected neighborhood
operator
5:
update A with x
6:
if fobj (x ) < fobj (x) then
7:
x = x
8:
end if
9:
depth = depth + 1
10: end while
11: return x
customers, a random number, rn(1 ≤ rn ≤ 7) is generated, and
then rn customers from this route are randomly selected to be
removed. These removed customers are finally reinserted into
the best position. The definition of the best position are the
same as N1 .
N3 exchanges a sequence of customers (i.e., subtour)
between two routes, preserving the orientation of the
sequences. In this operator, a sequence of customers (subtour)
includes all customers after a selected customer. It can
be handled as a single customer. For example, route 1,
r1 = <1, 2, 3, 4, 5>, is selected using the selectRoute function, and a customer, for example customer 2, is randomly
chosen from this route. The sequence of customers after
customer 2 includes customers 3–5. They can be treated as a
single customer S1. Then, we try to insert S1 into the best position. Suppose that the best position is found in the place after
customer 6 in route 2, r2 = <6, 7, 8, 9, 10>. The sequence
of customers after customer 6 including customers 7–9 is also
handled as a single customer S2. Then, S2 is reinserted into
the place after customer 2 in route 1 to complete the exchange
process. Finally, two new routes, r1 = <1, 2, 7, 8, 9, 10> and
r2 = <6, 3, 4, 5>, are obtained. More details can be found
in [21].
For a solution x, one of the neighborhood operators
(N1 , N2 , and N3 ) is randomly selected to generate a neighbor solution x of x. Then, x will be replaced by x if fobj (x )
is better than fobj (x). Once a neighbor solution x is generated, archive A is updated to avoid missing any nondominated
solution during the search process.
N1 only makes small changes to the current solution and
carries out the search within a restricted part of the solution
space, facilitating the algorithm’s convergence. In contrast,
N2 and N3 make larger changes to the current solution and
guide the search to different areas of the solution space. Thus,
MOLS acts as variable neighborhood search [45]. The simple
but powerful mixed neighborhood structures and the stochastic
elements of MOLS allow effective diversification for escaping
from local minima.
4) Feasibility Checking: In this paper, only feasible solutions are considered in all neighborhood operators. Efficient
feasibility checking is necessary for time and capacity constraints to reduce computational complexity.
For time constraints, a slack Sc(i,j) , which defines the maximum time shift allowed in route j after the ith customer,
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
WANG et al.: MO-VRPSDPTW: FORMULATION, INSTANCES, AND ALGORITHMS
can be precomputed as follows [46], [47]:
⎧
ec(i,j) − ac(i,j) ,
for i = Nj + 1
⎪
⎪
⎨
Sc(i,j) =
(20)
min(ec(i,j) + md − ac(i,j)
⎪
⎪
⎩
wc(i,j) + Sc(i+1,j) ),
for i = Nj , . . . , 1.
When a customer is inserted between the (i − 1)th and ith
customers in the jth route, the insertion is legal if the difference between the arrival time at the ith customer location
after and before the insertion is less or equal to Sc(i,j) . Hence,
the computational complexity of checking a feasible position
is reduced to O(1). Furthermore, only the routes which are
changed by neighborhood operators need to be reevaluated.
The feasibility checking process for capacity constraint in
MO-VRPSDPTW has additional complexity because of the
load fluctuation of vehicles. In order to speed up the feasibility
checking process, special metrics [48] are used to capture the
load fluctuation of vehicles along their routes. Two copies of
depot are numbered as 0 and N + 1. Each customer c(k, j) is
associated with the following quantities.
−
−
and δc(k,j)
indicate the amount of load picked
1) πc(k,j)
up and load to deliver, respectively, on board of vehicle j, when vehicle j reaches customer c(k, j) with
k = 1, . . . , Nj + 1.
+
+
2) πc(k,j)
and δc(k,j)
indicate the amount of load picked
up and load to deliver, respectively, on board of
vehicle j, when vehicle j leaves customer c(k, j) with
k = 0, . . . , Nj .
−
−
−
= maxi=1,...,k {πc(i,j)
+ δc(i,j)
} with k = 1, . . . ,
3) Mc(k,j)
Nj + 1 indicates the maximum total load of vehicle j
since its departure from the depot to customer c(k, j).
+
+
+
= maxi=k,...,Nj {πc(i,j)
+ δc(i,j)
} with k = 0, . . . , Nj
4) Mc(k,j)
c(k, j) indicates the maximum total load of vehicle j from
its departure from customer c(k, j) until to the depot.
These quantities are precomputed.
N1 and N2 both insert a customer into a route. The following
two conditions must hold to satisfy the capacity constraint
when customer i is inserted between c(k, j) and c(k, j + 1):
−
+ gi ≤ C
Mc(k+1,j)
+
Mc(k,j)
+ pi ≤ C.
(21)
N3 exchanges a sequence of customers. An exchange operator can be decomposed into two insert operators. Suppose
that the exchange operator involves c(k1 , j1 ) and c(k1 + 1, j1 )
in route j1 , and c(k2 , j2 ) and c(k2 + 1, j2 ) in route j2 , the feasibility checking for capacity constraint is made by checking
the following conditions:
−
−
−
− δc(k
+ δc(k
≤C
Mc(k
1 +1,j1 )
1 +1,j1 )
2 +1,j2 )
+
+
+
Mc(k
− πc(k
+ πc(k
≤C
2 ,j2 )
2 ,j2 )
1 ,j1 )
−
−
−
Mc(k
− δc(k
+ δc(k
≤C
2 +1,j2 )
2 +1,j2 )
1 +1,j1 )
+
+
+
Mc(k
− πc(k
+ πc(k
≤ C.
1 ,j1 )
1 ,j1 )
2 ,j2 )
(22)
Hence, the computational complexity of feasibility checking
for capacity constraint is also reduced to O(1).
During the process of the neighborhood operators, if a
removed customer cannot be reinserted into any existing route,
7
then a new route is created for this customer. Since the triangle
inequality may not hold for the travel time, once a customer is
removed from a route, feasibility checking for time constraints
should also be carried out on this route. More specifically, for
customers i, j, and k, tik ≥ tij + tjk may happen; a feasibility
checking for time constraints [see (16) and (17)] should be
implemented on this route once customer j is removed. If the
route is infeasible, then it is split into two or more feasible
routes.
5) Complexity Analysis: The complexity of local search for
f1 consists of reducing vehicles and updating the archive. The
worst case for reducing vehicles is to reinsert all customers,
which takes O(N 2 ). Thus, the total complexity of local search
procedure is max(O(N 2 ), O(m · |A|)) for f1 . The main complexity of the local search procedures for f2 , . . . , f5 consists of
generating a neighbor solution and updating the archive. The
complexity of generating a neighbor solution is O(N 2 ), and the
complexity of updating the archive is O(m · |A|), where |A|
is the size of the archive. Hence, the total complexity of local
search procedures is max(O(depth · N 2 ), O(depth · m · |A|))
for f2 , . . . , f5 .
D. MOMA for MO-VRPSDPTW
The MOEA/D framework [31] is chosen for solving the MO-VRPSDPTW defined by (18). The reason is
twofold [32], [40].
1) The research literature shows that MOEA/D seems to
be more suitable for tackling multiobjective combinatorial optimization problems because it can directly use
problem-specific (local search) techniques to intensify
the exploration of promising regions in solution space.
MOEA/D provides a very natural framework for using
single-objective search techniques.
2) Pareto dominance-based algorithms, for example,
NSGA-II, do not always work well on MOPs with
many objectives [21]. Decomposition (or scalarizing
function)-based fitness evaluation is a promising alternative to Pareto dominance-based fitness evaluation
especially for the many-objective problems and multiobjective combinatorial optimization problems [32], [40].
Instead of solving a MOP directly, MOEA/D explicitly decomposes it into Q scalar optimization subproblems.
MOEA/D solves these subproblems simultaneously by evolving a population of solutions. At each generation, the population is composed of the best solution found so far for
each subproblem. The neighborhood relations among these
subproblems are defined based on the distances among their
aggregation weight vectors. Each subproblem is optimized by
using information from its neighboring subproblems. In this
paper, MOEA/D has been adapted, and thus a MOMA is proposed to solve MO-VRPSDPTW. The procedure of MOMA
is given in Algorithm 4.
Population initialization is the same as the initialization
process of MOLS but without local search. At the same
time, the same -dominance archive as in MOLS is adopted.
The distinct components of the proposed MOMA include the
decomposition, crossover operator and local search, which are
described in detail as follows.
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
8
IEEE TRANSACTIONS ON CYBERNETICS
Algorithm 4 MOMA
1: Archive A = ∅
2: generate Q uniformly distributed weight vectors 1 , · · · , Q , where
i = (λi1 , · · · , λi5 ) /*Decomposition*/
3: for i = 1 to Q do
4:
compute the Euclidean distance between each pair of weight vectors
and get the T closest weight vectors to each weight vector. Set the
neighborhood B(i) = i1 , . . . , iT .
5:
initiate xi
6: end for
7: while stopping criterion is not met do
8:
for i = 1 to Q do
9:
choose p, q randomly from B(i)
10:
o = crossover(xp , xq ) /*Crossover Operator*/
11:
if ∃obj ∈ {1, . . . , 5}, λiobj == 1 then
12:
x = LSobj (o) and update archive A /*Objectivewise local
searches: Algorithms 2 and 3*/
13:
else
14:
x = LSi (o) and update archive A /*Decomposition-based local
search: Algorithm 5*/
15:
end if
16:
for each j ∈ B(i) do
17:
if gws (x |j ) ≤ gws (xj |j ) then
18:
xj = x
19:
end if
20:
end for
21:
end for
22: end while
23: return A
Algorithm 5 Local Search LSi (x)/*Decomposition-Based
Local Search*/
1: depth = 1
2: while depth < MaxDepth do
3:
randomly se lect a neighborhood operator from (N1 , N2 , N3 )
4:
x = generate a neighbor solution of x using the selected neighborhood
operator
5:
update A with x
6:
if gws (x |i ) < gws (x|i ) then
7:
x = x
8:
end if
9:
depth = depth + 1
10: end while
11: return x
1) Decomposition:
MOMA
decomposes
the
MO-VRPSDPTW defined by (18) into Q single-objective
subproblems using a weighted sum approach, thus the ith
subproblem is defined as
min g (x | ) =
ws
i
i
5
λik fk
i
x
Fig. 2.
(23)
k=1
where i = (λi1 , λi2 , . . . , λi5 ) is the weight vector of sub
5
i
ws i i
problem i(i = 1, . . . , Q), and
k=1 λk = 1. g (x | ) is
used to emphasize that i is a coefficient vector in (23)
while xi is a solution to be optimized. According to the
definitions in Section II-A, five objectives are of different
scales, and direct use of them will make the algorithm bias
toward objectives with larger scales. Therefore, normalization
is required. The maximal and minimal values of all feasible
solutions found in each objective are used to normalize each
objective [21], [40].
2) Crossover Operator: MOMA generates new solutions
using a crossover operator described as follows. First, a random number of routes are selected from the first parent, and
copied into the offspring. Then all the routes from the second
parent, which are not in conflict with customers already copied
from the first parent, are copied into the offspring. Thus, the
crossover operator makes offspring inherit routes from parents.
Once inherited routes are chosen, they can be regarded as seed
routes. All un-routed customers should be inserted into existing routes. If one customer cannot be inserted properly, a new
route should be created for this customer. The criterion of the
insertion is according to weighted sum of objectives described
in (23). An example has been given in Fig. 2 to illustrate the
crossover operator for two parents.
3) Local Search: The offspring produced by crossover is
further improved by local search described in Algorithm 5.
Example of crossover operator.
The acceptance rule of local search for multiobjective optimization is very important because it decides which solutions
generated in the path of local search should be accepted as
the improved ones. The acceptance rule in LSi (x) is to
choose the solution with the lower value of the aggregation
objective function [see (23)] in the path of local search as
the improved solution. This local search can be viewed as
decomposition-based local search.
The subproblem with λiobj = 1 corresponds to an extreme
solution, which focuses on optimizing objective obj. In this
case, objectivewise local search is called to optimize the
corresponding objective.
E. Characteristics of Proposed Algorithms
MOLS uses different objectivewise local search procedures
to optimize different objectives, which means more problemspecific knowledge can be used to guide the search. On the
other hand, the -dominance archive is adopted to maintain the
nondominated solutions with convergence and diversity properties. Although MOEA/D provides a very natural framework
for using single-objective search techniques, it is crucial to
know when and how a single-objective search technique can be
used [40]. The proposed MOMA combines the advanced features from both objectivewise and decomposition-based local
searches with MOEA/D, and proposes a hybrid MOEA/D for
MO-VRPSDPTW.
Both MOLS and MOMA are quite simple to implement,
and they can be applied to real-world problems effectively.
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
WANG et al.: MO-VRPSDPTW: FORMULATION, INSTANCES, AND ALGORITHMS
IV. E XPERIMENTAL R ESULTS
In order to assess the performance of the proposed algorithms, experiments are implemented in C on a PC (Pentium
2.70 GHz with 2 GB RAM). The parameter MaxDepth of
MOLS is set to 10. For parameter , we find that = 0.05
can obtain satisfactory solutions on all instances. More information about the effect of can be found in [44]. In MOMA,
T is set to 5. The settings of Q and 1 , . . . , Q are controlled
by a parameter H [31]. More specifically, 1 , . . . , Q are all
weight vectors in which each individual weight takes a value
from {0, 1/H, . . . , H/H}. Therefore, the number of such vecm−1
5−1
4
= CH+5−1
= CH+4
. For all instances,
tors is Q = CH+m−1
H = 8, and therefore Q = 495. More information about the
effect of these parameters can be found in [31]. The parameters of local searches and -dominance archive used in MOMA
are set as in MOLS.
MOMA evolves a population of 495 individuals for
200 generations. Average running time (seconds) of MOMA
over 30 runs for each instance is shown in the last column (T)
of Tables I and II. Unlike the population-based evolutionary algorithm, there is no explicit concept of generation in
MOLS. Thus, the running time of MOMA is set as the
stopping criterion for MOLS in our experiments for fair
comparison.
9
TABLE I
AVERAGE VALUES OF IGD, HV, AND C-M ETRIC OF MOLS
(D ENOTED AS L) AND MOMA (D ENOTED AS M) ON
O UR R EAL -W ORLD I NSTANCES
A. Benchmark Instances
The main aim of this paper is to propose new algorithms
for real-world MO-VRPSDPTW instances. In order to test the
stability of two algorithms across datasets, the algorithms are
also tested on the previous Solomon dataset-based instances.
Hence, the proposed algorithms are tested on two sets of
benchmark instances described as follows.
1) Our Real-World MO-VRPSDPTW Instances [29]:
Forty-five real-world instances [20] are created using the combinations of three sizes of customers, three types of vehicle
capacities and five time window profiles.
The opening time of the depot is 8 h. The time windows
are designed to imitate what the delivery company faces every
day and can be described by five kinds of profiles, denoted as
profiles 0–4.
In profile 0, all customers are available all day in 480 min.
In profile 1, three types of customers are considered: early
customers, midday customers, and late customers. In order
to cover the whole day with these customers, time windows are created with a length of 160 min/type of customer
(i.e., 480 min/three types of customers). Therefore, early customers, midday customers, and late customers will be served
in the time window [0, 160], [160, 320], and [320, 480] min,
respectively.
In profile 2, the length of each time window is set to
130 min: the opening hours will be [0, 130] for early customers, [175, 305] for midday customers, and [350, 480] for
late customers.
In profile 3, the length of each time window is set to
100 min. Therefore the time windows will be [0, 100] for early
customers, [190, 290] for midday customers, and [380, 480]
for late customers.
In profile 4, customers are associated with one of the time
windows from profiles aforementioned.
The capacity of each vehicle can be set to C = D +
δ/100(D − D), where δ ∈ [0, 100] is used to modulate
the slack
margin of the capacity [20], D = maxi {gi } and
D= N
i=1 gi . If δ takes value close to 0, the capacity of the
vehicle C will be very limited. If δ takes value close to 100,
the vehicle will have a capacity C close to the total demand.
MO-VRPSDPTW instances are created by using
the following combinations: 1) number of customers:
{50, 150, 250}; 2) three types of δ: {60, 20, 5}; and 3)
time windows: profiles {0, 1, 2, 3, 4}. Hence, a total of
45 MO-VRPSDPTW instances (3 ∗ 3 ∗ 5) are generated. As shown in the first column of Table I, the
format of the instance name is “num1 -num2 -num3 ,”
where num1 indicates the number of customers, num2
indicates the index of δ type, and num3 indicates the index
of time window profile. The demand for each customer is
set to 10, 20, or 30, each with probability 1/3, the pickup
quantity for each customer is also set to 10, 20, or 30, each
with probability 1/3, and the service time for each customer
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
10
IEEE TRANSACTIONS ON CYBERNETICS
TABLE II
AVERAGE VALUES OF IGD, HV, AND C-M ETRIC OF MOLS
(D ENOTED AS L) AND MOMA (D ENOTED AS M) ON
R EVISED S OLOMON I NSTANCES
no unary performance metric is able to give a comprehensive
measure on the performance of a multiobjective optimization
algorithm. In this paper, the following three metrics are used.
1) Inverted generational distance (IGD) [31] is used to
measure both the convergence and diversity of the
nondominated solutions obtained.
2) Hypervolume (HV) [50] is used to indicate the area in
the objective space that is dominated by at least one
of the nondominated solutions. The larger the HV, the
closer the corresponding nondominated solutions toward
the Pareto front. Therefore, it reflects the closeness of
the nondominated solutions to the Pareto front. HV also
measures both the convergence and diversity.
3) Coverage metric (C-metric) [51] is a commonly used
metric for comparing two sets of nondominated solutions (denoted as X and Y). C(X, Y) is defined as the
percentage of the solutions in Y that are Pareto dominated by at least one solution in X. C(X, Y) = 1 means
all nondominated solutions in Y are Pareto dominated
by solutions in X, and C(Y, X) = 1 means all nondominated solutions in X are Pareto dominated by solutions
in Y. The sum of C(X, Y) and C(Y, X) is not always
equal to 1, since some solutions in X and Y may not
Pareto dominate each other.
Higher (lower) value of HV (IGD) can be considered as
a better set of solutions approximating the true Pareto front
from the convergence and diversity viewpoints. To provide
additional information on convergence, C-metric is also used.
To compute the performance metrics, the objective vectors of
all nondominated solutions obtained by the considered algorithms are approximately seen as the Pareto front because the
true Pareto front of the problems is unknown. Furthermore,
all objective values are normalized because of the difference
among the ranges of objectives.
C. Experimental Results
is set to 10, 20, or 30 min, also each with probability 1/3.
A maximum delay of 30 min is allowed for each customer,
that is, md = 30 min [20].
2) Revised Solomon Instances by Wang and Chen [5]: They
are revised from Solomon instances, and can be downloaded
in [13]. This dataset includes 56 instances with 100 customers.
In this paper, the maximum delay allowed for each customer
i is set to 30% of the length of its time windows, that is,
mdi = 30% · (ei − bi ) as in [24], since the instances within
each problem class have different time window lengths.
B. Performance Metrics
The performance of an algorithm for MOPs is evaluated in
terms of both convergence and diversity. As discussed in [49],
A significant feature of the real-world MO-VRPSDPTW
dataset is that the triangle inequality may not hold
for the travel time. Thus, previous single-objective or
multiobjective optimization approaches for the Solomon
dataset [4], [5], [52], which assume that the triangle inequality always holds for the travel time, cannot be directly used
to solve real-world MO-VRPSDPTW instances. Furthermore,
most previous algorithms only consider VRPSDPTW/VRPTW
with two objectives [5], [52]; thus, they cannot effectively
solve the many-objective VRPSDPTW considered in this
paper. Two algorithms proposed in this paper are compared
with each other since there is no previous benchmark algorithm. Our algorithms can be seen as benchmark algorithms for
real-world MO-VRPSDPTW instances, and can be compared
by future research.
Average values of the metrics over 30 independent runs
of both algorithms on the test sets have been calculated, and
are shown in Tables I and II. Due to space limitations, standard deviations of the metrics are not presented in the tables.
To further show the difference between two algorithms, the
Wilcoxon signed-rank test [53]–[55] at 5% significance level
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
WANG et al.: MO-VRPSDPTW: FORMULATION, INSTANCES, AND ALGORITHMS
11
Fig. 3. Plots of nondominated solutions obtained by MOLS (denoted as ◦) and MOMA (denoted as ×) on the selected instance 250–0–0, respectively.
(a) f 1 − f2 plane. (b) f1 − f3 plane. (c) f2 − f3 plane.
is conducted. The significantly better values between MOLS
and MOMA are highlighted in boldface in tables. The result
of the test is summarized as w/t/l, which means that the performance of MOLS is better than, similar to, and worse than
that of MOMA in w, t, and l instances, respectively. From
Tables I and II, several observations can be found as follows.
1) For 45 real-world instances, MOLS significantly outperforms MOMA in most instances. Specifically, MOLS
significantly outperforms MOMA in 28 instances and is
outperformed by MOMA in ten instances in terms of
IGD. In terms of HV, MOLS significantly outperforms
MOMA in 26 instances and is outperformed by MOMA
in 12 instances. This means that MOLS can obtain better
results than MOMA in terms of both convergence and
diversity. MOLS also significantly outperforms MOMA
in 38 instances in terms of C-metric. It means that
MOLS shows better convergence performance.
2) For 56 revised Solomon instances, it is more obvious
that MOLS significantly outperforms MOMA in terms
of all metrics. Specifically, MOLS significantly outperforms MOMA in 46 instances and is outperformed
by MOMA in only one instance in terms of IGD. In
terms of HV, MOLS significantly outperforms MOMA
in 47 instances and is outperformed by MOMA in six
instances. MOLS significantly outperforms MOMA in
55 instances in terms of C-metric.
3) As shown above, for revised Solomon instances,
MOLS significantly outperforms MOMA in almost all
instances, while MOLS is significantly outperformed by
MOMA in some real-world instances. The superiority
of MOLS over MOMA in real-world instances is not so
obvious as in revised Solomon instances. Thus, future
researchers should test their algorithms on the proposed real-world MO-VRPSDPTW instances to prove
usefulness in real life environment for their algorithms.
To visually demonstrate the convergence and diversity properties of two algorithms, the projection of nondominated
solutions obtained by MOLS (denoted as ◦) and MOMA
(denoted as ×) in a selected instance 250–0–0, at f1 − f2 ,
f1 − f3 , and f2 − f3 planes are shown in Fig. 3. It is clearly
shown that most objective vectors generated by MOLS are
better than those generated by MOMA. Moreover, the distribution of solutions obtained by MOLS spreads much wider along
the obtained Pareto front. Hence, the fact that MOLS is better
than MOMA in this instance in terms of both convergence and
diversity is confirmed.
Finally, the relative benefits and limitations of the proposed algorithms are briefly discussed. In both MOLS and
MOMA, single-objective local search can be easily adopted
for MO-VRPSDPTW. In MOMA, diversity is naturally preserved by the diversity among subproblems, but some different
subproblems may have the same optimal solution for a combinatorial optimization problem [31]. Therefore, even a large
number of subproblems may not lead to a reasonably good
approximation to the Pareto front of the combinatorial optimization problem. In MOLS, no subproblem (weight) needs
to be defined in advance. Instead, objectivewise local search is
adopted to optimize each objective, which promotes the convergence and diversity of search. This may be the reason that
MOLS is better than MOMA in terms of both convergence
and diversity in most of instances.
V. C ONCLUSION
This paper has introduced a multiobjective variant of
VRPSDPTW and a set of realistic benchmark instances. Then
two algorithms have been designed for the MO-VRPSDPTW.
Extensive experiments have shown the effectiveness of the
proposed algorithms. The proposed algorithms can be seen
as benchmark algorithms for real-world MO-VRPSDPTW
instances, which can be used for comparison by future
research.
In the future, this paper can be extended in multiple
directions. Firstly, the proposed MO-VRPSDPTW model
can be extended to other green VRPs, for example, pollution VRP for reducing CO2 or greenhouse gas emissions,
by including broader objectives that reflect environmental
cost [2], [3]. Secondly, the proposed algorithms can also be
extended to solve other multiobjective VRPs in reverse logistics, including VRP with backhauls and VRP with mixed
pickup and delivery problems. Finally, from the perspective
of fundamental research, advanced multiobjective optimization algorithms for many-objective combinatorial optimization
problems should be further studied and developed since existing multiobjective optimization algorithms mainly focus on
continuous benchmark functions.
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
12
IEEE TRANSACTIONS ON CYBERNETICS
R EFERENCES
[1] C. Lin, K. Choy, G. Ho, S. Chung, and H. Lam, “Survey of green
vehicle routing problem: Past and future trends,” Expert Syst. Appl.,
vol. 41, no. 4, pp. 1118–1138, 2014.
[2] K. Devika, A. Jafarian, and V. Nourbakhsh, “Designing a sustainable closed-loop supply chain network based on triple bottom line
approach: A comparison of metaheuristics hybridization techniques,”
Eur. J. Oper. Res., vol. 235, no. 3, pp. 594–615, 2014.
[3] T. R. P. Ramos, M. I. Gomes, and A. P. Barbosa-Povoa, “Planning a
sustainable reverse logistics system: Balancing costs with environmental
and social concerns,” Omega, vol. 48, pp. 60–74, Oct. 2014.
[4] S. Kassem and M. Chen, “Solving reverse logistics vehicle routing problems with time windows,” Int. J. Adv. Manuf. Technol., vol. 68, nos. 1–4,
pp. 57–68, 2013.
[5] H.-F. Wang and Y.-Y. Chen, “A genetic algorithm for the simultaneous
delivery and pickup problems with time window,” Comput. Ind. Eng.,
vol. 62, no. 1, pp. 84–95, Feb. 2012.
[6] M. Lai and E. Cao, “An improved differential evolution algorithm for
vehicle routing problem with simultaneous pickups and deliveries and
time windows,” Eng. Appl. Artif. Intell., vol. 23, no. 2, pp. 188–195,
2010.
[7] B. Eksioglu, A. V. Vural, and A. Reisman, “The vehicle routing problem:
A taxonomic review,” Comput. Ind. Eng., vol. 57, no. 4, pp. 1472–1483,
2009.
[8] S.-C. Horng, “Combining artificial bee colony with ordinal optimization for stochastic economic lot scheduling problem,” IEEE Trans. Syst.,
Man, Cybern., Syst., vol. 45, no. 3, pp. 373–384, Mar. 2015.
[9] Y. Jin and J.-K. Hao, “Effective learning-based hybrid search for
bandwidth coloring,” IEEE Trans. Syst., Man, Cybern., Syst., to be
published.
[10] D. Li, M. Li, X. Meng, and Y. Tian, “A hyperheuristic approach for
intercell scheduling with single processing machines and batch processing machines,” IEEE Trans. Syst., Man, Cybern., Syst., vol. 45, no. 2,
pp. 315–325, Feb. 2015.
[11] E. Angelelli and R. Mansini, “The vehicle routing problem with
time windows and simultaneous pick-up and delivery,” in Quantitative
Approaches to Distribution Logistics and Supply Chain Management,
A. Klose, M. Speranza, and L. V. Wassenhove, Eds. Berlin, Germany:
Springer, 2002, pp. 249–267.
[12] L. Boubahri, S.-A. Addouche, and A. E. Mhamedi, “Multi-ant colonies
algorithms for the VRPSPDTW,” in Proc. Int. Conf. Commun. Comput.
Control Appl., Hammamet, Tunisia, 2011, pp. 1–6.
[13] H.-F. Wang and Y.-Y. Chen. (2014). VRPSDPTW Instances
Revised
From
Solomon
Dataset.
[Online].
Available:
http://oz.nthu.edu.tw/~d933810/test.htm
[14] M. M. Solomon, “Algorithms for the vehicle routing and scheduling
problems with time window constraints,” Oper. Res., vol. 35, no. 2,
pp. 254–265, 1987.
[15] C. Wang, F. Zhao, D. Mu, and J. W. Sutherland, “Simulated annealing for
a vehicle routing problem with simultaneous pickup-delivery and time
windows,” in Advances in Production Management Systems. Sustainable
Production and Service Supply Chains. Berlin, Germany: Springer, 2013,
pp. 170–177.
[16] A. Deng, C. Mao, and Y. Zhou, “Optimizing research of an improved
simulated annealing algorithm to soft time windows vehicle routing
problem with pick-up and delivery,” Syst. Eng. Theory Pract., vol. 29,
no. 5, pp. 186–192, 2009.
[17] R. Liu, X. Xie, V. Augusto, and C. Rodriguez, “Heuristic algorithms for
a vehicle routing problem with simultaneous delivery and pickup and
time windows in home health care,” Eur. J. Oper. Res., vol. 230, no. 3,
pp. 475–486, Nov. 2013.
[18] N. Jozefowiez, F. Semet, and E.-G. Talbi, “Multi-objective vehicle
routing problems,” Eur. J. Oper. Res., vol. 189, no. 2, pp. 293–309,
2008.
[19] N. Labadie and C. Prodhon, “A survey on multi-criteria analysis in
logistics: Focus on vehicle routing problems,” in Applications of MultiCriteria and Game Theory Approaches, L. Benyoucef, J.-C. Hennet, and
M. K. Tiwari, Eds. London, U.K.: Springer, 2014, pp. 3–29.
[20] J. Castro-Gutierrez, D. Landa-Silva, and J. M. P´erez, “Nature of realworld multi-objective vehicle routing with evolutionary algorithms,”
in Proc. IEEE Int. Conf. Syst. Man Cybern. (SMC), Anchorage, AK,
USA, 2011, pp. 257–264.
[21] Y. Zhou and J. Wang, “A local search-based multiobjective optimization algorithm for multiobjective vehicle routing problem with time
windows,” IEEE Syst. J., to be published.
[22] D. Tas, N. Dellaert, T. van Woensel, and T. de Kok, “Vehicle routing
problem with stochastic travel times including soft time windows and
service costs,” Comput. Oper. Res., vol. 40, no. 1, pp. 214–224, 2013.
[23] H. Hashimotoa, T. Ibaraki, S. Imahori, and M. Yagiura, “The vehicle routing problem with flexible time windows and traveling times,”
Discrete Appl. Math., vol. 154, no. 16, pp. 2271–2290, 2006.
[24] R. Eglese, Z. Fu, and L. Y. O. Li, “A unified tabu search algorithm for
vehicle routing problems with soft time windows,” J. Oper. Res. Soc.,
vol. 59, pp. 663–673, Feb. 2008.
[25] W. C. Chiang and R. A. Russell, “A metaheuristic for the vehiclerouting problem with soft time windows,” J. Oper. Res. Soc., vol. 55,
pp. 1298–1310, Jul. 2004.
[26] J. Castro-Gutierrez, “Multi-objective tools for the vehicle routing
problem with time windows,” Ph.D. dissertation, School Comput.
Sci., Univ. Nottingham, Nottingham, U.K., 2012. [Online]. Available:
http://etheses.nottingham.ac.uk/3713/
[27] A. Rodriguez and R. Ruiz, “A study on the effect of asymmetry on
real capacitated vehicle routing problems,” Comput. Oper. Res., vol. 39,
no. 7, pp. 2142–2151, 2012.
[28] C. L. Fleming, S. E. Griffis, and J. E. Bell, “The effects of triangle
inequality on the vehicle routing problem,” Eur. J. Oper. Res., vol. 224,
no. 1, pp. 1–7, 2013.
[29] J. Wang. (2014). MO-VRPTW and MO-VRPSDPTW Datasets. [Online].
Available: http://sist.sysu.edu.cn/~wangjiah/
[30] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist
multiobjective genetic algorithm: NSGA-II,” IEEE Trans. Evol. Comput.,
vol. 6, no. 2, pp. 182–197, Apr. 2002.
[31] Q. Zhang and H. Li, “MOEA/D: A multiobjective evolutionary algorithm
based on decomposition,” IEEE Trans. Evol. Comput., vol. 11, no. 6,
pp. 712–731, Dec. 2007.
[32] A. Zhou et al., “Multiobjective evolutionary algorithm: A survey of the
state of the art,” Swarm Evol. Comput., vol. 1, no. 1, pp. 32–49, 2011.
[33] K. Nag, T. Pal, and N. Pal, “ASMiGA: An archive-based steadystate micro genetic algorithm,” IEEE Trans. Cybern., vol. 45, no. 1,
pp. 40–52, Jan. 2015.
[34] M. Li, S. Yang, K. Li, and X. Liu, “Evolutionary algorithms with
segment-based search for multiobjective optimization problems,” IEEE
Trans. Cybern., vol. 44, no. 8, pp. 1295–1313, Aug. 2014.
[35] Z.-H. Zhan et al., “Multiple populations for multiple objectives:
A coevolutionary technique for solving multiobjective optimization
problems,” IEEE Trans. Cybern., vol. 43, no. 2, pp. 445–463, Apr. 2013.
[36] L. Paquete, T. Schiavinotto, and T. St¨utzle, “On local optima in multiobjective combinatorial optimization problems,” Ann. Oper. Res., vol. 156,
no. 1, pp. 83–97, 2007.
[37] F. Tricoire, “Multi-directional local search,” Comput. Oper. Res., vol. 39,
no. 12, pp. 3089–3101, 2012.
[38] E.-G. Talbi, M. Basseur, A. Nebro, and E. Alba, “Multi-objective optimization using metaheuristics: Non-standard algorithms,” Int. Trans.
Oper. Res., vol. 19, nos. 1–2, pp. 283–305, 2012.
[39] L. Ke, Q. Zhang, and R. Battiti, “Hybridization of decomposition and
local search for multiobjective optimization,” IEEE Trans. Cybern.,
vol. 44, no. 10, pp. 1808–1820, Oct. 2014.
[40] J. Wang and Y. Cai, “Multiobjective evolutionary algorithm for frequency assignment problem in satellite communications,” Soft Comput.,
to be published.
[41] J. Wang, C. Zhong, Y. Zhou, and Y. Zhou, “Multiobjective optimization algorithm with objective-wise learning for continuous multiobjective
problems,” J. Amb. Intell. Human. Comput., to be published.
[42] M. Laumanns, L. Thiele, K. Deb, and E. Zitzler, “Combining convergence and diversity in evolutionary multiobjective optimization,”
Evol. Comput., vol. 10, no. 3, pp. 263–282, 2002.
[43] K. Deb, M. Mohan, and S. Mishra, “Evaluating the -domination
based multi-objective evolutionary algorithm for a quick computation of
Pareto-optimal solutions,” Evol. Comput., vol. 13, no. 4, pp. 501–525,
2005.
[44] S. Bandyopadhyay, U. Maulik, and R. Chakraborty, “Incorporating
-dominance in AMOSA: Application to multiobjective 0/1 knapsack
problem and clustering gene expression data,” Appl. Soft Comput.,
vol. 13, no. 5, pp. 2405–2411, 2013.
[45] P. Hansen and N. Mladenovic, “Variable neighborhood search: Principles
and applications,” Eur. J. Oper. Res., vol. 130, no. 3, pp. 449–467, 2001.
[46] G. A. P. Kindervater and M. W. P. Savelsbergh, “Vehicle routing: Handling edge exchanges,” in Local Search in Combinatorial
Optimization, E. H. L. Aarts and J. K. Lenstra, Eds. Chichester, U.K.:
Wiley, 1997, pp. 337–360.
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
WANG et al.: MO-VRPSDPTW: FORMULATION, INSTANCES, AND ALGORITHMS
[47] N. Labadie and C. Prins, “Vehicle routing nowadays: Compact review
and emerging problems,” in Production Systems and Supply Chain
Management in Emerging Countries: Best Practices, G. Mej´ıa and
N. Velasco, Eds. Berlin, Germany: Springer, 2012, pp. 141–166.
[48] N. Bianchessi and G. Righini, “Heuristic algorithms for the vehicle routing problem with simultaneous pick-up and delivery,” Comput. Oper.
Res., vol. 34, no. 2, pp. 578–594, 2007.
[49] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and
V. G. da Fonseca, “Performance assessment of multiobjective
optimizers: An analysis and review,” IEEE Trans. Evol. Comput., vol. 7,
no. 2, pp. 117–132, Apr. 2003.
[50] M. Laumanns, E. Zitzler, and L. Thiele, “SPEA2: Improving the strength
Pareto evolutionary algorithm,” in Evolutionary Methods for Design,
Optimization and Control With Applications to Industrial Problems,
K. Giannakoglou, D. T. Tsahalis, J. Periaux, K. D. Papailiou, and
T. Fogarty, Eds. Barcelona, Spain: CIMNE, 2002, pp. 95–100.
[51] E. Zitzler and L. Thiele, “Multi-objective evolutionary algorithms:
A comparative study and the strength Pareto approach,” IEEE Trans.
Evol. Comput., vol. 3, no. 4, pp. 257–271, Nov. 1999.
[52] T.-C. Chiang and W.-H. Hsu, “A knowledge-based evolutionary algorithm for the multiobjective vehicle routing problem with time
windows,” Comput. Oper. Res., vol. 45, no. 5, pp. 25–37, 2014.
[53] F. Wilcoxon, “Individual comparisons by ranking methods,” Biometrics,
vol. 1, no. 6, pp. 80–83, 1945.
[54] J. Derrac, S. García, D. Molina, and F. Herrera, “A practical tutorial on
the use of nonparametric statistical tests as a methodology for comparing
evolutionary and swarm intelligence algorithms,” Swarm Evol. Comput.,
vol. 1, no. 1, pp. 3–18, 2011.
[55] J. Alcalá-Fdez et al., “KEEL: A software tool to assess evolutionary
algorithms to data mining problems,” Soft Comput., vol. 13, no. 3,
pp. 307–318, 2008.
Jiahai Wang (M’07) received the Ph.D. degree from
Toyama University, Toyama, Japan, in 2005.
In 2005, he joined Sun Yat-sen University,
Guangzhou, China, where he is currently an
Associate Professor with the Department of
Computer Science. His current research interests include computational intelligence and its
applications.
Ying Zhou received the Ph.D. degree from Sun
Yat-sen University, Guangzhou, China, in 2014.
In 2014, she joined the Shenzhen Institute of
Information Technology, Shenzhen, China. Her
current research interests include local search
algorithms and their applications, multiobjective
optimization, and other evolutionary computation
techniques.
Yong Wang (M’08) received the B.S. degree in
automation from the Wuhan Institute of Technology,
Wuhan, China, in 2003, and the M.S. degree in pattern recognition and intelligent systems and Ph.D.
degree in control science and engineering, both from
Central South University (CSU), Changsha, China,
in 2006 and 2011, respectively.
He is currently an Associate Professor with the
School of Information Science and Engineering,
CSU. His current research interests include evolutionary computation, single-objective optimization,
constrained optimization, multiobjective optimization, and their real-world
applications. He was a Reviewer for over 40 international journals.
Dr. Wang was a recipient of the Hong Kong Scholar from the
Mainland—Hong Kong Joint Post-Doctoral Fellows Program, China, in 2013,
the Excellent Doctoral Dissertation by Hunan Province, China, in 2013, the
New Century Excellent Talents in University by the Ministry of Education,
China, in 2013, and the 2015 IEEE Computational Intelligence Society
Outstanding Ph.D. Dissertation Award. He was a PC Member of over
20 international conferences.
13
Jun Zhang (M’02–SM’08) received the Ph.D.
degree in electrical engineering from the City
University of Hong Kong, Hong Kong, in 2002.
From 2003 to 2004, he was a Brain Korean
21 Post-Doctoral Fellow at the Department of
Electrical Engineering and Computer Science, Korea
Advanced Institute of Science and Technology,
Daejeon, Korea. Since 2004, he has been at Sun
Yat-sen University, Guangzhou, China, where he
is currently a Cheung Kong Professor. His current
research interests include computational intelligence,
cloud computing, big data, high performance computing, data mining, wireless sensor networks, operations research, and power electronic circuits. He
has authored seven research books and book chapters and over 100 technical
papers in the above areas.
Mr. Zhang was a recipient of the China National Funds for Distinguished
Young Scientists from the National Natural Science Foundation of China,
in 2011 and the First-Grade Award in Natural Science Research from
the Ministry of Education, China, in 2009. He is currently an Associate
Editor of the IEEE T RANSACTIONS ON E VOLUTIONARY C OMPUTATION,
the IEEE T RANSACTIONS ON I NDUSTRIAL E LECTRONICS, and the IEEE
T RANSACTIONS ON C YBERNETICS. He is the Founding and the Current Chair
of the IEEE Guangzhou Subsection and ACM Guangzhou Chapter.
C. L. Philip Chen (S’88–M’88–SM’94–F’07)
received the M.S. degree from the University of
Michigan, Ann Arbor, MI, USA, and the Ph.D.
degree from Purdue University, West Lafayette, IN,
USA, in 1985 and 1988, respectively, both in electrical engineering.
He was a tenured professor, a Department Head,
and an Associate Dean with two different universities in the U.S. for 23 years. He is currently the
Dean of the Faculty of Science and Technology
and a Chair Professor with the Department of
Computer and Information Science, University of Macau, Macau, China.
His current research interests include systems, cybernetics, and computational
intelligence.
Prof. Chen was the President of the IEEE Systems, Man, and Cybernetics
Society from 2012 to 2013. He has been an Editor-in-Chief of the IEEE
T RANSACTIONS ON S YSTEMS , M AN , AND C YBERNETICS : S YSTEMS since
2014 and an associate editor of several IEEE transactions. He is the
Chair of Technical Committee 9.1 Economic and Business Systems of
International Federation of Automatic Control. He is a Program Evaluator for
the Accreditation Board of Engineering and Technology Education in computer engineering, electrical engineering, and software engineering programs,
USA. He is a fellow of the American Association for the Advancement of
Science.
Zibin Zheng (S’05–M’11) received the Ph.D.
degree from the Department of Computer Science
and Engineering, Chinese University of Hong Kong
(CUHK), Hong Kong, in 2010.
He is an Associate Research Fellow with the
Shenzhen Research Institute, Chinese University of
Hong Kong. His current research interests include
cloud computing, service computing, and software
engineering.
Mr. Zheng was a recipient of the Outstanding
Thesis Award of CUHK in 2012, the ACM Special
Interest Group on Software Engineering Distinguished Paper Award in
International Conference on Software Engineering 2010, the IBM Ph.D.
Fellowship Award in 2010. He served as a Program Committee Member of
more than 20 conferences.