Information Sciences 406–407 (2017) 146–169

Contents lists available at ScienceDirect

Information Sciences
journal homepage: www.elsevier.com/locate/ins

Multi-vehicle selective pickup and delivery using
metaheuristic algorithms
Chuan-Kang Ting∗, Xin-Lan Liao, Yu-Hsuan Huang, Rung-Tzuo Liaw
Department of Computer Science and Information Engineering, National Chung Cheng University, Chiayi 62102, Taiwan

a r t i c l e

i n f oArticle history:
Received 17 December 2015
Revised 21 January 2017
Accepted 3 April 2017
Available online 8 April 2017
Keywords:
Selective pickup and delivery problem
Multiple vehicles
Tabu search
Genetic algorithm
Scatter searcha b s t r a c t
The pickup and delivery problem (PDP) addresses real-world problems in logistics andtransportation, and establishes a critical class of vehicle routing problems. This study
presents a novel variant of the PDP, called the multi-vehicle selective pickup and delivery problem (MVSPDP), and designs three metaheuristic algorithms for this problem. The
MVSPDP aims to ﬁnd the minimum-cost routes for a ﬂeet of vehicles collecting and supplying commodities, subject to the constraints on vehicle capacity and travel distance. The
problem formulation features relaxing the requirement of visiting all pickup nodes and enabling multiple vehicles for achieving transportation eﬃciency. To solve the MVSPDP, we
propose three metaheuristic algorithms: tabu search (TS), genetic algorithm (GA), and scatter search (SS). A ﬁxed-length representation is presented to indicate the varying numberof vehicles used and the selection of pickup nodes. Furthermore, we devise four operators
for TS, GA, and SS to handle the selection of pickup nodes, number of vehicles used, andtheir routes. The experimental results indicate that the three metaheuristic algorithms caneffectively solve the MVSPDP. In particular, TS outperforms GA and SS in solution quality
and convergence speed. In addition, the problem formulation produces substantially lower
transportation costs than the PDP does, thus validating the utility of the MVSPDP.
© 2017 Elsevier Inc. All rights reserved.

1. Introduction
The pickup and delivery problem (PDP) belongs to a critical class of vehicle routing problems and arises in several realworld logistic problems, such as the delivery of packages and letters. In the PDP, customers are classiﬁed into pickup nodes
and delivery nodes, which represent customers providing and demanding commodities, respectively. The PDP is used to ﬁnd
the shortest route that satisﬁes the requests of customers. Comprehensive surveys of the PDP are presented in [7,50,51,61].
Parragh et al. [50,51] divided PDP scenarios into transportation between the depot and customers as well as conveyanceamong customers. The ﬁrst type involves commodities picked from and delivered to the depot. This one-to-many-to-one(1-M-1) PDP is applicable to reverse logistics in simultaneously managing product distribution from the storehouse and material collection for remanufacture [33,34,55,56]. The second type implements the one-to-one (1-1) PDP structure, throughwhich commodities are transferred between paired pickup and delivery nodes (e.g., dial-a-ride system [18] and message
transmission in mobile networks [64]). This type can also be applied in many-to-many (M-M) transportation to supply a setof delivery nodes with commodities collected from numerous pickup nodes [1,8,9,24,37].

∗

Corresponding author.
E-mail address: ckting@cs.ccu.edu.tw (C.-K. Ting).

http://dx.doi.org/10.1016/j.ins.2017.04.001
0020-0255/© 2017 Elsevier Inc. All rights reserved.

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

147

In most PDPs, equal amounts of total supply and total demand are assumed, which implicitly imposes a constraint on
visiting all customers [7]. Ting and Liao [61] recently proposed the selective pickup and delivery problem (SPDP) by relaxing
the constraint that the vehicle must visit all nodes in the PDP. Instead of collecting commodities from all pickup nodes,
the SPDP enables the selection of pickup nodes from which to collect suﬃcient commodities for all delivery nodes. Sucha relaxation in visiting all pickup nodes can substantially reduce the transportation cost for cases that focus on satisfying
the demand of delivery nodes. The selectability of pickup nodes is particularly useful for the real-world applications dealing
with the stock and supply of commodities. For example, in the city bike rental service, arranging the routes for the trucks totransport bikes to the popular spots is a key issue for redistribution of bikes. The SPDP facilitates reducing the transportation
cost by picking up bikes from only some stations, rather than visiting all, and then delivering them to the demanded areas.
Compared with the PDP, the selectability of pickup nodes in the SPDP permits unequal amounts of total supply and total
demand and enables all delivery nodes to be supplied eﬃciently with an adequate number of commodities collected fromsome pickup nodes by using a single vehicle.
This study extends the SPDP in the number of vehicles to improve transportation eﬃciency. Speciﬁcally, we propose anovel problem formulation, called the multi-vehicle selective pickup and delivery problem (MVSPDP), aiming to minimize thetransportation cost for multiple vehicles to supply all delivery nodes with commodities obtained from some selected pickup
nodes, given that each vehicle exhibits identical limits on vehicle capacity and travel distance. Using multiple vehicles caneffectually enhance transportation eﬃciency; however, this extension induces additional constraints, objectives, and logistic
characteristics and, therefore, poses considerable challenges to solving the problem. First, two additional constraints must be
addressed in the MVSPDP: one constraint requires the vehicle load to be nonnegative and lower than the vehicle capacity
along the route; the other constraint imposes the maximum travel distance for workload balance or the refueling requirements of each vehicle. Second, the number of vehicles in use plays an essential role in the resultant routes. Determining theappropriate number of vehicles to use and arranging their routes are key issues in solving the MVSPDP. Third, the selectability of pickup nodes can substantially reduce the cost of the route but increase the problem complexity [61]. An increase in
the number of vehicles, nevertheless, compounds the diﬃculty in selecting suitable pickup nodes for optimal routes.
To solve the MVSPDP, we develop three metaheuristic algorithms: tabu search (TS) [29], genetic algorithm (GA) [39],and scatter search (SS) [28]. These three algorithms have been proved to be effective in resolving numerous search andoptimization problems [15,30,46,49]. The key challenge in these algorithms is to design effective operators that select pickup
nodes, organize visiting orders, and allocate multiple vehicles to minimize the total transportation cost in accordance withthe two constraints. In this study, we design the components and operators for the three metaheuristic algorithms to solve
the MVSPDP:• The representation for candidate solutions enables indicating the varying number of vehicles used and the selection of
pickup nodes by using a ﬁxed-length string.
• An evaluation function considering the solution feasibility is proposed to handle the constraints.
• An initialization method based on the sweep algorithm is developed to improve search eﬃciency.
• Four operators are devised to deal with route planning, the selection of pickup nodes, and the number of vehicles
adopted.
In this study, a series of experiments is conducted to examine and compare the performance of the three metaheuristic
algorithms in solving the MVSPDP. In addition, the experiments are conducted to investigate the effects of deterministic
and adaptive control on the probability of performing a speciﬁc operator. The remainder of this paper is organized as follows: Section 2 presents the formulation of the MVSPDP; Section 3 elucidates the proposed TS, GA, and SS methods; andSection 4 provides the experimental results. Finally, Section 5 draws the conclusions of this study.
2. The multi-vehicle selective pickup and delivery problem
This study presents a novel variant of the PDP, called the multi-vehicle selective pickup and delivery problem (MVSPDP).
The MVSPDP aims for the shortest routes for a ﬂeet of vehicles that transport commodities from some pickup nodes to alldelivery nodes, subject to the constraints on vehicle capacity and travel distance. The vehicles are regulated to start fromand return to the depot. The following subsection presents the problem formulation of the MVSPDP. In addition, we derivea bound of solution values for the MVSPDP.
2.1. Problem formulation

 Let  G = (V, E ) be an undirected complete graph with vertex set V = (v0 , . . . , vn ) and edge set E =vi , v j | vi , v j ∈ V, vi = v j , in which each edge (vi , v j ) has a cost cij > 0 and ci j = c ji . The MVSPDP involves one depot v0 and n customer nodes v1 , . . . , vn . Let di denote the demand of node vi . According to the demand, the vertex set V is
divided into a set of pickup nodes V + = {vi | vi ∈ V, di > 0}, a set of delivery nodes V − = {vi | vi ∈ V, di < 0}, and the depot v0
with d0 = 0. All m vehicles are assumed to exhibit the same capacity Q and maximal travel distance permitted R.
Resolving the MVSPDP requires selecting pickup nodes and arranging visiting orders such that multiple vehicles supply
all delivery nodes with suﬃcient commodities. Let binary decision variable xi jk = 1 if the kth vehicle travels from vi to v j in a

148

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

Fig. 1. An example representation for an MVSPDP instance with n = 16 customer nodes and m = 4 vehicles. The number inside circles denotes demand di .
The solid, dashed, and dotted lines represent the routes of three vehicles.

directed path; otherwise, xi jk = 0. Another binary variable yik is associated with the visiting of a vehicle to a particular node:
yik = 1 if the kth vehicle visits vi , and yik = 0 otherwise. Let V ± = V + ∪ V − , the objective and constraints of the MVSPDP areformulated as follows:

min

m



ci j xi jk(1)

k=1 vi ,v j ∈Vs.t.

v j ∈V
m

k=1
m

k=1
m




xi jk =

x jik = yik , ∀vi ∈ V, k ∈ {1, . . . , m}(2)v j ∈Vyik = 1, ∀vi ∈ V −(3)yik ≤ 1, ∀vi ∈ V +(4)

y0k ≤ m

(5)

k=1



vi ,v j ∈Vysk ≤ci j xi jk ≤ R, ∀k ∈ {1, . . . , m}(6)xi jk , ∀S ⊆ V ± , vs ∈ S, k ∈ {1, . . . , m}(7)vi ∈S v j ∈/ S
x0 jk , ∀S ⊆ V ± , k ∈ {1, . . . , m}y−1vi ∈S ikv j ∈S0 ≤ zSkdi yik ≤ Q, ∀S ⊆ V ± , k ∈ {1, . . . , m}v ,v j ∈S xi jkizSk = 

(8)
(9)

vi ∈S

