Scatter search for the stochastic travel-time vehicle routing problem with simultaneous pick-ups and deliveries

Computers & Operations Research 39 (2012) 2277–2290

Contents lists available at SciVerse ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

Scatter search for the stochastic travel-time vehicle routing problem with simultaneous pick-ups and deliveries Tao Zhang a, W.A. Chaovalitwongse b, Yuejie Zhang c,n a

School of Information Management and Engineering, Shanghai University of Finance and Economics, Shanghai 200433, PR China Departments of Industrial and Systems Engineering and Radiology, University of Washington, USA c School of Computer Science, Key Laboratory of Intelligent Information Processing, Fudan University, Shanghai 200433, PR China b

a r t i c l e i n f o

abstract

Available online 9 December 2011

In parallel with the growth of both domestic and international economies, there have been substantial efforts in making manufacturing and service industries more environmental friendly (i.e., promotion of environmental protection). Today manufacturers have become much more concerned with coordinating the operations of manufacturing (for new products) and recycling (for reuse of resources) together with scheduling the forward/reverse flows of goods over a supply chain network. The stochastic traveltime vehicle routing problem with simultaneous pick-ups and deliveries (STT-VRPSPD) is one of the major operations problems in bi-directional supply chain research. The STT-VRPSPD is a very challenging and difficult combinatorial optimization problem due to many reasons such as a nonmonotonic increase or decrease of vehicle capacity and the stochasticity of travel times. In this paper, we develop a new scatter search (SS) approach for the STT-VRPSPD by incorporating a new chanceconstrained programming method. A generic genetic algorithm (GA) approach for STT-VRPSPD is also developed and used as a reference for performance comparison. The Dethloff data will be used to evaluate the performance characteristics of both SS and GA approaches. The computational results suggest that the SS solutions are superior to the GA solutions. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Stochastic travel-time vehicle routing problem (STT-VRP) Integer programming VRP with simultaneous pick-ups and deliveries (VRPSPD) Genetic algorithm (GA) Scatter search (SS)

1. Introduction The STT-VRPSPD studied in this paper is motivated by an emerging real life application. With the realization of domestic and international environmental protection policies developed in recent years, enterprises have made every conscious effort to recycle their used products. Thus, nowadays manufacturers need to be concerned with coordinating the operations of manufacturing (for new products) and recycling (for reuse of resources) together with scheduling the forward and reverse flows of goods over a supply chain network. The forward flow delivers new products to customers, while the reverse flow picks up recyclable materials (e.g., empty soda bottles, waste papers, or unsold periodicals) as well as obsolete products that have to be processed for environmentally safe disposal (e.g., expiring pharmaceutical and chemicals products, or contaminated materials). To make this operations economically sound, logistics service providers have to consider simultaneous pick-up and delivery services in delivery/recycling systems to best utilize the vehicle space and reduce the delivery and recycling costs. In addition to

n

Corresponding author. Tel.: þ86 13 661473292; fax: þ 86 21 65904784. E-mail address: [email protected] (Y. Zhang).

0305-0548/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2011.11.021

the challenge of forward and reverse logistics, the transportation and delivery of commodities are vulnerable to changes in traffic density and travel speed that may be caused by accidents and weather changes, especially in large metropolitan areas. This, in turn, suggests that to solve such real life, challenging problems, the travel times should be modeled as stochastic rather than deterministic variables. Although the STT-VRPSPD considered in this paper is very practical, it is a very challenging and difficult problem. For instance, because customers require simultaneous pick-up and delivery service, the vehicle load capacity varies greatly and frequently as opposed to a monotonic decrease (delivery) or a monotonic increase (pick-up) of the capacity in the traditional vehicle routing problem with backhauls (VRPB), which is a special case of VRPSPD. Specifically in the VRPB, there are separate customer groups, ‘‘delivery only’’ and ‘‘pick-up only’’ while the VRPSPD must deal with customers with both delivery and pick-up demands simultaneously. Another difficulty is the limitation in the vehicle’s maximum travel time. Due to the stochasticity of travel times, when the vehicle’s travel time is close to the maximum, the vehicle may have to travel back to the depot even if its remaining load capacity can satisfy the demands of unvisited customers. These practical challenges impose several hard constraints to the problem, which will increase the number of service vehicles and the total travel time.

2278

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

In order to efficiently find practical solutions to the STTVRPSPD, we develop a chance-constrained programming model to transform the dynamic problem to a static problem. In addition, we construct a scatter search (SS) framework for the STTVRPSPD. Among various metaheuristics algorithms, SS is very appropriate for the STT-VRPSPD because of its flexibility and effectiveness in generating diverse solutions to combinatorial problems. This is even more important for the stochastic problem considered here. SS has also been shown to be integratable with the tabu search, and it has been proven more effective than genetic algorithm (GA). New contributions of our SS framework includes constructing new solution improvement and combination strategies as well as improving the Clarke and Wright (C–W) algorithm given in [63] to more effectively construct initial good solutions. These efforts make our SS framework suitable for the STT-VRPSPD as it can effectively deal with the stochasticity of travel times and the fluctuating capacity loads. The rest of this paper is organized as follows. Section 2 discusses the background and previous studies related to the STT-VRPSPD and the SS algorithm. Section 3 provides the problem description of the STT-VRPSPD chance-constrained programming model and its mathematical formulation. In Section 4, we explain the logic and procedure of our SS framework as well as provide a basic GA algorithm for STT-VRPSPD that we use as a benchmark algorithm. Section 5 discusses how the test instances are generated and provide computational results. Section 6 concludes our study and discusses future research directions.

2. Background and related work 2.1. Vehicle routing problem with simultaneous pick-up and delivery (VRPSPD) The VRPSPD is a variant of the capacitated vehicle routing problem where the vehicles are not only required to deliver goods but also to pick up some goods from the customers. VRPB and VRPSPD are among the reverse logistics classes of the traditional VRP, where the VRPB has been shown to be a special case of VRPSPD [1]. Both VRPB and VRPSPD have been proven to be NP-hard [1,2]. The VRPSPD was first introduced by Min in 1989 [3] for a book-distribution problem between 22 public libraries (customers) and a library administration center (depot) with a fixed number of service vehicles and a limited load capacity of homogeneous vehicles. The VRPSPD was then solved by using clustering and sorting algorithms, and subsequently optimizing the customer sequence in each cluster based on a standard traveling salesman problem (TSP). Over the years, there have been a number of studies on the VRPSPD. Halse [4] studied a variant of VRPSPD as well as many other special case problems in the VRP literature. Cases with a single depot and multiple vehicles and number of nodes varying between 22 and 150 are addressed. Lagrangian relaxation and a column generation approach is utilized. A cluster-first route-second algorithm is developed in which nodes are first distributed to vehicles and then the problem is solved using 3-opt algorithm. Fisher et al. [5] described a new algorithm based on a network flow relaxation which imposes necessary conditions on the flow of empty vehicles from order delivery points to order pick-up points. The new algorithm is fast and has performed well on a set of more than 430 test problems, which include a number of real problems obtained from the Shanghai Truck Transportation Corporation. Savelsbergh and Sol [6] studied a general pick-up and delivery problem and divided the class of pick-up and delivery problems into static and dynamic problems. Within either of these classes, the author distinguished single-vehicle and multiple-vehicle problems. The

paper has given an overview of the solution approaches to the various characteristics of pick-up and delivery problems. Gendreau et al. [7] studied a single vehicle VRPSPD to solve the TSP at first, and then arranged the order of pick-up and delivery on the TSP routing. In recent years, Dethloff did some deep research on this area. Dethloff [8] addressed the VRPSPD based on the point of view of reverse logistics, where he proposed a mathematical model and developed a solution algorithm that was based on an insertion criterion using the vehicle service degrees of freedom. He developed 40 instances to test his algorithm and compared his results with those of Salhi and Nagy [9] and reported an improvement on Min’s problem. Subsequently, Dethloff [1] discussed the relationship between VRPSPD and VRPB. The possibility of solving the VRPB by applying an insertion heuristic based on the concept of ‘residual capacities’ originally designed for the VRPSPD is investigated. Numerical results indicate that, for certain instances, this approach is more favorable than the application of a heuristic suggested for the VRPB in the literature. Tang and Galvao ¯ [10] proposed a mathematical model for the VRPSPD, and solve this problem through a tour partitioning heuristic developed for the TSP with simultaneous pick-up and delivery service. Subsequently, Tang and Galvao ¯ [11] constructed a mathematical model for the VRPSPD with a maximum distance constraint, and solved it by using a tabu search and mixed local search algorithms. Gronalt et al. [12] presented four different savings based heuristics for the VRPSPD of full truckloads under time window constraints to minimize empty vehicle movements. Salhi and Nagy [13] studied heuristic algorithms for a single and multiple depot vehicle routing problem with pick-ups and deliveries (VRPPD). The paper eschews the concept of insertion and proposes a method that treats pick-ups and deliveries in an integrated manner. They first found a solution to the VRP by allowing infeasibilities and then modified this solution to make it feasible for the single and multiple-depot VRPPD. New approaches for solving the VRPSPD are studied by a great number of researchers in recent years. Zachariadis et al. [14] proposed a hybrid solution approach incorporating the rationale of two well-known meta-heuristics which have been proven to be effective for routing problem variants, namely tabu search and guided local search. The intelligence of the proposed hybrid algorithm was designed to achieve a vast exploration of the search space, by escaping from local optima and intensifying at promising solution regions. Bianchessi and Righini [15] presented and compared constructive algorithms, local search algorithms, and tabu search algorithms, and address the VRPSPD issue of applying the tabu search paradigm to algorithms based on complex and variable neighborhoods. Arc-exchange-based and node-exchange-based neighborhoods were combined. Ai and Kachitvichyanukul [16] applied a real-value version of particle swarm optimization (PSO) algorithm for solving the VRPSPD. Gajpal and Abad [17] proposed an ant colony system (ACS) algorithm that employs a construction rule as well as two multi-route local search schemes to solve the VRPSPD problem. Subramanian et al. [18] presented a parallel approach for solving the VRPSPD problem. The parallel algorithm is embedded with a multi-start heuristic which consists of a variable neighborhood descent procedure, with a random neighborhood ordering, integrated in an iterated local search (ILS) framework. The experiments were performed in a cluster with a multi-core architecture using up to 256 cores. C - atay [19] proposed an ant colony algorithm employing a new saving-based visibility function and pheromone updating procedure. Despite recent advances made by the above mentioned studies, the VRPSPD remains to be a very challenging combinatorial optimization problem.

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

