- Email: [email protected]
Stochastic single vehicle routing with a predefined customer sequence and multiple depot returns
European Journal of Operational Research 197 (2009) 557–571
Contents lists available at ScienceDirect
European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor
Production, Manufacturing and Logistics
Stochastic single vehicle routing with a predefined customer sequence and multiple depot returns A. Tatarakis, I. Minis * Department of Financial and Management Engineering, University of the Aegean, 31 Fostini Street, 82100 Chios, Greece
a r t i c l e
i n f o
Article history: Received 18 December 2007 Accepted 3 July 2008 Available online 17 July 2008 Keywords: Logistics Stochastic vehicle routing Multiple product delivery Dynamic programming
a b s t r a c t We study the routing of a single vehicle that delivers multiple products under stochastic demand. Specifically, we investigate two practical variations of this problem: (i) The case in which each product type is stored in its dedicated compartment in the vehicle, and (ii) the case in which all products are stored together in the vehicle’s single compartment. Suitable dynamic programming algorithms are proposed to determine the minimum expected (routing) cost for each case. Furthermore, the optimal routing policy is derived by developing appropriate theorems. The efficiency of the algorithms is studied by solving large problem sets. Ó 2008 Published by Elsevier B.V.
1. Introduction Both industry and academia have long recognized the potential for optimizing operations in goods distribution, especially since the latter may contribute up to 20% or more to the total costs of a product (see [1]). Formally, most problems in goods distribution are related to the Vehicle Routing Problem (VRP), a generalization of the classic Travelling Salesman problem (TSP) (see [2–4]). The VRP seeks a set of efficient vehicle routes to serve a number of geographically dispersed customers. In this paper, we address a special case of delivery: We consider a single vehicle that serves customers at a predefined sequence; the vehicle is allowed to return to the depot for stock replenishment. The customer demands are random and are modelled by independent variables with known statistics (derived from historical data). The distances, or travel times, among all points in the network (depot and customer sites) are known. The product quantities loaded onto the vehicle are limited by appropriate capacities, which cannot be exceeded. The vehicle returns to the depot in order to refill as needed. Upon completion of service at each customer, the driver has to make a decision: (a) Proceed to the next customer, as long as the probability of being able to fully satisfy its demand (for all products) is acceptable and the corresponding distance (i.e. the distance from the next customer to the depot) is favourable; (b) return to the depot in order to refill, and resume the route by visiting the next customer in the predefined sequence. If the vehicle proceeds to the next customer but the actual demand of this customer turns out to be higher than the stock carried on board (for any product), the vehicle will deliver as much as feasible, will then go to the depot to refill, and return to the customer to fully satisfy its demand. We refer to this problem as the Stochastic Vehicle Routing with Depot Returns Problem (SVRDRP). Practical applications of the SVRDRP may arise in different settings. A prime example is the Ex-van deliveries model, in which during each customer visit the driver simultaneously sells and delivers high-demand commodities to retail outlets. The concept behind Ex-van deliveries is to facilitate regular stock replenishment of certain types of commodities, so that the outlet (super market, kiosk, etc.) can maximise its sales. The mission of each Ex-van vehicle is to visit all assigned customer sites, typically in predefined sequence, and replenish the stock of selected products. The demand of each customer point is not known in advance but it is revealed upon arrival. Therefore, the total demand of a scheduled customer sequence typically exceeds the total capacity of the vehicle for a particular item, forcing the vehicle to return to the depot in order to replenish its own stock, before resuming its route. In addition to the Ex-van environment, similar delivery characteristics may also arise in other practical settings. For example, material handling systems in a manufacturing shop often operate along fixed pathways that connect the material warehouse with workcenters located along these pathways. Consider the case of Automated Guided Vehicle Systems (AGVs), which are self-propelled vehicles typically
* Corresponding author. Tel.: +30 2271035454. E-mail addresses: [email protected] (A. Tatarakis), [email protected] (I. Minis). 0377-2217/$ - see front matter Ó 2008 Published by Elsevier B.V. doi:10.1016/j.ejor.2008.07.006
558
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
guided along a magnetic induction strip, or a painted strip on the shop floor, and transport multiple discrete parts to workcenters, obviously in a predefined sequence. Note that in addition to the main pathway connecting the workcenters, there are spurs connecting each workcenter with the material warehouse, allowing the return and reloading of the AGV (see [5]). In this case, demand may also be stochastic, depending on the state of the buffers of the workcenters been served. The SVRDRP for a single product has been studied by Yang et al. [22]. In this paper, we study the SVRDRP for two cases of multiple product deliveries and develop (a) algorithms to determine the minimum expected routing cost and (b) policies that allow the vehicle’s driver to make the optimal decisions after serving each customer based on the level of the stock on board. The rest of the paper is organized as follows. In the next section, we present relevant research to date. In Section 3, we present the compartmentalized load case and in Section 4, we present the unified load case. In both these Sections, we illustrate the results obtained through appropriate examples. In Section 5, we study the efficiency of the proposed algorithms. The conclusions of this work are summarized in Section 6. 2. Background The Vehicle Routing Problem (VRP) was introduced by Dantzig and Ramser, almost 50 years ago (see [6]). The objective of the VRP is to deliver goods to a set of customers with known demands following minimum-cost vehicle routes originating from and terminating to a depot (see [7–10]). A rich survey of significant research results in this fundamental problem of Operations Research is given by Toth and Vigo [11]. The Stochastic Vehicle Routing Problem (SVRP) refers to a family of problems that combine the characteristics of stochastic and integer programs, and are often regarded as computationally intractable (see [12]). Therefore, only relatively small instances can be solved to optimality and effective heuristics are hard to design and assess. A characteristic problem is the VRP with Stochastic Demands (see [13,14]). In this problem, the demands are usually independent random variables that may (or may not) follow a known distribution. This problem has significant practical applications. A typical example is garbage collection, in which it is impossible to know a priori the quantity to be collected from each site. Another example, is the delivery of petrol to petrol stations. In this case, it is not known how much petrol will be sold due to the time elapsed between customer order and delivery. In order to address the inherent uncertainty in this type of VRP, a recourse action (i.e. return to the depot in order to refill) is usually embedded into the formulation of the problem, and penalties are incurred in the case of a route failure (see [15]). Due to the stochastic nature of this problem, the objective function is the expected value of the total route cost; (see [16–19]). Bertsimas (in [16]) constructs an a priori sequence among all customers of minimal expected total length and proposes suitable heuristics. His approach is shown to be a strong alternative to the strategy of re-optimization. More recently, Birattari et al. (see [19]), employed five metaheuristics (simulated annealing, tabu search, local search, ant colony optimisation, and evolutionary algorithms) and tested the effect of hybridization (TSPapproximation and VRPSD-approximation). The Stochastic Vehicle Routing Problem (SVRP) for a single product has been addressed in [16–19]. Secomandi (see [20,21]) has proposed a re-optimization approach, in which after the demand of each customer has been revealed the remaining part of the problem is re-solved by using Neuro-dynamic programming. This approach may yield better solutions than the preventive restocking strategy (returning to the depot before a stock out actually occurs), but it is much more computationally expensive. Moreover, the initially planned route may be altered completely (therefore, resembling more the SVRP that the SVRDRP), and this situation may present a limitation in practice. Yang et al. (see [22]) investigate the single and multi-vehicle single product SVRP, which has strong similarities with our case, as it entails a predefined customer visit sequence. Instead of adopting a simple recourse action, an optimal restocking policy of the vehicle has been incorporated in the route design. In particular, the restocking points are deliberately planned along the route, such that the probability of the route failure, and the accompanying recourse cost (including any penalty) is reduced, and the expected total routing cost is minimized. The single vehicle, multiple product case with deterministic demands has recently been presented in [23], in which three practical cases are studied: (i) The case in which each product type is stored in its own compartment in the vehicle, (ii) the case in which all product types are stored together in the vehicle’s single compartment, and (iii) the case of pickup and delivery. For each case, suitable dynamic programming algorithms to determine the optimal routing of the vehicle were developed, and their efficiency was studied by solving large problem sets. In the present paper, we propose and investigate two practical variations of the VRDRP with stochastic demands (SVRDRP): (i) the case of multiple product deliveries when each product is stored in its own compartment in the vehicle, and (ii) the case of multiple product deliveries when all products are stored together in the vehicle’s single compartment. 3. The compartmentalized load case We consider a set of nodes V = {0, 1, . . . , n}, with node 0 denoting the depot and nodes 1, . . . , n corresponding to customers, and a set of arcs A = {(j, j + 1):j, j + 1 2 V} [ {(0, j):j 2 V, j – 0} [ {(j, 0):j 2 V, j – 0} that join the customers along the route 1 ? 2 ? ? n, as well as all customers to the depot. The travel cost (distance) of each arc (j, j0 ) is denoted by cjj0 > 0. The cjj0 values satisfy the triangular inequality. In this problem, we assume that the delivery vehicle is divided into K sections and each section is suitable for carrying one product type P only. Let Qi be the capacity of the vehicle for product i 2 {1, . . . , K}. Clearly, Ki¼1 Q i ¼ Q the total load capacity of the vehicle. Note that all 3 product quantities are calculated using the same unit of measure e.g. m or kg. We define nij to be the stochastic demand of customer j 2 {1, . . . , n} for product type i 2 f1; . . . ; Kg:nij follows a discrete distribution with mi possible values, ni1, ni2, . . ., nimi and probability mass function pijk ¼ Pðnij ¼ nik Þ. The objective of the problem is for the vehicle to serve all customers according to the predefined sequence 1, . . . , n and minimize the expected travel cost. The vehicle is at the depot initially and after serving all customers it returns to the depot. 3.1. The minimum expected cost For simplicity, but without loss of generality, we will present our approach for two product types. Let fj(z1, z2) denote the total minimum expected cost from customer j onward if the vehicle, after serving customer j, carries quantities z1 of product type-1 and z2 of product type-
559
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
2. If Sj represents the set of all possible loads that a vehicle can carry after serving customer j, then fj(z1, z2) for (z1, z2) 2 Sj satisfies the following dynamic programming recursion:
8 P P fjþ1 ðz1 n1k1 ; z2 n2k2 Þp1jþ1;k1 p2jþ1;k2 cj;jþ1 þ > > > 1k > k1 :n 1 6z1 k2 :n2k2 6z2 > > i > > P h P > 1k1 > þ 2c þ f ðz þ Q n ; Q Þ p1jþ1;k1 p2jþ1;k2 > 1 jþ1;0 jþ1 1 2 > > > k1 :z1 > i > > P h < þ P 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 k1 :n1k1 6z1 k2 :z2 fj ðz1 ; z2 Þ ¼ min > i > P P h > > > þ 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 > > > k1 :z1 > > > > m1 P m2 > P > > > fjþ1 ðQ 1 n1k1 ; Q 2 n2k2 Þp1jþ1;k1 p2jþ1;k2 : cj;0 þ c0;jþ1 þ
9 > > > > > > > > > > > > = > > > > ðaÞ > > > > > > > > ; )
ð1Þ
k1 ¼1 k2 ¼1
ðbÞ with the boundary condition fn(z1, z2) = cn,0 z1, z2 2 Sj. In Eq. (1) the minimum is between the two terms; the top term, (a), represents the expected cost of going directly to the next customer. The bottom term, (b), represents the expected cost of the restocking action. Part (a) consists of four summation terms (in the four rows), which correspond to the four areas of the z1, z2 space shown in Fig. 1. The first summation term of (a) represents the expected cost incurred if the stock of both items z1 and z2 is sufficient to fully satisfy the demands n1k1 and n2k2 of customer j + 1 (Area A in Fig. 1). In Fig. 2 the loads of the vehicle before and after serving customer j + 1 are schematically represented above the respective arcs (the demand of customer j + 1 is presented above the node). The stock left onboard after serving customer j + 1 is ðz1 n1k1 Þ P 0 and ðz2 n2k2 Þ P 0, respectively. The second summation of part (a) represents the expected cost incurred if the vehicle proceeds to the next customer directly and the stock z1 is not sufficient to fully satisfy the demand n1k1 of the next customer, while the stock z2 is sufficient to fully satisfy the demand n2k2 (Area B in Fig. 1). The path of the vehicle and the corresponding vehicle loads are shown in Fig. 3. In this case, the vehicle will visit customer j + 1, it will fully satisfy demand n2k2 for product type-2 but will partially satisfy demand n1k1 for product type-1; thus, it is required to return to the depot for stock replenishment. The stock left on board after serving customer j + 1 for
ξ 2k 2 Q2 C
D
A
B
z2
Q1
z1
ξ 1k1
Fig. 1. The solution space per customer point.
ξ 1k1 , ξ 2k 2
j
z1,z2
z1 - ξ 1k1 , z2 - ξ 2k 2
j+1
Fig. 2. The case that the demand is fully satisfied: z1 P n1k1 ; z2 P n2k2 .
ξ 1k
z1,z2 j
0, z2 - ξ 2k
, ξ 2k
Q1 – ( ξ 1k - z1), Q2
j+1 Q1, Q2
D
Fig. 3. The case in which z1 < n1k1 and n2k2 6 z2 .
560
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
ξ 1k
z1,z2 j
1k 2k , ξ 2k Q1 – ( ξ - z1), Q2 – ( ξ - z2)
j+1
0,0 Q1, Q2
D
Fig. 4. The case in both z1 < n1k1 and z2 < n2k2 .
ξ 1k
z1,z2
j
, ξ 2k
Q1 - ξ 1k , Q2 – ξ 2k
j+1 Q1, Q2
D
Fig. 5. The proactive depot return case.
the first time is 0 and ðz2 n2k2 Þ P 0, respectively. The stock left on board after serving customer j + 1 for the second time is Q 1 ðn1k1 z1 Þ for product type-1 and Q2 for product type-2, respectively. The third sum of part (a) represents the expected cost incurred if the stock z2 is not sufficient to fully satisfy the demand n2k2 of the next customer, while stock z1 is sufficient to satisfy demand n1k1 (Area C in Fig. 1). The path of the vehicle and the corresponding vehicle loads are analogous to those of the previous case. The fourth term represents the expected cost incurred if the stock of both items z1 and z2 is not sufficient to fully satisfy the demands n1k1 or n2k2 of the next customer (Area D in Fig. 1). The path of the vehicle and the corresponding vehicle loads are shown in Fig. 4. Part (b) of the minimization equation, represents the expected cost incurred if the vehicle proceeds to the next customer j + 1 via the depot and, therefore, the stock of both items on board z1 and z2 is replenished to Q1 and Q2, respectively. Due to the fact that this is a proactive depot return, the values of z1 and z2 do not affect the result; i.e. the latter is independent of z1 and z2. The path of the vehicle and the corresponding vehicle loads are shown in Fig. 5. 3.2. The optimal routing policy The value of the minimum expected cost can be estimated using Eq. (1). However, the decision required to guide the distribution operation in the SVRDRP environment is the following: Having served customer j, and if the vehicle remaining load is (z1, z2), which is the preferred action; (a) to proceed to customer j + 1, or (b) to visit the depot prior to serving customer j + 1. This decision can be modelled by the following variable: 0
if part (a) ≤ part (b) in Eq. (1); i.e. proceed to customer j+1
1
if part (a) > part (b) in Eq. (1); i.e. proceed to the depot
xj(z1, z2) =
Thus, the vehicle proceeds to the next customer if the corresponding expected cost (until the end of the route) is minimum; otherwise the route via the depot is preferred. This decision variable is also used in the deterministic counterpart of this problem (see [23]). In order to develop the policy to obtain this minimum expected cost (in an average sense) we will extend the threshold theorem presented by Yang et al. [22] to multiple dimensions. Lemma 1
fj ðz1 ; z2 Þ 6 fj ðQ 1 ; Q 2 Þ þ 2c0j
for all z1 ; z2 2 Sj :
ð2Þ
Proof. From Eq. (1) we obtain Q1 X
Q2 X
k1 :n1k1 ¼1
k2 :n2k2 ¼1
fj ðz1 ; z2 Þ 6 cj;0 þ c0;jþ1 þ
fjþ1 ðQ 1 n1k1 ; Q 2 n2k2 Þp1jþ1;k1 p2jþ1;k2 :
ð3Þ
For z1 = Q1 and z2 = Q2 part (a) of Eq. (1) is always less than or equal to part (b) since cj,j+1 6 cj,0 + c0,j+1; therefore
fj ðQ 1 ; Q 2 Þ ¼ cj;jþ1 þ
Q1 X
Q2 X
fjþ1 ðQ 1 n1k1 ; Q 2 n2k2 Þp1jþ1;k1 p2jþ1;k2
k1 :n1k1 ¼1 k2 :n2k2 ¼1
taking into consideration that the last three terms of part (a) of Eq. (1) become zero.
ð4Þ
561
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Substituting Eq. (4) in Eq. (3) results in
fj ðz1 ; z2 Þ 6 cj;0 þ c0;jþ1 þ ½fj ðQ 1 ; Q 2 Þ cj;jþ1 6 cj;0 þ ½c0;j þ cj;jþ1 þ ½fj ðQ 1 ; Q 2 Þ cj;jþ1 6 2cj;0 þ fj ðQ 1 ; Q 2 Þ
As before, consider z1 and z2, the item quantities on board after serving customer j. j
Theorem 1. For each customer j, there exists a threshold function h2 ðz1 Þ, such that the optimal decision, after serving customer j, is to continue to j customer j + 1 if z2 > h2 ðz1 Þ or return to the depot otherwise. Proof. We will first show by induction that for all (z1, z2) 2 Sj, fj(z1, z2) is a non-increasing function. That is, for z1, z2 2 Sj and d1,d2 P 0
fj ðz1 þ d1 ; z2 þ d2 Þ 6 fj ðz1 ; z2 Þ:
ð5Þ
This relationship is true for the last customer n, where fn(z1, z2) = cn0 is independent of (z1, z2). Hence fn(z1, z2) is non-increasing with respect to (z1, z2) 2 Sn. We will now prove that, if fj+1(z1, z2) is non-increasing with respect to (z1, z2) 2 Sj+1, then fj(z1, z2) is also non-increasing with respect to (z1, z2) 2 Sj. Let Hj(z1, z2) and H0j ðz1 ; z2 Þ denote the values of part (a) and part (b) inside the minimisation in Eq. (1). Furthermore, let
Hj ðz1 ; z2 Þ ¼ Haj ðz1 ; z2 Þ þ Hbj ðz1 ; z2 Þ þ Hcj ðz1 ; z2 Þ þ Hdj ðz1 ; z2 Þ;
ð6Þ
where (see also Fig. 1):
1k1
Hbj ðz1 ; z2 Þ
X
¼
X
X
Haj ðz1 ; z2 Þ ¼ cj;jþ1 þ
k1 :n
h
X
k1 :z1
Hcj ðz1 ; z2 Þ
X
¼
X
k1 :n1k1 6z1 k2 :z2
X
Hdj ðz1 ; z2 Þ ¼
fjþ1 ðz1 n1k1 ; z2 n2k2 Þp1jþ1;k1 p2jþ1;k2 ;
6z1 k2 :n2k2 6z2
X
h
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 ;
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 ;
h
ð6:1Þ
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :
k1 :z1
If we expand each one of the last three terms taking into account the regions of Fig. 6 we obtain
X
Hbj ðz1 ; z2 Þ ¼
h i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
X
k1 :z1
X
þ
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 ;
ð6:2Þ
k1 :z1 þd1
X
Hcj ðz1 ; z2 Þ ¼
h i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
k1 :n1k1 6z1 k2 :z2
þ
X
X
h
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 ;
ξ
2k
ð6:3Þ
k1 :n1k1 6z1 k2 :z2 þd2
2
Q2 (c2)
(d4)
(d3)
(c1)
(d1)
(d2)
z2+δ 2
z2 (b1)
(a)
z1
(b2)
z1+δ 1
Fig. 6. The definition space of Hj(z1, z2).
Q1
ξ
1k 1
562
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Hdj ðz1 ; z2 Þ ¼ þ
h i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
X
k1 :z1
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
X
k1 :z1 þd1
þ þ
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
2k2 k1 :z1
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :
X
ð6:4Þ
k1 :z1 þd1
Similarly Hj(z1 + d1, z2 + d2) can be written using Fig. 7 as follows: 0
0
0
0
Hj ðz1 þ d1 ; z2 þ d2 Þ ¼ Haj ðz1 þ d1 ; z2 þ d2 Þ þ Hbj ðz1 þ d1 ; z2 þ d2 Þ þ Hcj ðz1 þ d1 ; z2 þ d2 Þ þ Hdj ðz1 þ d1 ; z2 þ d2 Þ:
ð7Þ
Again, if we expand each one of the above terms taking into account the regions of Fig. 7 we obtain
X
X
k1 :n1k1 6z1
2k2
0
Haj ðz1 þ d1 ; z2 þ d2 Þ ¼ cj;jþ1 þ
X
þ
k2 :n
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2
X6z2
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2
k1 :z1
þ
X 1k1
X
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2
6z1 k2 :z2
k1 :n
X
þ
X
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2 ;
ð7:1Þ
k1 :z1
X
0
Hbj ðz1 þ d1 ; z2 þ d2 Þ ¼
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
k1 :n1k1 >z1 þd1 k2 :n2k2 6z2
X
þ
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 ;
ð7:2Þ
k1 :n1k1 >z1 þd1 k2 :z2
X
X
0
Hcj ðz1 þ d1 ; z2 þ d2 Þ ¼
k1 :n1k1 6z1 k2 :n2k2 >z2 þd2
þ
X
h
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 h
X
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 ;
ð7:3Þ
k1 :z1z2 þd2 0
Hdj ðz1 þ d1 ; z2 þ d2 Þ ¼
X
X
h i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; z2 þ d2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :
ð7:4Þ
k1 :n1k1 >z1 þd1 k2 :n2k2 >z2 þd2
Subtracting the terms that correspond to the same nine regions of Figs. 6 and 7 we obtain
h i h i a0 a0 Hj ðz1 ; z2 Þ Hj ðz1 þ d1 ; z2 þ d2 Þ ¼ Haj ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ þ Hbj 1 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ h i h i a0 b0 þ Hbj 2 ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ þ Hcj 1 ðz1 ; z2 Þ Hj 4 ðz1 þ d1 ; z2 þ d2 Þ h i h i a0 b0 þ Hdj 1 ðz1 ; z2 Þ Hj 3 ðz1 þ d1 ; z2 þ d2 Þ þ Hdj 2 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ h i h i c0 c0 þ Hcj 2 ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ þ Hdj 4 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ h i 0 þ Hdj 3 ðz1 ; z2 Þ Hdj ðz1 þ d1 ; z2 þ d2 Þ : ξ
2k
2
Q2 (c′1)
(c′2)
(d′)
(a′4)
(a′3)
(b′2)
(a′1)
(a′2)
(b′1)
z2+δ 2
z2
z1
z1+δ 1
Fig. 7. The definition space of Hj (z1 + d1, z2 + d2).
Q1
ξ
1k 1
ð8Þ
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
563
The detailed expression of Eq. (8) is given in Appendix A. In Eqs. (8) and (A.1), we can distinguish two types of terms. The first type includes the first, third, sixth, seventh, eighth, and ninth terms. Considering the first term we have a0
Haj ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ ¼
X
X
½fjþ1 ðz1 n1k1 ; z2 n2k2 Þ fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2 :
k1 :n1k1 6z1 k2 :n2k2 6z2
Since fj+1(z1, z2) is monotonically non-increasing, this term is non-negative. Similarly we can show that the differences corresponding to the third, sixth, seventh, eighth, and ninth terms are also non-negative. The second type of terms includes the second, forth and fifth terms. Considering the second term we have a0
Hbj 1 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ i X X h 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :1 ¼ k1 :z1
Recall that from Lemma 1:
fj ðz1 ; z2 Þ 6 fj ðQ 1 ; Q 2 Þ þ 2c0j
for all z 2 Sj ;
therefore,
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ 6 fjþ1 ðQ 1 ; Q 2 Þ þ 2cjþ1;0 : However,
fjþ1 ðQ 1 ; Q 2 Þ 6 fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ since n1k1 > z1 and fj+1 is non-increasing. Thus,
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ 6 fjþ1 ðQ 1 ; Q 2 Þ þ 2cjþ1;0 6 fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ þ 2cjþ1;0 : a0
Therefore, Hbj 1 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ P 0. Similarly we can show that the forth and fifth terms of Eq. (8) are also non-negative. Considering the non-negativity of all terms in Eq. (8) it is clear that Hj(z1, z2) Hj(z1 + d1,z2 + d2) P 0 and hence, Hj(z1, z2) is a nonincreasing function. Note now that part (b) of Eq. (1), H0j ðz1 ; z2 Þ, is independent of z1, z2 and, thus, a constant H0j ðz1 ; z2 Þ ¼ H0j in the z1, z2 space. Fig. 8 plots the terms Hj(z1, z2n) and H0j with o Hj(z1, z2) been non-increasing as shown above. From this Figure it is clear that fj ðz1 ; z2 Þ ¼ min Hj ðz1 ; z2 Þ; H0j ðz1 ; z2 Þ is a non-increasing function. Furthermore, it is also clear that due to intersection Hj ðz1 ; z2 Þ \ H0j for j j customer j and for each value of z1 there exists a value h2 , such that for z2 > h2 Hj ðz1 ; z2 Þ < H0j and, thus, the optimal decision after serving customer j is to continue to customer j + 1. Otherwise a return to the depot is preferable (in the mean value sense). This concludes the proof of Theorem 1. h
Fig. 8. The combined threshold graphical representation.
564
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
3.3. Complexity vs. the number of products The influence of the number of product types to the complexity of the problem is significant, and worth noting. As pointed out in Section 2, the formulation for the single product counterpart of the problem presented in this paper is given in Yang et al. (see [22]). According to this reference we have
fj ðzÞ ¼ min
8 P P fjþ1 ðz nk Þpjþ1;k þ ½b þ 2cjþ1;0 þ fjþ1 ðz þ Q nk Þpjþ1;k cj;jþ1 þ > > < k k:n 6z k:nk >z
ðaÞ
m P > > : cj;0 þ c0;jþ1 þ fjþ1 ðQ nk Þpjþ1;k
ðbÞ
k¼1
with the boundary condition
fn ðzÞ ¼ cn0 z 2 Sn : The complexity increase is evident by comparing the above equation to Eq. (1). In Eq. (1), the number of terms to be evaluated for each part, (a) or (b), is m1m2 while in the single product case, the corresponding number of terms is m. It is, thus, clear that the problem complexity increases as the space of possible loads increases. As a consequence, for example, in the three product case, the number of terms increases to m1m2m3 where m3 is the number of possible quantities for product i = 3. 3.4. Implementation and example In order to solve the compartmentalized case of multiple product delivery, we developed an appropriate algorithm that uses Dynamic Programming to derive the optimal solution in a reasonable amount of time. Based on the formulation presented in Section 3.1, the algorithm starts from the end of the route (last customer to be visited) and iterates towards the beginning of the route, calculating the remaining minimum expected cost from each customer site until the end of the route. This procedure computes the minimum expected cost of the route, given a distance matrix and the demand probability mass functions. Based on the result of the algorithm, the threshold function h2(z1) for each customer j can be obtained, in line with what has already been described in Section 3.2. Consider the five-customer network of Fig. 9. The vehicle capacity is Q = 10 units and is equally split between two products (Q1 = Q2 = 5); the demand nij for each product i (i = 1, 2) and customer j (j = 1, . . . , 5) is given in Appendix B, and the distances between the nodes cjj0 are given in Fig. 10.
Fig. 9. Five-customer network for the multiple product extension.
Customer1 5
4
3
z2
h12 (3)
2
1
No return Return
0 0
1
2
3
z1
4
5
Fig. 10. The driver decision based on the load combinations after serving customer 1.
565
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571 1
The problem is solved using the dynamic programming algorithm. Fig. 10 illustrates the threshold function h2 ðz1 Þ as obtained from the results for customer 1. In this case H01 ¼ 70; 2. For any (z1, z2) combination below the border defined by the switch between the red and green points, the optimal decision is to return to the depot in order to refill the vehicle. Conversely, for any combinations above this border the optimal decision would be to proceed directly to the next customer. This policy can be very simply and clearly communicated to the vehicle driver and result in significant cost savings for the fleet operator. 4. The unified load case In the unified load case the vehicle may carry any quantity of product i 2 {1, . . . , K}, provided that the total capacity Q of the vehicle is not exceeded. Upon completion of service at customer site j, the vehicle’s driver has to make the same decision as the one described in Section 3, i.e. proceed to customer j + 1, or return to the depot in order to refill the vehicle and resume the route to serve customer j + 1. The objective of the problem is to serve all customers (replenish their stock for all items) and minimize expected travel cost. As before, we declare nij the stochastic demand of customer j 2 {1, . . . , n} for product i. The demand per customer is no longer independent with respect to the product types, since for each customer j holds n1j þ n2j þ þ nkj 6 Q . 4.1. The minimum expected cost For simplicity, but without loss of generality, we will present our approach for the two product case. Let pj ðk1 ; k2 Þ ¼ Probðn1j ¼ n1k1 ; n2j ¼ n2k2 Þ represent the combined probability mass function. Also let fj(z1, z2) denote the minimum expected cost from customer j onward. Note that all product quantities are calculated using the same unit of measure e.g. m3 or kg. In this case there are additional issues to be considered. First, upon visiting the depot, an additional decision needs to be made regarding the quantities of stock to be loaded onto the vehicle. Let h be the quantity of product type-1 to be loaded onto the vehicle. Then in the 2product case, the quantity of product type-2 will be Q-h. Secondly, two sequential returns to the depot may occur in this problem. This may occur only if the vehicle does not proceed to client j + 1 directly but visits the depot first. In this case, if the quantity h for product type-1 (or Q-h for product type-2) loaded onto the vehicle at the depot is not adequate to fully satisfy the client’s demand, the vehicle will return to the depot once more, and make an additional, informed this time, decision. The stock s loaded to the vehicle for product type-1 (and Q-s for product type-2) will guarantee that the demand of customer j + 1 is fully satisfied. This is shown in Fig. 11 for the case where h < n1k1 . The dynamic programming formulation of the unified load problem is shown below
8 9 P P > cj;jþ1 þ fjþ1 ðz1 n1k1 ; z2 n2k2 Þpjþ1 ðk1 ; k2 Þ > > > > > 1k1 2k2 > > k1 :n 6z1 k2 :n 6z2 > > > > > " # > > > > > > P P > > 1k1 > > > > þ 2c þ min f ðh ðn z Þ; Q hÞ p ðk ; k Þ > 1 1 2 jþ1;0 jþ1 jþ1 > > > 1k > n 1 z1 6h6Q > = k1 :z1 > > " # > ðaÞ > > P P > > > > > > þ 2cjþ1;0 þ min fjþ1 ðh; Q h ðn2k2 z2 ÞÞ pjþ1 ðk1 ; k2 Þ > > > > > n2k2 z2 6h6Q > k1 :n1k1 6z1 k2 :z2 > > > > > " # > < > > P P > 1k1 2k2 fj ðz1 ; z2 Þ ¼ min > þ 2cjþ1;0 þ min fjþ1 ðh ðn z1 Þ; Q h ðn z2 ÞÞ pjþ1 ðk1 ; k2 Þ > ; > 1k1 2k2 > 1k 2k n z 6h6Q ðn z Þ > 1 2 k1 :z1 > > 9 8 P P > > cj;0 þ c0;jþ1 þ fjþ1 ðh n1k1 ; Q h n2k2 Þpjþ1 ðk1 ; k2 Þ > > > > > > > > > 1k1 2k2 > > > k1 :n 6h k2 :n 6Qh > > > > > > > > > > > > > P P = < > 1k1 > þ 2c þ min f ðs ðn hÞ; Q sÞ p ðk ; k Þ > 1 2 jþ1;0 jþ1 jþ1 > ðbÞ min 1k > n 1 h6s6Q > 06h6Q > > k1 :h > > > > > > > > > > > > P > > > > þ P > > > 2cjþ1;0 þ min fjþ1 fs; Q s ½n2k2 ðQ hÞg pjþ1 ðk1 ; k2 Þ > > ; : : k1 :n1k1 6h k2 :Qh
06s62Q hn2k2
with the boundary condition fn(z1, z2) = cn,0 z1, z2 2 Sj.
ξ 1k1
, ξ 2k 2 s - ( ξ 1k1 -θ), Q-s
j
j+1 z1,z2
0, Q-(θ - ξ 2k 2 ) θ , Q- θ
s, Q-s D
Fig. 11. The case that h < n1k1 z1 .
ð9Þ
566
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
In the above equation the first part in the minimization, term (a), represents the expected cost of going directly to the next customer, whereas the second part, term (b), represents the expected cost of the restocking action. Each term in parts (a) or (b) is determined following arguments, similar to those presented in Section 3.1. For example, consider the third term (in the second row) of part (a). This term represents the case in which upon arrival at customer j the stock on board for product z1 is not adequate to satisfy the customer’s demand, while the stock for product z2 is. In this case, the vehicle will have to fully satisfy the demand for z2 and then return to the depot, refill the capacity of both (z1, z2), and return to customer j in order to fully satisfy its demand for product z1. In order to identify the optimal refill loads for (z1, z2), an additional, internal this time, minimization is incorporated in this term. As an additional example, consider the second term (second row) of part (b). The path of the vehicle and the loads are shown in Fig. 11. This term represents the case in which after loading at the depot and upon arrival at customer j the stock on board for product-1, h, is not adequate to satisfy the customer’s demand, while the stock for product-2, Q-h, is. In this case, the vehicle will have to fully satisfy the demand for product-2, and then return to the depot, refill the capacity of both products, and return to customer j in order to fully satisfy its demand for product-1. As before, in order to identify the optimal refill loads for both products, an additional, internal minimization is incorporated in this term. 4.2. The optimal routing policy The value of the minimum expected cost could be estimated using Eq. (9). As in the previous case, the decision that guides the distribution operations is modelled by the variable xj(z1, z2). Furthermore, a threshold theorem analogous to Theorem 1 has been proven (see [24]). Proceeding in the same fashion as in Section 3.2 we first introduce a new lemma (Lemma 2) according to which
fj ðz1 ; z2 Þ 6 2c0j þ min fj ðh; Q hÞ for all z1 ; z2 2 Sj
ð10Þ
06h6Q
We have also shown that a theorem similar to Theorem 1 holds for the unified load case. j
Theorem 2. In the unified load case of Eq. (9), for each customer j and for each z1 6 Q, there exists a threshold h2 ðz1 Þ 6 Q z1 , such that the j optimal decision, after serving customer j, is to continue to customer j+1 if z2 > h2 ðz1 Þ, or return to the depot otherwise. The proof of the above theorem proceeds by showing that term (a) of Eq. (9) is a non-increasing function, and uses the fact that term (b) is independent of (z1, z2). The monotonicity of term (a) of Eq. (9) is shown by again segmenting the demand space appropriately and showing that for each segment of this space term (a) is non-increasing. Lemma 2 is employed to support showing the non-increasing property for some of these space segments. The full proof of Lemma 2 and Theorem 2 are given in reference [24]. Finally, for the complexity of the dynamic programming algorithm vs. the number of products, similar arguments with those discussed in Section 3.3 hold. 4.3. Implementation and example Similarly to the compartmentalized multiple product case, and based on the formulation presented in Section 4.1, we developed an appropriate algorithm that uses Dynamic Programming to derive the optimal solution of the unified load problem in a reasonable amount of time. Consider the five-customer network of Fig. 10. The total vehicle capacity (for both products) is Q = 10 units; the demand nij for each product i (i = 1, 2) per customer j (j = 1, . . . , 5) is given in Appendix B, and the distances between the nodes cjj0 are given in Fig. 10. Fig. 12 1 illustrates the threshold function h2 ðz1 Þ as obtained from the results for Customer 1. It has been calculated that H01 ¼ 67; 6. For any (z1, z2) combination below the border defined by the switch between the red and green points, the optimal decision is to return to the depot, x1(z1, z2) = 1, in order to refill the vehicle. Conversely, for any combination above this border the optimal decision would be to proceed directly to the next customer, x1(z1, z2) = 0. Again, this policy can be very simply and clearly communicated to the vehicle driver and result in significant cost savings for the fleet operator.
5. Algorithm performance 5.1. Compartmentalized load case The performance of the algorithm for the compartmentalized load case was found to be within acceptable levels. As an indication, for a test problem of 10 customer points, two products and vehicle capacity 10 units per product, the algorithm derived the minimum expected cost, as well as the threshold curves per customer, within 9 seconds. The number of combinations examined for each of the customer points 2–9 were approximately 1330. All experiments were run on a PC equipped with Intel Pentium IV, at 2.4 GHz and 512 MB of RAM. In order to further assess the performance of the algorithm for the two product compartmentalized load case a large number of problem test cases were created and run. We generated approximately 30,000 problems of appropriate characteristics and the results obtained are shown in Fig. 13. Three different problem test cases are presented in Fig. 13. The first test case (10,000 randomly generated problems) concerns a vehicle with total capacity Q = 10 equally split between the two compartments, and is shown in green1 colour. The second test case (another 10,000 randomly generated problems) concerns a vehicle with total capacity Q = 20 equally split between the two compartments, and is shown in blue colour. The third test case (another 10,000 randomly generated problems) concerns a vehicle with total capacity Q = 30 equally split between the two compartments, and is shown in red colour. Each point in these curves corresponds to the average solution time for 1000 ran-
1
For interpretation of colour in Figs. 13 and 14, the reader is referred to the web version of this article.
567
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Customer 1 11 Proceed to Customer 2
10
Proceed to Depot
9 8
Z2
7 6 5 4 3 2 1 0 1
0
2
3
4
5
6
7
8
9
10
11
Z1 Fig. 12. The driver decision based on the load combinations after serving customer 1.
domly generated problems. The demand distributions were generated randomly. From Fig. 13 it can be clearly seen that for a given number of customers, the increase of the vehicle capacity (Q = 10, 20, 30) results in a significant (almost exponential) increase in the computational time of the algorithm (note that the scale of the y-axis is logarithmic). On the other hand, if the capacity of the vehicle is kept constant, and the number of customers is increased, the computational time increases, almost linearly. 5.2. Unified load case The performance of the algorithm for the unified load case was also found to be within acceptable levels. As an indication, for a testproblem of 10 customer points, two products and total vehicle capacity of 10 units, the algorithm derived the minimum expected cost, as well as the threshold curves per customer within 1420 seconds (=23.6 minutes). The number of combinations examined for each of the customer points 2–9 were approximately 8600. One needs to take into consideration that this problem is significantly more complex than its compartmentalized counterpart (approximately 1330 combinations per customer) due to the additional minimization steps required to identify the optimal h and s values for each customer point. In order to further assess the performance of the algorithm, a large number of problem test cases were created and run (see Fig. 14). Three different problem test cases are plotted in this figure. The first test case (10,000 randomly generated problems) concerns a vehicle with a total capacity (sum for both products) Q = 5 and is shown in green colour. Each one of the 10 points of the curve is the average of 1000 randomly generated problems. The second test case (again 10,000 randomly generated problems) concerns a vehicle with a total capacity Q = 10 and is shown in blue colour. The third test case (10,000 randomly generated problems) included a vehicle with a total capacity Q = 15, and is shown in red colour. The demand distributions were generated randomly (uniform distributions). The trends obtained are similar to the ones corresponding to the compartmentalized case. For example, from Fig. 14 it can be clearly seen that for a given number of customers, the increase of the vehicle capacity resulted in an almost exponential increase in the compu-
2 Product case
Computation time (sec)
10000
1000
100
10 Q= 30 Q= 20 Q= 10
1 5
6
7
8
9
10
11
Customers Fig. 13. Performance results of the algorithm.
12
13
14
15
568
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
100000
Computation time (sec)
10000
1000
100
10
Q=15 Q =10 Q=5
1 5
6
7
9
8
10
11
12
13
14
15
Customers Fig. 14. Performance results of the algorithm for the unified load case.
tation time of the algorithm (note that the y-axis of the graph is in a logarithmic scale). On the other hand, if the capacity of the vehicle is kept constant, and the number of customers is increased, the computational time also increases. 6. Conclusions In this paper, we investigated two practical cases of the Stochastic Vehicle Routing with Depot Returns Problem (SVRDRP); the compartmentalized and the unified load case. The objective of both these cases was to serve all customers with a single vehicle and minimize the expected value of the travel cost. For both cases we presented the characteristics of each problem, a method to determine the minimum expected cost, and the theoretical results that permit us to determine the optimal decision after serving each customer. Both cases were formulated through dynamic programming, and for both it was proven that there exists an appropriate threshold function for each customer that distinguishes two regions in the space of possible loads after serving the customer: The region of loads for which the optimal decision (after serving the customer) is to return to the depot, and the region for which the optimal decision is to continue to the next customer. For both cases, through the execution of a large number of randomly generated problems, it was found that the increase of the vehicle capacity results in an almost exponential increase in the computational time of the algorithm. On the other hand, if the capacity of the vehicle is kept constant, and the number of customers is increased, the computational time increases (almost) linearly with the number of customers. The unified load case proved to be significantly more complex than the compartmentalized load one, as expected. This is mainly due to the increased complexity of the formulation, which includes additional minimization steps in order to identify the vehicle replenishment stock for each customer point.
Appendix A. Detailed expression of Eq. (8) Eq. (8) is expressed as follows: Hj ðz1 ; z2 Þ Hj ðz1 þ d1 ; z2 þ d2 Þ ¼
2 4cj;jþ1 þ
X
X
1k1
fjþ1 ðz1 n
; z2 n
2k2
Þp1jþ1;k1 p2jþ1;k2
cj;jþ1
k1 :n1k1 6z1 k2 :n2k2 6z2
2 þ4
h
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1
p1jþ1;k1 p2jþ1;k2
i
2
X
þ4
3
X
1k1
fjþ1 ðz1 þ d1 n
; z2 þ d2 n
2k2
Þp1jþ1;k1 p2jþ1;k2 5
k1 :n1k1 6z1 k2 :n2k2 6z2
X
X
X
X
X
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ
k1 :z1
X
h
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1 þd1
X
h
X
2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n
1k1
; Q 2Þ
i
3
2
p1jþ1;k1 p2jþ1;k2 5
þ4
k1 :n1k1 >z1 þd1 k2 :n2k2 6z2
p1jþ1;k1 p2jþ1;k2
X
X
X
k1 :n1k1 6z1 k2 :z2
X
k1 :n1k1 6z1 k2 :z2
3 fjþ1 ðz1 þ d1 n
1k1
2k2
; z2 þ d2 n
Þp1jþ1;k1 p2jþ1;k2 5
h
2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ
i
569
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
2
X
þ4
h i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
k1 :z1
X
X
3 1k1
fjþ1 ðz1 þ d1 n
; z2 þ d2 n
2k2
Þp1jþ1;k1 p2jþ1;k2 5
k1 :z1
2 þ4
X
h
X
k1 :z1 þd1
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
3 i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 5
k1 :n1k1 >z1 þd1 k2 :z2
2 þ4
X
h
X
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
k1 :n1k1 6z1 k2 :z2 þd2
X
h
X
2k2
2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n
Þ
i
3 p1jþ1;k1 p2jþ1;k2 5
k1 :n1k1 6z1 k2 :n2k2 >z2 þd2
2 þ4
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1
X
h
X
2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n
2k2
3 i 1 2 Þ pjþ1;k1 pjþ1;k2 5
k1 :z1z2 þd2
2 þ4
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1 þd1
X
X
h
2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n
1k1
2k2
; z2 þ d 2 þ Q 2 n
Þ
i
3 p1jþ1;k1 p2jþ1;k2 5
ðA:1Þ
k1 :n1k1 >z1 þd1 k2 :n2k2 >z2 þd2
Appendix B. The probability mass functions of Examples 3.3 and 4.3 Tables B.1 and B.2 provide the probability mass functions for clients 1–5 in the example of Section 3.3. Table B.3 provides the probability mass functions for clients 1–5 in the example of Section 4.3.
Table B.1 The probability mass function for product 1 Customer
1 2 3 4 5
Demand 0
1
2
3
4
5
0.01 0.01 0.01 0.01 0.01
0.7 0.6 0.7 0.5 0.2
0.1 0.2 0.2 0.2 0.6
0.1 0.1 0.05 0.2 0.1
0.08 0.08 0.03 0.08 0.08
0.01 0.01 0.01 0.01 0.01
Table B.2 The probability mass function for product 2 Customer
1 2 3 4 5 Vehicle capacity = (5.5).
Demand 0
1
2
3
4
5
0.01 0.01 0.01 0.01 0.01
0.1 0.2 0.03 0.8 0.08
0.1 0.6 0.05 0.1 0.1
0.7 0.1 0.2 0.05 0.2
0.08 0.08 0.7 0.03 0.6
0.01 0.01 0.01 0.01 0.01
570
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Table B.3 The probability mass functions for the unified load case Demand for product 1
Demand for product 2
Customer 1
Customer 2
Customer 3
Customer 4
Customer 5
10 9 9 8 8 8 7 7 7 7 6 6 6 6 6 5 5 5 5 5 5 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
0 1 0 2 1 0 3 2 1 0 4 3 2 1 0 5 4 3 2 1 0 6 5 4 3 2 1 0 7 6 5 4 3 2 1 0 8 7 6 5 4 3 2 1 0 9 8 7 6 5 4 3 2 1 0 10 9 8 7 6 5 4 3 2 1 0
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.7 0.01 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.1 0.21 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.01 0.6 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.001 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.7 0.01 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.1 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.002 0.01 0.01 0.6 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.6 0.001 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.101 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
Vehicle capacity = 10.
References [1] M. Reimann, K. Doerner, R.F. Hartl, D-ants: Savings based ants divide and conquer the vehicle routing problem, Computers & Operations Research 31 (4) (2003) 563–591. [2] N. Christofides, The traveling salesman problem, in: N. Christofides, R. Mingozzi, P. Toth, C. Sandi (Eds.), Combinatorial Optimization, Wiley, New York, 1979. [3] G. Cornuejols, G.L. Nemhauser, Tight bounds for Christofides traveling salesman heuristic, Mathematical Programming 14 (1) (1978) 116–121.
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571 [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
571
M. Gendreau, G. Laporte, A. Hertz, An approximation algorithm for the traveling salesman problem with Backhauls, Operations Research 45 (4) (1997) 639–641. U. Rembold, C. Blume, R. Dillmann, Computer Integrated Manufacturing Technology and Systems, Marcel Dekker, New York, 1985. G.B. Dantzig, J.H. Ramser, The truck dispatching problem, Management Science 6 (1) (1959) 80–91. G. Clarke, J. Wright, Scheduling of vehicles from a central depot to a number of delivery points, Operations Research 12 (4) (1964) 568–581. A.A. Assad, Modeling and implementation issues, in: B.L. Golden, A.A. Assad (Eds.), Vehicle Routing: Methods and Studies, North-Holland, Amsterdam, 1988. B.L. Golden, A.A. Assad (Eds.), Vehicle Routing: Methods and Studies, Studies in Management Science and Systems, vol. 16, Elsevier Science Publications, Amsterdam, 1988, pp. 7–45. G. Laporte, I.H. Osman, Routing problems: A bibliography, Annals of Operation Research 61 (1) (1995) 227–262. P. Toth, D. Vigo, The Vehicle Routing Problem, SIAM, Philadelphia, 2002. M. Gendreau, G. Laporte, R. Seguin, Stochastic vehicle routing, European Journal of Operational Research 88 (1996) 3–12. M. Dror, G. Laporte, P. Trudeau, Vehicle routing with stochastic demands: Properties and solution frameworks, Transportation Science 23 (3) (1989) 166–176. G. Laporte, F. Louveaux, H. Mercure, Models and exact solutions for a class of stochastic location-routing problems, European Journal of Operational Research 39 (1989) 71–78. W.R. Stewart, B.L. Golden, Stochastic vehicle routing: A comprehensive approach, European Journal of Operation Research 14 (1983) 371–385. D.J. Bertsimas, A vehicle routing problem with stochastic demand, Operations Research 40 (3) (1992) 574–585. P. Trudeau, M. Dror, Stochastic inventory routing: Route design with stockouts and route failures, Transportation Science 26 (3) (1992) 171–184. M. Dror, Modeling vehicle routing with uncertain demands as a stochastic program: Properties of the corresponding solution, European Journal of Operational Research 64 (3) (1993) 432–441. M. Birattari, C. Blum, M. Sampels, O. Rossi-Doria, M. Chiarandini, L. Pauqete, T. Schiavinotto, L. Bianchi, M. Mastrolilli, T. Bousonville, H. Juillie, D. Huber, Hybrid metaheuristics for the vehicle routing problem with stochastic demands, Technical Report (2005), IDSIA/USI-SUPSI, Dalle Molle Institute for Artificial Intelligence. N. Secomandi, Comparing neuro-dynamic programming algorithms for the vehicle routing problem with stochastic demands, Computers & Operations Research 27 (2000) 1201–1225. N. Secomandi, A rollout policy for the vehicle routing problem with stochastic demands, Operations Research 49 (5) (2001) 796–802. W.-H. Yang, K. Mathur, R.H. Ballou, Stochastic vehicle routing problem with restocking, Transportation Science 34 (2000) 99–112. P. Tsirimpas, A. Tatarakis, I. Minis, E.G. Kyriakidis, Single vehicle routing with a predefined customer sequence and multiple depot returns, European Journal of Operations Research 187 (2) (2008) 483–495. A. Tatarakis, A class of single vehicle routing problems with predefined customer sequence and depot returns, Ph.D. thesis, Department of Financial and Management Engineering, University of the Aegean, Chios, Greece, 2007.
Contents lists available at ScienceDirect
European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor
Production, Manufacturing and Logistics
Stochastic single vehicle routing with a predefined customer sequence and multiple depot returns A. Tatarakis, I. Minis * Department of Financial and Management Engineering, University of the Aegean, 31 Fostini Street, 82100 Chios, Greece
a r t i c l e
i n f o
Article history: Received 18 December 2007 Accepted 3 July 2008 Available online 17 July 2008 Keywords: Logistics Stochastic vehicle routing Multiple product delivery Dynamic programming
a b s t r a c t We study the routing of a single vehicle that delivers multiple products under stochastic demand. Specifically, we investigate two practical variations of this problem: (i) The case in which each product type is stored in its dedicated compartment in the vehicle, and (ii) the case in which all products are stored together in the vehicle’s single compartment. Suitable dynamic programming algorithms are proposed to determine the minimum expected (routing) cost for each case. Furthermore, the optimal routing policy is derived by developing appropriate theorems. The efficiency of the algorithms is studied by solving large problem sets. Ó 2008 Published by Elsevier B.V.
1. Introduction Both industry and academia have long recognized the potential for optimizing operations in goods distribution, especially since the latter may contribute up to 20% or more to the total costs of a product (see [1]). Formally, most problems in goods distribution are related to the Vehicle Routing Problem (VRP), a generalization of the classic Travelling Salesman problem (TSP) (see [2–4]). The VRP seeks a set of efficient vehicle routes to serve a number of geographically dispersed customers. In this paper, we address a special case of delivery: We consider a single vehicle that serves customers at a predefined sequence; the vehicle is allowed to return to the depot for stock replenishment. The customer demands are random and are modelled by independent variables with known statistics (derived from historical data). The distances, or travel times, among all points in the network (depot and customer sites) are known. The product quantities loaded onto the vehicle are limited by appropriate capacities, which cannot be exceeded. The vehicle returns to the depot in order to refill as needed. Upon completion of service at each customer, the driver has to make a decision: (a) Proceed to the next customer, as long as the probability of being able to fully satisfy its demand (for all products) is acceptable and the corresponding distance (i.e. the distance from the next customer to the depot) is favourable; (b) return to the depot in order to refill, and resume the route by visiting the next customer in the predefined sequence. If the vehicle proceeds to the next customer but the actual demand of this customer turns out to be higher than the stock carried on board (for any product), the vehicle will deliver as much as feasible, will then go to the depot to refill, and return to the customer to fully satisfy its demand. We refer to this problem as the Stochastic Vehicle Routing with Depot Returns Problem (SVRDRP). Practical applications of the SVRDRP may arise in different settings. A prime example is the Ex-van deliveries model, in which during each customer visit the driver simultaneously sells and delivers high-demand commodities to retail outlets. The concept behind Ex-van deliveries is to facilitate regular stock replenishment of certain types of commodities, so that the outlet (super market, kiosk, etc.) can maximise its sales. The mission of each Ex-van vehicle is to visit all assigned customer sites, typically in predefined sequence, and replenish the stock of selected products. The demand of each customer point is not known in advance but it is revealed upon arrival. Therefore, the total demand of a scheduled customer sequence typically exceeds the total capacity of the vehicle for a particular item, forcing the vehicle to return to the depot in order to replenish its own stock, before resuming its route. In addition to the Ex-van environment, similar delivery characteristics may also arise in other practical settings. For example, material handling systems in a manufacturing shop often operate along fixed pathways that connect the material warehouse with workcenters located along these pathways. Consider the case of Automated Guided Vehicle Systems (AGVs), which are self-propelled vehicles typically
* Corresponding author. Tel.: +30 2271035454. E-mail addresses: [email protected] (A. Tatarakis), [email protected] (I. Minis). 0377-2217/$ - see front matter Ó 2008 Published by Elsevier B.V. doi:10.1016/j.ejor.2008.07.006
558
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
guided along a magnetic induction strip, or a painted strip on the shop floor, and transport multiple discrete parts to workcenters, obviously in a predefined sequence. Note that in addition to the main pathway connecting the workcenters, there are spurs connecting each workcenter with the material warehouse, allowing the return and reloading of the AGV (see [5]). In this case, demand may also be stochastic, depending on the state of the buffers of the workcenters been served. The SVRDRP for a single product has been studied by Yang et al. [22]. In this paper, we study the SVRDRP for two cases of multiple product deliveries and develop (a) algorithms to determine the minimum expected routing cost and (b) policies that allow the vehicle’s driver to make the optimal decisions after serving each customer based on the level of the stock on board. The rest of the paper is organized as follows. In the next section, we present relevant research to date. In Section 3, we present the compartmentalized load case and in Section 4, we present the unified load case. In both these Sections, we illustrate the results obtained through appropriate examples. In Section 5, we study the efficiency of the proposed algorithms. The conclusions of this work are summarized in Section 6. 2. Background The Vehicle Routing Problem (VRP) was introduced by Dantzig and Ramser, almost 50 years ago (see [6]). The objective of the VRP is to deliver goods to a set of customers with known demands following minimum-cost vehicle routes originating from and terminating to a depot (see [7–10]). A rich survey of significant research results in this fundamental problem of Operations Research is given by Toth and Vigo [11]. The Stochastic Vehicle Routing Problem (SVRP) refers to a family of problems that combine the characteristics of stochastic and integer programs, and are often regarded as computationally intractable (see [12]). Therefore, only relatively small instances can be solved to optimality and effective heuristics are hard to design and assess. A characteristic problem is the VRP with Stochastic Demands (see [13,14]). In this problem, the demands are usually independent random variables that may (or may not) follow a known distribution. This problem has significant practical applications. A typical example is garbage collection, in which it is impossible to know a priori the quantity to be collected from each site. Another example, is the delivery of petrol to petrol stations. In this case, it is not known how much petrol will be sold due to the time elapsed between customer order and delivery. In order to address the inherent uncertainty in this type of VRP, a recourse action (i.e. return to the depot in order to refill) is usually embedded into the formulation of the problem, and penalties are incurred in the case of a route failure (see [15]). Due to the stochastic nature of this problem, the objective function is the expected value of the total route cost; (see [16–19]). Bertsimas (in [16]) constructs an a priori sequence among all customers of minimal expected total length and proposes suitable heuristics. His approach is shown to be a strong alternative to the strategy of re-optimization. More recently, Birattari et al. (see [19]), employed five metaheuristics (simulated annealing, tabu search, local search, ant colony optimisation, and evolutionary algorithms) and tested the effect of hybridization (TSPapproximation and VRPSD-approximation). The Stochastic Vehicle Routing Problem (SVRP) for a single product has been addressed in [16–19]. Secomandi (see [20,21]) has proposed a re-optimization approach, in which after the demand of each customer has been revealed the remaining part of the problem is re-solved by using Neuro-dynamic programming. This approach may yield better solutions than the preventive restocking strategy (returning to the depot before a stock out actually occurs), but it is much more computationally expensive. Moreover, the initially planned route may be altered completely (therefore, resembling more the SVRP that the SVRDRP), and this situation may present a limitation in practice. Yang et al. (see [22]) investigate the single and multi-vehicle single product SVRP, which has strong similarities with our case, as it entails a predefined customer visit sequence. Instead of adopting a simple recourse action, an optimal restocking policy of the vehicle has been incorporated in the route design. In particular, the restocking points are deliberately planned along the route, such that the probability of the route failure, and the accompanying recourse cost (including any penalty) is reduced, and the expected total routing cost is minimized. The single vehicle, multiple product case with deterministic demands has recently been presented in [23], in which three practical cases are studied: (i) The case in which each product type is stored in its own compartment in the vehicle, (ii) the case in which all product types are stored together in the vehicle’s single compartment, and (iii) the case of pickup and delivery. For each case, suitable dynamic programming algorithms to determine the optimal routing of the vehicle were developed, and their efficiency was studied by solving large problem sets. In the present paper, we propose and investigate two practical variations of the VRDRP with stochastic demands (SVRDRP): (i) the case of multiple product deliveries when each product is stored in its own compartment in the vehicle, and (ii) the case of multiple product deliveries when all products are stored together in the vehicle’s single compartment. 3. The compartmentalized load case We consider a set of nodes V = {0, 1, . . . , n}, with node 0 denoting the depot and nodes 1, . . . , n corresponding to customers, and a set of arcs A = {(j, j + 1):j, j + 1 2 V} [ {(0, j):j 2 V, j – 0} [ {(j, 0):j 2 V, j – 0} that join the customers along the route 1 ? 2 ? ? n, as well as all customers to the depot. The travel cost (distance) of each arc (j, j0 ) is denoted by cjj0 > 0. The cjj0 values satisfy the triangular inequality. In this problem, we assume that the delivery vehicle is divided into K sections and each section is suitable for carrying one product type P only. Let Qi be the capacity of the vehicle for product i 2 {1, . . . , K}. Clearly, Ki¼1 Q i ¼ Q the total load capacity of the vehicle. Note that all 3 product quantities are calculated using the same unit of measure e.g. m or kg. We define nij to be the stochastic demand of customer j 2 {1, . . . , n} for product type i 2 f1; . . . ; Kg:nij follows a discrete distribution with mi possible values, ni1, ni2, . . ., nimi and probability mass function pijk ¼ Pðnij ¼ nik Þ. The objective of the problem is for the vehicle to serve all customers according to the predefined sequence 1, . . . , n and minimize the expected travel cost. The vehicle is at the depot initially and after serving all customers it returns to the depot. 3.1. The minimum expected cost For simplicity, but without loss of generality, we will present our approach for two product types. Let fj(z1, z2) denote the total minimum expected cost from customer j onward if the vehicle, after serving customer j, carries quantities z1 of product type-1 and z2 of product type-
559
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
2. If Sj represents the set of all possible loads that a vehicle can carry after serving customer j, then fj(z1, z2) for (z1, z2) 2 Sj satisfies the following dynamic programming recursion:
8 P P fjþ1 ðz1 n1k1 ; z2 n2k2 Þp1jþ1;k1 p2jþ1;k2 cj;jþ1 þ > > > 1k > k1 :n 1 6z1 k2 :n2k2 6z2 > > i > > P h P > 1k1 > þ 2c þ f ðz þ Q n ; Q Þ p1jþ1;k1 p2jþ1;k2 > 1 jþ1;0 jþ1 1 2 > > > k1 :z1
9 > > > > > > > > > > > > = > > > > ðaÞ > > > > > > > > ; )
ð1Þ
k1 ¼1 k2 ¼1
ðbÞ with the boundary condition fn(z1, z2) = cn,0 z1, z2 2 Sj. In Eq. (1) the minimum is between the two terms; the top term, (a), represents the expected cost of going directly to the next customer. The bottom term, (b), represents the expected cost of the restocking action. Part (a) consists of four summation terms (in the four rows), which correspond to the four areas of the z1, z2 space shown in Fig. 1. The first summation term of (a) represents the expected cost incurred if the stock of both items z1 and z2 is sufficient to fully satisfy the demands n1k1 and n2k2 of customer j + 1 (Area A in Fig. 1). In Fig. 2 the loads of the vehicle before and after serving customer j + 1 are schematically represented above the respective arcs (the demand of customer j + 1 is presented above the node). The stock left onboard after serving customer j + 1 is ðz1 n1k1 Þ P 0 and ðz2 n2k2 Þ P 0, respectively. The second summation of part (a) represents the expected cost incurred if the vehicle proceeds to the next customer directly and the stock z1 is not sufficient to fully satisfy the demand n1k1 of the next customer, while the stock z2 is sufficient to fully satisfy the demand n2k2 (Area B in Fig. 1). The path of the vehicle and the corresponding vehicle loads are shown in Fig. 3. In this case, the vehicle will visit customer j + 1, it will fully satisfy demand n2k2 for product type-2 but will partially satisfy demand n1k1 for product type-1; thus, it is required to return to the depot for stock replenishment. The stock left on board after serving customer j + 1 for
ξ 2k 2 Q2 C
D
A
B
z2
Q1
z1
ξ 1k1
Fig. 1. The solution space per customer point.
ξ 1k1 , ξ 2k 2
j
z1,z2
z1 - ξ 1k1 , z2 - ξ 2k 2
j+1
Fig. 2. The case that the demand is fully satisfied: z1 P n1k1 ; z2 P n2k2 .
ξ 1k
z1,z2 j
0, z2 - ξ 2k
, ξ 2k
Q1 – ( ξ 1k - z1), Q2
j+1 Q1, Q2
D
Fig. 3. The case in which z1 < n1k1 and n2k2 6 z2 .
560
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
ξ 1k
z1,z2 j
1k 2k , ξ 2k Q1 – ( ξ - z1), Q2 – ( ξ - z2)
j+1
0,0 Q1, Q2
D
Fig. 4. The case in both z1 < n1k1 and z2 < n2k2 .
ξ 1k
z1,z2
j
, ξ 2k
Q1 - ξ 1k , Q2 – ξ 2k
j+1 Q1, Q2
D
Fig. 5. The proactive depot return case.
the first time is 0 and ðz2 n2k2 Þ P 0, respectively. The stock left on board after serving customer j + 1 for the second time is Q 1 ðn1k1 z1 Þ for product type-1 and Q2 for product type-2, respectively. The third sum of part (a) represents the expected cost incurred if the stock z2 is not sufficient to fully satisfy the demand n2k2 of the next customer, while stock z1 is sufficient to satisfy demand n1k1 (Area C in Fig. 1). The path of the vehicle and the corresponding vehicle loads are analogous to those of the previous case. The fourth term represents the expected cost incurred if the stock of both items z1 and z2 is not sufficient to fully satisfy the demands n1k1 or n2k2 of the next customer (Area D in Fig. 1). The path of the vehicle and the corresponding vehicle loads are shown in Fig. 4. Part (b) of the minimization equation, represents the expected cost incurred if the vehicle proceeds to the next customer j + 1 via the depot and, therefore, the stock of both items on board z1 and z2 is replenished to Q1 and Q2, respectively. Due to the fact that this is a proactive depot return, the values of z1 and z2 do not affect the result; i.e. the latter is independent of z1 and z2. The path of the vehicle and the corresponding vehicle loads are shown in Fig. 5. 3.2. The optimal routing policy The value of the minimum expected cost can be estimated using Eq. (1). However, the decision required to guide the distribution operation in the SVRDRP environment is the following: Having served customer j, and if the vehicle remaining load is (z1, z2), which is the preferred action; (a) to proceed to customer j + 1, or (b) to visit the depot prior to serving customer j + 1. This decision can be modelled by the following variable: 0
if part (a) ≤ part (b) in Eq. (1); i.e. proceed to customer j+1
1
if part (a) > part (b) in Eq. (1); i.e. proceed to the depot
xj(z1, z2) =
Thus, the vehicle proceeds to the next customer if the corresponding expected cost (until the end of the route) is minimum; otherwise the route via the depot is preferred. This decision variable is also used in the deterministic counterpart of this problem (see [23]). In order to develop the policy to obtain this minimum expected cost (in an average sense) we will extend the threshold theorem presented by Yang et al. [22] to multiple dimensions. Lemma 1
fj ðz1 ; z2 Þ 6 fj ðQ 1 ; Q 2 Þ þ 2c0j
for all z1 ; z2 2 Sj :
ð2Þ
Proof. From Eq. (1) we obtain Q1 X
Q2 X
k1 :n1k1 ¼1
k2 :n2k2 ¼1
fj ðz1 ; z2 Þ 6 cj;0 þ c0;jþ1 þ
fjþ1 ðQ 1 n1k1 ; Q 2 n2k2 Þp1jþ1;k1 p2jþ1;k2 :
ð3Þ
For z1 = Q1 and z2 = Q2 part (a) of Eq. (1) is always less than or equal to part (b) since cj,j+1 6 cj,0 + c0,j+1; therefore
fj ðQ 1 ; Q 2 Þ ¼ cj;jþ1 þ
Q1 X
Q2 X
fjþ1 ðQ 1 n1k1 ; Q 2 n2k2 Þp1jþ1;k1 p2jþ1;k2
k1 :n1k1 ¼1 k2 :n2k2 ¼1
taking into consideration that the last three terms of part (a) of Eq. (1) become zero.
ð4Þ
561
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Substituting Eq. (4) in Eq. (3) results in
fj ðz1 ; z2 Þ 6 cj;0 þ c0;jþ1 þ ½fj ðQ 1 ; Q 2 Þ cj;jþ1 6 cj;0 þ ½c0;j þ cj;jþ1 þ ½fj ðQ 1 ; Q 2 Þ cj;jþ1 6 2cj;0 þ fj ðQ 1 ; Q 2 Þ
As before, consider z1 and z2, the item quantities on board after serving customer j. j
Theorem 1. For each customer j, there exists a threshold function h2 ðz1 Þ, such that the optimal decision, after serving customer j, is to continue to j customer j + 1 if z2 > h2 ðz1 Þ or return to the depot otherwise. Proof. We will first show by induction that for all (z1, z2) 2 Sj, fj(z1, z2) is a non-increasing function. That is, for z1, z2 2 Sj and d1,d2 P 0
fj ðz1 þ d1 ; z2 þ d2 Þ 6 fj ðz1 ; z2 Þ:
ð5Þ
This relationship is true for the last customer n, where fn(z1, z2) = cn0 is independent of (z1, z2). Hence fn(z1, z2) is non-increasing with respect to (z1, z2) 2 Sn. We will now prove that, if fj+1(z1, z2) is non-increasing with respect to (z1, z2) 2 Sj+1, then fj(z1, z2) is also non-increasing with respect to (z1, z2) 2 Sj. Let Hj(z1, z2) and H0j ðz1 ; z2 Þ denote the values of part (a) and part (b) inside the minimisation in Eq. (1). Furthermore, let
Hj ðz1 ; z2 Þ ¼ Haj ðz1 ; z2 Þ þ Hbj ðz1 ; z2 Þ þ Hcj ðz1 ; z2 Þ þ Hdj ðz1 ; z2 Þ;
ð6Þ
where (see also Fig. 1):
1k1
Hbj ðz1 ; z2 Þ
X
¼
X
X
Haj ðz1 ; z2 Þ ¼ cj;jþ1 þ
k1 :n
h
X
k1 :z1
Hcj ðz1 ; z2 Þ
X
¼
X
k1 :n1k1 6z1 k2 :z2
X
Hdj ðz1 ; z2 Þ ¼
fjþ1 ðz1 n1k1 ; z2 n2k2 Þp1jþ1;k1 p2jþ1;k2 ;
6z1 k2 :n2k2 6z2
X
h
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 ;
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 ;
h
ð6:1Þ
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :
k1 :z1
If we expand each one of the last three terms taking into account the regions of Fig. 6 we obtain
X
Hbj ðz1 ; z2 Þ ¼
h i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
X
k1 :z1
X
þ
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 ;
ð6:2Þ
k1 :z1 þd1
X
Hcj ðz1 ; z2 Þ ¼
h i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
k1 :n1k1 6z1 k2 :z2
þ
X
X
h
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 ;
ξ
2k
ð6:3Þ
k1 :n1k1 6z1 k2 :z2 þd2
2
Q2 (c2)
(d4)
(d3)
(c1)
(d1)
(d2)
z2+δ 2
z2 (b1)
(a)
z1
(b2)
z1+δ 1
Fig. 6. The definition space of Hj(z1, z2).
Q1
ξ
1k 1
562
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Hdj ðz1 ; z2 Þ ¼ þ
h i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
X
k1 :z1
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
X
k1 :z1 þd1
þ þ
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
2k2 k1 :z1
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :
X
ð6:4Þ
k1 :z1 þd1
Similarly Hj(z1 + d1, z2 + d2) can be written using Fig. 7 as follows: 0
0
0
0
Hj ðz1 þ d1 ; z2 þ d2 Þ ¼ Haj ðz1 þ d1 ; z2 þ d2 Þ þ Hbj ðz1 þ d1 ; z2 þ d2 Þ þ Hcj ðz1 þ d1 ; z2 þ d2 Þ þ Hdj ðz1 þ d1 ; z2 þ d2 Þ:
ð7Þ
Again, if we expand each one of the above terms taking into account the regions of Fig. 7 we obtain
X
X
k1 :n1k1 6z1
2k2
0
Haj ðz1 þ d1 ; z2 þ d2 Þ ¼ cj;jþ1 þ
X
þ
k2 :n
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2
X6z2
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2
k1 :z1
þ
X 1k1
X
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2
6z1 k2 :z2
k1 :n
X
þ
X
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2 ;
ð7:1Þ
k1 :z1
X
0
Hbj ðz1 þ d1 ; z2 þ d2 Þ ¼
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
k1 :n1k1 >z1 þd1 k2 :n2k2 6z2
X
þ
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 ;
ð7:2Þ
k1 :n1k1 >z1 þd1 k2 :z2
X
X
0
Hcj ðz1 þ d1 ; z2 þ d2 Þ ¼
k1 :n1k1 6z1 k2 :n2k2 >z2 þd2
þ
X
h
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 h
X
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 ;
ð7:3Þ
k1 :z1
Hdj ðz1 þ d1 ; z2 þ d2 Þ ¼
X
X
h i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; z2 þ d2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :
ð7:4Þ
k1 :n1k1 >z1 þd1 k2 :n2k2 >z2 þd2
Subtracting the terms that correspond to the same nine regions of Figs. 6 and 7 we obtain
h i h i a0 a0 Hj ðz1 ; z2 Þ Hj ðz1 þ d1 ; z2 þ d2 Þ ¼ Haj ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ þ Hbj 1 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ h i h i a0 b0 þ Hbj 2 ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ þ Hcj 1 ðz1 ; z2 Þ Hj 4 ðz1 þ d1 ; z2 þ d2 Þ h i h i a0 b0 þ Hdj 1 ðz1 ; z2 Þ Hj 3 ðz1 þ d1 ; z2 þ d2 Þ þ Hdj 2 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ h i h i c0 c0 þ Hcj 2 ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ þ Hdj 4 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ h i 0 þ Hdj 3 ðz1 ; z2 Þ Hdj ðz1 þ d1 ; z2 þ d2 Þ : ξ
2k
2
Q2 (c′1)
(c′2)
(d′)
(a′4)
(a′3)
(b′2)
(a′1)
(a′2)
(b′1)
z2+δ 2
z2
z1
z1+δ 1
Fig. 7. The definition space of Hj (z1 + d1, z2 + d2).
Q1
ξ
1k 1
ð8Þ
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
563
The detailed expression of Eq. (8) is given in Appendix A. In Eqs. (8) and (A.1), we can distinguish two types of terms. The first type includes the first, third, sixth, seventh, eighth, and ninth terms. Considering the first term we have a0
Haj ðz1 ; z2 Þ Hj 1 ðz1 þ d1 ; z2 þ d2 Þ ¼
X
X
½fjþ1 ðz1 n1k1 ; z2 n2k2 Þ fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þp1jþ1;k1 p2jþ1;k2 :
k1 :n1k1 6z1 k2 :n2k2 6z2
Since fj+1(z1, z2) is monotonically non-increasing, this term is non-negative. Similarly we can show that the differences corresponding to the third, sixth, seventh, eighth, and ninth terms are also non-negative. The second type of terms includes the second, forth and fifth terms. Considering the second term we have a0
Hbj 1 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ i X X h 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ p1jþ1;k1 p2jþ1;k2 :1 ¼ k1 :z1
Recall that from Lemma 1:
fj ðz1 ; z2 Þ 6 fj ðQ 1 ; Q 2 Þ þ 2c0j
for all z 2 Sj ;
therefore,
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ 6 fjþ1 ðQ 1 ; Q 2 Þ þ 2cjþ1;0 : However,
fjþ1 ðQ 1 ; Q 2 Þ 6 fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ since n1k1 > z1 and fj+1 is non-increasing. Thus,
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ 6 fjþ1 ðQ 1 ; Q 2 Þ þ 2cjþ1;0 6 fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ þ 2cjþ1;0 : a0
Therefore, Hbj 1 ðz1 ; z2 Þ Hj 2 ðz1 þ d1 ; z2 þ d2 Þ P 0. Similarly we can show that the forth and fifth terms of Eq. (8) are also non-negative. Considering the non-negativity of all terms in Eq. (8) it is clear that Hj(z1, z2) Hj(z1 + d1,z2 + d2) P 0 and hence, Hj(z1, z2) is a nonincreasing function. Note now that part (b) of Eq. (1), H0j ðz1 ; z2 Þ, is independent of z1, z2 and, thus, a constant H0j ðz1 ; z2 Þ ¼ H0j in the z1, z2 space. Fig. 8 plots the terms Hj(z1, z2n) and H0j with o Hj(z1, z2) been non-increasing as shown above. From this Figure it is clear that fj ðz1 ; z2 Þ ¼ min Hj ðz1 ; z2 Þ; H0j ðz1 ; z2 Þ is a non-increasing function. Furthermore, it is also clear that due to intersection Hj ðz1 ; z2 Þ \ H0j for j j customer j and for each value of z1 there exists a value h2 , such that for z2 > h2 Hj ðz1 ; z2 Þ < H0j and, thus, the optimal decision after serving customer j is to continue to customer j + 1. Otherwise a return to the depot is preferable (in the mean value sense). This concludes the proof of Theorem 1. h
Fig. 8. The combined threshold graphical representation.
564
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
3.3. Complexity vs. the number of products The influence of the number of product types to the complexity of the problem is significant, and worth noting. As pointed out in Section 2, the formulation for the single product counterpart of the problem presented in this paper is given in Yang et al. (see [22]). According to this reference we have
fj ðzÞ ¼ min
8 P P fjþ1 ðz nk Þpjþ1;k þ ½b þ 2cjþ1;0 þ fjþ1 ðz þ Q nk Þpjþ1;k cj;jþ1 þ > > < k k:n 6z k:nk >z
ðaÞ
m P > > : cj;0 þ c0;jþ1 þ fjþ1 ðQ nk Þpjþ1;k
ðbÞ
k¼1
with the boundary condition
fn ðzÞ ¼ cn0 z 2 Sn : The complexity increase is evident by comparing the above equation to Eq. (1). In Eq. (1), the number of terms to be evaluated for each part, (a) or (b), is m1m2 while in the single product case, the corresponding number of terms is m. It is, thus, clear that the problem complexity increases as the space of possible loads increases. As a consequence, for example, in the three product case, the number of terms increases to m1m2m3 where m3 is the number of possible quantities for product i = 3. 3.4. Implementation and example In order to solve the compartmentalized case of multiple product delivery, we developed an appropriate algorithm that uses Dynamic Programming to derive the optimal solution in a reasonable amount of time. Based on the formulation presented in Section 3.1, the algorithm starts from the end of the route (last customer to be visited) and iterates towards the beginning of the route, calculating the remaining minimum expected cost from each customer site until the end of the route. This procedure computes the minimum expected cost of the route, given a distance matrix and the demand probability mass functions. Based on the result of the algorithm, the threshold function h2(z1) for each customer j can be obtained, in line with what has already been described in Section 3.2. Consider the five-customer network of Fig. 9. The vehicle capacity is Q = 10 units and is equally split between two products (Q1 = Q2 = 5); the demand nij for each product i (i = 1, 2) and customer j (j = 1, . . . , 5) is given in Appendix B, and the distances between the nodes cjj0 are given in Fig. 10.
Fig. 9. Five-customer network for the multiple product extension.
Customer1 5
4
3
z2
h12 (3)
2
1
No return Return
0 0
1
2
3
z1
4
5
Fig. 10. The driver decision based on the load combinations after serving customer 1.
565
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571 1
The problem is solved using the dynamic programming algorithm. Fig. 10 illustrates the threshold function h2 ðz1 Þ as obtained from the results for customer 1. In this case H01 ¼ 70; 2. For any (z1, z2) combination below the border defined by the switch between the red and green points, the optimal decision is to return to the depot in order to refill the vehicle. Conversely, for any combinations above this border the optimal decision would be to proceed directly to the next customer. This policy can be very simply and clearly communicated to the vehicle driver and result in significant cost savings for the fleet operator. 4. The unified load case In the unified load case the vehicle may carry any quantity of product i 2 {1, . . . , K}, provided that the total capacity Q of the vehicle is not exceeded. Upon completion of service at customer site j, the vehicle’s driver has to make the same decision as the one described in Section 3, i.e. proceed to customer j + 1, or return to the depot in order to refill the vehicle and resume the route to serve customer j + 1. The objective of the problem is to serve all customers (replenish their stock for all items) and minimize expected travel cost. As before, we declare nij the stochastic demand of customer j 2 {1, . . . , n} for product i. The demand per customer is no longer independent with respect to the product types, since for each customer j holds n1j þ n2j þ þ nkj 6 Q . 4.1. The minimum expected cost For simplicity, but without loss of generality, we will present our approach for the two product case. Let pj ðk1 ; k2 Þ ¼ Probðn1j ¼ n1k1 ; n2j ¼ n2k2 Þ represent the combined probability mass function. Also let fj(z1, z2) denote the minimum expected cost from customer j onward. Note that all product quantities are calculated using the same unit of measure e.g. m3 or kg. In this case there are additional issues to be considered. First, upon visiting the depot, an additional decision needs to be made regarding the quantities of stock to be loaded onto the vehicle. Let h be the quantity of product type-1 to be loaded onto the vehicle. Then in the 2product case, the quantity of product type-2 will be Q-h. Secondly, two sequential returns to the depot may occur in this problem. This may occur only if the vehicle does not proceed to client j + 1 directly but visits the depot first. In this case, if the quantity h for product type-1 (or Q-h for product type-2) loaded onto the vehicle at the depot is not adequate to fully satisfy the client’s demand, the vehicle will return to the depot once more, and make an additional, informed this time, decision. The stock s loaded to the vehicle for product type-1 (and Q-s for product type-2) will guarantee that the demand of customer j + 1 is fully satisfied. This is shown in Fig. 11 for the case where h < n1k1 . The dynamic programming formulation of the unified load problem is shown below
8 9 P P > cj;jþ1 þ fjþ1 ðz1 n1k1 ; z2 n2k2 Þpjþ1 ðk1 ; k2 Þ > > > > > 1k1 2k2 > > k1 :n 6z1 k2 :n 6z2 > > > > > " # > > > > > > P P > > 1k1 > > > > þ 2c þ min f ðh ðn z Þ; Q hÞ p ðk ; k Þ > 1 1 2 jþ1;0 jþ1 jþ1 > > > 1k > n 1 z1 6h6Q > = k1 :z1
06s62Q hn2k2
with the boundary condition fn(z1, z2) = cn,0 z1, z2 2 Sj.
ξ 1k1
, ξ 2k 2 s - ( ξ 1k1 -θ), Q-s
j
j+1 z1,z2
0, Q-(θ - ξ 2k 2 ) θ , Q- θ
s, Q-s D
Fig. 11. The case that h < n1k1 z1 .
ð9Þ
566
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
In the above equation the first part in the minimization, term (a), represents the expected cost of going directly to the next customer, whereas the second part, term (b), represents the expected cost of the restocking action. Each term in parts (a) or (b) is determined following arguments, similar to those presented in Section 3.1. For example, consider the third term (in the second row) of part (a). This term represents the case in which upon arrival at customer j the stock on board for product z1 is not adequate to satisfy the customer’s demand, while the stock for product z2 is. In this case, the vehicle will have to fully satisfy the demand for z2 and then return to the depot, refill the capacity of both (z1, z2), and return to customer j in order to fully satisfy its demand for product z1. In order to identify the optimal refill loads for (z1, z2), an additional, internal this time, minimization is incorporated in this term. As an additional example, consider the second term (second row) of part (b). The path of the vehicle and the loads are shown in Fig. 11. This term represents the case in which after loading at the depot and upon arrival at customer j the stock on board for product-1, h, is not adequate to satisfy the customer’s demand, while the stock for product-2, Q-h, is. In this case, the vehicle will have to fully satisfy the demand for product-2, and then return to the depot, refill the capacity of both products, and return to customer j in order to fully satisfy its demand for product-1. As before, in order to identify the optimal refill loads for both products, an additional, internal minimization is incorporated in this term. 4.2. The optimal routing policy The value of the minimum expected cost could be estimated using Eq. (9). As in the previous case, the decision that guides the distribution operations is modelled by the variable xj(z1, z2). Furthermore, a threshold theorem analogous to Theorem 1 has been proven (see [24]). Proceeding in the same fashion as in Section 3.2 we first introduce a new lemma (Lemma 2) according to which
fj ðz1 ; z2 Þ 6 2c0j þ min fj ðh; Q hÞ for all z1 ; z2 2 Sj
ð10Þ
06h6Q
We have also shown that a theorem similar to Theorem 1 holds for the unified load case. j
Theorem 2. In the unified load case of Eq. (9), for each customer j and for each z1 6 Q, there exists a threshold h2 ðz1 Þ 6 Q z1 , such that the j optimal decision, after serving customer j, is to continue to customer j+1 if z2 > h2 ðz1 Þ, or return to the depot otherwise. The proof of the above theorem proceeds by showing that term (a) of Eq. (9) is a non-increasing function, and uses the fact that term (b) is independent of (z1, z2). The monotonicity of term (a) of Eq. (9) is shown by again segmenting the demand space appropriately and showing that for each segment of this space term (a) is non-increasing. Lemma 2 is employed to support showing the non-increasing property for some of these space segments. The full proof of Lemma 2 and Theorem 2 are given in reference [24]. Finally, for the complexity of the dynamic programming algorithm vs. the number of products, similar arguments with those discussed in Section 3.3 hold. 4.3. Implementation and example Similarly to the compartmentalized multiple product case, and based on the formulation presented in Section 4.1, we developed an appropriate algorithm that uses Dynamic Programming to derive the optimal solution of the unified load problem in a reasonable amount of time. Consider the five-customer network of Fig. 10. The total vehicle capacity (for both products) is Q = 10 units; the demand nij for each product i (i = 1, 2) per customer j (j = 1, . . . , 5) is given in Appendix B, and the distances between the nodes cjj0 are given in Fig. 10. Fig. 12 1 illustrates the threshold function h2 ðz1 Þ as obtained from the results for Customer 1. It has been calculated that H01 ¼ 67; 6. For any (z1, z2) combination below the border defined by the switch between the red and green points, the optimal decision is to return to the depot, x1(z1, z2) = 1, in order to refill the vehicle. Conversely, for any combination above this border the optimal decision would be to proceed directly to the next customer, x1(z1, z2) = 0. Again, this policy can be very simply and clearly communicated to the vehicle driver and result in significant cost savings for the fleet operator.
5. Algorithm performance 5.1. Compartmentalized load case The performance of the algorithm for the compartmentalized load case was found to be within acceptable levels. As an indication, for a test problem of 10 customer points, two products and vehicle capacity 10 units per product, the algorithm derived the minimum expected cost, as well as the threshold curves per customer, within 9 seconds. The number of combinations examined for each of the customer points 2–9 were approximately 1330. All experiments were run on a PC equipped with Intel Pentium IV, at 2.4 GHz and 512 MB of RAM. In order to further assess the performance of the algorithm for the two product compartmentalized load case a large number of problem test cases were created and run. We generated approximately 30,000 problems of appropriate characteristics and the results obtained are shown in Fig. 13. Three different problem test cases are presented in Fig. 13. The first test case (10,000 randomly generated problems) concerns a vehicle with total capacity Q = 10 equally split between the two compartments, and is shown in green1 colour. The second test case (another 10,000 randomly generated problems) concerns a vehicle with total capacity Q = 20 equally split between the two compartments, and is shown in blue colour. The third test case (another 10,000 randomly generated problems) concerns a vehicle with total capacity Q = 30 equally split between the two compartments, and is shown in red colour. Each point in these curves corresponds to the average solution time for 1000 ran-
1
For interpretation of colour in Figs. 13 and 14, the reader is referred to the web version of this article.
567
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Customer 1 11 Proceed to Customer 2
10
Proceed to Depot
9 8
Z2
7 6 5 4 3 2 1 0 1
0
2
3
4
5
6
7
8
9
10
11
Z1 Fig. 12. The driver decision based on the load combinations after serving customer 1.
domly generated problems. The demand distributions were generated randomly. From Fig. 13 it can be clearly seen that for a given number of customers, the increase of the vehicle capacity (Q = 10, 20, 30) results in a significant (almost exponential) increase in the computational time of the algorithm (note that the scale of the y-axis is logarithmic). On the other hand, if the capacity of the vehicle is kept constant, and the number of customers is increased, the computational time increases, almost linearly. 5.2. Unified load case The performance of the algorithm for the unified load case was also found to be within acceptable levels. As an indication, for a testproblem of 10 customer points, two products and total vehicle capacity of 10 units, the algorithm derived the minimum expected cost, as well as the threshold curves per customer within 1420 seconds (=23.6 minutes). The number of combinations examined for each of the customer points 2–9 were approximately 8600. One needs to take into consideration that this problem is significantly more complex than its compartmentalized counterpart (approximately 1330 combinations per customer) due to the additional minimization steps required to identify the optimal h and s values for each customer point. In order to further assess the performance of the algorithm, a large number of problem test cases were created and run (see Fig. 14). Three different problem test cases are plotted in this figure. The first test case (10,000 randomly generated problems) concerns a vehicle with a total capacity (sum for both products) Q = 5 and is shown in green colour. Each one of the 10 points of the curve is the average of 1000 randomly generated problems. The second test case (again 10,000 randomly generated problems) concerns a vehicle with a total capacity Q = 10 and is shown in blue colour. The third test case (10,000 randomly generated problems) included a vehicle with a total capacity Q = 15, and is shown in red colour. The demand distributions were generated randomly (uniform distributions). The trends obtained are similar to the ones corresponding to the compartmentalized case. For example, from Fig. 14 it can be clearly seen that for a given number of customers, the increase of the vehicle capacity resulted in an almost exponential increase in the compu-
2 Product case
Computation time (sec)
10000
1000
100
10 Q= 30 Q= 20 Q= 10
1 5
6
7
8
9
10
11
Customers Fig. 13. Performance results of the algorithm.
12
13
14
15
568
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
100000
Computation time (sec)
10000
1000
100
10
Q=15 Q =10 Q=5
1 5
6
7
9
8
10
11
12
13
14
15
Customers Fig. 14. Performance results of the algorithm for the unified load case.
tation time of the algorithm (note that the y-axis of the graph is in a logarithmic scale). On the other hand, if the capacity of the vehicle is kept constant, and the number of customers is increased, the computational time also increases. 6. Conclusions In this paper, we investigated two practical cases of the Stochastic Vehicle Routing with Depot Returns Problem (SVRDRP); the compartmentalized and the unified load case. The objective of both these cases was to serve all customers with a single vehicle and minimize the expected value of the travel cost. For both cases we presented the characteristics of each problem, a method to determine the minimum expected cost, and the theoretical results that permit us to determine the optimal decision after serving each customer. Both cases were formulated through dynamic programming, and for both it was proven that there exists an appropriate threshold function for each customer that distinguishes two regions in the space of possible loads after serving the customer: The region of loads for which the optimal decision (after serving the customer) is to return to the depot, and the region for which the optimal decision is to continue to the next customer. For both cases, through the execution of a large number of randomly generated problems, it was found that the increase of the vehicle capacity results in an almost exponential increase in the computational time of the algorithm. On the other hand, if the capacity of the vehicle is kept constant, and the number of customers is increased, the computational time increases (almost) linearly with the number of customers. The unified load case proved to be significantly more complex than the compartmentalized load one, as expected. This is mainly due to the increased complexity of the formulation, which includes additional minimization steps in order to identify the vehicle replenishment stock for each customer point.
Appendix A. Detailed expression of Eq. (8) Eq. (8) is expressed as follows: Hj ðz1 ; z2 Þ Hj ðz1 þ d1 ; z2 þ d2 Þ ¼
2 4cj;jþ1 þ
X
X
1k1
fjþ1 ðz1 n
; z2 n
2k2
Þp1jþ1;k1 p2jþ1;k2
cj;jþ1
k1 :n1k1 6z1 k2 :n2k2 6z2
2 þ4
h
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1
p1jþ1;k1 p2jþ1;k2
i
2
X
þ4
3
X
1k1
fjþ1 ðz1 þ d1 n
; z2 þ d2 n
2k2
Þp1jþ1;k1 p2jþ1;k2 5
k1 :n1k1 6z1 k2 :n2k2 6z2
X
X
X
X
X
fjþ1 ðz1 þ d1 n1k1 ; z2 þ d2 n2k2 Þ
k1 :z1
X
h
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1 þd1
X
h
X
2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n
1k1
; Q 2Þ
i
3
2
p1jþ1;k1 p2jþ1;k2 5
þ4
k1 :n1k1 >z1 þd1 k2 :n2k2 6z2
p1jþ1;k1 p2jþ1;k2
X
X
X
k1 :n1k1 6z1 k2 :z2
X
k1 :n1k1 6z1 k2 :z2
3 fjþ1 ðz1 þ d1 n
1k1
2k2
; z2 þ d2 n
Þp1jþ1;k1 p2jþ1;k2 5
h
2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ
i
569
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
2
X
þ4
h i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
X
k1 :z1
X
X
3 1k1
fjþ1 ðz1 þ d1 n
; z2 þ d2 n
2k2
Þp1jþ1;k1 p2jþ1;k2 5
k1 :z1
2 þ4
X
h
X
k1 :z1 þd1
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
3 i 2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n1k1 ; Q 2 Þ p1jþ1;k1 p2jþ1;k2 5
k1 :n1k1 >z1 þd1 k2 :z2
2 þ4
X
h
X
i 2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
k1 :n1k1 6z1 k2 :z2 þd2
X
h
X
2k2
2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n
Þ
i
3 p1jþ1;k1 p2jþ1;k2 5
k1 :n1k1 6z1 k2 :n2k2 >z2 þd2
2 þ4
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1
X
h
X
2cjþ1;0 þ fjþ1 ðQ 1 ; z2 þ d2 þ Q 2 n
2k2
3 i 1 2 Þ pjþ1;k1 pjþ1;k2 5
k1 :z1
2 þ4
X
h
X
i 2cjþ1;0 þ fjþ1 ðz1 þ Q 1 n1k1 ; z2 þ Q 2 n2k2 Þ p1jþ1;k1 p2jþ1;k2
k1 :z1 þd1
X
X
h
2cjþ1;0 þ fjþ1 ðz1 þ d1 þ Q 1 n
1k1
2k2
; z2 þ d 2 þ Q 2 n
Þ
i
3 p1jþ1;k1 p2jþ1;k2 5
ðA:1Þ
k1 :n1k1 >z1 þd1 k2 :n2k2 >z2 þd2
Appendix B. The probability mass functions of Examples 3.3 and 4.3 Tables B.1 and B.2 provide the probability mass functions for clients 1–5 in the example of Section 3.3. Table B.3 provides the probability mass functions for clients 1–5 in the example of Section 4.3.
Table B.1 The probability mass function for product 1 Customer
1 2 3 4 5
Demand 0
1
2
3
4
5
0.01 0.01 0.01 0.01 0.01
0.7 0.6 0.7 0.5 0.2
0.1 0.2 0.2 0.2 0.6
0.1 0.1 0.05 0.2 0.1
0.08 0.08 0.03 0.08 0.08
0.01 0.01 0.01 0.01 0.01
Table B.2 The probability mass function for product 2 Customer
1 2 3 4 5 Vehicle capacity = (5.5).
Demand 0
1
2
3
4
5
0.01 0.01 0.01 0.01 0.01
0.1 0.2 0.03 0.8 0.08
0.1 0.6 0.05 0.1 0.1
0.7 0.1 0.2 0.05 0.2
0.08 0.08 0.7 0.03 0.6
0.01 0.01 0.01 0.01 0.01
570
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571
Table B.3 The probability mass functions for the unified load case Demand for product 1
Demand for product 2
Customer 1
Customer 2
Customer 3
Customer 4
Customer 5
10 9 9 8 8 8 7 7 7 7 6 6 6 6 6 5 5 5 5 5 5 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
0 1 0 2 1 0 3 2 1 0 4 3 2 1 0 5 4 3 2 1 0 6 5 4 3 2 1 0 7 6 5 4 3 2 1 0 8 7 6 5 4 3 2 1 0 9 8 7 6 5 4 3 2 1 0 10 9 8 7 6 5 4 3 2 1 0
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.7 0.01 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.1 0.21 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.01 0.6 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.001 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.7 0.01 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.1 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.002 0.01 0.01 0.6 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.6 0.001 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.2 0.101 0.01 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
Vehicle capacity = 10.
References [1] M. Reimann, K. Doerner, R.F. Hartl, D-ants: Savings based ants divide and conquer the vehicle routing problem, Computers & Operations Research 31 (4) (2003) 563–591. [2] N. Christofides, The traveling salesman problem, in: N. Christofides, R. Mingozzi, P. Toth, C. Sandi (Eds.), Combinatorial Optimization, Wiley, New York, 1979. [3] G. Cornuejols, G.L. Nemhauser, Tight bounds for Christofides traveling salesman heuristic, Mathematical Programming 14 (1) (1978) 116–121.
A. Tatarakis, I. Minis / European Journal of Operational Research 197 (2009) 557–571 [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
571
M. Gendreau, G. Laporte, A. Hertz, An approximation algorithm for the traveling salesman problem with Backhauls, Operations Research 45 (4) (1997) 639–641. U. Rembold, C. Blume, R. Dillmann, Computer Integrated Manufacturing Technology and Systems, Marcel Dekker, New York, 1985. G.B. Dantzig, J.H. Ramser, The truck dispatching problem, Management Science 6 (1) (1959) 80–91. G. Clarke, J. Wright, Scheduling of vehicles from a central depot to a number of delivery points, Operations Research 12 (4) (1964) 568–581. A.A. Assad, Modeling and implementation issues, in: B.L. Golden, A.A. Assad (Eds.), Vehicle Routing: Methods and Studies, North-Holland, Amsterdam, 1988. B.L. Golden, A.A. Assad (Eds.), Vehicle Routing: Methods and Studies, Studies in Management Science and Systems, vol. 16, Elsevier Science Publications, Amsterdam, 1988, pp. 7–45. G. Laporte, I.H. Osman, Routing problems: A bibliography, Annals of Operation Research 61 (1) (1995) 227–262. P. Toth, D. Vigo, The Vehicle Routing Problem, SIAM, Philadelphia, 2002. M. Gendreau, G. Laporte, R. Seguin, Stochastic vehicle routing, European Journal of Operational Research 88 (1996) 3–12. M. Dror, G. Laporte, P. Trudeau, Vehicle routing with stochastic demands: Properties and solution frameworks, Transportation Science 23 (3) (1989) 166–176. G. Laporte, F. Louveaux, H. Mercure, Models and exact solutions for a class of stochastic location-routing problems, European Journal of Operational Research 39 (1989) 71–78. W.R. Stewart, B.L. Golden, Stochastic vehicle routing: A comprehensive approach, European Journal of Operation Research 14 (1983) 371–385. D.J. Bertsimas, A vehicle routing problem with stochastic demand, Operations Research 40 (3) (1992) 574–585. P. Trudeau, M. Dror, Stochastic inventory routing: Route design with stockouts and route failures, Transportation Science 26 (3) (1992) 171–184. M. Dror, Modeling vehicle routing with uncertain demands as a stochastic program: Properties of the corresponding solution, European Journal of Operational Research 64 (3) (1993) 432–441. M. Birattari, C. Blum, M. Sampels, O. Rossi-Doria, M. Chiarandini, L. Pauqete, T. Schiavinotto, L. Bianchi, M. Mastrolilli, T. Bousonville, H. Juillie, D. Huber, Hybrid metaheuristics for the vehicle routing problem with stochastic demands, Technical Report (2005), IDSIA/USI-SUPSI, Dalle Molle Institute for Artificial Intelligence. N. Secomandi, Comparing neuro-dynamic programming algorithms for the vehicle routing problem with stochastic demands, Computers & Operations Research 27 (2000) 1201–1225. N. Secomandi, A rollout policy for the vehicle routing problem with stochastic demands, Operations Research 49 (5) (2001) 796–802. W.-H. Yang, K. Mathur, R.H. Ballou, Stochastic vehicle routing problem with restocking, Transportation Science 34 (2000) 99–112. P. Tsirimpas, A. Tatarakis, I. Minis, E.G. Kyriakidis, Single vehicle routing with a predefined customer sequence and multiple depot returns, European Journal of Operations Research 187 (2) (2008) 483–495. A. Tatarakis, A class of single vehicle routing problems with predefined customer sequence and depot returns, Ph.D. thesis, Department of Financial and Management Engineering, University of the Aegean, Chios, Greece, 2007.