Constraint (2) restricts incoming and outgoing ﬂows to be equal and gives the value of decision variable yik . Constraint(3) guarantees that each delivery node can be visited exactly once, and constraint (4) enables selectability of pickup nodes.
The number of available vehicles is limited to m as speciﬁed in (5). Constraint (6) further limits the maximal travel distanceR permitted for each vehicle. The subtour elimination for each route is presented in constraint (7), which requires that eachsubset of customers assigned to a vehicle has at least one ﬂow going out of that subset. Regarding (8), the ﬂoor part considers the connection of nodes in the subset S; that is, the number of traversed edges in a route equals the number ofvisited nodes minus one, i.e., v ,v ∈S xi jk = v ∈S yik − 1. Hence, the variable zSk = 1 if there exists an incoming ﬂow issuedijifrom v0 to the subsequent connected route of the kth vehicle; otherwise, zSk = 0. For example, considering the route of thethird vehicle in Fig. 1, i.e., v0→v6→v10→v13→v16→v0 , the connected nodes S = {v6 , v10 }, {v6 , v10 , v13 }, and {v6 , v10 , v13 , v16 }have v ,v ∈S xi j3 = v ∈S yi3 − 1 = 1, 2, and 3, respectively, which give zS3 = 1 in that x063 = 1. By contrast, the unconnectediji

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

nodes in the route, e.g. S = {v10 , v16 }, result in zS3 = 0 due tovi ,v j ∈S xi j3 <vi ∈S yi3 − 1 as well as149



v j ∈S x0 j3 = 0. Accord-

ingly, constraint (9) conﬁnes the vehicle load within capacity Q.
The MVSPDP relaxes the requirement for visiting all pickup nodes and uses multiple vehicles. This problem formulation is pertinent to real-world logistic applications, particularly for enabling some suppliers (pickup nodes) to satisfy thedemands of all customers (delivery nodes) by using multiple vehicles. Furthermore, the selectability of pickup nodes cansubstantially reduce the transportation cost in the PDP. Some problems are related to but differ from the SPDP. The singlecommodity vehicle routing problem with pickup and delivery service (1-VRPPD) [43] uses multiple vehicles to solve theone-commodity pickup-and-delivery problem [37]. The 1-VRPPD involves locating minimum-cost routes to transport homogeneous commodities among unpaired customers, according to the constraints that each node must be visited exactly once
and the depot can balance the total demand. Falcon et al. [25] applied this problem to manage carrier-based coverage repair
in wireless sensor networks, in which pickup nodes are optional. The principal difference between these problems and the
MVSPDP lies in the depot capability and unitary demand. The linehaul-feeder vehicle routing problem with virtual depots[12] enables selectability of reloading spots; virtual depots refer to public parking lots in which a large vehicle dispatches
commodities to small vehicles for relaying delivery in narrow streets. By contrast, in the MVSPDP, each node is visited atmost once to collect a suﬃcient number of commodities and serve all of the delivery customers.
Furthermore, the MVSPDP can be viewed as a combination of the SPDP and the vehicle routing problem (VRP); therefore,
the properties of both the SPDP and VRP must be considered in resolving the MVSPDP. For example, the VRP accounts forthe number of vehicles and depots [23,40], time windows [6,32,42,59,66], workload balance [36], and dynamics [44]. In 1-1transportation, the precedence constraint guarantees that a requested pickup node is visited before reaching the designateddestination [4,47]. Additional constraints on the duration of a single route and the number of passengers disembarking at
a site exist in some industries; for example, oil companies are subject to refueling requirements and security concerns on
the production platform [63]. The multiple requests in emergency transportation operations and the control of automatedguided vehicle dispatching in manufacturing create additional constraints or objectives regarding the passenger ride time,
occupancy rate, and resource required for service quality [26,48,52]. In addition to load splitting, which allows multiple vehicles to serve mutual customers in the VRP, pairing transportation with transfer opportunity enables goods to be stored
temporarily in transshipment nodes, which are visited several times by multiple vehicles to satisfy requests [20,53]. The
1-M-1 PDP is a generalization of the capacitated vehicle routing problem and presents various limitations on solution types[7], including Hamiltonian [5,13,31], delivery-ﬁrst pickup-second [62,65], and lasso tours [38]. Zhu et al. [67] recently considered the 1-M-1 PDP with three objectives and dynamic requests. Regarding intermodal shipments using containers, individual trips are permitted when containers are fully loaded, and a delivery trip and a pickup trip can be merged if thetime and container size are compatible [58]. Noteworthily, compensation for the travel cost is considered when selecting
pickups in the 1-M-1 structure [35], whereas selection of pickups is executed in the MVSPDP to collect a suﬃcient numberof commodities for supplying delivery nodes.
2.2. Bound from the minimum k-degree center tree
This subsection attempts to bound the values of solutions to the MVSPDP. Christoﬁdes et al. [16] derived a bound forthe multiple traveling salesman problem (M-TSP), which can also serve as a lower bound for the VRP because the VRP canbe viewed as an M-TSP with additional constraints. In light of the fact that MVSPDP is a generalization of the VRP, we use
their method to compute the bound for the MVSPDP. Speciﬁcally, a solution to the M-TSP is decomposed into three sets,i.e., St , S0 , and S1 , where St consists of the edges forming a k-degree center tree (k-DCT) in the solution, and S0 and S1
comprise the remaining edges connected to v0 and those between customer nodes, respectively. A k-DCT is a tree whosebranches connect to all vertices in G and the degree of v0 is k. The number of available vehicles limits the maximal numberof edges in S0 , i.e., |S0 | ≤ m, leading to k = 2m − |S0 |. Let binary decision variable ξitj = 1 if (vi , v j ) is in the edge set St , andξitj = 0 otherwise; in addition, let ξi0j and ξi1j denote such binary decision variables for S0 and S1 , respectively. According toChristoﬁdes et al. [16], the M-TSP can be formulated as follows:min

ci j (ξitj + ξi0j + ξi1j )(10)ξitj ≥ 1, ∀S ⊂ V, S = ∅(11)

ξ0t j = 2m − y(12)ξitj = n(13)ξ00j = y, y ∈ Z∗

(14)

vi ,v j ∈Vs.t.

vi ∈S,v j ∈/ S

v j ∈Vvi ,v j ∈Vv j ∈V150

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

ˆ  = 3, and δ = 7.67. The number inside circles denotes demand di . The grounds ofFig. 2. An example initialization for an MVSPDP instance with n = 16, m
the circles indicate respective partitions. The solid, dashed, and dotted lines represent the routes of three vehicles.

Fig. 3. Example of relocation. The dotted circles indicate the node and target place chosen for relocation.

Fig. 4. Example of inversion. The dotted circle indicates the part to be inverted.

ξi1j = m − y(15)

vi ,v j ∈V ±(ξitj + ξi0j + ξi1j ) = 2,∀vi ∈ V ±(16)v j ∈V

Constraint (11) guarantees the connectivity of the k-DCT, in which the degree of v0 is limited by constraint (12). Constraints(13), (14), and (15) specify the numbers of edges required for St , S0 , and S1 , respectively. Constraint (16) conﬁnes the degreeof each vertex in an M-TSP to two, except for v0 .
The Lagrangian relaxation produces three decomposed problems Pt , P0 , and P1 for each candidate value of latent variable
y and introduces constraint (16) into a general objective:gκ ( p, y ) =vi ,v j ∈V

(ci j + pi + p j )ξiκj − 2



pi ,

(17)

vi ∈V ±

where p = ( p1 , . . . , pn ) represents the non-negative penalties associated with the n customer nodes, and κ ∈ {t, 0, 1} denotes
the index of ξitj , ξi0j , and ξi1j . For a given value of y, let gt (p, y) be the value of optimal solution to Pt deﬁned by (17), (11),(12), and (13); let g0 (p, y) be the value of optimal solution to P0 deﬁned by (17) and (14); and let g1 (p, y) be the value of

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

151

Fig. 5. Example of selection: a) selecting a pickup node; b) unselecting a pickup node.

optimal solution to P1 deﬁned by (17) and (15). A lower bound to the M-TSP can be computed by

max max gt ( p, y ) + g0 ( p, y ) + g1 ( p, y ) .
y≤m

p≥0

(18)

Note that this lower bound is valid for the VRP and the MVSPDP, albeit omitting the constraints on the vehicle load andtravel distance. Restated, the values of feasible solutions to these two problems must be no smaller than this lower bound;
otherwise, contradiction will occur during the approximation of the decomposed problems.
3. Metaheuristics for the MVSPDP
This study develops three metaheuristic algorithms for the MVSPDP: tabu search (TS), genetic algorithm (GA), and scattersearch (SS). A ﬁxed-length representation is presented for the three algorithms to accommodate the MVSPDP. TS is initializedusing the modiﬁed sweep algorithm, whereas GA and SS are initialized by combining the modiﬁed sweep algorithm withrandom initialization to increase the population diversity.
In addition, we propose four operators serving as the neighborhood functions in TS, mutation operators in GA, and improvement operators in SS. Using multiple search operators in metaheuristic algorithms has shown to be helpful for exploring different structures in a problem and achieved considerable successes [14,54,57,60]. In view of the beneﬁts, this study
devises four operators for the three algorithms to deal with the visiting order, selection of pickup nodes, and the number ofvehicles adopted.
The designs of the proposed algorithms are detailed in the subsequent sections.
3.1. Tabu search
Tabu search is a metaheuristic algorithm known for its effectiveness in combinatorial optimization problems. This metaheuristic algorithm uses an explicit memory structure to record the search trajectory. Based on recorded information, TS
guides the search considering both intensiﬁcation and diversiﬁcation. In TS, the search is processed using a series of moves
from a solution to its neighbor, where the neighborhood is generally deﬁned by a distance function, such as the Hamming distance. The tabu-versus-aspiration strategy controls the moves: On the one hand, TS uses the tabu list to record theforbidden moves and prevents the search from being mired in the local optima by indicating certain moves as tabu. The
tabu tenure affects the level of restriction on choosing the forbidden moves; long tabu tenure guides the search to explore

152

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

Fig. 6. Example of addition/deletion.