2.2. Vehicle routing problem with stochastic travel times Research on the stochastic vehicle routing problem (SVRP) can be divided into five categories: VRP with stochastic demands (VRPSD), VRP with stochastic customers (VRPSC), VRP with stochastic travel times (VRPSTT), VRP with stochastic serving times (VRPSST) and dynamic stochastic VRP (DSVRP). The VRPSTT is arguably one of the most challenging and practical variants of the VRP, in which the travel time between every node pair is a random variable related to traffic jam, road maintenance or weather conditions. Kao [20] proposed two solution approaches, based on dynamic programming and implicit enumeration, to solve the TSPSTT. In a subsequent study by Sniedovich [21], the dynamic programming approach was shown effective in obtaining close-to-optimal solutions when the problem is monotonic. In practice, it is extremely difficult to prove that a TSPSTT is monotonic. Nevertheless, Carraway et al. [22] later developed a generalized dynamic programming approach that overcame the limitations of Kao’s algorithm for TSPSTT. Laporte et al. [23] proposed a variant of VRPSTT chance-constrained programming model and compensation model, where a penalty function was introduced to pro-rate the delay time. By exploring the dual problem of VRPSTT, they used an L-Shaped algorithm to successfully solve the dual problem with 10, 15, and 20 vehicles, respectively. Park and Song [24] subsequently constructed three new heuristic algorithms based on an extension of VRP algorithms. Most recently, Xie [25] extended the VRPSTT chanceconstrained programming model by considering the same capacity of vehicles, and used a generation algorithm with some degree of success. Shao et al. [26] provided a hybrid particle swarm optimization algorithm to solve chance-constraint model with consideration of capacity. The above papers captured the stochastic travel times by chance-constraint programming and compensation model. Woensel et al. [27] created an innovative modeling scheme, which introduced the traffic congestion component based on queueing theory, to capture the stochastic behavior of travel times. In this model, the expected travel times were obtained from the distance and the speed in each time slice and variance of the travel time was calculated to evaluate routes on the risk involved. Based on Woensel et al. [27], Lecluyse et al. [28] introduced the variability in traffic flow into the model, which was used to evaluate the routes based on the uncertainty involved. At the same time, the objective function was extended by standard deviation of the travel times. Connors and Sumalee [29] and Chen and Zhou [30] studied the stochasticity of travel times from the view of travelers’ equilibrium. The former study considered a modeling framework where the perceived value and perceived probabilities of travel time outcomes are obtained via nonlinear transformations of the actual travel times and their probabilities. Based on maximizing the perceived values, the analysis proved that there was only one existence of equilibrium. The latter modeled a novel a-reliable mean-excess traffic equilibrium in which both reliability and unreliability aspects of travel time variability in the route choice decision process are considered. A route-based traffic assignment algorithm via the selfadaptive alternating direction method was used to solve this inequality problem. Although there have been a number of studies on VRPSTT in the literature, practical algorithms and solution approaches for such problem are still needed. This paper proposes a new SS framework to find very good solutions to the STT-VRPSPD effectively and efficiently. 2.3. Scatter search (SS) The SS framework was proposed by Glover in 1997 [31], and the first SS algorithm was subsequently developed by Glover in

2279

1998 [32]. The SS framework has been widely studied because it is very flexible and can be easily applied to many diverse optimization problems. For instance, some steps of the SS framework can be omitted or the procedures can be rearranged. In the recent years, the SS framework has been widely used to solve many complex combinatorial optimization problems, and it was shown to be more effective than other metaheuristics in many instances. Campos et al. [33] developed a SS approach to solve the linear ordering problem (LOP) by using the concept of adaptive memory to modify the greedy function in the diversification generation method. A nearest neighbor search algorithm, proposed by Laguna et al. [34], was then used in the improvement method, and the SS approach was shown to outperform the tabu search and a greedy algorithm. In a subsequent study, Campos et al. [35] conducted a more thorough analysis and obtained the same conclusion. Later, Campos et al. [36] extended the SS framework by designing a tabu rule for permutation problems. The SS framework was compared with the OptQuest and Frontline software packages, and the results suggested that the SS framework was more effective in solving the LOP but less effective in solving the TSP while the SS was the most efficient among the three approaches. Glover et al. [37] used the SS framework to solve nonlinear optimization problems, which were based on [38], and showed that the SS framework was much more effective and efficient than the GA. Habiba and Nabila [39] proposed a SS approach for the satisfiability problem and the weighted maximum satisfiability problem, and the results showed that the SS approach outperformed the greedy randomize adaptive search procedure (GRASP). Hung et al. [40,41] used a SS approach to solve the binary decision diagram minimum problem and showed that the SS approach outperformed the GA. In addition, Hung and Song [42] employed a SS approach to effectively solve several telecommunication problems. In addition to general single-objective optimization problems, the SS framework can be extended to solve multi-objective optimization problems, which often arise in a more complex environment. The key idea of SS in multi-objective optimization is to maintain a good search time to balance the solutions’ concentration and diversity so that the solutions on the Pareto border has certain dispersion. Nebro et al. [43] was the first to introduce the SS framework in the area of multi-objective optimization. Marti et al. [44] used a multi-objective integer programming problem as a test instance to compare the efficiency of SS with CPLEX 6.5, and demonstrated that the SS was 17 times faster than CPLEX but its solution was inferior to that of CPLEX. Subsequently, Laguna and Marti [45] designed numerical tests to investigate the SS for solving global optimization problems with multimodal function. It showed that the SS’ solution was better than that of GA with the same computational time, and the probability of obtaining optimal solutions was two times more than that of GA. Beausoleil [46] proposed a SS algorithm for solving nonlinear multi-objective problems, using a number of starting points of the tabu search method as the diversification generation method, and applying Cramer choice function to divide the reference set. Molina et al. [47] proposed a nonlinear multi-objective SS algorithm based on global estimation method, and compared the solutions with those of Beausoleil [46]. Other applications of the SS framework in multi-objective optimization can be found in [48,49]. The SS framework can also be seamlessly integrated with other evolutionary algorithms to form hybrid algorithms for complex combinatorial problems. Hybrid algorithm can make good use of the characteristics of different algorithms to achieve complementary advantages to improve the algorithm optimal performance and efficiency. Blazewicz et al. [50] proposed a framework integrating SS and the tabu search, in which the tabu table was employed to control the direction of solution movement to speed up the process.

2280

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

Trafalis and Kasap [51] proposed a hybrid algorithm of combined genetic algorithm, tabu search, and SS algorithms by embedding the genetic operator and the tabu table in the SS framework to build better solution sets to guide the search direction. Santana-Quintero et al. [52] employed the SS algorithm in a particle swarm optimization framework as a local search method to improve the concentration of particle swarm optimization algorithm. The SS framework has been widely used in several diverse fields, especially in logistics and supply chain. Russell and Chiang [53] developed a SS framework with a tabu-search-based improvement solution method to solve the VRPTW. Greistorfer [54] combined the SS and the tabu search to solve the routing arrangement problem. Keskin and Uster [55] used a path relinking (PR) subroutine as a subset of the SS to solve the middle point of the two-stage layout problem with capacity constraints. Alvarez et al. [56] employed the SS to solve a commodity network design problem. Gutierrez et al. [57] compared the SS with the RAND method in solving a multi-product inventory problem. In addition, the SS algorithm has also been used in production management problems. Gonzalez and Adenso-Diaz [58] used the SS to find the best split sequence after the end of complex product life cycles. Nowicki and Smutnicki [59] combined the SS with the PR to solve flow shop scheduling problems. Alvarez-Valdes et al. [60] used the SS to solve the project scheduling problem subject to the renewable resources and non-renewable resources constraints. Rahimi-Vahed et al. [61] used the SS to solve a mixed-model assembly line sequencing problem, and compared its performance with other evolutionary algorithms. Although there has been a body of work in applying the SS to routing and logistics problems, the application SS for STT-VRPSPD is still new. To our knowledge, this work is among the pioneering studies in applying the SS to solve the STT-VRPSPD. The STT-VRPSPD is closely related to the stochastic travel time multi-TSP (STT-MTSP) chance-constrained programming model proposed in Laporte et al. [23]; however, the vehicle load capacity constraint was not included in the STT-MTSP model. Subsequently, Xie [25] investigated the STT-MTSP chance-constrained programming model with the load capacity constraint and use the GA framework to solve the problem. Based on [25,62], this study extends the chance-constrained programming model and the load capacity constraint of the STT-MTSP to construct the STT-VRPSPD chance-constrained model with the vehicle load capacity constraint.