unvisited territory, encouraging diversiﬁcation. On the other hand, the aspiration criterion provides the opportunity to override the tabu restriction. That is, the aspiration criterion enables superior neighboring solutions to be chosen, despite therestriction of tabu moves. This criterion supports the ability of intensiﬁcation for the search.
Designing TS involves solution representation, an evaluation function, initialization, a neighborhood function, and an
adaptive search phase. The following subsections describe our designs of TS for the MVSPDP .
3.1.1. Representation and evaluation function
This study presents a representation for candidate solutions to the MVSPDP based on the chromosome representationof Tan et al. [57]. In addition to indicating the visiting orders of nodes for vehicles, the proposed representation considers
selecting pickup nodes in the MVSPDP. Speciﬁcally, the representation consists of two parts: The ﬁrst part comprises m + 1
nonnegative integers, of which the ﬁrst m integers indicate the number of nodes for m vehicles to visit, and the ﬁnal integer
provides the number of unselected pickup nodes. In the example representation in Fig. 1, the ﬁrst value “5” indicates thatthe ﬁrst vehicle must visit ﬁve nodes, while the fourth value “0” signiﬁes that the fourth vehicle is unused and, therefore,
has no corresponding route. The second part of the representation indicates the visiting order of nodes for each vehicle.
The ﬁnal sequence records the unselected pickup nodes, which are not visited by any vehicle. Noteworthily, although them + 1 sequences are variable-length, the sum of their lengths is ﬁxed to the number of nodes n. Hence, the proposed
representation applies a ﬁxed-length encoding, which is capable of indicating the varying number of vehicles used and theselection of pickup nodes.
In TS, the evaluation function determines the promising directions and guides the search toward the optima. The proposed evaluation function considers the objective of the MVSPDP and deals with infeasible solutions. Formally, the evaluation function for a solution x is deﬁned asf (x ) =m

ci j xi jk + α pd + β c̄ pl ,(19)

k=1 vi ,v j ∈V

where c̄ is the average cost among all edges, and α and β denote the weights for penalizing the constraint violation on
loading and distance; both weights are set to 2 in this study. The penalty for violating load constraint pl counts the numberof times that vehicles are overloaded or have insuﬃcient commodities on board, and penalty pd accounts for the travel

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

153

Fig. 7. Crossover for the proposed GA: whole arithmetic crossover for the ﬁrst part and order crossover for the second part. The gray squares denotethe unselected pickup nodes, and the dashed lines indicate the cut points for order crossover. The green numbers in the resultant offspring highlight thechanges in the number of visited nodes. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version ofthis article.)

Fig. 8. Example of SS creating a solution from a and b. The dotted line indicates the cut point.

154

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

Table 1MVSPDP test instances.
Instance

mvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

Modiﬁed from

vrpnc01
vrpnc06
vrpnc02
vrpnc07
vrpnc03
vrpnc08
vrpnc12
vrpnc14vrpnc04
vrpnc09
vrpnc05
vrpnc10vrpnc11
vrpnc13kelly05
kelly01
kelly06
kelly02
kelly07
kelly03kelly08
kelly04

m

6
7
12
13
10
10
14
11
16
16
22
22
10
13
6
13
11
13
12
13
15
14

Q

85
85
85
85
105
95
140
115
130
100
110
105
110
65
450
366
450
350
450
450
600
500

R

180
260
190
230
280
290
300
1100
280
300
300
300
350
800
30 0 0
1083
2500
1500
2166
2400
20 0 02666#customers
Pickup |V + |Delivery |V − |

26
27
39
41
49
44
46
46
69
73
96
91
62
65
98
128
142
172
188
202
231
238

24
23
36
34
51
56
54
54
81
77
103
108
58
55
102
112
138
148
172
198
209
242

Table 2
Parameter setting.
Parameter

Value(a) Tabu search.
Representation
Initialization
Neighborhood function
Neighborhood search
Tabu tenure
Neighborhood Size
Termination#Evaluations per iteration
Total #evaluations

Integer + order
Modiﬁed sweep algorithm
Relocation, inversion, selection, addition/deletion
RandomAdaptive
200
1500 iterations
200
30 0,0 0 0(b) Genetic algorithm.
Representation
Initialization
Recombination (pc )Mutation (pm )Parent selection
Survival selection
Population size
Termination#Evaluations per generation
Total #evaluations

Integer + order
Modiﬁed sweep algorithm (50%), random (50%)
Whole arithmetic (0.9), order crossover (0.9)
Relocation (0.6), inversion (0.6), selection (1/|V + |), addition/deletion (1/m)
2-tournamentμ+λ
100
30 0 0 generations
100
30 0,0 0 0(c) Scatter search.
Representation
Initialization
Population size
Reference size
Improvement method
HC neighborhood function
HC neighborhood size
HC termination
Termination#Evaluations per generation
Total #evaluations

Integer + order
Modiﬁed sweep algorithm (50%), random (50%)
25
5Hill climbing (HC)
Relocation, inversion, selection, addition/deletion100
10 iterations
60 generations
5010
325,625

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

155

Table 3
The average route length (avg.), feasible rate (feas.), and average running time (time) for the TS, GA, and SS using theﬁxed probability of applying the four operators on the MVSPDP instances with γ = 0. Boldface marks the best result ineach instance.
InstanceTS-ﬁxedmvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

GA-ﬁxed

SS-ﬁxed

Avg.

Feas. (%)

Time (s)Avg.

Feas. (%)

Time (s)Avg.

Feas. (%)

Time (s)

473
434
622
568
723
716
656
529
833
853
938
978
1186
693
8440
6584
10,51810,35912,745
14,95415,751
17,798

90
93
93
97
83
93
87
100
87
97
97
90
97
93
17
90
100
100
97
97
93
90

0.16
0.21
0.19
0.23
0.28
0.30
0.34
0.40
0.34
0.35
0.40
0.40
0.25
0.50
0.77
0.40
0.63
0.52
0.59
0.70
0.58
0.64

475
444
642
572
757
739
629
500
921
886
1090
1142
1623
825
11,109
–
–
–
–
–
–
–

93
93
77
90
53
93
90
100
80
87
83
87
50
80
3
0
0
0
0
0
0
0

2.15
3.12
4.95
6.15
8.35
8.53
9.34
8.51
17.51
17.40
28.85
29.21
11.98
12.12
24.16
–
–
–
–
–
–
–

473
444
668
601
810
745
625
522
995
985
1169
1149
1372
800
–
–
–
–
–
–
–
24,690

73
87
83
97
80
90
93
100
87
83
50
87
37
90
0
0
0
0
0
0
0
3

0.19
0.24
0.22
0.25
0.29
0.33
0.58
0.61
0.47
0.47
0.84
0.85
0.46
0.76
–
–
–
–
–
–
–
5.43

Table 4
The average route length (avg.), feasible rate (feas.), and average running time (time) for the TS, GA, and SS using the ﬁxed
probability of applying the four operators on the MVSPDP instances with γ = 32. Boldface marks the best result in each
instance.
InstanceTS-ﬁxedmvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

GA-ﬁxed

SS-ﬁxed

Avg.

Feas. (%)

Time (s)Avg.

Feas. (%)

Time (s)Avg.

Feas. (%)

Time (s)

404
365
496
442
564
603
515
464
666
664
781
812
775
563
6190
4529
7560
6512
8643
9299
10,066
12,320

97
100
93
93
100
100
97
100
100
100
97
100
100
100
100
97
97
97
93
100
97
93

0.15
0.20
0.19
0.23
0.27
0.28
0.43
0.49
0.34
0.37
0.39
0.39
0.25
0.55
0.70
0.38
0.82
0.57
0.54
0.76
0.64
0.82

406
363
500
453
558
602
514
452
690
690
814
862
807
579
7112
5804
10,0749354
13,725
15,611
17,57921,865

93
100
93
93
100
93
97
100
97
90
90
93
93
100
93
83
90
73
70
87
57
57

2.36
3.08
6.07
6.07
8.61
7.94
8.51
7.90
17.33
17.02
28.96
27.93
12.00
12.62
21.64
30.53
37.50
50.43
61.40
73.16
89.72
99.08

407
373
510
448
570
611
503
461
673
679
797
837
793
625
6982
6943
9841
10,825
13,602
24,56911,913
15,179

97
100
87
87
93
100
90
100
100
93
97
90
93
100
97
63
90
53
10
23
93
93

0.17
0.21
0.22
0.25
0.31
0.33
0.45
0.44
0.43
0.43
0.55
0.58
0.34
0.49
0.84
1.18
1.30
1.95
2.48
3.24
3.37
3.42

distance exceeding the limitation R. That is,

pl =

m 

k=1 S⊆V ±

pd =

m

k=1

Sk ,



max 0,(20)
vi ,v j ∈V

ci j xi jk − R ,

(21)

156

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169
Table 5The average route length (avg.), feasible rate (feas.) in percentage, and average running time (time) in seconds, for theTS, GA, and SS using the adaptive probability of applying the four operators, in comparison with CPLEX, on the MVSPDP
instances with γ = 0. Boldface marks the best result in each instance.
Instance

mvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

TS-adaptive

GA-adaptive

SS-adaptive

CPLEX

Avg.

Feas.

Time

Avg.

Feas.

Time

Avg.

Feas.

Time

Length

Time

477
433
623
562
720
715
621
508
827
808
898
927
1152
629
8288
6490
10,367
10,112
12,373
14,296
15,22817,357

73
100
77
93
73
93
87
100
97
100
100
97
100
100
17
93
100
93
100
100
100
97

0.15
0.20
0.17
0.21
0.24
0.25
0.46
0.65
0.29
0.31
0.35
0.36
0.30
0.76
0.46
0.40
0.45
0.49
0.53
0.55
0.55
0.65

468
437
688
635
815
778
714
499
1207
1258
1636
1713
1680
1336
–
–
–
–
–
–
–
–

97
100
73
80
17
100
37
100
23
10
10
13
10
30
0
0
0
0
0
0
0
0

1.84
2.99
3.57
5.50
8.15
7.99
7.99
8.16
16.19
16.54
26.96
26.33
11.66
11.84
–
–
–
–
–
–
–
–

482
447
667
590
806
747
614
519
1080
911
1232
1291
1565
767
–
–
–
–
–
–
–
22,260

87
97
83
93
80
93
90
100
77
73
43
57
57
100
0
0
0
0
0
0
0
3

0.17
0.23
0.22
0.24
0.29
0.35
0.46
0.70
0.48
0.68
0.85
1.03
0.54
0.77
–
–
–
–
–
–
–
4.07

480
428
760
767
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–

289408321872
12183421345694
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–

Table 6
The average route length (avg.), feasible rate (feas.) in percentage, and average running time (time) in seconds, for the TS, GA,
and SS using the adaptive probability of applying the four operators, in comparison with CPLEX, on the MVSPDP instanceswith γ = 32. Boldface marks the best result in each instance.
InstanceTS-adaptive

GA-adaptive

SS-adaptive

CPLEX

Avg.

Feas.

Time

Avg.

Feas.

Time

Avg.

Feas.

Time

LengthTime

mvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

405
369
496
453
560
595
515
459
654
645
761
793
785
528
5937
4481
7297
6341
8391
8938
9677
12,017

90
100
90
100
93
97
93
100
93
90
90
93
100
100
100
97
90
97
97
97
97
100

0.14
0.18
0.18
0.18
0.24
0.24
0.31
0.44
0.29
0.31
0.34
0.36
0.24
0.43
0.54
0.37
0.60
0.39
0.50
0.63
0.54
0.70

399
363
536
471
596
629
561
456
824
804
1078
1157
1041
847
7164
5849
10,0348674
13,333
14,025
17,557
20,471

97
100
97
97
97
93
83
100
100
100
100
97
93
100
100
93
100
100
100
97
67
73

3.43
3.45
6.11
6.33
9.36
9.05
9.10
9.50
18.41
19.04
31.78
30.83
14.82
15.42
30.26
45.25
57.91
78.93
98.19
112.96147.43134.87

407
373
509
458
574
609
505
458
685
679
801
829
784
605
6758
6433
8730
9713
16,78024,18010,976
13,348

90
100
90
97
93
100
90
100
97
97
100
100
100
100
97
63
80
47
3
30
93
90

0.18
0.22
0.23
0.25
0.31
0.33
0.61
0.64
0.41
0.41
0.55
0.61
0.45
0.61
0.76
1.29
1.21
1.88
2.64
3.23
3.03
2.75

361
343
458
392
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–

29831932783013766431485938
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–
–

where


Sk =

0

0 ≤ zSk


vi ∈S

1

di yik ≤ Q,

(22)otherwise.

The evaluation function (19) considers route length as well as the penalties for violating the constraints on vehicle load andmaximal travel distance. Accordingly, TS tends to locate the shortest feasible routes for a ﬂeet.

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

157

3.1.2. Initialization
Beyond random initialization, this study proposes an improved initialization method based on the sweep algorithm. The
sweep algorithm has been widely used for dispatching vehicles [27]. It assigns disjoint sets of customers to vehicles according to the polar-coordinate angle, and the routes are subsequently improved by the Lin-Kernighan algorithm [41], yielding anear-optimal solution to the single-depot vehicle-dispatch problem. Given that the number of pickup nodes is variable dueto their selectability in the MVSPDP, we perform such an assignment on only the delivery nodes. In addition, we estimate
the delivery requests served for a vehicle. First, the number sˆ of selected pickup nodes can be approximated as follows:



1− 

sˆ =


vi ∈V div ∈V + d j

where vi ∈V div j ∈V + d j| V + |,

(23)

ˆ can then beis the ratio of surplus commodities to the total supply. The number of required vehicles mjestimated by(sˆ + |V − | ) · c̄ˆ =m

R − c̄

+ 0.5 ,

(24)

where the numerator (sˆ + |V − | ) · c̄ provides an approximate travel distance for visiting the selected pickup nodes and alldelivery nodes, and the denominator (R − c̄ ) represents the maximal travel distance permitted for a vehicle, where c̄ issubtracted because the route in the numerator omits one edge (i.e., the edge from the ﬁnal node to the depot). Accordingly,
the estimated number of delivery requests served by a vehicle can be determined usingδ=vi ∈V − |di |
,
ˆ
m

(25)

ˆ  ∈ [min{mˆ , m}, m] is randomly chosen to diversify the initial state.
where mTo generate an initial solution, we sort the vertices according to the polar-coordinate angle and randomly select one
vertex as the starter. Following the sorting order, delivery nodes are sequentially assigned to a vehicle until the number ofrequired commodities exceeds δ . The ﬁnal node is then randomly determined to be included in the route or a new route
served by the next vehicle. This assignment separates the sorted vertices into partitions. Accordingly, each vehicle begins at
the ﬁrst pickup node in its partition and adds subsequent nodes (pickup or delivery) in ascending order of polar-coordinate
angle. The pickup nodes are stacked in the waiting list according to the sorting order. Given that a nonnegative vehicle
load must be maintained, a pickup node is popped from the waiting list whenever a negative vehicle load occurs. Fig. 2
illustrates the proposed initialization, where v4 in Fig. 2a is selected as the starter, producing sorted vertices as shown in
Fig. 2b. Fig. 2c presents three partitions, namely {v1 , . . . , v5 }, {v6 , . . . , v10 }, and {v11 , . . . , v16 }, in which v5 and v11 are theﬁnal nodes assigned: v5 is included in the ﬁrst route, but v11 is allocated to the subsequent vehicle. Additionally, v12 is
the ﬁrst pickup node of the third partition, but its supply is insuﬃcient before visiting v13 ; therefore, we pop v8 from thestack to keep the vehicle load nonnegative, causing crosses in the initial routes. Although it may bring about infeasible
solutions for the travel distance and vehicle load, using δ and polar-coordinate sorting can effectively generate promising
initial solutions.
3.1.3. Neighborhood functions
In this study, we devise four operators to construct the neighbors of the current solution. The four operators address
different aspects of the MVSPDP: The relocation operator and the inversion operator account for inter-route and intra-routevariations, respectively. The selection operator processes the selection of pickup nodes in the MVSPDP. Furthermore, theaddition/deletion operator adjusts the number of vehicles used.
For generating the neighbors, the proposed TS randomly performs one of the four operators on the current solution andrepeats this procedure until the predetermined neighborhood size is achieved. Two methods are proposed for determining
the probability of an operator being selected: ﬁxed probability and adaptive probability control. More details on the four
operators are presented below.
Relocation. This operator produces an inter-route variation by moving a randomly chosen node to a previous or subsequentroute. The destination is limited to keep the sequence of the polar angles of nodes as complete as possible. Fig. 3 illustrates
an application of the relocation operator in which node v9 is randomly chosen and relocated to the subsequent route.
Inversion. The inversion operator produces an intra-route variation that inverts a random partial route to alter the visitingorder. An example of this operator is presented in Fig. 4. The chosen partial route v5 → v9 → v12 is inverted into v12 → v9 →v5 .
Selection. The selection operator addresses the key factor in the MVSPDP: the selection of pickup nodes. The operator isdesigned to randomly move a pickup node among the m + 1 sequences. The selection operator randomly chooses a pickup
node from the ﬁnal sequence and then inserts this node into one of the ﬁrst m routes as a selected pickup node. In addition,

158

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

Fig. 9. Comparison of average route length, feasible rate, and average running time for the TS, GA, and SS using the ﬁxed probability of applying the four
operators on the MVSPDP instances with γ = 0 (left) and 32 (right).

this operator may unselect a pickup node by moving it from the ﬁrst m routes to the ﬁnal sequence. Note that in the
representation, the ﬁrst m sequences indicate the routes for m vehicles while the last sequence consists of the unselected
pickup nodes. As illustrated in Fig. 5, node v1 is selected and inserted between v2 and v7 ; by contrast, pickup node v13 isunselected as it is moved to the ﬁnal sequence.
Addition/Deletion. The number of vehicles used affects the routes and total distance traveled by all vehicles in the MVSPDP.
This study develops the addition/deletion operator to adjust the number of vehicles used; whether vehicles are added ordeleted is determined at uniform random. When adding a vehicle, the operator randomly selects a vehicle and then splits
its route into two parts: the ﬁrst part is the route for the selected vehicle, and the second part is separated as a route for a
new vehicle. For example, the new route in Fig. 6a is separated from the second route.
When deleting a vehicle, the operator merges a randomly selected route with its adjacent route in the ﬁrst part of
the representation. However, the merger may be infeasible because of vehicle overload. To address this issue, we propose aheuristic to trim unnecessary pickup nodes. More speciﬁcally, this method iteratively removes the last visited pickup node bymoving it to the sequence of unselection if this elimination can keep the number of remaining commodities nonnegative. Inother words, the demand of removed pickup nodes must be less than or equal to v ∈V di yik . Fig. 6b illustrates the deletioni
operation. Since the merged route has ﬁve remaining commodities, the deletion operator removes v13 with d13 = 3 to reduce

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

159

Fig. 10. Comparison of average route length, feasible rate, and average running time for the TS, GA, and SS using the adaptive probability of applying thefour operators, in comparison with CPLEX, on the MVSPDP instances with γ = 0 (left) and 32 (right).

the number of unnecessary pickup nodes. This procedure terminates at v6 because d6 = 3, which exceeds v ∈V di yik = 2 atithat time. The proposed method reduces the probability of overload and implicitly shortens the merged route.
3.1.4. Adaptive search phase
The proposed TS comprises two search phases: normal and intensiﬁcation phases. Each phase manipulates the recencybased short-term memory to achieve a particular strategy. In the TS, the tabu list records only the edges added or removed
by neighborhood functions instead of the entire solution. The list is used to forbid the change of these edges during tabu
tenure to increase diversiﬁcation. The tabu tenure is randomly initialized within {5, . . . , 15} in this study. Additionally, theintensiﬁcation phase halves the tabu tenure for the TS to exploit the neighborhood of elite solutions.
The search phase is adaptively controlled according to variable φ . The TS conducts intensiﬁcation only if φ is below thethreshold θ ; otherwise, it performs a normal search. The intensiﬁcation phase empties tabu list and halves the tabu tenure,and then restarts the search with the best-so-far solution. The value of φ is adjusted according to counter r, which increases
as the best solution improves, and decreases otherwise. Speciﬁcally, φ is varied by
φ · κ if r = −15 and φ > 0.2
φ = φ /κ if r = 20 and φ < 0.8
φotherwise