3. Problem description and mathematical formulation

(4) The stochastic travel time for each arc follows a normal distribution, which is known. (5) All vehicles are homogeneous, and their capacity is known. (6) Each client can only be visited once. In order to model this problem mathematically, we define the following notation: m C n V V0

maximum number of vehicles utilized; set of homogeneous vehicles, C ¼ f1,2,3,. . .,mg; number of customers; set of clients, V ¼ f1,2,3,. . .,ng; set of nodes, V0 ¼{0}[V, and 0 represents the depot; set of stochastic vectors of travel times on arcs; x’s sample space, which is finite sample space; travel time of vehicle l on arc (i,j) in the case of x; vehicle travel speed on arc (i,j); distance between nodes i and j, Dij ¼ vij t xijl ; travel cost of vehicle l traversing arc (i,j) – here we assume cijl ¼Dij ; unit cost of vehicle departure; upper limit of the total travel time of vehicle l in a day; delivery demand of customer i; pick-up demand of customer i; maximum load capacity of vehicle l; load of vehicle l after visiting customer i, before visiting customer j;

x X t xijl vij Dij cijl

l Bl di pi Ql qijl ( xijl ¼

1,

if vehicle l travels from node i to j in its tour

0,

otherwise

( zil ¼

1,

if vehicle l visits node i

0,

otherwise

(1) There is only one depot. (2) The locations of the depot and customer nodes are known. (3) The customers’ demands (both pick-up and delivery) are known.

iA V 0 ,

i,jA V 0 ,

lAC

lAC

The mathematical programming formulation of STT-VRPSPD is then given by min

m X n X n X

cijl xijl þ l

l¼1i¼0j¼0

s:t: m X n X xijl ¼ 1,

m X n X

x0jl

ð1Þ

l¼1j¼1

j ¼ 1,2,. . .,n

ð2Þ

l ¼ 1i ¼ 0 n X

xijl 

i¼0

The STT-VRPSPD can be formally defined as follows. Given a complete weighted digraph G ¼ ðV 0 ,AÞ, where V 0 ¼ f0,1,. . .,ng is a set of nodes, in which node 0 represents the depot, and V ¼ V 0 =f0g is the client node set, in which each and every node has simultaneous pick-up and delivery demands. The travel time on each arc in A is a stochastic variable with known distribution, which can be continuous or discrete. The vehicle load capacity is limited, and the total travel time of each vehicle is limited. Each vehicle must start its tour from the depot, transverse along the demand nodes (both pick-up and delivery) in the tour, and travel back to the depot to finish its tour. Each client must be served by only one vehicle. The objective of STT-VRPSPD is to find a feasible set of routes with the minimum total expected travel costs that satisfy all customers’ demands subject to the vehicle load capacity limitation and the maximum travel time for each vehicle. The STT-VRPSPD studied here is also subject to the following assumptions:

,

,

n X

n X

x0jl þ

j¼1 n X

xjil ¼ 0,

j ¼ 0,1,2,. . .,n,

l ¼ 1,2,. . .,m

ð3Þ

i¼0 n X

xj0l ¼ 2z0l ,

l ¼ 1,2,. . .,m

ð4Þ

j¼1

x0il r1,

l ¼ 1,2,. . .,m

ð5Þ

i¼1

0 P@

n X n X

1 x

t ijl rBl A Z a,

l ¼ 1,2,. . .,m

ð6Þ

i¼0j¼0

XX xijl r 9S91,

S D V,

l ¼ 1,2,. . .,m

ð7Þ

iAS jAS n X

q0jl ¼

j¼1 m X n X l¼1i¼0

n X n X

xijl dj ,

l ¼ 1,2,. . .,m

ð8Þ

qjil pj ,

ð9Þ

i¼0j¼1

qijl dj ¼

m X n X l¼1i¼0

j ¼ 0,1,2,. . .,n

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

0 r qijl rQ l ,

i ¼ 0,1,2,. . .n,

j ¼ 0,1,2,. . .n,

zil A f0,1g,

i ¼ 0,1,2,. . .n,

l ¼ 1,2,. . .,m

xijl A f0,1g,

i ¼ 0,1,2,. . .n,

j ¼ 0,1,2,. . .n,

l ¼ 1,2,. . .,m

ð10Þ

2281

procedure is illustrated in Fig. 1. Our approach can be outlined in the follow steps (for our pseudo-code, please see Appendix A).

ð11Þ l ¼ 1,2,. . .,m

ð12Þ

The objective function of STT-VRPSPD in Eq. (1) minimizes the total expected cost of all vehicle tours and the cost of vehicle departure. The constraints in Eq. (2) ensure that each customer is visited only once while the constraints in Eq. (3) are the flow preservation constraints; that is, vehicle l must depart from customer j after visiting customer j. The constraints in Eq. (4) set a limit that each vehicle must start and end its tour at the depot (node 0). Constraints (5) assure that each vehicle is only used once. We use a chance-constrained condition to ensure that the total travel time of vehicle l is less than the upper bound Bl with probability a in Formula (6). In order to avoid sub-tours, we use the oddity sub-tour elimination constraints in Eq. (7). The constraints in Eq. (8) enforce that the load capacity of each vehicle departing from the depot is equal to its total delivery demand. The load-preservation constraints in Eq. (9) express the proper movement of load delivery and/or pick-up. The capacity constraints in Eq. (10) ensure that the vehicle load at any given time does not exceed the vehicle capacity. The constraints in Eqs. (11) and (12) define the integral variables xijl and zil . We note that the a value in Eq. (6) can be viewed as a confidence level determined a priori. These chance constraints can be transformed into fixed constraints that can be solved by traditional methods for the nonstochastic version. The transformation can be described as follows. Recall that the vehicle’s travel time on each arc is assumed to be independent, and normally distributed, which can be represented by NðDij =vijl , s2ijl Þ, where vijl is the average travel speed of vehicle l on arc (i,j), and sijl is the standard deviation of vehicle’s travel time on arc (i,j). Thus we can define the slack Pn Pn travel time for vehicle l as t x x B , which follows Pni ¼ 0Pnj ¼ 0 ijl ijl l P P the normal distribution Nð i ¼ 0 j ¼ 0 Dij xijl =vijl Bl , ni¼ 0 nj¼ 0 2 2 sijl xijl Þ, and the chance constraints in Eq. (6) can be rewritten by 0 1 Pn Pn D x =v B ijl lC i¼0 j ¼ 0 ij ijl B ð13Þ P @Z r  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AZa Pn Pn 2 2 i¼0 j ¼ 0 sijl xijl It is very important to note that Eq. (13) exists if and only if the following expression holds: Pn Pn ¼0 j ¼ 0 Dij xijl =vijl Bl ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F1 ðaÞ r iq ð14Þ Pn Pn 2 2 i¼0 j ¼ 0 sijl xijl

4.1. Preprocessing: coding The initial step of our approach is to preprocess the parameter input of the STT-VRPSPD by encoding it to be compatible with the SS framework. We extend and apply the coding structure by using the integral number coding of all customer nodes in a full array with the depot as 0. For instance, when there are 9 customers, we use natural numbers 1–9 to represent the customer nodes where 0 denotes the depot node. If solution code ‘‘0123046780950’’ is a solution to the STT-VRPSPD, we can conclude that three vehicles are utilized – the first one visiting customers 1, 2, and 3, respectively, the second one visiting customers 4, 6, 7, and 8, respectively, and the third one visiting customers 9 and 5, respectively. 4.2. Create initial solutions After we encode the problem into the integral number coding, we extend the C–W algorithm to create initial solutions to the STT-VRPSPD. Specifically, we generalize the saving function in the traditional C–W algorithm to include not only the distance costs but also the travel times, pick-ups, and deliveries. Specifically, the saving function in the traditional C–W algorithm is defined by Sði,jÞ ¼ D0i þ Dj0 Dij , where i and j are customer nodes, and 0 is the depot. Our modified saving function when we add customer j in the route is given by S’ði,jÞ ¼ tð0,iÞ þ tðj,0Þtði,jÞ, ( Sn ði,jÞ ¼

S0 ði,jÞ þ bðdj pj Þ, 0

S ði,jÞgðpj dj Þ,

ð15Þ dj 4 pj , dj rpj ,

ð16Þ

where S0 ði,jÞ is the travel time saving, bUðdj pj Þ is the benefit of increased space from adding customer j because the delivery demand is larger than the pick-up demand; gUðpj dj Þ is the benefit of a decreased space after adding customer j so that customer j should be served toward the end of the route and the penalty cost gðpj dj Þ will be incurred. Here, let the values of greedy parameters b and g range between 0 and 1, i.e., b A ½0,1 and g A ½0,1. The use of these greedy parameters makes our approach more flexible and applicable to controlling the diversification of the solution. Specifically, to make our approach more greedy, we can increase the values of b and g. The flowchart of our modified C–W algorithm for STT-VRPSPD is illustrated in Fig. 2. The algorithm can be described as follows:

4. Scatter search framework for STT-VRPSPD

 The first step is to construct N routes for N customers, i.e., each From the previous section, STT-VRPSPD’s formulation is stochastic, and its combinatorial nature makes it extremely difficult to solve. In this work, we propose to employ a SS framework, first introduced in the late 1970s [31], to solve the STT-VRPSPD. The SS is a global search strategy based on evolutions. It first generates an initial population, selects high quality solutions, and some comparative scattered solutions as a reference set, uses ‘‘scatterconstringency’’ intelligent iterations to increase the quality and diversity of the solutions in the reference set, and integrate the solutions in the reference set with the current solution to obtain the best solutions. In short, there are five main components in the SS framework: diversification generation method, improvement method, reference set update method, subset generation method, and solution combination method. One of the main contributions of this work is to develop a solution approach based on the basic SS framework, and apply it to solve the STT-VRPSPD. The overall

route serves only one customer.

 The second step is to calculate the saving of grouping custo-





mers i and j together in the same route, Sn ði,jÞ. If the saving is positive, then we add ði,jÞ in the potential combinable set M ¼ fði,jÞ9Sn ði,jÞ 4 0g. The third step is to check whether set M is empty. If it is empty, we terminate the algorithm and return the current route list. Otherwise, we check if customers i and j satisfy the following conditions: i and j are not in the same formed routing path, i and j both are directly connected with the depot, one is the start node and another is the end node. If the conditions are satisfied, we proceed to the fourth step; otherwise, we go to the fifth step. The fourth step is to examine whether the routing is feasible (based on the vehicle load capacity and travel time limitation constraints) after connecting customers i and j according to the

2282

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

Create initial solution

Stop

Y Diversification Generation Method Improvement Method

Stop?

N

Initial solution set

Reference Set Update Method

Improvement Method

Refset Subset

Solution Combination Method Fig. 1. Flowchart of our scatter search framework.

Create initial routing

j

i

Calculate S* (i, j) and sorting

j

i

0

0

j Check each item in S* (i, j) orderly, judge whether the item satisfying connection condition.

m

i

j i

n

n

0 Y

0

i

i j

Connect i and j, update routing

m

j

N l

m

l n

N

m

n

M = M − {(i, j)} 0

0 Fig. 3. Connection strategies for nodes i and j.

Is set M empty?

Y Output initial solution Fig. 2. Flowchart of modified C–W algorithm for STT-VRPSPD.



direction of the paths. If it is feasible, then connect customers i and j in the route (connection strategies as shown in Fig. 3) and go to the fifth step. Otherwise, go to the fifth step directly. The fifth step is to update the set M by letting M ¼ Mfði,jÞg and go to the third step. Through this method, we create Psize initial solutions.

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

4.3. Diversification generation method (DGM) The DGM is used to exchange segments of the routing solutions to diversify the solutions. The proposed DGM will first reproduce this solution by randomly segmenting the solution and exchanging the solutions at the segments. For example, we assume that 0123046780950 is a solution generated by our modified C–W algorithm. We segment the solution as 012930467980950, where ‘‘9’’ indicates an exchange position. The resulting solution after applying DGM is 017930462980950 and it will be added to the initial solution population. We repeat the random selection, reproduction, and exchange steps until there are Psizesolutions in the initial solution population. 4.4. Improvement method (IM) It is very important to note that the DGM does not guarantee the solution’s feasibility after the exchange. Thus, the IM is a local search used to make infeasible solutions in the initial solution population feasible as well as improve the solution quality. 4.4.1. Generating feasible solutions The first part of IM is to make infeasible solutions feasible, i.e., satisfy all the constraints. We develop the following steps for the STT-VRPSPD and our genetic coding: Step 1: Check vehicle’s load capability. A) For vehicle l, if the initial load is larger than Q l (vehicle’s load limitation), then we calculate the cji =di ratio, which defines the cost-demand efficiency of customer i. Nodes i and j are visited by vehicle l, and node j denotes the immediate predecessor of customer i. We remove the customer with the largest ratio and repeat this step until the initial load capacity is less than or equal to Q l . For the customers removed during this step, we introduce new vehicles to serve them. Next we will go to Step B. B) If the initial load capacity is not larger than Q l , we examine if the vehicle will be overloaded after serving any customer in its route. We then remove the first customer in the route that will cause the first overload, and put the customer in a new route. We repeat this step until the vehicle will be not overloaded after serving any customer in its route. Next we go to Step 2. Step 2: Check vehicle’s travel times. If the time constraints in Eq. (8) are violated for any given route, we remove the customer with the largest value of t ji =di ratio, and put it into a new route. We repeat this step until there are no violations. 4.4.2. Solution improvement The second part of IM is to improve the quality of feasible solutions found in the previous subsection. We develop the following steps to improve both routing time and load capacity utilization. Step 1: Improving the travel times. First we identify the vehicle with the longest travel time. We then use the C–W rule to make an attempt to insert the customers of that route to other routes. Specifically, assuming customer i is a node between nodes h and s, which are currently served by the vehicle with the longest travel time. For any node pair j and k of another route, we then calculate tði,j,k,h,sÞ ¼ tðj,iÞ þtði,kÞ þ tðh,sÞtðh,iÞtði,sÞtðj,kÞ to quantify the savings in travel time. Specifically, if tði,j,k,h,sÞ Z0, we continue to select an insert position. On the other hand, if tði,j,k,h,sÞ o0, then we check if the insertion violates the time windows and load capacity constraints. If there is no violation,

2283

insert customer i between nodes j and k. We continue this procedure until there are no possible legitimate insertions and check other routes. Finally if there is no improvement on the travel times, then we go to Step 2; otherwise, we move to other feasible solutions. Step 2: Improving load capacity. After we are finished with Step 1, we attempt to improve the utilization of vehicle’s capacity. We first identify a vehicle with the maximum load capacity, defined by load capacity ¼d(deliver goods) þp(pick-up goods). Then we move the customers served by that vehicle to other vehicles while checking the feasibility of the vehicle capacity constraints. We repeat this step until there is no customers left or there is no other possible insertion, and go back to Step 1. 4.5. Reference set update method (RSUM) The RSUM is used to store and update high-quality (greedy) and diversified (semi-random) solutions. We form the reference set of size 9b9, and the set is consisted of two mutually exclusive subsets: a subset of high-quality solutions, called b1, and a subset of diversified solutions, called b2. The solutions in b are sorted by the solutions’ objective values in an ascending order. We choose the top 9b19 solutions to be in b1. The solutions in b2 are chosen by selecting the top 9b29 solutions with the largest diversity distances from bb1 . For instance, we can calculate the diversity distance as follows. Consider the solution 098745630210 in bb1 , which is not a high quality solution, there are 11 arcs (0,9), (9,8), (8,7), (7,4), (4,5), (5,6), (6,3), (3,0), (0,2), (2,1), (1,0) in this solution. Consider another solution 0123046780950 in b1, which is a high quality solution, there are 12 arcs (0,1), (1,2), (2,3), (3,0), (0,4), (4,6), (6,7), (7,8), (8,0), (0,9), (9,5), (5,0) in this solution. It appears that the arcs (3,0) and (0,9) are common in both solutions. Thus the diversity distance of the two solutions is maxð12,11Þ2¼10. Assume that we calculate the diversity distances of this non-high quality solution with two other high quality solutions and the diversity distances are 5 and 6. The diversity distance of this nonhigh quality solution is minð5,6,10Þ ¼ 5: 4.6. Subset generation method (SGM) There have been a number of studies on the SGM [37] to generate better solutions in the next solution combination. Here we use every pair of two solutions in the reference set to form a subset. We also limit the total number of subsets to be not more than the number of initial solutions, Psize. In the SGM, we can obtain C 29b9 subsets by combining every pair of two solutions from b. If C 29b9 is larger than Psize, then select Psizesubsets randomly. Subsequently we will use the solution combination method (in the next subsection) for every subset to generate a new solution. 4.7. Solution combination method (SCM) The SCM is a procedure to maintain characteristics of highquality solutions (e.g., optimal subtour in the routing) while promote the diversity in the solution set to overcome the situation when we are trapped in local optimal solutions. The SCM requires the use of the IM because combined solutions may or may not be feasible. In this paper, the SCM works as follows. First we randomly select segments (i.e., part of the solution) for solution exchange and calculate an average routing cost c/n before solution exchange. We keep the segment that has a lower average cost, exchange the solution segments and sort all customer nodes. Consider the following example. Assume we have two solutions, 0987945630210 and 01239046780950, where ‘‘9’’ represents a randomly selected exchange position. We then calculate average routing costs of segment 0987 and segment 0123. Assuming the

2284

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

average cost of segment 0987 is lower, the resulting exchanged solution is 09870123046050.

5. Generic genetic algorithm for STT-VRPSPD The genetic algorithm (GA), proposed by Holland [64], is among the most commonly used heuristic search algorithms. It uses a probability search method that is widely applicable to many real life problems. The GA was established on the basis of natural selection and colony genetics, in which it simulates the law of natural in biology evolutionism. In this paper, we develop a solution approach that is based on a generic GA framework, and use it to compare with our SS. The philosophies of GA and SS are quite different, although both are evolutionary algorithms. Table 1 lists the main differences between GA and SS. In this paper, the procedure of our GA for STT-VRPSPD can be described as follows.

5.1. Chromosome coding We employ the integer coding for chromosomes, where integer numbers from 1 to n are customers. A routing solution can be constructed based on the array of these integer numbers. For individual chromosome, we insert node 0, which denotes the depot, into the array based on the vehicle load capacity and travel time limitation constraints so that we can obtain a feasible solution. Consider an example of 9 customers. We can generate a chromosome (2, 4, 1, 8, 7, 5, 3, 6, 9). If we insert node 0 and get (0, 2, 4, 1, 0, 8, 7, 5, 3, 0, 6, 9, 0) based on the vehicle load capacity and travel time limitation constraints, it means we need three vehicles to serve the customers. Node 0 does not have to be inserted during the selection, crossover and mutation operations; however, node 0 is only inserted as we obtain the solutions and calculate the fitness function. The GA algorithm can be operated easily and we can process the cluster and order simultaneously by this method.

5.2. Create an initial population We create an initial population of solutions by setting the maximum population size maxpop¼100. Each solution can be constructed by randomly creating a full array of n customers.

5.3. Fitness function setting We use a weight scheme to design the fitness function that is defined by f i ¼ z0 =zi , where f i is the fitness of chromosome i, z0 is the minimum travel cost among all chromosomes in the current population, zi is the travel cost of chromosome i. The relative P fitness pi of chromosome i is then given by pi ¼ f i = f i . Table 1 Comparison between scatter search and genetic algorithm.

Selection strategy Combination strategy Evolution strategy

GA

SS

Random selection of parents generation Crossover and mutation

Selection from the reference set Structural combination

Survival of the fittest

Maintain the solution’s quality and diversity

5.4. Genetic operation design There are three main operators in the GA: selection, crossover, and mutation. For the selection operator, we employ a proportion selection approach (roulette wheel selection) to probabilistically select individuals in the population based on their fitness function values. That is, the higher the fitness function value, the more likely the individual will be selected. In particular, we start the selection operator by calculating the total fitness function value for each individual, the relative fitness of the population, and probabilistically selecting an individual from the population. For the crossover operator, we employ the PMX crossover operator. First we perform random pairing of individuals by randomly grouping N individuals into N/2 groups. Then we perform the crossover operation with the probability pc for each group to generate two new individuals (offsprings), and subsequently check the solution quality of the offsprings. If any one of the offspring solutions is better, we perform the crossover. Otherwise, we copy the parents’ chromosome to be a new offspring. For the mutation operator, we employ a reverse mutation operator by randomly creating two mutation points and reversing the mutation segments to generate a new individual with the mutation probability pm . For example, consider a solution 321|4675|8, where ‘‘9’’ indicates the random mutation point. A new individual 32157648 will be generated after the mutation operator. 5.5. Stopping criterion We terminate the GA by pre-determining the maximum number of iterations. Here we choose the maximum iterations¼150. 5.6. Overall procedure The flowchart of our generic GA algorithm for STT-VRPSPD is shown in Fig. 4, and can be described as follows. Step 1: Initialization—given a population size maxpop, current generation (k¼0), randomly create an initial population pop(k) and calculate the fitness values for all solutions in the initial population; Step 2: Selection—Probabilistically select individuals as parents for the next generation and create a new population of the size maxpop; Step 3: Crossover—Pair chromosomes of individuals according to the crossover probability pc and PMX, and keep the offsprings if their fitness values are better than that of their parents to form a new population; Step 4: Mutation—Perform the chromosome reverse mutation according to mutation probability pm and check the stopping criteria. If they are not met, go back to Step 2; otherwise, go to Step 5; Step 5: Terminate the algorithm and output the solution. 6. Computational results In this study, we test our SS and GA algorithms by using the Dethloff [8] instances and the instances we generate. All the programming and implementation of our computational studies were done using the C language on a personal computer running Windows XP with Intel Pentium Dual 1.87 GHz CPU and 2 G of RAM. 6.1. Test data The size of individual random Dethloff test instances is 50 (customers). There are two distribution patterns of customers

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

6.2. Scatter search parameter selection

Begin

Read GA, initial variable

Create initial solution randomly

Stop criterion ?

2285

Y

Output solution

N Calculate each chromosome’s fitness

Select chromosome for next generation by roulette wheel selection

Cross pairs of chromosomes according to cross over probability pm and PMX

Chromosome reverse mutation according to mutation probability pm Fig. 4. Flowchart of generic genetic algorithm framework.

used in this study. The first pattern is SCA, in which the coordinates of customers are distributed in the interval [0,100]. The second pattern is CON, in which the coordinates of half of the customers are distributed in the same way as SCA while the other half’s coordinates are distributed highly in the interval [100/3,200/3], which is similar to the city distribution. The travel distance information of the test instances was generated based on the Euclidean distance. Referring to the generation method of the Dethloff instances, we generate other 4 groups of test instances with 200 or 400 customers. Familiar with the Dethloff instances, in SCA pattern, the coordinates of customers are distributed in the interval [0,200]. In CON pattern, the coordinates of half of the customers are distributed in the same way as SCA while the other half’s coordinates are distributed highly in the interval [200/3,400/3]. The customers’ delivery demands (di) are randomly generated, to be an integer between 0 and 100. The pick-up demands (pi) were calculated by using the function pi ¼ ð0:5 þ r i Þdi , where r i is a random number uniformly distributed between 0 and 1. The vehicle’s maximum load capacity was P determined by Q ¼ ds =m, where m is a pre-determined parameter and ds is customer s’s delivery demand. In this study, we chose the value of m to be 3. Thus we generated 6 groups of test instances, namely SCA3, CON3, S200, C200, S400, and C400. Specifically, test instances SCA3 and CON3 have the same vehicle’s maximum load capacity, test instances S200 and C200 have the same vehicle’s maximum load capacity while test instances S400 and C400 have the same vehicle’s maximum load capacity. We also make the following assumptions for the STT-VRPSPD considered here. The vehicle’s travel speed on every arc is fixed to n ¼ 40. The travel cost cij along arc (i,j) is equal to the travel distance on arc Dij. The vehicle travel time t xijl on arc (i,j) is stochastic and follows a standard distribution with the mean Dij =v and the standard deviation s ¼ 1. There is a upper bound B of the vehicle’s travel time on each arc.

In the SS framework, we note that there are several parameters be to pre-determined: the size of initial population Psize, the size of reference set b, the size of high quality solution reference set b1, and the size of diversity solution reference set b2. Psize is a very influential parameter in the SS framework. If Psizeis large, the efficiency will be decreased while the quality of the initial reference set will be increased, vice versa. Similarly, 9b9 also has the same effect on the algorithm. The sizes of b1 and b2 also play a role in balancing between the algorithm’s efficiency and solution quality. We varied the values of P size to 30, 50, and 100, 9b9 to 15 and 20, and 9 b19 and 9 b29 to 5 and 10. To select the best parameter setting, we generated 10 instances of SCA3-0 and tested our SS approach on those 10 instances with B ¼ 15, a ¼ 0:9, v ¼ 40, n ¼ 50. The performance characteristics of different parameter settings of our SS approach are shown in Table 2. From the table, it appears that the setting Psize ¼100, 9b9 ¼ 15, 9b1 9 ¼ 10, 9b2 9 ¼ 5 appears to provide the best solution quality. Thus, we shall use this setting for the rest of this study. 6.3. Computational results As mentioned earlier, 60 test instances with the parameter setting a ¼ 0:9, B ¼ 15, l ¼ 50 were used to test the performance of our SS approach in this study. The parameters of GA were set as follows:pc ¼ 0:6, pm ¼ 0:1, the number of iterations¼150, the population size¼ 100. The parameters of SS were set as follows: Psize ¼100, 9b9 ¼ 15, 9b1 9 ¼ 10, 9b2 9 ¼ 5. To test the SS and the GA algorithms, we used 6 groups of instances with different customer numbers to test the algorithms. Table 3 shows the performance characteristics of the SS approach in comparison with the GA approach on instances SCA3 and CON3. Table 4 shows the performance characteristics of the SS approach in comparison with the GA approach on instances S200 and C200. Table 5 shows the performance characteristics of the SS approach in comparison with the GA approach on instances S400 and C400. For each test instance, both SS and GA approaches were run 10 times and the results were reported based on an average of the 10 runs. In the tables, the first column represents the names of randomly generated test instances, and the reported performance characteristics include the Best solution (objective function value), Worst solution, Average solution, Average CPU time, and Average number of vehicles of the ten replications. The last column in the table shows the improvement of the SS and GA solutions, which is calculated by improvement ¼ ðav:ObjGA av:ObjSS Þ= av:ObjGA  100%. From Tables 3–5, the SS approach outperforms the GA approach in terms of the solution quality in almost all test instances. For the smallest problems with n¼50 the average solution improvement of the SS over the GA is about 13%, whereas in larger problems n¼200 and 400 the improvements vary around 10.7% and 6.8%, respectively. In addition, it can be observed from Table 1 that the difference between the average and worst solutions of the SS is less than that of Table 2 Performance characteristics of SS tested on SCA3 with different parameter settings. Psize b

b1 b2 Total cost

vehicles Average computing time for optimal (s)

Average iterative times for optimal

100 100 100 50 30

10 5 10 10 10

4 5 4 4 4

24 8 19 35 25

15 10 20 15 15

5 5 10 5 5

1002 1268 1005 1102 1091

22.4 8.2 20.2 11.1 7.6

2286

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

Table 3 Performance characteristics of SS and GA tested using SCA3 and CON3 instances. Instance GA

Sca3-0 Sca3-1 Sca3-2 Sca3-3 Sca3-4 Sca3-5 Sca3-6 Sca3-7 Sca3-8 Sca3-9 Con3-0 Con3-1 Con3-2 Con3-3 Con3-4 Con3-5 Con3-6 Con3-7 Con3-8 Con3-9 Avg

SS

Improvement (%)

Best solution

Worst solution

Average solution

CPU time (s)

# of vehicles

Best solution

Worst solution

Average solution

CPU time (s)

# of vehicles

1060 1068 1034 1032 1112 1058 1029 1001 1073 1035 1077 823 857 819 851 876 798 819 752 838 951

1475 1426 1407 1320 1359 1346 1218 1238 1249 1196 1248 990 1021 1091 1017 1031 1050 1017 913 1097 1185

1245 1201 1185 1197 1214 1198 1102 1131 1091 1131 1169 898 962 952 932 966 915 926 822 961 1060

2.0 3.1 3.1 2.8 3.1 3.1 3.0 3.8 3.3 3.0 2.0 2.1 2.5 2.3 2.3 2.6 2.7 2.3 2.3 3.3

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

966 994 976 954 916 948 898 931 882 935 955 791 848 799 799 830 731 799 710 846 877

1013 1024 1055 1086 1079 1013 1004 1002 1138 1032 1086 883 1029 1052 877 875 773 874 732 961 979

1002 999 1017 1004 1017 981 959 971 951 993 1018 828 886 883 844 855 751 830 720 902 921

22.4 16.5 19.2 19.1 21.2 21.2 18.1 21.3 13.6 21.5 18.2 17.5 15.3 12.9 16.5 16.6 16.3 15.6 15.8 10.6 Minimum Mean Maximum

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

19.5 16.8 17.2 16.1 16.2 18.1 13.0 14.1 12.8 12.2 12.9 7.8 7.9 7.2 9.4 11.5 17.9 10.4 12.4 6.1 6.1 13.0 19.5

Table 4 Performance characteristics of SS and GA tested using S200 and C200 instances. Instance GA

S200-0 S200-1 S200-2 S200-3 S200-4 S200-5 S200-6 S200-7 S200-8 S200-9 C200-0 C200-1 C200-2 C200-3 C200-4 C200-5 C200-6 C200-7 C200-8 C200-9 Avg

SS

Improvement (%)

Best solution

Worst solution

Average solution

CPU time (s)

# of vehicles

Best solution

Worst solution

Average solution

CPU time (s)

# of vehicles

7873 7845 8132 8062 7844 7662 8213 7809 8114 8225 6412 6197 6189 6004 6311 5892 6128 6075 6109 6080 7059

8227 8352 8480 8375 8681 8368 8440 8727 8208 8448 6701 6648 6831 6454 6629 6513 6707 6567 6770 6562 7534

8058 8144 8303 8209 8408 8120 8286 8293 8151 8348 6531 6371 6503 6189 6418 6180 6347 6324 6332 6403 7296

362 346 388 355 329 340 373 314 355 350 371 335 338 399 355 440 362 382 402 388

20 20 20 20 20 20 20 20 20 20 16 16 16 16 16 16 16 16 16 16

6556 7107 6603 6752 6977 6784 7588 6306 6545 6693 5939 5555 5008 4658 5151 5212 5422 5207 5423 5356 6042

8124 8313 7686 8816 7784 7154 8309 8615 7944 7613 7017 6019 6496 6014 5879 6174 5859 5986 5971 6018 7090

7504 7494 6964 7518 7473 7031 7982 6960 7281 7167 6409 5761 5785 5396 5457 5355 5705 5642 5613 5755 6513

431 433 403 343 271 391 378 372 490 454 450 558 533 485 494 468 515 565 496 511 Minimum Mean Maximum

20 20 20 20 20 20 20 20 20 20 16 16 16 16 16 16 16 16 16 16

the GA. This implies that the SS approach is more stable than the GA approach. In Tables 4 and 5, the average solution improvement becomes smaller, and the difference between the average and worst solutions of the SS increases. However, the best solutions of the SS are still more superior to that of the GA. For the problems with n¼50 and 200, both SS and GA approaches yielded the same number of required vehicles. On the other hand, for the problems with n¼400 the number of required vehicles of the SS is less than that of the GA. The main advantage of the SS approach over the GA approach is

6.9 8.0 16.1 8.4 11.1 13.4 3.7 16.1 10.7 14.1 1.8 9.6 11.0 12.8 14.9 13.3 10.1 10.8 11.4 10.1 1.8 10.7 16.1

mainly due to its special feature to efficiently explore a broader search space. Table 6 shows the scalability of both SS and GA approaches when the problem size increases. The results are shown based on an average of 10 replications. As the problem gets larger, it is also noted that the required computational time of the SS also increases at a lower rate than that of the GA. The SS requires more computational time than the GA when the problem is small, while the SS requires less computational time than the GA when

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

2287

Table 5 Performance characteristics of SS and GA tested using S400 and C400 instances. Instance GA

S400-0 S400-1 S400-2 S400-3 S400-4 S400-5 S400-6 S400-7 S400-8 S400-9 C400-0 C400-1 C400-2 C400-3 C400-4 C400-5 C400-6 C400-7 C400-8 C400-9 Avg

SS

Improvement (%)

Best solution

Worst solution

Average solution

CPU time (s)

# of vehicles

Best solution

Worst solution

Average solution

CPU time (s)

# of vehicles

15609 15780 15092 15650 15932 15745 15514 15223 15941 15217 12756 12132 12081 11515 12356 12230 11991 12128 12426 11637 13848

16420 16416 16796 16367 17136 17116 16804 16560 17067 16443 13389 13279 12514 12566 12734 12857 12891 12966 12871 12549 14787

16148 16152 16084 16124 16434 16575 15980 16102 16403 16100 12925 12616 12282 12156 12553 12513 12333 12534 12702 12170 14344

1738 1795 1722 1852 1728 1423 1715 1574 1541 1566 1656 1386 1519 1609 1671 1462 1418 1612 1623 1845

40 40 40 40 40 40 40 40 41 40 32 32 32 32 32 32 32 32 32 31

14123 14088 14244 13281 13691 14274 14250 12587 14161 14178 11863 11036 11429 10272 11497 10985 11398 11296 11017 11474 12557

14526 15278 15059 14859 15001 15193 15323 15681 15053 15212 13163 12946 12643 12191 12726 12305 12960 12832 13695 12103 13937

14332 14180 14403 14233 14366 14634 14893 14402 14850 14763 12554 12219 12096 11264 12156 12068 12077 12503 12099 11840 13297

1343 1233 1117 1378 1122 1172 1359 1310 1203 1128 1305 1241 1238 1257 1193 1707 1121 1290 1203 1282 Minimum Mean Maximum

37 37 37 37 37 37 37 38 37 37 32 31 31 31 31 31 31 31 31 31

50 200 400

GA

SS

Objective value

Computational time (s)

Objective value

Computational time (s)

1060 7296 14787

2.7 364 1623

921 6513 13297

17.5 452 1260

12 expectation time

# of customers

expectation travel time (GA) upper bound of ETT (GA) expectation travel time (SS) upper bound of ETT (SS)

14

Table 6 Performance characteristics of SS and GA with increasing numbers of customers.

11.3 12.2 10.5 11.7 12.6 11.7 6.8 10.6 9.5 8.3 2.9 3.2 1.5 7.3 3.2 3.6 2.1 0.2 4.8 2.7 0.2 6.8 12.6

10 8 6 4 2 0 0

2500

GA SS

20

30 customer node

40

50

60

Fig. 6. Graphical representation of expected travel time of vehicles.

1500 1000 500 0 0

30

60

90

120 150 180 210 240 270 300 330 iterative times

Fig. 5. Evolution of solution convergence by SS and GA.

the problem is large. This clearly indicates that as the problem scale grows, the SS approach appears to be more scalable. Fig. 5 illustrates the evolution of best solutions found by the SS and GA approaches. The solutions are based on the best solution of 10 replications of SCA3-6. The x-axis represents the number of iterations and the y-axis represents the objective function value. From the figure, the convergence speed to high-quality solutions of the SS approach is very prominent. The SS approach found good solutions much faster than the GA approach. This suggests that the SS approach is much better than the GA approach in terms of both efficacy and efficiency.

load capacity

objective value

2000

10

1000 900 800 700 600 500 400 300 200 100 0

delivery load pick-up load total load

0

5

10

15

20

25 30 35 40 customer node

45

50

55

60

Fig. 7. Graphical representation of expected vehicle load of GA solutions.

Fig. 6 illustrates the upper bound of expected travel time, which were calculated based on 10 replications. The figure is plotted based on the best solution of 10 SCA3-6 instances. From the figure, the expected travel time of each vehicle in the SS solution is very close to the upper bound of expected travel time, which means each vehicle make full use of limitation expected

load capacity

2288

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

1000 900 800 700 600

Appendix A. Pseudo-code for scatter search algorithm

delivery load pick-up load total load

500 400 300 200 100 0 0

5

10

15

20

25 30 35 40 customer node

45

50

55

60

Fig. 8. Graphical representation of expected vehicle load of SS solutions.

travel time. Thus, it suggests that each vehicle in the SS solution utilizes the work-time equally. Figs. 7 and 8 illustrate the corresponding vehicle load graphs for the GA and SS approaches tested using the same SCA3-6 instance, respectively. In the figures, the load of each vehicle is different when it departs from the depot. We can observe an evolution of the load of each vehicle. We note that the vehicle’s capacity for this instance is 839. In Fig. 7, we observe that the GA solutions yield the loads of the vehicles in the range of 700 and 790, which represents about 83–94% capacity utilization. On the other hand, in Fig. 8, the SS solutions yield the loads in the range of 780 and 810, which represents about 93–97% capacity utilization. Thus, the SS solutions appear to be more effective.

7. Conclusion In this paper we studied the STT-VRPSPD and modeled the problem using the basis of VRPSPD and VRPSTT, i.e., constructing a chance-constrained programming model of vehicle’s load capacity constraint. We developed a SS approach solve the STTVRPSPD and constructed a generic GA approach to compare and reference the performance of our SS approach. In the computational experiment, we employed the Dethloff data to test both approaches. We first performed a preliminary analysis to finetune the parameters of our SS approach. Subsequently, we tested and compared our SS and GA approaches on 60 random Dethloff instances. The computational results suggest that the SS approach is more effective and superior to the GA approach in obtaining high-quality solutions. Although the GA approach runs faster than the SS approach when the problem is small, the SS approach runs faster than the GA approach when the problem is large. However, our SS framework still needs some improvement in terms of chance-constrained programming because it appears that the load capacity constraint and the time constraint influence the solutions. The best solution we obtained mainly depends on the tighter one of the two constraints. We plan to investigate this issue in our future study.

Acknowledgments This work is supported by the National Natural Science Fund of China under Grant nos. 71171126 and 61170095, the Shanghai Philosophy and Social Science Planning Fund under Grant no. 2011BGL015, the Natural Science Fund of Shanghai under Grant nos. 09ZR1420400 and 09ZR1403000.

Construct N routes for N demand nodes and one route has only one node Calculate any two nodes’ saving value Sn ði,jÞ Sort the value in descending order within the set M while (set Mis not empty) do {if theði,jÞ mapping the first item of Sn ði,jÞ satisfies the conditions then Examine the routing after connecting node i and j if the routing is feasible then Connect node i and j else Let M ¼ MfSn ði,jÞg} for p¼1 to Psize do Select a solution created above randomly and duplicate it Select two exchange bits in the duplicated solution randomly Exchange the selected bits and create a new solution end for while (the whole algorithm does not satisfy the end condition) do { Sorting the solutions according to the objective values on ascending. for each infeasible solutions do while (the load capacity of the vehicle is larger than Q l ) do {Calculate cji =di to discard the highest-cost and lowestefficiency node Update the load capacity of the vehicle} Put the discarded nodes the new vehicle if the load capacity of vehicle is lesser than Q l then while (overloads nodes exist) do { Take out the first overloaded node Put it into the new vehicle routing Update the overloads nodes} while (there is vehicle out of the time limitation constraint) do {Take out the node with the largest value of t ji =di Put it into the new vehicle routing Update vehicle} end for for each feasible solutions do Find out the vehicle with the longest total travel time for each client node i in the vehicle above between the nodes h and S Select randomly nodes j and k in another vehicle of the same solution Calculate: tði,j,k,h,sÞ ¼ tðj,iÞ þtði,kÞ þ tðh,sÞtðh,iÞtði,sÞtðj,kÞ if (tði,j,k,h,sÞ o0 and all constraints are satisfied after inserting) then Insert node i between nodes j and k end for if the travel time of the solution is not improved then Find out maximum-load-capacity vehicle for each client node in the route with the maximumload-capacity vehicle Select randomly two nodes j and k in another vehicle if constraints are satisfied then Insert the node s between nodes j and k end for end for Let top 9b1 9solutions consist of high quality solution reference set. Let top 9b2 9solutions with largest difference consist of diversity solution reference set.

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

Combine any two solutions in reference set to form subsets C 29b9 if number of subsets C 29b9 is larger than Psize then Select Psize subsets to combine to subsets else combine any two elements in all subsets C 29b9 Select segment to exchange on two original solutions randomly Calculates average routing cost c=n before exchange node Keep the segment which has lower average cost Discard kept clients in the other solution and sort the rest client nodes after kept routing orderly }

References [1] Dethloff J. Relation between vehicle routing problems: an insertion heuristic for the vehicle routing problem with simultaneous delivery and pick-up applied to the vehicle routing problem with backhauls. Journal of the Operational Research Society 2002;53(1):115–8. [2] Anily. S. The vehicle-routing problem with delivery and back-haul options. Naval Research Logistics 1996;43(3):415–34. [3] Min H. The multiple vehicle routing problems with simultaneous delivery and pick-up points. Transportation Research A 1989;23(5):377–86. [4] Halse, K. Modeling and solving complex vehicle routing problems. Ph.D. thesis. Institute of Mathematical Statistics and Operations Research, Technical University of Denmark, Lyngby, 1992. [5] Fisher ML, Tang B, Zheng Z. A network flow based heuristic for bulk pickup and delivery routing. Transportation Science 1994;29:45–55. [6] Savelsbergh MWP, Sol M. The general pickup and delivery problem. Transportation Science 1995;29(1):17–29. [7] Gendreau M, Laporte G, Vigo D. Heuristics for the traveling salesman problem with pickup and delivery. Computers and Operations Research 1999;26(7): 699–714. [8] Dethloff J. Vehicle routing and reverse logistics: the vehicle routing problem with simultaneous delivery and pick-up. OR Spektrum 2001;23(1):79–96. [9] Salhi S, Nagy G. A cluster insertion heuristic for single and multiple depot vehicle routing problems with backhauling. Journal of the Operational Research Society 1999;50:1034–42. [10] Tang Montane´ FA, Galva~ o RD. Vehicle routing problems with simultaneous pick-up and delivery service. Journal of the Operations Research Society of India (OPSEARCH) 2002;39(1):19–33. [11] Tang Montane´ FA, Galva~ o RD. A tabu search algorithm for the vehicle routing problem with simultaneous pick-up and delivery service. Computer and Operations Research 2006;33(3):595–619. [12] Gronalt M, Hartl RF, Reimann M. New saving based algorithms for time constrained pickup and delivery of full truckloads. Euorpean Jounral of Operational Research 2003;15:520–35. [13] Salhi S, Nagy G. Heuristic algorithms for single and multiple depot vehicle routing problems with pickups and deliveries. European Journal of Operational Research 2005;162(1):126–41. [14] Zachariadis EE, Tarantilis CD, Kiranoudis. CT. A hybrid metaheuristic algorithm for the vehicle routing problem with simultaneous delivery and pickup service. Expert Systems with Applications 2009;36(2):1070–81. [15] Bianchessi N, Righini G. Heuristic algorithms for the vehicle routing problem with simultaneous pick-up and delivery. Computers & Operations Research 2007;34(2):578–94. [16] Ai TJ, Kachitvichyanukul. V. A particle swarm optimization for the vehicle routing problem with simultaneous pickup and delivery. Computers & Operations Research 2009;36(5):1693–702. [17] Gajpal Y, Abad P. An ant colony system (ACS) for vehicle routing problem with simultaneous delivery and pickup. Computers & Operations Research 2009;36(12):3215–23. [18] Subramanian A, Drummond LMA, Bentes C, Ochi LS, Farias R. A parallel heuristic for the vehicle routing problem with simultaneous pickup and delivery. Computers & Operations Research 2010;37(11):1899–911. [19] C - atay B. A new saving-based ant algorithm for the vehicle routing problem with simultaneous pickup and delivery. Expert Systems with Applications 2010;37:6809–17. [20] Kao EPC. A preference order dynamic program for a stochastic traveling salesman problem. Operations Research 1978;26(6):1033–45. [21] Sniedovich M. Technical note—analysis of a preference order traveling salesman problem. Operations Research 1981;29(6):1234–7. [22] Carraway RL, Morin TL, Moskovitz H. Generalized dynamic programming for stochastic combinatorial optimization. Operations Research 1989;37(5):819–29. [23] Laporte G, Louveaux F, Mercure. H. The vehicle routing problem with stochastic travel times. Transportation Science 1992;26(3):161–70.

2289

[24] Park YB, Song S. Vehicle scheduling problems with time-varying speed. Computers & Industrial Engineering 1997;33(3–4):853–6. [25] Xie, B. Research on stochastic vehicle routing problems. Ph.D. thesis, Xinan Jiaotong University, China, 2003. [26] Shao, Z, Gao, S, Wang, S. A hybrid particle swarm optimization algorithm for vehicle routing problem with stochastic travel time. Fuzzy information and engineering, vol. 54. Berlin/Heidelberg: Springer; 2009.p. 566–74. [27] Woensel TV, Kerbache L, Peremans H, Vandaele N. A queueing framework for routing problems with time-dependent travel times. Journal of Mathematical Modelling and Algorithms 2007;6(1):151–73. [28] Lecluyse C, Van Wosensel T, Peremans H. Vehicle routing with stochastic time-dependent travel times. 4OR: A Quarterly Journal of Operations Research 2009;7(4):363–77. [29] Connors RD, Sumalee A. A network equilibrium model with travellers’ perception of stochastic travel times. Transportation Research, Part B 2009;43:614–24. [30] Chen A, Zhou Z. The a-reliable mean-excess traffic equilibrium model with stochastic travel times. Transportation Research, Part B 2010;44:493–513. [31] Glover F. Heuristics for integer programming using surrogate constraints. Decision Sciences (S0011-7315) 1977;1:156–66. [32] Glover F. A template for scatter search and path relinking. In: Hao JK, Lutton E, Ronald E, Schoenauer M, Snyers D, editors. Artificial evolution, lecture notes in computer science, 1363. Germany: Springer; 1998. p. 13–54. [33] Campos V, Laguna M, Marti R. Scatter search for the linear ordering problem. In: Corne D, Dorigo M, Glover F, editors. New ideas in optimization. UK: McGraw-Hill; 1999. p. 331–41. [34] Laguna M, Marti R, Campos V. Intensification and diversification with elite tabu search solutions for the linear ordering problem. Computers and Operations Research 1998;26(12):1217–30. [35] Campos V, Glover F, Laguna M, Marti R. An experimental evaluation of a scatter search for the linear ordering problem. Journal of Global Optimization 2001;21(4):397–414. [36] Campos V, Laguna M, Marti R. Context-independent scatter and tabu search for permutation problems. INFORMS Journal on Computing 2005;17(1): 111–22. [37] Glover F, Laguna M, Marti R. Scatter search. Advances in evolutionary computing: theory and applications. In: Ghosh A, Tsutsui S, editors. New York: Springer-Verlag; 2003. p. 519–37. ISBN:3-540-43330-9. [38] Michalewicz Z, Logan TD, Swaminathan S. Evolutionary operators for continuous convex parameter spaces. In: The 3rd annual conference on evolutionary programming. River Edge, NJ: World Scientific Publishing; 1994. p. 84–97. [39] Habiba, D, Nabila, A. Scatter search for SAT and W-MAX-SAT problems. In: The 33rd southeastern symposium on system theory, IEEE, Athens, OH, USA, March, 2001, p. 105–9. [40] Hung WNN, Song X. BDD variable ordering by scatter search. Computer design. ICCD proceedings of the international conference on computer design. Austin, Texas: IEEE Computer Society; 2001. p. 368–73. [41] Hung WNN, Song X, Aboulhamid EM, Driscoll. MA. BDD minimization by scatter search. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, IEEE 2002;21(8):974–9. [42] Hung, WNN, Song, X. On optimal cell assignments in PCS networks. In: The 21st IEEE International Proceedings of the Performance, Computing, and Communications Conference, IEEE, Phoenix, AZ, USA; 2002. p. 225–232. [43] Nebro, AJ, Luna, F, Alba, E. New ideas in applying scatter search to multiobjective optimization. In: Evolutionary multi-criterion optimization (S03029743). Lecture Notes in Computer Science. Berlin/Heidelberg: Springer. 2005, vol. 3410. p. 443–58. [44] Marti R, Lourenco H, Laguna M. Assigning proctors to exams with scatter search. In: Laguna M, Gonza´lez Velarde JL, editors. Computing tools for modeling, optimization and simulation. Kluwer Academic Publishers; 2000. p. 215–27. [45] Laguna M, Marti R. Experimental testing of advanced scatter search designs for global optimization of multimodal functions. Journal of Global Optimization 2005;33(2):235–55. [46] Beausoleil RP. MOSS multiobjective scatter search applied to non-linear multiple criteria optimization. European Journal of Operational Research (S0377-2217) 2005;169(2):426–49. [47] Molina J, Laguna M, Marti R, Caballero. R. SSPMO: A scatter tabu search procedure for non-linear multiobjective optimization. Informs Journal on Computing (S0899-1499) 2007;19(1):91–100. [48] Francois JL, Martin-del-Campo C, Morales LB. Development of a scatter search optimization algorithm for boiling water reactor fuel lattice design. Nuclear Science and Engineering (S0029-5639) 2007;155(3):367–77. [49] Krishna AG, Rao KM. Multi-objective optimization of surface grinding operations using scatter search approach. International Journal of Advanced Manufacturing Technology (S0268-3768) 2006;29(5):475–80. [50] Blazewicz J, Glover F, Kasprzak M. DNA sequencing-tabu and scatter search combined. Informs Journal on Computing (S0899-1499) 2004;16(3):232–40. [51] Trafalis TB, Kasap S. A novel metaheuristics approach for continuous global optimization. Journal of Global Optimization (S0925-5001) 2002;23(2):171–90. [52] Santana-Quintero, LV, Ramirez, N, Coello., CC. A multi-objective particle swarm optimizer hybridized with scatter search. In: The 5th Mexican international conference on artificial intelligence. Apizaco, Mexico, November 13–17, Lecture Notes in Artificial Intelligence. 2006, vol. 4293. p. 294–304.

2290

T. Zhang et al. / Computers & Operations Research 39 (2012) 2277–2290

[53] Russell R, Chiang W. Scatter search for the vehicle routing problem with time windows. European Journal of Operational Research (S0377-2217) 2006;169(2): 606–22. [54] Greistorfer P. A tabu scatter search metaheuristic for the arc routing problem. Computers & Industrial Engineering (S0360-8352) 2003;44(2):249–66. [55] Keskin BB, Uster H. A scatter search-based heuristic to locate capacitated transshipment points. Computers & Operations Research (S0305-0548) 2007;34(10):3112–25. [56] Alvarez AM, Gonza´lez-Velarde JL, De-Alba K. Scatter search for network design problem. Annals of Operations Research (S0254-5330) 2005;138(1):159–78. [57] Gutierrez, MA de-los-Cobos, AS, Goddard, J. A comparison between scatter search and the RAND method for solving the joint replenishment problem. In: The 5th Mexican international conference on artificial intelligence. Leipzig, Apizaco, Mexico, Lecture Notes in Artificial Intelligence, November 13–November 17, 2005. p. 287–95. [58] Gonzalez B, Adenso-Diaz B. A scatter search approach to the optimum disassembly sequence problem. Computers & Operations Research (S03050548) 2006;33(6):1776–93.

[59] Nowicki E, Smutnicki C. Some aspects of scatter search in the flow-shop problem. European Journal of Operational Research (S0377-2217) 2005;169(2): 654–66. [60] Alvarez-Valdes R, Crespo E, Tamarit. JM. A scatter search algorithm for project scheduling under partially renewable resources. Journal of Heuristics (S1381-1231) 2006;12(1–2):95–113. [61] Rahimi-Vahed AR, Rabbani M, Tavakkoli-Moghaddam. R. A multi-objective scatter search for a mixed-model assembly line sequencing problem. Advanced Engineering Informatics (S1474-0346) 2007;21(1):85–99. [62] Zhang T, Tian W, Zhang Y, Liu S. The improved ant colony algorithm for the VRPSPD with maximum distance constraint. Systems engineering-theory & Practice 2008;28(1):132–40. [63] Clarke G, Wright JW. Scheduling of vehicles from a central depot to a number of delivery points. Operations Research 1964;12(4):568–81. [64] Holland JH. Adaptation in natural and artificial systems. Cambridge, Massachusetts: MIT Press; 1992.