(26)

160

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

Fig. 11. Anytime behavior of the average route length over 30 runs of TS, GA, and SS using the adaptive probability, in comparison with CPLEX and its
bound, on four MVSPDP instances with γ = 0.

As the search stagnates so long that r reaches −15, the value of φ decreases due to the positive coeﬃcient κ < 1. If thissituation consecutively occurs τ times, then the search is regarded as frozen, and thus we inactivate r and φ for tf generations. In addition, the tabu tenure is doubled in subsequent iterations to encourage exploring the unvisited regions andescape from the local optima. In this study, the coeﬃcients κ and φ are empirically set to 0.8 and 0.4, respectively; θ is
ﬁxed to 0.2; τ = max{2, n/50 + 0.5 }; and the frozen time tf is 50 generations.
In summary, the adaptive search phase adjusts tabu tenure in response to the search progress for a balance between
exploitation and exploration. Algorithm 1 shows the procedure of the proposed TS.
3.2. Genetic algorithm
Genetic algorithm is well known for its effectiveness in global optimization. Generally, GA uses a set (population) of
candidate solutions (individuals or chromosomes) that can exchange information with each other by using genetic operators to explore the solution space. The operations in GA include selection, crossover, and mutation. The selection operator
chooses individuals from the population to generate new candidate solutions (offspring) through crossover and mutation.
The crossover operator exchanges and recombines the genetic information of selected individuals to produce offspring, andthe mutation slightly changes offspring to provide diversity. In GA, the ‘survival of the ﬁttest’ principle is applied in selecting
individuals with highest ﬁtness values to survive into the next generation.
This study involves developing a GA for the MVSPDP. The designs focus on representation, ﬁtness function, and genetic
operators. The details on the proposed GA are described in the subsequent sections.
3.2.1. Representation and ﬁtness function
Individuals in the GA are represented using the representation proposed for the TS. Restated, an individual is composedof two parts: the ﬁrst part indicates the use of vehicles, and the second part represents their respective routes, i.e., visit-

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

161

Fig. 12. Anytime behavior of the average route length over 30 runs of TS, GA, and SS using the adaptive probability, in comparison with CPLEX and its
bound, on four MVSPDP instances with γ = 32.

ing orders of nodes (cf. Fig. 1). The half-half strategy is adopted for population initialization: the TS initialization operator
generates half of the initial population; the other half is generated at random to provide diversity. The ﬁtness function uses(19), which accounts for the route length and the constraints on vehicle capacity and travel distance in the MVSPDP.

3.2.2. Genetic operators
The genetic operators in the GA include selection, crossover, and mutation. The selection operators, including parent
selection and survival selection, are associated with the ﬁtness values of individuals. This study adopts 2-tournament parent
selection because of its recognized performance [2,22]. The 2-tournament selection randomly picks two individuals from the
population and selects the better individual as a parent. Performing this selection twice yields a pair of parents to be usedin subsequent reproduction (i.e., crossover and mutation).
The crossover operator generates offspring by exchanging the information of two selected candidate solutions. In thelight of the representation of individuals, we adopt the whole arithmetic crossover [45] with α = 0.8 on the integer stringof the ﬁrst part and conduct the order crossover [21] on the merged routes in the second part. Note that the offspring
inherit the selection of pickup nodes from the second part of their parents. However, the sum of resultant integer strings
may exceed or fall short of the total number of customers. To address this issue, the proposed crossover operator entails
selecting a random gene and reducing or increasing its value by 1 until reaching the requirement. The merged route isthen divided among vehicles according to the designated numbers in the integer string. Fig. 7 demonstrates the proposedcrossover operator with which the offspring obtain new routes as well as a diverse selection of pickup nodes.
Moreover, this study utilizes the four neighborhood functions of TS as the mutation operators for the GA. The relocationand inversion operators are used to produce inter-route and intra-route variations, respectively. The selection operator alters
the selection of pickup nodes in the MVSPDP. The addition/deletion operator can increase or decrease the number of vehicles

162

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

Algorithm 1: Tabu search for the MVSPDP.
Input: MVSPDP instance
Output: solution x∗
initialize(x) ;x∗ ← x ;tenure ← rand(5, 15 );phase ← Normal;
r, c, t ← 0, φ ← 0.4;
repeatN ← neighborhood(x );repeatx ← best of N;
N ← N\{x };until f (x ) < f (x∗ ) or (x x ) ∈/ T abuList;
adjust tenure according to phase;
update T abuList with x x and tenure;x ← x ;if f (x ) < f (x∗ ) then x∗ ← x, r ← r + 1;
else r ← r − 1;adjust φ according to r;
if c = τ thent ← tf ;tenure ← 2 · tenure;t ← t − 1;
if φ ≤ θ and t < 0 thenif phase = Normal thenx ← x∗ ;T abuList ← ∅;phase ← Intensiﬁcation;c ← c + 1;
else
phase ← Normal;c ← 0;
if t = 0 then r ← 0, φ ← 0.4;// modified sweep algorithm// x∗ :best-so-far solution//:symmetric difference// Eq. (26)// intensification phase// normal phase

until terminated;

used. Each operator has a speciﬁc mutation rate, which is predetermined or adapted during the run. Section 4 examines theeffects of these two strategies.
The selection-crossover-mutation process is repeated until the offspring population is generated. The survival selection then picks out the ﬁtter individuals for the subsequent generation. This study uses (μ + λ) survival selection, in
which the ﬁttest μ individuals are selected from the union of μ parents and λ offspring to survive into the next
generation.

3.3. Scatter search
Scatter search manipulates the population of solutions and uses a combination method to create new solutions [30]. The
key components of SS are listed as follows.
• Diversiﬁcation generation: A collection of diverse solutions is created using arbitrary or seed solutions.
• Improvement: A solution is transformed into one or more enhanced solutions.
• Reference set update: This procedure develops and maintains a reference set containing numerous best solutions
found.
• Subset generation: A subset of solutions is generated from the reference set as the basis for solution combination.
• Solution combination: The subset of generated solutions is combined into one or more solutions.

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

163

Table 7
The average route length (avg.) and feasible rate (feas.) for the GA using the ﬁxed and adaptive probability of applying the four operators.
Column “time inc.” lists the increment (or decrement) in the running time of GA-adaptive over GA-ﬁxed. Boldface indicates the higher feasible
rate and signiﬁcantly shorter average route length according to a one-tailed t-test with conﬁdence level α = 0.05. The dash indicates that no
feasible solution is found in 30 runs.

γ =0Instance

γ = 32GA-ﬁxedmvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

GA-adaptive

GA-ﬁxed

GA-adaptive

Ave.

Feas. (%)

Ave.

Feas. (%)

Time inc. (%)

Ave.

Feas. (%)

Ave.

Feas. (%)

Time inc. (%)

475
444
642
572
757
739
629
500
921
886
1090
1142
1623
825
11,109
–
–
–
–
–
–
–

93
93
77
90
53
93
90
100
80
87
83
87
50
80
3
0
0
0
0
0
0
0

468
437
688
635
815
778
714
499
1207
1258
1636
1713
1680
1336
–
–
–
–
–
–
–
–

97
100
73
80
17
100
37
100
23
10
10
13
10
30
0
0
0
0
0
0
0
0

−14.42−4.17−27.88−10.57−2.40−6.33−14.45
−4.11
−7.54
−4.94
−6.55
−9.86
−2.67
−2.31

406
363
500
453
558
602
514
452
690
690
814
862
807
579
7112
5804
10,0749354
13,725
15,611
17,57921,865

93
100
93
93
100
93
97
100
97
90
90
93
93
100
93
83
90
73
70
87
57
57

399
363
536
471
596
629
561
456
824
804
1078
1157
1041
847
7164
5849
10,0348674
13,333
14,025
17,557
20,471

97
100
97
97
97
93
83
100
100
100
100
97
93
100
100
93
100
100
100
97
67
73

45.34
12.01
0.66
4.28
8.71
13.98
6.93
20.25
6.23
11.87
9.74
10.38
23.50
22.19
39.83
48.21
54.43
56.51
59.92
54.40
64.32
36.12

3.3.1. Diversiﬁcation generation and improvement
The representation and evaluation function of SS follow those of the proposed TS and GA. In addition, the diversiﬁcation
generation method uses the half-half strategy, in which both the modiﬁed sweep algorithm and random initialization are
adopted, as the proposed GA does. The improvement method is based on hill climbing (HC), for which the neighborhood
function employs the four operators of TS presented in Section 3.1.3.
3.3.2. Solution combination
The combination method exchanges the information collected during the search. The proposed combination method isused to create two solutions from a pair of routes in the reference set to retain their characteristics. Fig. 8 illustrates theprocedure for creating a solution. The solution c1 inherits the ﬁrst segment (i.e., the route of the ﬁrst vehicle) from solution
a and sequentially ﬁlls the remaining genes without duplicates from solution b. The split procedure is then applied to cut
the solution and determine the number of nodes that each vehicle serves. The other solution is generated similarly but inan opposite direction, that is, by copying the ﬁrst segment of b and ﬁlling the remainder from a.
In this study, the split procedure is based on the saving heuristic [17]. The cost (distance) before linking two nodes vi
and v j is c0i + ci0 + c0 j + c j0 , whereas the cost after linking them is c0i + ci j + c j0 . Thus, the beneﬁt s from aggregation of vi
and v j issi j = ci0 + c0 j − ci j .

(27)

The edges are linked according to the amount of saved cost. However, the constraints on travel distance (6) and vehicle
capacity (9) may cause an individual to be cut into too many segments when using the split procedure; consequently, thegenerated solutions become too similar. To address this problem, the split procedure uses two alternative constraints,



2



vi ∈V


vi ,v j ∈V

di
−|
|
V
−



≤ zSk

ci j xi jk ≤ R −

di yik ≤ Q

(28)

vi ∈S

n · c̄
50(29)

to relax the lower bound of the vehicle load and tighten the upper bound of the travel distance, respectively.
3.3.3. Reference set
The reference set is updated through an operation similar to the selection operation of GA, which chooses solutions togenerate new solutions. In SS, the reference set Ref is updated half from the original solutions and half from the newly

164

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169
Table 8The average route length (avg.) and feasible rate (feas.) for the TS using the ﬁxed and adaptive probability of applying the four operators. Column“time inc.” lists the increment (or decrement) in the running time of TS-adaptive over TS-ﬁxed. Boldface indicates the higher feasible rate andsigniﬁcantly shorter average route length according to a one-tailed t-test with conﬁdence level α = 0.05. The dash indicates that no feasible
solution is found in 30 runs.
Instanceγ =0

γ = 32TS-ﬁxedmvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

TS-adaptive

TS-ﬁxedTS-adaptive

Ave.

Feas. (%)

Ave.

Feas. (%)

Time inc. (%)

Ave.

Feas. (%)

Ave.

Feas. (%)

Time inc. (%)

473
434
622
568
723
716
656
529
833
853
938
978
1186
693
8440
6584
10,51810,35912,745
14,95415,751
17,798

90
93
93
97
83
93
87
100
87
97
97
90
97
93
17
90
100
100
97
97
93
90

477
433
623
562
720
715
621
508
827
808
898
927
1152
629
8288
6490
10,367
10,112
12,373
14,296
15,22817,357

73
100
77
93
73
93
87
100
97
100
100
97
100
100
17
93
100
93
100
100
100
97

−6.25
−4.76
−10.53−8.70−14.29−16.6735.29
62.50−14.71−11.43−12.50−10.0020.00
52.00−40.260.00−28.57−5.77−10.17−21.43
−5.17
1.56

404
365
496
442
564
603
515
464
666
664
781
812
775
563
6190
4529
7560
6512
8643
9299
10,066
12,320

97
100
93
93
100
100
97
100
100
100
97
100
100
100
100
97
97
97
93
100
97
93

405
369
496
453
560
595
515
459
654
645
761
793
785
528
5937
4481
7297
6341
8391
8938
9677
12,017

90
100
90
100
93
97
93
100
93
90
90
93
100
100
100
97
90
97
97
97
97
100

−6.67
−10.00−5.26−21.74−11.11−14.29−27.91−10.20−14.71−16.22−12.82−7.69−4.00−21.82−22.86−2.63−26.83−31.58−7.41−17.11−15.63−14.63

generated solutions. Speciﬁcally, from the original solutions, the b1 best and |Re2 f | − b1 randomly-selected solutions are included in the reference set. Similarly, the remaining half of the reference set is composed of the b2 best and |Re2 f | − b2
randomly-selected solutions generated by the combination method.
4. Experimental results
This study conducts a series of experiments to evaluate the proposed TS, GA, and SS for the MVSPDP. The test instances
are modiﬁed from the Capacitated VRP Instances [11], ranging from 50 to 480 customer nodes. In modifying these instances
for the MVSPDP, we change the demands of some random nodes to be negative, subject to the total demand being nonnegative. Table 1 lists the information of instances used in the experiments, where the total number of customers is indicated in the instance name. In addition, we introduce parameter γ as a gain in supply for each pickup node, i.e., di = di + γ
for all vi ∈ V + , to investigate the inﬂuence of the selectability of pickup nodes, where the vehicle capacity is increased by 2γ
as well. The test instances can be downloaded via http://cilab.cs.ccu.edu.tw/MVSPDP.zip. Table 2 lists the parameter setting
for TS, GA, and SS in the experiments. Each algorithm is tested with 30 independent runs, considering its stochastic nature.
The experiments are performed on Windows 7 platform and Intel i7-920 machines.
First, we investigate the performance of each algorithm using ﬁxed and adaptive control strategies for the probability of executing four operations, namely relocation, inversion, selection, and addition/deletion. In the experiments of ﬁxed
probability, the neighborhood function of TS and SS randomly conducts one of the four proposed operations with equal
probability, whereas GA uses commonly suggested mutation rates for the four operators and enables more than one operator to be used in a mutation operation. Tables 3 and 4 list the results of TS, GA, and SS on the MVSPDP instances with
γ = {0, 32}, where the dash indicates that no feasible solution is found in 30 runs. The feasible rate represents the percentage of runs achieving feasible solutions among the 30 runs. Boldface on the feasible rate marks the best result for the
instance, whilst boldface on the average route length denotes that the algorithm signiﬁcantly outperforms the second best
method according to a one-tailed t-test with conﬁdence level α = 0.05, and multiple boldfaced results infer a nonsigniﬁcantdifference between them. Fig. 9 further compares the results of the three metaheuristic algorithms. The experimental results
show that, as the number of nodes increases, TS outperforms both GA and SS in the average route length and feasible rate.
Precisely, TS excels with statistical signiﬁcance in 17 and 11 out of 22 instances for γ = 0 and 32, respectively. Comparedwith TS, GA and SS hardly yield feasible solutions for the MVSPDP instances with > 200 nodes and γ = 0. Furthermore, TS
uses lowest computation time, whereas GA requires more than 10 times computation time of TS and SS. These outcomesindicate the effectiveness and eﬃciency of the proposed TS for resolving the MVSPDP.
Second, this study proposes adaptive control over the probability of performing each operator. For TS and SS, this probability is determined according to the proportion of times that a particular operator is selected as the neighborhood function

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

165

Table 9The route length (avg.) and feasible rate (feas.) for the SS using the ﬁxed and adaptive probability of applying the four operators. Column “timeinc.” lists the increment (or decrement) in the running time of SS-adaptive over SS-ﬁxed. Boldface indicates the higher feasible rate and signiﬁcantly shorter average route length according to a one-tailed t-test with conﬁdence level α = 0.05. The dash indicates that no feasible solution isfound in 30 runs.
Instanceγ =0

γ = 32SS-ﬁxedmvspdp50a
mvspdp50b
mvspdp75a
mvspdp75b
mvspdp100a
mvspdp100b
mvspdp100c
mvspdp100d
mvspdp150a
mvspdp150b
mvspdp199a
mvspdp199b
mvspdp120a
mvspdp120b
mvspdp200
mvspdp240
mvspdp280
mvspdp320
mvspdp360
mvspdp400
mvspdp440
mvspdp480

SS-adaptive

SS-ﬁxedSS-adaptive

Ave.

Feas. (%)

Ave.

Feas. (%)

Time inc. (%)

Ave.

Feas. (%)

Ave.

Feas. (%)

Time inc. (%)

473
444
668
601
810
745
625
522
995
985
1169
1149
1372
800
–
–
–
–
–
–
–
24,690

73
87
83
97
80
90
93
100
87
83
50
87
37
90
0
0
0
0
0
0
0
3

482
447
667
590
806
747
614
519
1080
911
1232
1291
1565
767
–
–
–
–
–
–
–
22,260

87
97
83
93
80
93
90
100
77
73
43
57
57
100
0
0
0
0
0
0
0
3

−10.53
−4.17
0.00
−4.00
0.00
6.06
−20.69
14.75
2.13
44.68
1.19
21.18
17.39
1.32

407
373
510
448
570
611
503
461
673
679
797
837
793
625
6982
6943
9841
10,825
13,602
24,56911,913
15,179

97
100
87
87
93
100
90
100
100
93
97
90
93
100
97
63
90
53
10
23
93
93

407
373
509
458
574
609
505
458
685
679
801
829
784
605
6758
6433
8730
9713
16,78024,18010,976
13,348

90
100
90
97
93
100
90
100
97
97
100
100
100
100
97
63
80
47
3
30
93
90

5.88
4.76
4.55
0.00
0.00
0.00
35.56
45.45
−4.65
−4.65
0.00
5.17
32.35
24.49
−9.52
9.32
−6.92
−3.59
6.45
−0.31
−10.09−19.59−25.05during the search process. The lower bound of the probability is set at 0.05 to prevent stagnation, and the probability foreach operator is 0.25 during the frozen time in TS. As for GA, after an offspring is generated, we adaptively adjust the probability by multiplying the mutation rate by the ratio of its ﬁtness to the less ﬁt parent. Additionally, the mutation rates ofrelocation and inversion operators are limited within [0.3, 0.9], and the rates for the selection and addition/deletion operators are [1/2|V + |, 2/|V + |] and [1/2m, 2/m], respectively. Tables 5 and 6 summarize the performance of TS, GA, and SS withthe adaptive probability of applying each operator on the MVSPDP instances with γ = {0, 32}.
Fig. 10 compares the route length, feasible rate, and average running time of the three metaheuristic algorithms. Similar to the tendency observed regarding ﬁxed probability, TS signiﬁcantly outperforms GA and SS in 14 and 17 out of 22
instances for γ = 0 and 32, respectively. Moreover, the experimental results indicate the considerable advantages of TS overGA and SS in running time and feasible rate for solving the MVSPDP. Based on the superior performance, the TS with adaptive probability is a favorable method for the MVSPDP. Considering the search behavior of metaheuristic algorithms, thelocal search methods like TS are good at locating the optima in a region, whereas the global search methods such as GA
excel in identifying the promising regions for global optima. In this regard, the preferable performance of TS implies the advantage of local search methods to solve the MVSPDP. In addition, the comprehensive survey of Cordeau et al. [19] indicatedthat a considerable number of tabu search algorithms have been proposed for the VRP, while the use of adaptive memory
procedure initiated the study of population search heuristics for the VRP [10]. Baker and Ayechew [3] stated that, although
GA performs well, it is inferior to tabu search in terms of solution quality on the VRP. As the proposed MVSPDP is a variantof the VRP, these empirical validations are referable for the MVSPDP; in fact, our experimental results are consistent withthese ﬁndings.
Next, we compare the anytime behavior of the proposed TS, GA, and SS in terms of route length against time. The
comparison further includes the mixed integer programming solver IBM ILOG CPLEX, which presents a performance baseline
for the MVSPDP. A termination criterion of six million nodes traversed is set for the CPLEX due to long running time. Figs. 11
and 12 depict the progress of route length against running time for the test algorithms on four MVSPDP instances with
γ = {0, 32}. The results show that TS achieves the fastest convergence, and SS converges faster than GA does. In addition,
the three proposed methods, i.e., TS, SS, and GA, are much more eﬃcient than CPLEX. The ﬁgures show the closeness ofthe obtained route lengths to the gradually increasing lower bound of CPLEX. However, the long running time of CPLEX(> 400 h) hinders the experiments from getting the optimal solution and lower bound. Tables 5 and 6 compare the route
length obtained from the four test algorithms. According to the results, CPLEX can mostly gain shorter routes than the three
proposed algorithms on the small instances; nonetheless, its running time is at least 106 times longer than TS and SS, and
105 times longer than GA. By contrast, the three proposed algorithms, especially the TS, can yield competitive results inreasonable time.

166

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

Fig. 13. Anytime behavior of the average route length over 30 runs of TS and SS using the ﬁxed probability (f) and adaptive probability (a) of applying thefour operators on mvspdp199a with γ = 32.

Fig. 14. The average route lengths obtained from TS, GA, and SS using the adaptive probability of applying the four operators at different gain values γ .

Tables 7–9 compare the results of ﬁxed and adaptive probability to investigate their effects on the three algorithms.
According to Table 7, the GA with ﬁxed probability surpasses that with adaptive probability in route length, but is inferior
in feasible rate. By contrast, Tables 8 and 9 show that adaptively tuning the probability of executing operations for TS andSS can signiﬁcantly shorten the route length on instances involving substantial numbers of nodes and high gains. Notably,
TS with adaptive probability is computationally less time-consuming than that with ﬁxed probability in most cases. BothTS and SS with adaptive probability achieve similar feasible rates to their corresponding versions with ﬁxed probability. In
addition, the anytime behavior of TS and SS, shown in Fig. 13, demonstrates that adaptive probability effectively improves
the convergence speed of TS. These results validate the beneﬁts of employing adaptive probability in TS and SS, in which
the adaptive control mechanism facilitates selecting promising operators in the course of search process.
Finally, this study experiments with different gain values to examine the inﬂuence of the selectability of pickup nodes.
Fig. 14 presents only the results on mvspdp199a and mvspdp360 due to space limitation. The results indicate that GA and
SS infrequently ﬁnd feasible solutions on mvspdp360 when γ is low, whereas the proposed TS achieves the shortest routes
and yields the lowest standard deviation among the test algorithms. The route length generally decreases as γ increases,
conﬁrming that the ﬂexibility in selecting pickup nodes intensiﬁes when the commodities outnumber those required. The
selectability of pickup nodes in the MVSPDP expands the search space of VRP. In other words, the sparsity of feasible solutions and the low gain γ raise the diﬃculty of searching for the optima; however, as γ increases, shorter routes could beachieved by selecting the required pickup nodes.
Fig. 15 illustrates the routes obtained from TS on mvspdp199a with γ = {0, 32}. The resultant routes validate the beneﬁts of the MVSPDP: the route length and the number of selected pickup nodes decrease as γ increases. Additionally, theselectability reduces the number of vehicles used (Fig. 15b). These results indicate that relaxing the constraint on visiting
all nodes provides economic solutions for real-world applications that focus on supplying all of the customers demanding

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

167

Fig. 15. The route obtained from TS on mvspdp199a with γ = {0, 32}. The square denotes the depot, solid circles denote pickup nodes, and hollow circles
denote delivery nodes.

commodities. The solutions may contain “crosses” in routes, which seem to be undesirable in the sense of shortest path.
Nevertheless, the crosses are necessary for the optimal solutions in the SPDP and the MVSPDP because of the constraints
on vehicle capacity and nonnegative vehicle loads.

5. Conclusions
This study presents the multi-vehicle selective pickup and delivery problem (MVSPDP), which aims for the shortest routes
for a ﬂeet of vehicles to collect and supply commodities in accordance with constraints regarding vehicle capacity and travel
distance. The problem formulation relaxes the requirement for visiting all nodes and uses multiple vehicles. The MVSPDP is
pertinent to real-world logistic applications that focus on supplying the demands of all customers from a certain number ofproviders by using multiple vehicles.
To resolve the MVSPDP, this study proposed three metaheuristic algorithms: TS, GA, and SS. The ﬁxed-length representation enables indicating the varying number of vehicles and the selection of pickup nodes. Four neighborhood functions are
introduced to TS, each of which focuses on one critical aspect of solving the MVSPDP. In addition, these functions are usedas mutation operators in GA and as improvement methods for SS. Moreover, we proposed an initialization method based on
the sweep algorithm for the MVSPDP.
A series of experiments is conducted to examine the proposed TS, GA, and SS. The experimental results indicate theutility of the four operators in selecting pickup nodes and arranging the visiting order for the MVSPDP, which leads toshorter routes than does the PDP. In addition, the results validate that the three metaheuristic algorithms can effectively
solve the MVSPDP; in particular, TS outperforms both GA and SS in route length and feasible rate on most of the instances
at various gain values. Moreover, TS requires less running time to achieve the solutions than the other two algorithms do.
According to these outcomes, TS with the adaptive probability of applying the four operators is a preferable method for
resolving the MVSPDP.
Future studies include some directions. First, this study focuses on design of TS, GA, and SS; nonetheless, other established algorithms, such as ant colony optimization (ACO) and particle swarm optimization (PSO), are also promising for the
MVSPDP. Improving the selection of appropriate pickup nodes and arrangement of routes can also be considered. Second,
the MVSPDP can be extended to various transportation scenarios and requirements, such as partial pickup and delivery,time windows, and multiple objectives. Partial pickup and delivery enables each customer to be served by multiple vehicles.
The MVSPDP with time windows requires customers to be served at their preferred periods. The multi-objective MVSPDP
involves additional objectives, such as the number of vehicles and waiting time. In addition, the costs of unselecting pickup
customers and adding vehicles can be further taken into account. Third, formal analysis on the bound and behavior ofmetaheuristic algorithms for the MVSPDP is a challenging yet fruitful topic for future work.

Acknowledgment
The authors would like to thank the editor and reviewers for their valuable comments and suggestions. This work was
supported by the National Science Council of Taiwan, under contract NSC10 0-2628-E-194-0 09-MY3.

168

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

References[1] S. Anily, J. Bramel, Approximation algorithms for the capacitated traveling salesman problem with pickups and deliveries, Nav. Res. Logist. 46 (6)(1999) 654–670.
[2] T. Back, D.B. Fogel, Z. Michalewicz (Eds.), Handbook of Evolutionary Computation, IOP Publishing Ltd, 1997.
[3] B.M. Baker, M.A. Ayechew, A genetic algorithm for the vehicle routing problem, Comput. Oper. Res. 30 (5) (2003) 787–800.
[4] R. Baldacci, E. Bartolini, A. Mingozzi, An exact algorithm for the pickup and delivery problem with time windows, Oper. Res. 59 (2) (2011) 414–426.
[5] M. Battarra, G. Erdoǧan, G. Laporte, D. Vigo, The traveling salesman problem with pickups, deliveries, and handling costs, Transp. Sci. 44 (3) (2010)
383–399.
[6] A.K. Beheshti, S.R. Hejazi, A novel hybrid column generation-metaheuristic approach for the vehicle routing problem with general soft time window,
Inf. Sci. 316 (2015) 598–615.
[7] G. Berbeglia, J.-F. Cordeau, I. Gribkovskaia, G. Laporte, Static pickup and delivery problems: A classiﬁcation scheme and survey, TOP 15 (1) (2007) 1–31.
[8] C. Bordenave, M. Gendreau, G. Laporte, Heuristics for the mixed swapping problem, Comput. Oper. Res. 37 (1) (2010) 108–114.
[9] C. Bordenave, M. Gendreau, G. Laporte, A branch-and-cut algorithm for the preemptive swapping problem, Networks 59 (4) (2012) 387–399.
[10] B. Bozkaya, E. Erkut, G. Laporte, A tabu search heuristic and adaptive memory procedure for political districting, Eur. J. Oper. Res. 144 (1) (2003) 12–26.
[11] Capacitated VRP Instances, http://neo.lcc.uma.es/vrp/vrp- instances/capacitated- vrp- instances/, 2015.
[12] H.K. Chen, H.W. Chou, C.F. Hsueh, T.Y. Ho, The linehaul-feeder vehicle routing problem with virtual depots, IEEE Trans. Autom. Sci. Eng. 8 (4) (2011)694–704.
[13] P. Chen, H.K. Huang, X.Y. Dong, An ant colony system based heuristic algorithm for the vehicle routing problem with simultaneous delivery and pickup,
in: Proceedings of the 2007 IEEE Conference on Industrial Electronics and Applications (ICIEA 2007), 2007, pp. 136–141.
[14] X. Chen, Y.S. Ong, A conceptual modeling of meme complexes in stochastic search, IEEE Trans. Syst. Man Cybern. – Part C 42 (5) (2012) 612–625.
[15] R. Chiong, T. Weise, Z. Michalewicz, Variants of Evolutionary Algorithms for Real-World Applications, Springer, 2012.
[16] N. Christoﬁdes, A. Mingozzi, P. Toth, Exact algorithms for the vehicle routing problem, based on spanning tree and shortest path relaxations, Math.
Program. 20 (1) (1981) 255–282.
[17] G. Clarke, J.W. Wright, Scheduling of vehicles from a central depot to a number of delivery points, Oper. Res. 12 (4) (1964) 568–581.
[18] J.-F. Cordeau, G. Laporte, S. Ropke, Recent models and algorithms for one-to-one pickup and delivery problems, in: B. Golden, S. Raghavan, E. Wasil(Eds.), The Vehicle Routing Problem: Latest Advances and New Challenges, Operations Research/Computer Science Interfaces Series, vol.43, Springer
US, 2008, pp. 327–357.
[19] J.-F. Cordeau, G. Laporte, M.W. Savelsbergh, D. Vigo, Vehicle routing, in: Handbooks in Operations Research and Management Science, vol.14, Elsevier,2007, pp. 367–428.
[20] C.E. Cortés, M. Matamala, C. Contardo, The pickup and delivery problem with transfers: Formulation and a branch-and-cut solution method, Eur. J.
Oper. Res. 200 (3) (2010) 711–724.
[21] L. Davis, Applying adaptive algorithms to epistatic domains., in: Proceedings of the International Joint Conference on Artiﬁcial Intelligence, vol.85,
1985, pp. 162–164.
[22] A.E. Eiben, J.E. Smith, Introduction to Evolutionary Computing, Natural Computing, Springer-Verlag, 2003.
[23] B. Eksioglu, A.V. Vural, A. Reisman, The vehicle routing problem: A taxonomic review, Comput. Ind. Eng. 57 (4) (2009) 1472–1483.
[24] G. Erdoğan, J.-F. Cordeau, G. Laporte, A branch-and-cut algorithm for solving the non-preemptive capacitated swapping problem, Discrete Appl. Math.
158 (15) (2010) 1599–1614.
[25] R. Falcon, X. Li, A. Nayak, I. Stojmenovic, The one-commodity traveling salesman problem with selective pickup and delivery: An ant colony approach,in: Proceedings of the 2010 IEEE Congress on Evolutionary Computation (CEC 2010), 2010, pp. 1–8.
[26] T. Garaix, C. Artigues, D. Feillet, D. Josselin, Optimization of occupancy rate in dial-a-ride problems via linear fractional column generation, Comput.
Oper. Res. 38 (10) (2011) 1435–1442.
[27] B.E. Gillett, L.R. Miller, A heuristic algorithm for the vehicle-dispatch problem, Oper. Res. 22 (2) (1974) 340–349.
[28] F. Glover, Heuristics for integer programming using surrogate constraints, Dec. Sci. 8 (1) (1977) 156–166.
[29] F. Glover, M. Laguna, Tabu Search, Kluwer Academic Publishers, 1997.
[30] F. Glover, M. Laguna, R. Martí, Fundamentals of scatter search and path relinking, Control cybern. 39 (3) (20 0 0) 653–684.
[31] F.P. Goksal, I. Karaoglan, F. Altiparmak, A hybrid discrete particle swarm optimization for vehicle routing problem with simultaneous pickup anddelivery, Comput. Ind. Eng. 65 (1) (2013) 39–53.
[32] Y.-J. Gong, J. Zhang, O. Liu, R.-Z. Huang, H.-H. Chung, Y.-H. Shi, Optimizing the vehicle routing problem with time windows: A discrete particle swarm
optimization approach, IEEE Trans. Syst. Man Cybern. Part C 42 (2) (2012) 254–267.
[33] I. Gribkovskaia, G. Laporte, One-to-many-to-one single vehicle pickup and delivery problems, in: B. Golden, S. Raghavan, E. Wasil (Eds.), The Vehicle
Routing Problem: Latest Advances and New Challenges, Operations Research/Computer Science Interfaces, vol.43, Springer US, 2008, pp. 359–377.
[34] I. Gribkovskaia, G. Laporte, A. Shyshou, The single vehicle routing problem with deliveries and selective pickups, Comput. Oper. Res. 35 (9) (2008)
2908–2924.
[35] G. Gutiérrez-Jarpa, G. Desaulniers, G. Laporte, V. Marianov, A branch-and-price algorithm for the vehicle routing problem with deliveries, selective
pickups and time windows, Eur. J. Oper. Res. 206 (2) (2010) 341–349.
[36] K. Haridass, J. Valenzuela, A.D. Yucekaya, T. McDonald, Scheduling a log transport system using simulated annealing, Inf. Sci. 264 (2014) 302–316.
[37] H. Hernández-Pérez, J.J. Salazar-González, The one-commodity pickup-and-delivery travelling salesman problem, in: M. Jünger, G. Reinelt, G. Rinaldi(Eds.), Combinatorial Optimization — Eureka, You Shrink!, Lecture Notes in Computer Science, vol.2570, Springer Berlin Heidelberg, 2003, pp. 89–104.
[38] A. Hoff, I. Gribkovskaia, G. Laporte, A. Løkketangen, Lasso solution strategies for the vehicle routing problem with pickups and deliveries, Eur. J. Oper.
Res. 192 (3) (2009) 755–766.
[39] J. Holland, Adaptation in Natural and Artiﬁcial Systems., University of Michigan Press, 1975.
[40] H.C.W. Lau, T.M. Chan, W.T. Tsui, W.K. Pang, Application of genetic algorithms to solve the multidepot vehicle routing problem, IEEE Trans. Autom. Sci.
Eng. 7 (2) (2010) 383–392.
[41] S. Lin, B.W. Kernighan, An effective heuristic algorithm for the traveling-salesman problem, Oper. Res. 21 (2) (1973) 498–516.
[42] J. Luo, X. Li, M.-R. Chen, H. Liu, A novel hybrid shuﬄed frog leaping algorithm for vehicle routing problem with time windows, Inf. Sci. 316 (2015)
266–292.
[43] G. Martinovic, I. Aleksi, A. Baumgartner, Single-commodity vehicle routing problem with pickup and delivery service, Math. Problems Eng. 2008
(2008), doi:10.1155/2008/697981. Article ID 697981.
[44] M. Mavrovouniotis, S. Yang, Ant algorithms with immigrants schemes for the dynamic vehicle routing problem, Inf. Sci. 294 (2015) 456–477.
[45] Z. Michalewicz, Genetic Algorithms+ Data Structures= Evolution Programs, Springer, 1996.
[46] Z. Michalewicz, D.B. Fogel, How to Solve It: Modern Heuristics, Springer, 2004.
[47] Y. Nagata, S. Kobayashi, Guided ejection search for the pickup and delivery problem with time windows, in: P. Cowling, P. Merz (Eds.), Evolutionary
Computation in Combinatorial Optimization, Lecture Notes in Computer Science, vol.6022, Springer Berlin Heidelberg, 2010, pp. 202–213.
[48] D. Naso, B. Turchiano, Multicriteria meta-heuristics for AGV dispatching control based on computational intelligence, IEEE Trans. Syst. Man Cybern.
Part B 35 (2) (2005) 208–226.
[49] I.H. Osman, J.P. Kelly, Meta-heuristics: Theory and applications, Springer, 1996.
[50] S.N. Parragh, K.F. Doerner, R.F. Hartl, A survey on pickup and delivery problems—Part I: Transportation between customers and depot, Journal für
Betriebswirtschaft 58 (2008) 21–51.

C.-K. Ting et al. / Information Sciences 406–407 (2017) 146–169

169[51] S.N. Parragh, K.F. Doerner, R.F. Hartl, A survey on pickup and delivery problems—Part II: Transportation between pickup and delivery locations, Journal
für Betriebswirtschaft 58 (2008) 81–117.
[52] S.N. Parragh, K.F. Doerner, R.F. Hartl, X. Gandibleux, A heuristic two-phase solution approach for the multi-objective dial-a-ride problem, Networks 54(4) (2009) 227–242.
[53] H.L. Petersen, S. Ropke, The pickup and delivery problem with cross-docking opportunity, in: J.W. Böse, H. Hu, C. Jahn, X. Shi, R. Stahlbock, S. Voß(Eds.), Computational Logistics, Lecture Notes in Computer Science, vol.6971, Springer Berlin Heidelberg, 2011, pp. 101–113.
[54] A.K. Qin, V.L. Huang, P.N. Suganthan, Differential evolution algorithm with strategy adaptation for global numerical optimization, IEEE Trans. Evol.
Comput. 13 (2) (2009) 398–417.
[55] P. Sasikumar, G. Kannan, Issues in reverse supply chains, part II: Reverse distribution issues — an overview, Int. J. Sustainable Eng. 1 (4) (2008)
234–249.
[56] P. Sasikumar, G. Kannan, Issues in reverse supply chain, part III: Classiﬁcation and simple analysis, Int. J. Sustainable Eng. 2 (1) (2009) 2–27.
[57] K.C. Tan, C.Y. Cheong, C.K. Goh, Solving multiobjective vehicle routing problem with stochastic demand via evolutionary computation, Eur. J. Oper. Res.
177 (2) (2007) 813–839.
[58] K.C. Tan, Y.H. Chew, L.H. Lee, A hybrid multi-objective evolutionary algorithm for solving truck and trailer vehicle routing porblems, Eur. J. Oper. Res.
172 (3) (2006) 855–885.
[59] K.C. Tan, Y.H. Chew, L.H. Lee, A hybrid multiobjective evolutionary algorithm for solving vehicle routing problem with time windows, Comput. Optim.
Appl. 34 (1) (2006) 115–151.
[60] E. Teymourian, V. Kayvanfar, G.M. Komaki, M. Zandieh, Enhanced intelligent water drops and cuckoo search algorithms for solving the capacitatedvehicle routing problem, Inf. Sci. (2015).
[61] C.K. Ting, X.L. Liao, The selective pickup and delivery problem: Formulation and a memetic algorithm, Int. J. Prod. Econ. 141 (1) (2013) 199–211.
[62] P. Toth, D. Vigo, VRP with backhauls, in: P. Toth, D. Vigo (Eds.), The vehicle routing problem, SIAM Monographs on Discrete Mathematics and Applications 9, Society for Industrial and Applied Mathematics, 2002, pp. 195–224.
[63] N. Velasco, P. Dejax, C. Guéret, C. Prins, A non-dominated sorting genetic algorithm for a bi-objective pick-up and delivery problem, Eng. Optim. 44(3) (2012) 305–325.
[64] H.A. Waisanen, D. Shah, M.A. Dahleh, A dynamic pickup and delivery problem in mobile networks under information constraints, IEEE Trans. Automat.
Control 53 (6) (2008) 1419–1433.
[65] X.P. Wang, H.Y. Cao, A dynamic vehicle routing problem with backhaul and time window, in: Proceedings of the 2008 IEEE International Conference
on Service Operations and Logistics, and Informatics (SOLI 2008), 1, 2008, pp. 1256–1261.
[66] E.T. Yassen, M. Ayob, M.Z.A. Nazri, N.R. Sabar, Meta-harmony search algorithm for the vehicle routing problem with time windows, Inf. Sci. 325 (2015)
140–158.
[67] Z. Zhu, J. Xiao, S. He, Z. Ji, Y. Sun, A multi-objective memetic algorithm based on locality-sensitive hashing for one-to-many-to-one dynamic pickup-and-delivery problem, Inf. Sci. 329 (2016) 73–89.