- Email: [email protected]
A tabu search heuristic for the split delivery vehicle routing problem with production and demand calendars
European Journal of Operational Research 202 (2010) 122–130
Contents lists available at ScienceDirect
European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor
Production, Manufacturing and Logistics
A tabu search heuristic for the split delivery vehicle routing problem with production and demand calendars Marie-Claude Bolduc a,b, Gilbert Laporte a,c, Jacques Renaud a,b,*, Fayez F. Boctor a,b a
Interuniversity Research Center on Enterprise Networks, Logistics and Transportation (CIRRELT), Université Laval, Québec, Canada G1K 7P4 Faculté des Sciences de l’Administration, Université Laval, Québec, Canada G1K 7P4 c Canada Research Chair in Distribution Management, HEC Montréal, 3000, chemin de la Côte-Sainte-Catherine, Montréal, Canada H3T 2A7 b
a r t i c l e
i n f o
Article history: Received 17 July 2008 Accepted 6 May 2009 Available online 18 May 2009 Keywords: Vehicle routing Distribution calendar Production calendar Inventory Internal fleet Common carrier Split delivery Tabu
a b s t r a c t The purpose of this article is to propose a tabu search heuristic for the split delivery Vehicle Routing Problem with Production and Demand Calendars (VRPPDC). This new problem consists of determining which customers will be served by a common carrier, as well as the delivery routes for those served by the private fleet, in order to minimize the overall transportation and inventory costs. We first model this problem and then propose a simple decomposition procedure that can be used to provide a starting solution. Next, we introduce a new tabu search heuristic and we describe two new neighbor reduction strategies. Finally, we present the results of our extensive computational tests. According to these tests, our reduction strategies are efficient not only at reducing computing time but also at improving the overall solution quality. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction This paper proposes a tabu search heuristic for the split delivery Vehicle Routing Problem with Production and Demand Calendars (VRPPDC), defined as follows. Let G ¼ ðV; AÞ be a graph, where V ¼ f0; . . . ; ng is the vertex set and A ¼ fði; jÞ : i; j 2 V; i – jg is the arc set. Vertex 0 is a distribution centre (DC), and the remaining vertices represent customers. The DC must deal with an a priori production calendar provided by the factory. This calendar fixes the availability of each product over a time horizon of T periods. The production calendar ðg pt Þ defines the number of units of product p ¼ f1; . . . ; Pg made available at the beginning of period t ¼ f1; . . . ; Tg; g pt is not cumulative and may vary from period to period. Finished products p can be stored at the DC, incurring a periodic unit inventory holding cost hp . Each customer i has a demand calendar dpit , which specifies the number of new units of product p that must be received in period t (i.e., by the delivery date). A private fleet of m identical vehicles is based at the DC. At each period t, each customer i may be served by the private fleet vehicles ðk ¼ 1; . . . ; mÞ, by common carriers ðk ¼ 0Þ, or by any com* Corresponding author. Address: Interuniversity Research Center on Enterprise Networks, Logistics and Transportation (CIRRELT), Université Laval, Québec, Canada G1K 7P4. E-mail addresses: [email protected] (M.-C. Bolduc), [email protected] (G. Laporte), [email protected] (J. Renaud), [email protected] (F.F. Boctor). 0377-2217/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2009.05.008
bination of the two. Each vehicle in the private fleet can perform at most one route of L time units per period, has a capacity Q k expressed in terms of volume, and each product p has a volume bp . For the private fleet, cij and t ij represent the cost and the travel time associated with arc ði; jÞ. The common carriers are assumed to have an unlimited capacity Q 0 . They manage their routes and charge for each delivery as a function of distance and quantity. In practice, the common carriers charge a fixed cost f and a variable cost which is a step function of the distance li between the depot and customer i, and of the total quantity qit transported to i in period t. The lengths of the discret, and the varization intervals for distance and quantity are l and q iable costs are cl and cq . Generally, for a given origin-destination pair, the cost function related to the quantity delivered has up to seven step intervals (see for example the Accu-Rate 2008 rating software available for download on the USF Corporation web site http://www.usfc.com). The total cost charged by the common carrier for delivering quantity qit to customer i is therefore:
e: eit ¼ f þ cl dli =le þ cq dqit =q
ð1Þ
The VRPPDC consists of determining the routes of the private fleet, as well as the customers that must be served by a common carrier, in order to minimize the overall transportation and inventory costs over a given time horizon. The resulting delivery calendar must satisfy the following constraints: (i) for each period, each private fleet vehicle performs only one route, starting and ending at the
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
DC; (ii) the capacity Q of each private fleet vehicle must be respected; (iii) the maximal length L of each route must be respected; (iv) product availability for shipment or inventory is fixed according to the production calendar; and (v) customers must receive their products no later than their contracted delivery dates without penalty for the DC. The detailed formulation is provided in the appendix. The VRPPDC is a real-life problem encountered in the snowmobile industry where all retailer (customer) orders are contracted prior to the selling season, and before the manufacturers have launched their production and stored their supplies in the DC. The contract between the manufacturer and the retailers indicate the quantity of each product by the due date, i.e., the demand schedule. The DC is thus bound by contract to ship each product by the date specified in the demand schedule. Because delivered products become property of the retailer, there is no penalty for the DC if the products are shipped before their contracted delivery dates. Since the DC usually owns a limited fleet, common carriers are often used to increase transportation capacity in peak periods. The literature on the VRPPDC is rather limited. To our knowledge, most of the articles that focus on problems involving production calendars have dealt with third-party logistics for transportation. There are, however, some exceptions. Fumero and Vercellis (1999) have developed a mathematical VRPPDC model in which only a private fleet is used. Their objective was to determine both distribution and production calendars, with the latter being adapted to feed the former. This study used ‘‘on time” deliveries only. Boudia et al. (2007) have proposed a GRASP with a reactive mechanism for solving the single product problem. Inventory management is closely related to the classical, mono-product, Inventory Routing Problem – IRP (see Campbell et al., 2002; Moin and Salhi, 2007; Zhao et al., 2008). However, in the VRPPDC, quantities requested by customers are fixed over time by the demand calendars and not all products are available at the depot at the beginning of the planning horizon. This paper makes the following contributions. We first introduce a new complex VRP variant, which, to our knowledge, has never been studied before. We then model this problem as a mixed integer program, and we solve the model using a tabu search heuristic with split deliveries. Last, we present two new neighbor reduction strategies which, when tested, proved efficient not only at reducing computing time but also at improving overall solution quality. The remainder of this article is organized as follows. In Section 2, we present our VRPPDC local search heuristic. Section 3 describes our computational experiments and our results, and Section 4 offers our conclusions.
2. VRPPDC local search heuristic The VRPPDC can be formulated as a mixed integer model in which the objective function is non-linear because of common carrier cost function eit . The model, presented in the appendix, cannot be solved optimally for any realistic size instance. Since tabu search is a very powerful tool for solving the classical VRP (Tarantilis and Kiranoudis, 2002) and its variations (Osman and Salhi, 1996), we have developed a composite local search heuristic made up of three components: (i) construction of an initial feasible solution, (ii) a tabu search phase which allows split deliveries and the use of common carriers and (iii) a final improvement phase. The tabu search considers intermediate infeasible solutions by using penalty functions like those developed by Gendreau et al. (1994) and Cordeau et al. (1997). More specifically, vehicle capacity and route length may be violated during the search (soft constraints) but they must be respected in the final solution. However, the production calendar (i.e., inventory availability) and the demand calendar (i.e., inventory
123
quantities and the subsequent delivery dates) are always enforced during the search (hard constraints). These hard and soft constraints lead to an important difference with respect to the other tabu searches proposed in the literature: the double exploration of each neighbor, respectively called ‘‘complete switch” and ‘‘partial switch”. The complete switch evaluates the possibility of transferring all the product quantities for one customer from the original route (corresponding to a vehicle-period combination) to another one, while still respecting the hard constraints. Such a move may or may not generate feasible solutions if some soft constraints are not respected. The partial switch is a split delivery strategy, consisting of moving only the quantities that respect the target vehicle’s availability in terms of space and time. Such a partial transfer will always lead to a feasible route as it respects both the soft and hard constraints. The following subsections present the initial feasible solution, the tabu search, and the improvement phases. 2.1. Initial feasible solution If the quantities to be delivered are determined for each product and each period, the routing problem for each period corresponds to a multi-product version of the Vehicle Routing Problem with a Private fleet and a Common carrier (VRPPC) as defined by Bolduc et al. (2008). The VRPPC consists of minimizing the total cost of serving once a set of customers using either the limited private fleet vehicles or the common carrier and, for the private fleet, to perform the routes. Two heuristics have recently been developed for the VRPPC: the Selection, Routing and Improvement heuristic (SRI – Bolduc et al., 2007) and the Randomized construction, Improvement, Perturbation metaheuristic (RIP – Bolduc et al., 2008). SRI is a three phase heuristic: (1) selection of customers to be served by the common carrier, (2) construction of an initial solution, and (3) improvement procedure. This heuristic is quite fast and generates better results than a previous algorithm. RIP is a five step metaheuristic: (1) randomized savings construction phase, (2) a 4-opt* route improvement procedure, (3) 2*-interchange inter-route improvement procedure, (4) 2-add-drop improvement procedure, and (5) switch procedure. The 4-opt* heuristic is a route improvement procedure which may exchange up to four vertices (Lin, 1965; Renaud et al., 1996) and the 2*-interchange is a restricted version of the k-interchange procedure (Osman, 1993; Bolduc et al., 2008). These two procedures use sophisticated decision rules to reduce the neighborhood size and accelerate the search without decreasing solution quality. In the VRPPC computational results, RIP generated better results than SRI, but required much longer computational processing times. The heuristics developed for the VRPPC can be adapted to generate solutions to the single-period subproblems. Indeed, adaptations are necessary since these VRPPC heuristics were designed to handle only one product. The rest of this subsection presents a two-phase method designed to allow SRI and RIP to be used for the solution of singleperiod multi-product subproblems. The two-phases solve the problem for each period independently, and the best result of each phase is retained for each period. In each phase, the common carrier cost for each customer is calculated using (1). 2.1.1. Two-phase method using SRI and RIP Phase 1: Based on the demand calendars, generate a list of customers with a demand due in the current period. For each customer, generate a new temporary customer for each product p ordered (each customer may be split into, at most, P customers) and evaluate the common carrier cost
124
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
Table 1 Notation used in the tabu search algorithm. BðsÞ NðsÞ cðsÞ hðsÞ qðsÞ dðsÞ f ðsÞ gðsÞ s; s s0 s
a b
qi rikt g h d k
mp sikt v l j
attribute set for solutions ðBðsÞ ¼ fði; k; tÞg; i ¼ f1; . . . ; ng; k ¼ f1; . . . ; mg; t ¼ f1; . . . ; TgÞ neighborhood for solution s transportation cost for solution s holding inventory costs for solution s total excess quantity for solution s total excess duration for solution s total cost for solution s; f ðsÞ ¼ cðsÞ þ hðsÞ þ aqðsÞ þ bdðsÞ total penalized cost for solution s (i.e., f ðsÞ, to which is added a penalty term for the most frequently visited customers in order to diversify the search) solutions best solution identified in the current neighborhood NðsÞ best feasible solution identified penalty factor for capacity violation penalty factor for tour duration violation number of times a customer i has been added to the solution aspiration criterion of attribute ði; k; tÞ total number of iterations to be performed tabu duration parameter used to update a and b iteration counter movable quantities of product p last iteration for which attribute ði; k; tÞ is declared tabu number of iterations without improving s average number of iterations needed to pass from one feasible solution s0 to the next number of iterations a solution s0 remains infeasible
for this scenario. This cost is needed to solve the VRPPC. Solve this subproblem for each period, using either the SRI or the RIP heuristic. This strategy may produce large subproblems since each customer may be split a number of times. Phase 2: Generate a list of customers with a demand due in the current period. For each customer, create a single unified demand, corresponding to the sum of all product demands, and evaluate the common carrier cost for this scenario. Solve this subproblem for each period, using either the SRI or the RIP heuristic.
correspond to the objective function (2) in the appendix. Violations of the vehicle capacity or route duration constraints are penalized according to the following functions:
qðsÞ ¼
k¼1 t¼1
dðsÞ ¼
2
p¼1 i¼1
T 6X n X
6 4
#þ bp qpikt Q k 3þ
n X
i¼0
;
j¼0 j–i
7 t ij xijkt L7 5 ;
where ½xþ ¼ maxf0; xg. When s is feasible, qðsÞ and dðsÞ both equal to zero. 2.2.1.2. Attribute set management and movable quantities. Each solution s has an associated attribute set BðsÞ ¼ fði; k; tÞg, in which customer i is served by vehicle k (0: common carriers; 1 to m: private fleet) at period t. Moving from a solution s to a neighbor solution s0 2 NðsÞ modifies BðsÞ. The neighborhood NðsÞ of a solution s is composed of all solutions obtained by moving a quantity v piktk0 t0 0 associated to customer i from route ðk; tÞ to ðk ; t0 Þ, where 0 0 0 0 k; k ¼ f0; . . . ; mg; t; t ¼ f1; . . . ; Tg and k – k or t – t. For a customer i and the initial route ðk; tÞ under consideration, the movable quantities v piktk0 t0 depend on the type of switch. In the complete switch, when the target period t 0 is earlier than the original period t; v piktk0 t0 is the minimum movable quantity possible between the already planned delivered quantity ðqpikt Þ and the available inventory at the DC between periods t0 and tðspt Þ. However, when the period t0 is later than the original period t; v piktk0 t0 is simply set equal to already planned delivered quantities. In the partial switch, v piktk0 t0 also takes the available vehicle capacity into account. Note 0 that quantities may be split between vehicles k and k in the same period or over different periods t and t0 . To summarize, v piktk0 t0 may be formulated as follows: complete switch:
v piktk t
¼
8 > < min q
; min fs g ; if t0 < t ph pikt 0
> :q
h¼t ;...;t
otherwise;
pikt ;
partial switch:
v piktk t
0 0
2.2. Tabu search phase
m X
k¼1 t¼1
0 0
Repeating these strategies for each period leads to one single feasible initial solution, composed of all the routes obtained for each period. This solution can be enhanced through a tabu exploration phase and an improvement phase. Table 1 provides the notation used in the tabu search algorithm.
" m X T P X n X X
8 n P > > 0 ; if t 0 < t ; min fs g; Q q min q 0 > ph k pikt phk t < h¼t 0 ;...;t h¼1 ¼ n > P > > min qpikt ; Q k qphk0 t0 ; otherwise: : h¼1
The following subsections present the tabu search components, the complete VRPPDC tabu search algorithm and the neighbor reduction strategies. 2.2.1. Tabu search components The tabu search phase is the core of the VRPPDC heuristic (see Section 2.2.2). It is based on the following components which are described below: penalized cost function, attribute set management and movable quantities and finally, the exploration phase itself. 2.2.1.1. Penalized cost function. During the exploration phase a penalized cost function f ðsÞ ¼ cðsÞ þ hðsÞ þ aqðsÞ þ bdðsÞ is used, where
cðsÞ ¼
n X n X m X T X i¼0
hðsÞ ¼
j¼0 k¼1 t¼1 j–i
P X T X p¼1 t¼1
hp spt ;
cij xijkt þ
n X n X T X i¼0
j¼0 t¼1 j–i
eit xij0t ;
An aspiration criterion rikt is associated with each attribute ði; k; tÞ. It is the initial cost solution if ði; k; tÞ belongs to the attributes included in the initial solution and it is set equal to 1 for all the other attributes. The aspiration criteria are updated during the search. 0
2.2.1.3. Exploration phase. For given attributes ði; k; tÞ and ði; k ; t0 Þ and for a quantity v piktk0 t0 , the exploration phase performs the following steps: 1. Set the current solution s equal to the initial solution in the 0 neighborhood ðs ¼ sÞ. (In the following steps, ðk ; t0 Þ 2 s and ðk; tÞ 2 sÞ. 0 2. If customer i is not already in route ðk ; t0 Þ, insert this customer into the route and apply the 3-opt procedure (Lin, 1965). 0 3. Move v piktk0 t0 for all p from route ðk; tÞ to ðk ; t0 Þ. If v piktk0 t0 < qpikt (i.e., only some of the original delivery quantity is moved), prioritize the product that has the highest inventory cost. 4. If v piktk0 t0 ¼ qpikt , remove the customer from this route and reoptimize it using 3-opt.
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
5. For periods t to t0 , update the inventories spt and reevaluate the global inventory cost. 6. Reevaluate the routing cost for both routes. 7. If the solution s is non-tabu ðsik0 t0 < kÞ or if a best feasible solu0 tion for attribute ði; k ; t0 Þ is obtained (f ðsÞ < rik0 t0 and s is feasible), then apply the following operations: (a) If the current solution s is worse than the initial solution s in the neighborhood ðf ðsÞ > f ðsÞÞ; f ðsÞ is penalized pffiffiffiffiffiffiffiffiffiffi for diversification gðsÞ ¼ f ðsÞ þ 0:02 mnT ðcðsÞ þ hðsÞÞqi =k; otherwise, gðsÞ ¼ f ðsÞ. (b) If the current solution s is the new best neighbor ðgðsÞ < f ðs0 ÞÞ, set s0 ¼ s. When a solution s0 is retained as the best solution found in NðsÞ, customer i cannot be reinserted into route ðk; tÞ for the next h iterations. During these h iterations, a tabu move involving ði; k; tÞ is allowed only if its global cost is lower than the aspiration criterion rikt . This criterion is updated whenever a new best feasible solution is obtained. To diversify the search, the attributes repeatedly inserted into the solution are penalized by a term proportional to qi (see Step 7(a)). Because of the split quantities, penalizing the attribute ði; k; tÞ is not appropriate since using a penalty qikt may lead to moving the same customer to several different routes. To diversify the solution, different customers are instead moved at each iteration and their moves to other routes are penalized using qi . The other variables used in the determination of gðsÞ are the same as those mentioned by Cordeau et al. (1997), except for the constant factor, which is set to 0.02 instead of 0.015. After testing many values, we found that 0.02 generated the best results. 2.2.2. Complete VRPPDC Tabu search algorithm A step-by-step description of the VRPPDC heuristic is given below. 1. Obtain the initial feasible solution s using the two-phase method using SRI and RIP outlined in Section 2.1. 2. Set qi ¼ 0, for all i. 3. For each attribute ði; k; tÞ included in the initial best feasible solution s , set rikt ¼ f ðs Þ; otherwise, set rikt ¼ 1. 4. For each tabu search iteration ðk ¼ 1; . . . ; gÞ, perform the following operations: (a) Set the best neighborhood solution value f ðs0 Þ ¼ 1. (b) Set the initial solution in the neighborhood s ¼ s . 0 0 (c) For each attributes ði; k; tÞ and ði; k ; t 0 Þ, where k – k or 0 0 t – t , do the following. If t > t and the current solution s do respect the delivery dates, go to Step 4(c1); otherwise, go to Step 4(c2). (c1) COMPLETE SWITCH. 1. For each product p, calculate the movable quantity P v piktk0 t0 . If Pp¼1 v piktk0 t0 ¼ 0 (nothing can be moved), restart Step 4(c) using the next neighbor. 2. Apply the Exploration phase described in Section 2.2.1. (c2) PARTIAL SWITCH. 1. For each product p, calculate the movable quantity P v piktk0 t0 . If Pp¼1 v piktk0 t0 ¼ 0, restart Step 4(c) using the next neighbor. 2. Apply the Exploration phase described in Section 2.2.1 for all p, where v piktk0 t0 – 0. (d) Set the attribute ði; k; tÞ of the best neighbor s0 tabu ðsikt ¼ k þ hÞ and adjust the diversification criterion ðqi ¼ qi þ 1Þ. (e) If the best neighbor s0 is a new best feasible solution (i.e., f ðs0 Þ < f ðs Þ and s0 is feasible), perform the following operations:
125
(e1) Set s ¼ s0 . (e2) For each attribute ði; k; tÞ included in the solution s0 , update the aspiration criteria ðrikt ¼ minfrikt ; f ðs0 ÞgÞ. (f) If qðs0 Þ ¼ 0 (no quantity violation), then a ¼ a=ð1 þ dÞ. Otherwise, a ¼ ð1 þ dÞa. (g) If dðs0 Þ ¼ 0 (no duration violation), then b ¼ b=ð1 þ dÞ. Otherwise, b ¼ ð1 þ dÞb.
2.2.3. Neighbor reduction strategies The tabu search can be very time consuming due to the large size to the neighborhood NðsÞ and also to the non-linear cost functions that must be constantly reevaluated. Thus, we propose two neighbor reduction strategies designed to reduce the computing time. The first strategy, called Random, involves randomly selecting some neighbors from NðsÞ for evaluation, thus allowing only a predefined proportion of the total neighborhood to be considered. This strategy is quite simple and may be used on any kind of problem. By forcing a decrease in the number of neighbors, the strategy is able to decrease the time needed to find a solution. However, because the strategy does not necessarily select the best neighbor, it could lead to a different, and maybe worse, solution. The second strategy, called Distance, takes the problem configuration into account. The idea is simple: a given neighbor is evaluated only if the customer under consideration is within a fixed distance from the nearest customer on the route into which it will be inserted. The fixed distance is predefined as being equal to a fraction of the longest distance on the geographical map. This strategy makes it possible not to evaluate the insertion of a customer into a distant route. 2.3. Improvement phase When the tabu iterations are completed, the following improvement heuristic is applied using a first improvement strategy and until there is no improvement: 1. advance the delivery of complete blocks of product demand; 2. advance the delivery of partial blocks of product demand; 3. switch an internal customer (i.e., one served by a private vehicle) with an external customer (i.e., one served by a common carrier). 0
For each route ðk ; t 0 Þ, the first two steps verify whether it is possible to deliver, at t0 , a delivery planned for a later time t. Step 1 evaluates whether, in terms of inventory availability and the remaining vehicle capacity, a given customer’s deliveries, planned for a later period (i.e., in period t, where t > t0 Þ can be completely 0 integrated into route ðk ; t 0 Þ at an earlier period. After all current customers have been evaluated, this procedure is repeated for the partial blocks (Step 2). For a product that has been planned for delivery in quantity qpikt , the delivery of any quantity between 1 and qpikt may be moved to an earlier period. These two steps allow the global inventory cost to be reduced by filling the vehicles scheduled as much as possible. Step 3 is similar to the Osman 1–1 exchange (1993). The difference between Step 3 and the original method is that the switches between customers served by a private fleet vehicle and those served by common carriers are the only ones evaluated. When a customer switch is performed, the entire product quantities for these two customers are moved. In other words, there are no partial transfers. A switch can be conducted for any vehicle-periodcustomer combination, as long as the problem constraints are respected and the global cost (i.e., transportation and inventory costs) is reduced. Then, the 3-opt heuristic is applied to each modified route.
126
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
ing the demand for five periods of any product in any one period, the production of each period was assigned entirely to the most urgent product remaining to be delivered, chosen with respect to the minimum ratio between the current stock level and the demand to be filled. The parameters of the common carrier cost function were set ¼ 0:25Q . as: f ¼ 60; cl ¼ 50; l ¼ 25; cq ¼ 15; q When determining the number and capacity of the private fleet vehicles, we have assumed that, on average, about 75% of the demand would be filled by this fleet, while the other 25% would be filled by a common carrier. The number of vehicles was initially set equal to 2, and this number was increased until the demand could be covered by the number of vehicles selected, each with a capacity between 100 and 150 units. The exact vehicle capacity was determined with respect to the demands. For example, a demand of 500 units per period would require a total fleet capacity of at least 375 (75% of 500). Since each vehicle is assumed to have a capacity between 100 and 150 units, two vehicles would have a total capacity within the interval [200, 300], while three vehicles would have a capacity of [300, 450]. Thus, for a demand of 375, three vehicles would be necessary, each with a capacity of 125 (i.e., 375 divided by 3). The volume of each product was set to 1, and the vehicle speed was set to 80 km/h, the maximum route length to 13 time units (hours), and the variable distance cost to $2/km. Table 2 presents the number of vehicles, their capacity and the total demand for each instance.
3. Computational results We first describe our test instances, then present the tabu search parameters, and finally report our results. 3.1. Test instances For each of the 100 instances, 50 customer coordinates were randomly generated on a plane measuring 240 km by 200 km. The number of products was set equal to 3; the inventory costs per product unit were $0.15, $0.25 and $0.35, respectively; and the planning horizon was equal to 10 periods. To determine the demand calendar for each customer, the following rules were applied. First, each customer demand level was classified accordingly to the following probabilities: 25% had a low demand per period, 50% had an average demand per period (i.e., twice the low demand), and 25% had a high demand per period (i.e., four times the low demand). The upper bounds of low demand products were set equal to 2, 4 and 3 units for product 1, 2 and 3, respectively. For each period, demands were randomly generated according to a continuous uniform distribution between 0 and an upper bound based on the customer category. For each customer, the ordering interval was randomly selected between one and four periods; the day of the first order was then selected. The quantity of each order was the sum of the customer’s daily demand over the ordering interval. The production calendar was designed according to the customer demand calendars. For each product, the initial inventory at the DC was equal to the average demand over five periods, and the production lot size was the average demand divided by the production frequency, which was set randomly between one and five periods. Since the factory was assumed capable of produc-
3.2. Tabu search parameters The number of iterations of the tabu search algorithm in the VRPPDC heuristic was set g ¼ 100; 000 iterations. Two parameter sets (A and B) were also tested.
Table 2 Some instance characteristics. Instance
m
Q
Total demand
Instance
m
Q
Total demand
Instance
m
Q
Total demand
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
3 3 2 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
110 130 140 120 100 110 120 110 110 110 110 110 120 150 120 120 110 120 100 100 100 130 120 120 120 100 110 110 120 110 100 100 110 130
4230 5089 3692 4550 4015 4242 4789 4146 4090 4188 4292 4373 4601 3774 4601 4683 4404 4459 3989 3977 4063 5121 4606 4656 4704 3924 4439 4161 4457 4391 3974 3888 4253 4982
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
3 3 3 3 2 2 3 3 3 3 2 2 3 3 3 3 3 3 3 3 3 2 3 3 2 3 3 3 3 3 3 3 3
110 110 100 100 140 150 110 140 110 120 130 140 100 140 130 120 120 110 110 110 120 140 110 120 150 110 110 120 110 120 120 100 110
4216 4190 3998 4064 3643 3815 4226 5291 4191 4850 3324 3570 3887 5483 4896 4502 4452 4211 4249 4254 4515 3716 4290 4590 3867 4103 4153 4575 4212 4581 4483 3900 4392
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3
110 120 110 120 110 110 110 120 120 100 110 110 110 120 130 110 130 110 110 110 120 110 120 100 120 120 100 120 110 150 100 100 110
4422 4475 4287 4528 4300 4146 4343 4690 4787 3991 4178 4266 4220 4541 5015 4420 3283 4099 4257 4251 4490 4284 4534 4004 4517 4615 3922 4608 4253 3827 4027 4070 4214
127
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
3.2.1. Parameter set A Parameter set A used the VRPPDC heuristic as described in Section 2.2.2, with the following parameter values: a ¼ 1; b ¼ 1; h ¼ ½7:5log10 n and d ¼ 1:5. 3.2.2. Parameter Set B This second set was designed to calibrate the search parameter values. First, the penalty factors were set at a ¼ 0:5 and b ¼ 0:5. Because the VRPPDC heuristic starts from an initial feasible solution, lower penalty factors can be used at first. In addition, the problem considered in this paper can lead to many tabu iterations with infeasible solutions. Since many infeasible solutions may be generated in a neighborhood before a new feasible solution is obtained, we have added a step in our algorithm before Step 4(f) to update the penalty factors a and b: ðf 0 Þ If s0 is feasible, or if j 6 l, then d ¼ 1:5; otherwise, d ¼ 1 until s0 becomes feasible. To prevent the tabu search algorithm from being stuck in a local optimum with only a few neighbors, we also adjusted the criterion in Step 4(d) (Section 2.2.2) to read: h ¼ ½7:5log10 n þ dn=100e. This variable tabu duration allows the search to visit other neighbors. 3.3. Results The heuristic was coded in Microsoft VB.Net 2005 and run under Windows XP on a personal computer with a Xeon 3.6 GHz processor and 1.00 Gb of RAM. All reported times are expressed in seconds, and all statistics represent the average results over a set of 100 test instances. The following results are compared with the best-known solutions obtained for all instances (see Table 3). For each instance, we have recorded the best solution found during all our tests which include the computational experiments used in the paper but also some other tests conducted with parameters that have been tried during the development of the algorithms but not presented in the paper. This section presents our computational experiments with our tabu search heuristic. We first present the results of the initial feasible solution and improvement phase.
The neighbor reduction strategies are then analyzed. Finally, results for the complete VRPPDC heuristic are presented. 3.3.1. Evaluation of the initial feasible solution and improvement phase Our two-phase method (Section 2.1.1) used to obtain the initial feasible solution is clearly inefficient, since, as shown in Table 4, the average deviation above the best-known solution is 57.23% for the RIP-based heuristic and 64.01% for SRI-based heuristic. The improvement phase (Section 2.3) yields slightly better results: 49.15% for RIP and 51.28% for SRI. However, RIP solution times are much longer than those of SRI, with an average solution time of 301 s compared to only two seconds for SRI. The results reported in this section for the VRPPDC heuristic, called SRITabu, were initiated from the SRI-based solution which was much faster. 3.3.2. Neighbor reduction strategies To evaluate the proposed neighbor reduction strategies, all 100 instances were solved with SRITabu and 500 iterations, using parameter set A and the two neighbor reduction strategies. The goal of the evaluation was to demonstrate that, using our neighbor reduction strategies, solutions of the same quality could be obtained in less time. Table 5 present the results of our tests. The first column in the table presents the average deviation for the original SRITabu after 500 iterations. Columns 2, 3 and 4 contain the results obtained with the Random reduction strategy that eliminates 25, 50 and 75% of the neighbors. The best quality result for this strategy was obtained with a random 50% reduction, producing an average deviation of 11.57%. This result is close to that of the original SRITabu, which produced an average deviation of 10.88%, but it is much faster, with an average of 16 s compared to 29. Columns 5, 6 and 7 contain the results using the Distance reduction strategy. The first test discarded all neighbors that were not closer than 15% of the longest distance on the map (column 5), the second, those that were not closer than 25% of the longest distance (column 6), and the last, those that were not closer than 35% (column 7). Overly restricting the distance produced a significant deterioration in solution quality: the deviation obtained with
Table 3 Best-known solutions. Instance
Best
Instance
Best
Instance
Best
Instance
Best
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
23,113.57 25,391.41 18,268.49 23,585.96 25,173.42 27,788.49 27,581.27 24,156.73 21,783.57 25,171.74 28,200.50 28,937.65 26,087.67 18,151.88 25,250.77 27,917.18 27,729.59 22,731.64 22,629.94 25,401.72 26,403.15 26,359.05 27,477.55 24,901.55 26,618.87
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
24,745.46 27,181.70 25,466.19 23,228.56 29,169.52 24,635.79 25,759.79 22,809.43 23,604.40 26,281.08 24,837.63 26,247.15 25,163.93 19,983.13 20,232.92 24,333.88 24,578.88 24,127.50 24,602.39 20,382.55 20,443.80 22,683.97 26,214.30 24,993.88 24,243.11
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
27,856.98 24,289.10 20,583.28 24,700.93 27,918.29 19,058.70 25,576.13 27,075.20 20,262.76 23,327.91 24,512.32 24,346.14 25,515.31 27,200.36 25,266.38 26,227.86 25,801.10 28,676.45 25,742.31 24,900.09 25,947.90 25,470.19 23,799.70 23,840.44 26,594.19
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
26,828.62 27,672.02 22,771.41 26,817.43 24,737.76 26,776.09 23,887.33 28,674.86 19,024.32 24,307.21 23,319.70 25,379.44 25,831.23 23,950.72 25,659.14 25,024.39 22,336.13 25,274.56 25,586.53 25,879.00 23,210.87 19,520.18 22,771.17 26,181.31 24,832.96
128
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
Table 4 Results for the initial feasible solution and the improvement phase heuristics. Initial feasible solution
RIP
SRI
Improvement phase
Without
With
Without
With
Average deviation from the best-known solution (%) Average time (s) Number of best solutions
57.23
49.15
64.01
51.28
301 0
301 0
2 0
2 0
Table 5 Comparison of the neighbor reduction strategies for SRITabu with 500 iterations. Original SRITabu
Neighbor reduction strategy Random
% of reduction Average deviation from the best-known solution (%) Average time (s) Number of best solutions
Distance
0 10.88
25 11.68
50 11.57
75 12.16
15 18.48
25 10.36
35 10.57
29 0
23 0
16 0
9 0
18 0
21 0
23 0
distances closer than 15% is 18.48%, while 10.36% was obtained with distances closer than 25%. Both the 25% and 35% distance reduction strategies produced better results than the SRITabu by itself and were also faster, although the distance reduction strategy is, on average, slower than the random strategy. The next subsection presents the results of our tests using the best parameter combination for each neighbor reduction strategy: ‘‘closer than 25%” for the distance strategy and ‘‘50%” for the random reduction strategy. 3.3.3. VRPPDC heuristic: SRITabu results Tables 6 and 7 show the average results obtained by the SRITabu heuristic over the 100 instances, using 500 to 100,000 iterations. The tables are separated horizontally into three main sections: the results obtained with the SRITabu heuristic only, the results obtained with the ‘‘closer than 25%” distance reduction strategy, and then those obtained with the ‘‘50%” random reduction strategy. Each horizontal section shows the worst, best and average deviation from the best-known solution, the average time in
seconds, the number of best solutions found, and the number of solutions found that deviated less than 0.5% from the best solutions. Table 6 uses parameter set A, while Table 7 employs the parameter set B. The best results obtained by SRITabu heuristic with parameter set A are shown in the ‘‘100,000 iterations” column of Table 6. For the SRITabu only (i.e., no neighbor reduction strategy), the average deviation is 2.50%, and the solutions were obtained in an average of 5725 s. This ‘‘no strategy” method produced six best solutions and 24 results with less than 0.5% of deviation. The distance strategy, on the other hand, generated an average deviation of 2.62% and used half the computing time (i.e., an average of 2719 s). This strategy produced four best solutions and 11 results with less than 0.5% of deviation. However, the random strategy generated the best overall results. With only 10,000 iterations, this strategy yielded a 2.51% deviation, a result similar to those obtained by the two others methods with 100,000 iterations. Moreover, this strategy produced these results much faster, within only 300 s on average. With 100,000 iterations, the random strategy obtained the lowest average deviation (1.12%) almost as fast as the distance strategy obtained a much higher average deviation (i.e., an average of 2879 s) and obtained four best solutions and 29 solutions with less than 0.5% of deviation. As these results show, the random reduction strategy halved both the average deviation and the computing time with respect to the original SRITabu application, using parameter set A. From an empirical point of view, these results are significant. Table 7 presents the results for parameter set B. Using this strategy, almost no improvement was obtained over 50,000 iterations. For parameter set B, the best quality results were obtained by the original SRITabu with an average deviation of 1.77%. This strategy generated six best solutions and 14 results with less than 0.5% of deviation. However, this strategy was also the most time consuming of the three, taking 2702 s on average. In fact, both the random and distance reduction strategies obtained results faster, with an average of 1433 and 1420 s, respectively. Unfortunately, these solutions have larger deviation of 2.42% and 2.59%, respectively. With parameter set B, only the original SRITabu performed better than the various scenarios using parameter set A. Examining both tables, strategy by strategy allows parameter set A to be compared to B. In terms of quality, the ‘‘no strategy” results obtained with A at 100,000 iterations were obtained with B at
Table 6 Results of the SRITabu heuristic with parameter set A.
No strategy Deviation
Time Number of solutions Distance 25% Deviation
Time Number of solutions Random 50% Deviation
Time Number of solutions
Iterations
500
1000
2000
5000
10,000
15,000
20,000
50,000
100,000
Worst Best Average Average (s) Best Under 0.5%
18.11% 4.65% 10.88% 29 0 0
17. 03% 2.13% 7.36% 57 0 0
17.03% 0.33% 5.24% 115 0 1
17.03% 0.14% 3.73% 290 0 2
17.03% 0.00% 3.18% 579 1 5
17.03% 0.00% 2.96% 871 1 7
17.03% 0.00% 2.83% 1163 1 11
17.03% 0.00% 2.56% 2892 5 21
17.03% 0.00% 2.50% 5725 6 24
Worst Best Average Average (s) Best Under 0.5%
18.69% 4.36% 10.36% 21 0 0
15.19% 3.53% 7.29% 39 0 0
14.10% 1.86% 5.28% 72 0 0
14.10% 0.80% 3.95% 162 0 0
14.10% 0.05% 3.42% 304 0 1
14.10% 0.00% 3.22% 444 1 2
14.10% 0.00% 3.08% 582 1 3
14.10% 0.00% 2.75% 1.385 4 10
14.10% 0.00% 2.62% 2.719 4 11
Worst Best Average Average (s) Best Under 0.5%
20.34% 5.94% 11.57% 16 0 0
16.09.% 2.87% 7.60% 32 0 0
10.37% 1.29% 5.24% 65 0 0
8.45% 0.13% 3.25% 153 0 2
8.11% 0.13% 2.51% 300 0 3
7.33% 0.13% 2.14% 443 0 4
7.33% 0.13% 1.88% 589 0 4
7.12% 0.00% 1.40% 1.475 2 11
7.12 0.00% 1.12% 2.879 4 29
129
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130 Table 7 Results of the SRITabu heuristic with parameter set B.
No strategy Deviation
Time Number of solutions Distance 25% Deviation
Time Number of solutions Random 50% Deviation
Time Number of solutions
Iterations
500
1000
2000
5000
10,000
15,000
20,000
50,000
100,000
Worst Best Average Average (s) Best Under 0.5%
36.81% 4.15% 11.90% 30 0 0
22.80% 2.92% 8.09% 59 0 0
22.80% 0.97% 5.50% 117 0 0
14.39% 0.00% 3.15% 282 1 3
9.15% 0.00% 2.33% 550 2 6
8.32% 0.00% 2.04% 817 3 9
8.32% 0.00% 1.88% 1095 4 12
8.32% 0.00% 1.77% 2702 6 14
8.32% 0.00% 1.77% 5165 6 14
Worst Best Average Average (s) Best Under 0.5%
30.13% 4.84% 11.16% 18 0 0
16.80% 2.96% 7.82% 34 0 0
10.86% 1.89% 5.32% 64 0 0
7.70% 0.39% 3.49% 145 0 1
7.70% 0.53% 2.84% 274 0 0
7.70% 0.53% 2.70% 409 0 0
7.70% 0.53% 2.63% 550 0 0
7.70% 0.28% 2.59% 1.420 0 1
7.70% 0.28% 2.59% 2.874 0 1
Worst Best Average Average (s) Best Under 0.5%
21.93% 5.42% 11.92% 16 0 0
16.52% 2.87% 8.03% 31 0 0
15.17% 1.37% 5.69% 60 0 0
10.07% 0.48% 3.58% 144 0 1
9.57% 0.48% 2.70% 285 0 1
9.57% 0.48% 2.55% 427 0 1
9.57% 0.48% 2.53% 572 0 1
9.57% 0.48% 2.42% 1.433 0 1
9.57% 0.48% 2.42% 2.800 0 1
10,000 iterations and with only one tenth of the computing time. However, the same improvement is not true for the random and the distance strategies. Overall, the best combination appears to be the random neighbor reduction strategy applied with parameter set A. This is the only combination that produces average results under 2% within less than 600 s, and it leads to the best overall performance, with an average deviation of 1.12%. 4. Conclusion
uikt : upper bound for the load of vehicle k upon leaving customer i in period t; spt : inventory of product p at the DC at the end of period t; qpikt : quantity of product p leaving the DC to be delivered to cusP tomer i with vehicle k in period t. Note that qit ¼ Pp¼1 qpi0t ; t 2 f1; . . . ; Tg. The model is expressed as follows:
Minimize
n X n X m X T X i¼0
We have introduced the split delivery Vehicle Routing Problem with Production and Demand Calendars, a rich variant of the VRP which, to our knowledge, has never previously been studied. The problem was modeled as a mixed integer program and solved by means of a tabu search heuristic with split deliveries. Our results demonstrate that it is possible to obtain an excellent solution within a reasonable amount of time. We have also introduced several neighbor reduction strategies which were proved efficient not only at reducing computing time but also at improving the overall solution quality.
þ
j¼0 k¼1 t¼1 j–i
P X T X
This work was partially supported by the Canadian Natural Sciences and Engineering Research Council under Grants OPG0036509, 39682-05 and OPG0172633. This support is gratefully acknowledged. Thanks are due to three referees for their valuable comments.
n X n X T X i¼0
eit xij0t
j¼0 t¼1 j–i
hp spt
ð2Þ
p¼1 t¼1
subject to spt ¼ sp;t1 þ g pt
n X m X
qpikt
i¼1 k¼0
ðp 2 f1; . . . ; Pg; t 2 f1; . . . ; TgÞ m X t X
qpikh P
k¼0 h¼1
Acknowledgements
cij xijkt þ
t X
dpih
h¼1
ðp 2 f1; . . . ; Pg; i 2 f1; . . . ; ng; t 2 f1; . . . ; TgÞ n X
m X
x0jkt 6 m ðt 2 f1; . . . ; TgÞ
ð4Þ ð5Þ
j¼1 k¼1 n X
xhjkt ¼
j¼0 j–h
n X
xihkt
i¼0 i–h
Appendix A. Problem formulation
ðh 2 f0; . . . ; ng; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
In addition to the notation in introduction, define the following variables:
n X
n X
i¼0
j¼0 j–i
8 > < 1 if vehicle k visits a vertex j immediately xijkt ¼ after vertex i at period t > : 0 otherwise; 1 if vehicle k is used in period t ykt ¼ 0 otherwise;
ð3Þ
P X p¼1
t ij xijkt 6 L ðk 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
qpikt 6 Q k
n X
ð6Þ ð7Þ
xjikt
j¼0 j–i
ði 2 f1; . . . ; ng; k 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ
ð8Þ
130
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130 n X n X i¼0
(10) ensure that the vehicle capacity is never exceeded. Constraints (11) were introduced by Miller et al. (1960) to eliminate subtours. Kara et al. (2004) describe a lifting procedure for (11).
xijkt 6 ðn þ 1Þykt
j¼0 j–i
ðk 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ P X n X
ð9Þ
bp qpikt 6 Q k ykt
p¼1 i¼1
ðk 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
ð10Þ
uikt ujkt þ ðn 1Þxijkt 6 n 2 ði; j 2 f1; . . . ; ng; i – j; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ ð11Þ xijkt 2 f0; 1g ði; j 2 f0; . . . ; ng; i – j; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
ð12Þ
ykt 2 f0; 1g ðk 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ
ð13Þ
uikt P 0 ði 2 f1; . . . ; ng; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
ð14Þ
qpikt P 0 ðp 2 f1; . . . ; Pg; i 2 f1; . . . ; ng; k 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ
ð15Þ
spt P 0 ðp 2 f1; . . . ; Pg; t 2 f0; . . . ; TgÞ:
ð16Þ
The objective function minimizes the routing costs of the private fleet, the non-linear costs of the common carrier (see (1)), and the DC’s inventory costs. Constraints (3) maintain the product balance by adjusting the DC’s inventory level according to incoming production and outgoing shipments. Constraints (4) ensure that the contracted customer demand calendars are satisfied. Constraints (5) specify that at most m private fleet vehicles can be used each day, while constraints (6) indicate that the same vehicle k must arrive at and leave from customer h. Constraints (7) require that each private fleet route lasts at most L time units. Constraints (8) and (9) respectively specify that deliveries can only be made if a tour exists already and that a tour exists only if a vehicle is used. Constraints
References Bolduc, M.-C., Renaud, J., Boctor, F.F., 2007. A heuristic for the routing and carrier selection problem. European Journal of Operational Research 183, 926–932. Bolduc, M.-C., Renaud, J., Boctor, F.F., Laporte, G., 2008. A perturbation metaheuristic for the vehicle routing problem with private fleet and common carriers. Journal of Operational Research Society 59, 776–787. Boudia, M., Louly, M.A.O., Prins, C., 2007. A reactive GRASP and path relinking for a combined production–distribution problem. Computers and Operations Research 34, 3402–3419. Campbell, A.M., Clarke, L.W., Savelsbergh, M.W.P., 2002. Inventory routing in practice. In: Toth, P., Vigo, D. (Eds.), The Vehicle Routing Problem, 2002, SIAM Monographs on Discrete Mathematics and Applications, Philadelphia, pp. 309– 330. Cordeau, J.-F., Gendreau, M., Laporte, G., 1997. A tabu search heuristic for periodic and multi-depot vehicle routing problem. Networks 30, 105–119. Fumero, F., Vercellis, C., 1999. Synchronized development of production inventory, and distribution schedules. Transportation Science 33, 330–340. Gendreau, M., Hertz, A., Laporte, G., 1994. A tabu search heuristic for vehicle routing problem. Management Science 40, 1276–1290. Kara, I., Laporte, G., Bektas, T., 2004. A note on the Lifted Miller–Tucker–Zemlin subtour elimination constraints for the capacitated vehicle routing problem. European Journal of Operational Research 158, 793–795. Lin, S., 1965. Computer solutions of the traveling salesman problem. Bell System Technical Journal, 2245–2269. Miller, C.E., Tucker, A.W., Zemlin, R.A., 1960. Integer programming formulation of traveling salesman problems. Journal of the Association for Computing Machinery 7, 326–329. Moin, N.H., Salhi, S., 2007. Inventory routing problems: a logistical overview. Journal of the Operational Research Society 58, 1185–1194. Osman, I.H., 1993. Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem. Annals of Operations Research 41, 421–451. Osman, I.H., Salhi, S., 1996. Local search strategies for the vehicle fleet mix problem. In: Rayward-Smith, V.J., Osman, I.H., Reeves, C.R., Smith, G.D. (Eds.), Modern Heuristic Search Methods. Wiley, Chichester, pp. 131–153. Renaud, J., Boctor, F.F., Laporte, G., 1996. A fast composite heuristic for the symmetric traveling salesman problem. INFORMS Journal on Computing 8, 134–143. Tarantilis, C.D., Kiranoudis, C.T., 2002. BoneRoute: An adaptive memory-based method for effective fleet management. Annals of Operations Research 115, 227–241. Zhao, Q.-H., Chen, S., Zang, C.-X., 2008. Model and algorithm for inventory/routing decision in a three-echelon logistics system. European Journal of Operational Research 191, 623–635.
Contents lists available at ScienceDirect
European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor
Production, Manufacturing and Logistics
A tabu search heuristic for the split delivery vehicle routing problem with production and demand calendars Marie-Claude Bolduc a,b, Gilbert Laporte a,c, Jacques Renaud a,b,*, Fayez F. Boctor a,b a
Interuniversity Research Center on Enterprise Networks, Logistics and Transportation (CIRRELT), Université Laval, Québec, Canada G1K 7P4 Faculté des Sciences de l’Administration, Université Laval, Québec, Canada G1K 7P4 c Canada Research Chair in Distribution Management, HEC Montréal, 3000, chemin de la Côte-Sainte-Catherine, Montréal, Canada H3T 2A7 b
a r t i c l e
i n f o
Article history: Received 17 July 2008 Accepted 6 May 2009 Available online 18 May 2009 Keywords: Vehicle routing Distribution calendar Production calendar Inventory Internal fleet Common carrier Split delivery Tabu
a b s t r a c t The purpose of this article is to propose a tabu search heuristic for the split delivery Vehicle Routing Problem with Production and Demand Calendars (VRPPDC). This new problem consists of determining which customers will be served by a common carrier, as well as the delivery routes for those served by the private fleet, in order to minimize the overall transportation and inventory costs. We first model this problem and then propose a simple decomposition procedure that can be used to provide a starting solution. Next, we introduce a new tabu search heuristic and we describe two new neighbor reduction strategies. Finally, we present the results of our extensive computational tests. According to these tests, our reduction strategies are efficient not only at reducing computing time but also at improving the overall solution quality. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction This paper proposes a tabu search heuristic for the split delivery Vehicle Routing Problem with Production and Demand Calendars (VRPPDC), defined as follows. Let G ¼ ðV; AÞ be a graph, where V ¼ f0; . . . ; ng is the vertex set and A ¼ fði; jÞ : i; j 2 V; i – jg is the arc set. Vertex 0 is a distribution centre (DC), and the remaining vertices represent customers. The DC must deal with an a priori production calendar provided by the factory. This calendar fixes the availability of each product over a time horizon of T periods. The production calendar ðg pt Þ defines the number of units of product p ¼ f1; . . . ; Pg made available at the beginning of period t ¼ f1; . . . ; Tg; g pt is not cumulative and may vary from period to period. Finished products p can be stored at the DC, incurring a periodic unit inventory holding cost hp . Each customer i has a demand calendar dpit , which specifies the number of new units of product p that must be received in period t (i.e., by the delivery date). A private fleet of m identical vehicles is based at the DC. At each period t, each customer i may be served by the private fleet vehicles ðk ¼ 1; . . . ; mÞ, by common carriers ðk ¼ 0Þ, or by any com* Corresponding author. Address: Interuniversity Research Center on Enterprise Networks, Logistics and Transportation (CIRRELT), Université Laval, Québec, Canada G1K 7P4. E-mail addresses: [email protected] (M.-C. Bolduc), [email protected] (G. Laporte), [email protected] (J. Renaud), [email protected] (F.F. Boctor). 0377-2217/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2009.05.008
bination of the two. Each vehicle in the private fleet can perform at most one route of L time units per period, has a capacity Q k expressed in terms of volume, and each product p has a volume bp . For the private fleet, cij and t ij represent the cost and the travel time associated with arc ði; jÞ. The common carriers are assumed to have an unlimited capacity Q 0 . They manage their routes and charge for each delivery as a function of distance and quantity. In practice, the common carriers charge a fixed cost f and a variable cost which is a step function of the distance li between the depot and customer i, and of the total quantity qit transported to i in period t. The lengths of the discret, and the varization intervals for distance and quantity are l and q iable costs are cl and cq . Generally, for a given origin-destination pair, the cost function related to the quantity delivered has up to seven step intervals (see for example the Accu-Rate 2008 rating software available for download on the USF Corporation web site http://www.usfc.com). The total cost charged by the common carrier for delivering quantity qit to customer i is therefore:
e: eit ¼ f þ cl dli =le þ cq dqit =q
ð1Þ
The VRPPDC consists of determining the routes of the private fleet, as well as the customers that must be served by a common carrier, in order to minimize the overall transportation and inventory costs over a given time horizon. The resulting delivery calendar must satisfy the following constraints: (i) for each period, each private fleet vehicle performs only one route, starting and ending at the
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
DC; (ii) the capacity Q of each private fleet vehicle must be respected; (iii) the maximal length L of each route must be respected; (iv) product availability for shipment or inventory is fixed according to the production calendar; and (v) customers must receive their products no later than their contracted delivery dates without penalty for the DC. The detailed formulation is provided in the appendix. The VRPPDC is a real-life problem encountered in the snowmobile industry where all retailer (customer) orders are contracted prior to the selling season, and before the manufacturers have launched their production and stored their supplies in the DC. The contract between the manufacturer and the retailers indicate the quantity of each product by the due date, i.e., the demand schedule. The DC is thus bound by contract to ship each product by the date specified in the demand schedule. Because delivered products become property of the retailer, there is no penalty for the DC if the products are shipped before their contracted delivery dates. Since the DC usually owns a limited fleet, common carriers are often used to increase transportation capacity in peak periods. The literature on the VRPPDC is rather limited. To our knowledge, most of the articles that focus on problems involving production calendars have dealt with third-party logistics for transportation. There are, however, some exceptions. Fumero and Vercellis (1999) have developed a mathematical VRPPDC model in which only a private fleet is used. Their objective was to determine both distribution and production calendars, with the latter being adapted to feed the former. This study used ‘‘on time” deliveries only. Boudia et al. (2007) have proposed a GRASP with a reactive mechanism for solving the single product problem. Inventory management is closely related to the classical, mono-product, Inventory Routing Problem – IRP (see Campbell et al., 2002; Moin and Salhi, 2007; Zhao et al., 2008). However, in the VRPPDC, quantities requested by customers are fixed over time by the demand calendars and not all products are available at the depot at the beginning of the planning horizon. This paper makes the following contributions. We first introduce a new complex VRP variant, which, to our knowledge, has never been studied before. We then model this problem as a mixed integer program, and we solve the model using a tabu search heuristic with split deliveries. Last, we present two new neighbor reduction strategies which, when tested, proved efficient not only at reducing computing time but also at improving overall solution quality. The remainder of this article is organized as follows. In Section 2, we present our VRPPDC local search heuristic. Section 3 describes our computational experiments and our results, and Section 4 offers our conclusions.
2. VRPPDC local search heuristic The VRPPDC can be formulated as a mixed integer model in which the objective function is non-linear because of common carrier cost function eit . The model, presented in the appendix, cannot be solved optimally for any realistic size instance. Since tabu search is a very powerful tool for solving the classical VRP (Tarantilis and Kiranoudis, 2002) and its variations (Osman and Salhi, 1996), we have developed a composite local search heuristic made up of three components: (i) construction of an initial feasible solution, (ii) a tabu search phase which allows split deliveries and the use of common carriers and (iii) a final improvement phase. The tabu search considers intermediate infeasible solutions by using penalty functions like those developed by Gendreau et al. (1994) and Cordeau et al. (1997). More specifically, vehicle capacity and route length may be violated during the search (soft constraints) but they must be respected in the final solution. However, the production calendar (i.e., inventory availability) and the demand calendar (i.e., inventory
123
quantities and the subsequent delivery dates) are always enforced during the search (hard constraints). These hard and soft constraints lead to an important difference with respect to the other tabu searches proposed in the literature: the double exploration of each neighbor, respectively called ‘‘complete switch” and ‘‘partial switch”. The complete switch evaluates the possibility of transferring all the product quantities for one customer from the original route (corresponding to a vehicle-period combination) to another one, while still respecting the hard constraints. Such a move may or may not generate feasible solutions if some soft constraints are not respected. The partial switch is a split delivery strategy, consisting of moving only the quantities that respect the target vehicle’s availability in terms of space and time. Such a partial transfer will always lead to a feasible route as it respects both the soft and hard constraints. The following subsections present the initial feasible solution, the tabu search, and the improvement phases. 2.1. Initial feasible solution If the quantities to be delivered are determined for each product and each period, the routing problem for each period corresponds to a multi-product version of the Vehicle Routing Problem with a Private fleet and a Common carrier (VRPPC) as defined by Bolduc et al. (2008). The VRPPC consists of minimizing the total cost of serving once a set of customers using either the limited private fleet vehicles or the common carrier and, for the private fleet, to perform the routes. Two heuristics have recently been developed for the VRPPC: the Selection, Routing and Improvement heuristic (SRI – Bolduc et al., 2007) and the Randomized construction, Improvement, Perturbation metaheuristic (RIP – Bolduc et al., 2008). SRI is a three phase heuristic: (1) selection of customers to be served by the common carrier, (2) construction of an initial solution, and (3) improvement procedure. This heuristic is quite fast and generates better results than a previous algorithm. RIP is a five step metaheuristic: (1) randomized savings construction phase, (2) a 4-opt* route improvement procedure, (3) 2*-interchange inter-route improvement procedure, (4) 2-add-drop improvement procedure, and (5) switch procedure. The 4-opt* heuristic is a route improvement procedure which may exchange up to four vertices (Lin, 1965; Renaud et al., 1996) and the 2*-interchange is a restricted version of the k-interchange procedure (Osman, 1993; Bolduc et al., 2008). These two procedures use sophisticated decision rules to reduce the neighborhood size and accelerate the search without decreasing solution quality. In the VRPPC computational results, RIP generated better results than SRI, but required much longer computational processing times. The heuristics developed for the VRPPC can be adapted to generate solutions to the single-period subproblems. Indeed, adaptations are necessary since these VRPPC heuristics were designed to handle only one product. The rest of this subsection presents a two-phase method designed to allow SRI and RIP to be used for the solution of singleperiod multi-product subproblems. The two-phases solve the problem for each period independently, and the best result of each phase is retained for each period. In each phase, the common carrier cost for each customer is calculated using (1). 2.1.1. Two-phase method using SRI and RIP Phase 1: Based on the demand calendars, generate a list of customers with a demand due in the current period. For each customer, generate a new temporary customer for each product p ordered (each customer may be split into, at most, P customers) and evaluate the common carrier cost
124
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
Table 1 Notation used in the tabu search algorithm. BðsÞ NðsÞ cðsÞ hðsÞ qðsÞ dðsÞ f ðsÞ gðsÞ s; s s0 s
a b
qi rikt g h d k
mp sikt v l j
attribute set for solutions ðBðsÞ ¼ fði; k; tÞg; i ¼ f1; . . . ; ng; k ¼ f1; . . . ; mg; t ¼ f1; . . . ; TgÞ neighborhood for solution s transportation cost for solution s holding inventory costs for solution s total excess quantity for solution s total excess duration for solution s total cost for solution s; f ðsÞ ¼ cðsÞ þ hðsÞ þ aqðsÞ þ bdðsÞ total penalized cost for solution s (i.e., f ðsÞ, to which is added a penalty term for the most frequently visited customers in order to diversify the search) solutions best solution identified in the current neighborhood NðsÞ best feasible solution identified penalty factor for capacity violation penalty factor for tour duration violation number of times a customer i has been added to the solution aspiration criterion of attribute ði; k; tÞ total number of iterations to be performed tabu duration parameter used to update a and b iteration counter movable quantities of product p last iteration for which attribute ði; k; tÞ is declared tabu number of iterations without improving s average number of iterations needed to pass from one feasible solution s0 to the next number of iterations a solution s0 remains infeasible
for this scenario. This cost is needed to solve the VRPPC. Solve this subproblem for each period, using either the SRI or the RIP heuristic. This strategy may produce large subproblems since each customer may be split a number of times. Phase 2: Generate a list of customers with a demand due in the current period. For each customer, create a single unified demand, corresponding to the sum of all product demands, and evaluate the common carrier cost for this scenario. Solve this subproblem for each period, using either the SRI or the RIP heuristic.
correspond to the objective function (2) in the appendix. Violations of the vehicle capacity or route duration constraints are penalized according to the following functions:
qðsÞ ¼
k¼1 t¼1
dðsÞ ¼
2
p¼1 i¼1
T 6X n X
6 4
#þ bp qpikt Q k 3þ
n X
i¼0
;
j¼0 j–i
7 t ij xijkt L7 5 ;
where ½xþ ¼ maxf0; xg. When s is feasible, qðsÞ and dðsÞ both equal to zero. 2.2.1.2. Attribute set management and movable quantities. Each solution s has an associated attribute set BðsÞ ¼ fði; k; tÞg, in which customer i is served by vehicle k (0: common carriers; 1 to m: private fleet) at period t. Moving from a solution s to a neighbor solution s0 2 NðsÞ modifies BðsÞ. The neighborhood NðsÞ of a solution s is composed of all solutions obtained by moving a quantity v piktk0 t0 0 associated to customer i from route ðk; tÞ to ðk ; t0 Þ, where 0 0 0 0 k; k ¼ f0; . . . ; mg; t; t ¼ f1; . . . ; Tg and k – k or t – t. For a customer i and the initial route ðk; tÞ under consideration, the movable quantities v piktk0 t0 depend on the type of switch. In the complete switch, when the target period t 0 is earlier than the original period t; v piktk0 t0 is the minimum movable quantity possible between the already planned delivered quantity ðqpikt Þ and the available inventory at the DC between periods t0 and tðspt Þ. However, when the period t0 is later than the original period t; v piktk0 t0 is simply set equal to already planned delivered quantities. In the partial switch, v piktk0 t0 also takes the available vehicle capacity into account. Note 0 that quantities may be split between vehicles k and k in the same period or over different periods t and t0 . To summarize, v piktk0 t0 may be formulated as follows: complete switch:
v piktk t
¼
8 > < min q
; min fs g ; if t0 < t ph pikt 0
> :q
h¼t ;...;t
otherwise;
pikt ;
partial switch:
v piktk t
0 0
2.2. Tabu search phase
m X
k¼1 t¼1
0 0
Repeating these strategies for each period leads to one single feasible initial solution, composed of all the routes obtained for each period. This solution can be enhanced through a tabu exploration phase and an improvement phase. Table 1 provides the notation used in the tabu search algorithm.
" m X T P X n X X
8 n P > > 0 ; if t 0 < t ; min fs g; Q q min q 0 > ph k pikt phk t < h¼t 0 ;...;t h¼1 ¼ n > P > > min qpikt ; Q k qphk0 t0 ; otherwise: : h¼1
The following subsections present the tabu search components, the complete VRPPDC tabu search algorithm and the neighbor reduction strategies. 2.2.1. Tabu search components The tabu search phase is the core of the VRPPDC heuristic (see Section 2.2.2). It is based on the following components which are described below: penalized cost function, attribute set management and movable quantities and finally, the exploration phase itself. 2.2.1.1. Penalized cost function. During the exploration phase a penalized cost function f ðsÞ ¼ cðsÞ þ hðsÞ þ aqðsÞ þ bdðsÞ is used, where
cðsÞ ¼
n X n X m X T X i¼0
hðsÞ ¼
j¼0 k¼1 t¼1 j–i
P X T X p¼1 t¼1
hp spt ;
cij xijkt þ
n X n X T X i¼0
j¼0 t¼1 j–i
eit xij0t ;
An aspiration criterion rikt is associated with each attribute ði; k; tÞ. It is the initial cost solution if ði; k; tÞ belongs to the attributes included in the initial solution and it is set equal to 1 for all the other attributes. The aspiration criteria are updated during the search. 0
2.2.1.3. Exploration phase. For given attributes ði; k; tÞ and ði; k ; t0 Þ and for a quantity v piktk0 t0 , the exploration phase performs the following steps: 1. Set the current solution s equal to the initial solution in the 0 neighborhood ðs ¼ sÞ. (In the following steps, ðk ; t0 Þ 2 s and ðk; tÞ 2 sÞ. 0 2. If customer i is not already in route ðk ; t0 Þ, insert this customer into the route and apply the 3-opt procedure (Lin, 1965). 0 3. Move v piktk0 t0 for all p from route ðk; tÞ to ðk ; t0 Þ. If v piktk0 t0 < qpikt (i.e., only some of the original delivery quantity is moved), prioritize the product that has the highest inventory cost. 4. If v piktk0 t0 ¼ qpikt , remove the customer from this route and reoptimize it using 3-opt.
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
5. For periods t to t0 , update the inventories spt and reevaluate the global inventory cost. 6. Reevaluate the routing cost for both routes. 7. If the solution s is non-tabu ðsik0 t0 < kÞ or if a best feasible solu0 tion for attribute ði; k ; t0 Þ is obtained (f ðsÞ < rik0 t0 and s is feasible), then apply the following operations: (a) If the current solution s is worse than the initial solution s in the neighborhood ðf ðsÞ > f ðsÞÞ; f ðsÞ is penalized pffiffiffiffiffiffiffiffiffiffi for diversification gðsÞ ¼ f ðsÞ þ 0:02 mnT ðcðsÞ þ hðsÞÞqi =k; otherwise, gðsÞ ¼ f ðsÞ. (b) If the current solution s is the new best neighbor ðgðsÞ < f ðs0 ÞÞ, set s0 ¼ s. When a solution s0 is retained as the best solution found in NðsÞ, customer i cannot be reinserted into route ðk; tÞ for the next h iterations. During these h iterations, a tabu move involving ði; k; tÞ is allowed only if its global cost is lower than the aspiration criterion rikt . This criterion is updated whenever a new best feasible solution is obtained. To diversify the search, the attributes repeatedly inserted into the solution are penalized by a term proportional to qi (see Step 7(a)). Because of the split quantities, penalizing the attribute ði; k; tÞ is not appropriate since using a penalty qikt may lead to moving the same customer to several different routes. To diversify the solution, different customers are instead moved at each iteration and their moves to other routes are penalized using qi . The other variables used in the determination of gðsÞ are the same as those mentioned by Cordeau et al. (1997), except for the constant factor, which is set to 0.02 instead of 0.015. After testing many values, we found that 0.02 generated the best results. 2.2.2. Complete VRPPDC Tabu search algorithm A step-by-step description of the VRPPDC heuristic is given below. 1. Obtain the initial feasible solution s using the two-phase method using SRI and RIP outlined in Section 2.1. 2. Set qi ¼ 0, for all i. 3. For each attribute ði; k; tÞ included in the initial best feasible solution s , set rikt ¼ f ðs Þ; otherwise, set rikt ¼ 1. 4. For each tabu search iteration ðk ¼ 1; . . . ; gÞ, perform the following operations: (a) Set the best neighborhood solution value f ðs0 Þ ¼ 1. (b) Set the initial solution in the neighborhood s ¼ s . 0 0 (c) For each attributes ði; k; tÞ and ði; k ; t 0 Þ, where k – k or 0 0 t – t , do the following. If t > t and the current solution s do respect the delivery dates, go to Step 4(c1); otherwise, go to Step 4(c2). (c1) COMPLETE SWITCH. 1. For each product p, calculate the movable quantity P v piktk0 t0 . If Pp¼1 v piktk0 t0 ¼ 0 (nothing can be moved), restart Step 4(c) using the next neighbor. 2. Apply the Exploration phase described in Section 2.2.1. (c2) PARTIAL SWITCH. 1. For each product p, calculate the movable quantity P v piktk0 t0 . If Pp¼1 v piktk0 t0 ¼ 0, restart Step 4(c) using the next neighbor. 2. Apply the Exploration phase described in Section 2.2.1 for all p, where v piktk0 t0 – 0. (d) Set the attribute ði; k; tÞ of the best neighbor s0 tabu ðsikt ¼ k þ hÞ and adjust the diversification criterion ðqi ¼ qi þ 1Þ. (e) If the best neighbor s0 is a new best feasible solution (i.e., f ðs0 Þ < f ðs Þ and s0 is feasible), perform the following operations:
125
(e1) Set s ¼ s0 . (e2) For each attribute ði; k; tÞ included in the solution s0 , update the aspiration criteria ðrikt ¼ minfrikt ; f ðs0 ÞgÞ. (f) If qðs0 Þ ¼ 0 (no quantity violation), then a ¼ a=ð1 þ dÞ. Otherwise, a ¼ ð1 þ dÞa. (g) If dðs0 Þ ¼ 0 (no duration violation), then b ¼ b=ð1 þ dÞ. Otherwise, b ¼ ð1 þ dÞb.
2.2.3. Neighbor reduction strategies The tabu search can be very time consuming due to the large size to the neighborhood NðsÞ and also to the non-linear cost functions that must be constantly reevaluated. Thus, we propose two neighbor reduction strategies designed to reduce the computing time. The first strategy, called Random, involves randomly selecting some neighbors from NðsÞ for evaluation, thus allowing only a predefined proportion of the total neighborhood to be considered. This strategy is quite simple and may be used on any kind of problem. By forcing a decrease in the number of neighbors, the strategy is able to decrease the time needed to find a solution. However, because the strategy does not necessarily select the best neighbor, it could lead to a different, and maybe worse, solution. The second strategy, called Distance, takes the problem configuration into account. The idea is simple: a given neighbor is evaluated only if the customer under consideration is within a fixed distance from the nearest customer on the route into which it will be inserted. The fixed distance is predefined as being equal to a fraction of the longest distance on the geographical map. This strategy makes it possible not to evaluate the insertion of a customer into a distant route. 2.3. Improvement phase When the tabu iterations are completed, the following improvement heuristic is applied using a first improvement strategy and until there is no improvement: 1. advance the delivery of complete blocks of product demand; 2. advance the delivery of partial blocks of product demand; 3. switch an internal customer (i.e., one served by a private vehicle) with an external customer (i.e., one served by a common carrier). 0
For each route ðk ; t 0 Þ, the first two steps verify whether it is possible to deliver, at t0 , a delivery planned for a later time t. Step 1 evaluates whether, in terms of inventory availability and the remaining vehicle capacity, a given customer’s deliveries, planned for a later period (i.e., in period t, where t > t0 Þ can be completely 0 integrated into route ðk ; t 0 Þ at an earlier period. After all current customers have been evaluated, this procedure is repeated for the partial blocks (Step 2). For a product that has been planned for delivery in quantity qpikt , the delivery of any quantity between 1 and qpikt may be moved to an earlier period. These two steps allow the global inventory cost to be reduced by filling the vehicles scheduled as much as possible. Step 3 is similar to the Osman 1–1 exchange (1993). The difference between Step 3 and the original method is that the switches between customers served by a private fleet vehicle and those served by common carriers are the only ones evaluated. When a customer switch is performed, the entire product quantities for these two customers are moved. In other words, there are no partial transfers. A switch can be conducted for any vehicle-periodcustomer combination, as long as the problem constraints are respected and the global cost (i.e., transportation and inventory costs) is reduced. Then, the 3-opt heuristic is applied to each modified route.
126
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
ing the demand for five periods of any product in any one period, the production of each period was assigned entirely to the most urgent product remaining to be delivered, chosen with respect to the minimum ratio between the current stock level and the demand to be filled. The parameters of the common carrier cost function were set ¼ 0:25Q . as: f ¼ 60; cl ¼ 50; l ¼ 25; cq ¼ 15; q When determining the number and capacity of the private fleet vehicles, we have assumed that, on average, about 75% of the demand would be filled by this fleet, while the other 25% would be filled by a common carrier. The number of vehicles was initially set equal to 2, and this number was increased until the demand could be covered by the number of vehicles selected, each with a capacity between 100 and 150 units. The exact vehicle capacity was determined with respect to the demands. For example, a demand of 500 units per period would require a total fleet capacity of at least 375 (75% of 500). Since each vehicle is assumed to have a capacity between 100 and 150 units, two vehicles would have a total capacity within the interval [200, 300], while three vehicles would have a capacity of [300, 450]. Thus, for a demand of 375, three vehicles would be necessary, each with a capacity of 125 (i.e., 375 divided by 3). The volume of each product was set to 1, and the vehicle speed was set to 80 km/h, the maximum route length to 13 time units (hours), and the variable distance cost to $2/km. Table 2 presents the number of vehicles, their capacity and the total demand for each instance.
3. Computational results We first describe our test instances, then present the tabu search parameters, and finally report our results. 3.1. Test instances For each of the 100 instances, 50 customer coordinates were randomly generated on a plane measuring 240 km by 200 km. The number of products was set equal to 3; the inventory costs per product unit were $0.15, $0.25 and $0.35, respectively; and the planning horizon was equal to 10 periods. To determine the demand calendar for each customer, the following rules were applied. First, each customer demand level was classified accordingly to the following probabilities: 25% had a low demand per period, 50% had an average demand per period (i.e., twice the low demand), and 25% had a high demand per period (i.e., four times the low demand). The upper bounds of low demand products were set equal to 2, 4 and 3 units for product 1, 2 and 3, respectively. For each period, demands were randomly generated according to a continuous uniform distribution between 0 and an upper bound based on the customer category. For each customer, the ordering interval was randomly selected between one and four periods; the day of the first order was then selected. The quantity of each order was the sum of the customer’s daily demand over the ordering interval. The production calendar was designed according to the customer demand calendars. For each product, the initial inventory at the DC was equal to the average demand over five periods, and the production lot size was the average demand divided by the production frequency, which was set randomly between one and five periods. Since the factory was assumed capable of produc-
3.2. Tabu search parameters The number of iterations of the tabu search algorithm in the VRPPDC heuristic was set g ¼ 100; 000 iterations. Two parameter sets (A and B) were also tested.
Table 2 Some instance characteristics. Instance
m
Q
Total demand
Instance
m
Q
Total demand
Instance
m
Q
Total demand
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
3 3 2 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
110 130 140 120 100 110 120 110 110 110 110 110 120 150 120 120 110 120 100 100 100 130 120 120 120 100 110 110 120 110 100 100 110 130
4230 5089 3692 4550 4015 4242 4789 4146 4090 4188 4292 4373 4601 3774 4601 4683 4404 4459 3989 3977 4063 5121 4606 4656 4704 3924 4439 4161 4457 4391 3974 3888 4253 4982
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
3 3 3 3 2 2 3 3 3 3 2 2 3 3 3 3 3 3 3 3 3 2 3 3 2 3 3 3 3 3 3 3 3
110 110 100 100 140 150 110 140 110 120 130 140 100 140 130 120 120 110 110 110 120 140 110 120 150 110 110 120 110 120 120 100 110
4216 4190 3998 4064 3643 3815 4226 5291 4191 4850 3324 3570 3887 5483 4896 4502 4452 4211 4249 4254 4515 3716 4290 4590 3867 4103 4153 4575 4212 4581 4483 3900 4392
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3
110 120 110 120 110 110 110 120 120 100 110 110 110 120 130 110 130 110 110 110 120 110 120 100 120 120 100 120 110 150 100 100 110
4422 4475 4287 4528 4300 4146 4343 4690 4787 3991 4178 4266 4220 4541 5015 4420 3283 4099 4257 4251 4490 4284 4534 4004 4517 4615 3922 4608 4253 3827 4027 4070 4214
127
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
3.2.1. Parameter set A Parameter set A used the VRPPDC heuristic as described in Section 2.2.2, with the following parameter values: a ¼ 1; b ¼ 1; h ¼ ½7:5log10 n and d ¼ 1:5. 3.2.2. Parameter Set B This second set was designed to calibrate the search parameter values. First, the penalty factors were set at a ¼ 0:5 and b ¼ 0:5. Because the VRPPDC heuristic starts from an initial feasible solution, lower penalty factors can be used at first. In addition, the problem considered in this paper can lead to many tabu iterations with infeasible solutions. Since many infeasible solutions may be generated in a neighborhood before a new feasible solution is obtained, we have added a step in our algorithm before Step 4(f) to update the penalty factors a and b: ðf 0 Þ If s0 is feasible, or if j 6 l, then d ¼ 1:5; otherwise, d ¼ 1 until s0 becomes feasible. To prevent the tabu search algorithm from being stuck in a local optimum with only a few neighbors, we also adjusted the criterion in Step 4(d) (Section 2.2.2) to read: h ¼ ½7:5log10 n þ dn=100e. This variable tabu duration allows the search to visit other neighbors. 3.3. Results The heuristic was coded in Microsoft VB.Net 2005 and run under Windows XP on a personal computer with a Xeon 3.6 GHz processor and 1.00 Gb of RAM. All reported times are expressed in seconds, and all statistics represent the average results over a set of 100 test instances. The following results are compared with the best-known solutions obtained for all instances (see Table 3). For each instance, we have recorded the best solution found during all our tests which include the computational experiments used in the paper but also some other tests conducted with parameters that have been tried during the development of the algorithms but not presented in the paper. This section presents our computational experiments with our tabu search heuristic. We first present the results of the initial feasible solution and improvement phase.
The neighbor reduction strategies are then analyzed. Finally, results for the complete VRPPDC heuristic are presented. 3.3.1. Evaluation of the initial feasible solution and improvement phase Our two-phase method (Section 2.1.1) used to obtain the initial feasible solution is clearly inefficient, since, as shown in Table 4, the average deviation above the best-known solution is 57.23% for the RIP-based heuristic and 64.01% for SRI-based heuristic. The improvement phase (Section 2.3) yields slightly better results: 49.15% for RIP and 51.28% for SRI. However, RIP solution times are much longer than those of SRI, with an average solution time of 301 s compared to only two seconds for SRI. The results reported in this section for the VRPPDC heuristic, called SRITabu, were initiated from the SRI-based solution which was much faster. 3.3.2. Neighbor reduction strategies To evaluate the proposed neighbor reduction strategies, all 100 instances were solved with SRITabu and 500 iterations, using parameter set A and the two neighbor reduction strategies. The goal of the evaluation was to demonstrate that, using our neighbor reduction strategies, solutions of the same quality could be obtained in less time. Table 5 present the results of our tests. The first column in the table presents the average deviation for the original SRITabu after 500 iterations. Columns 2, 3 and 4 contain the results obtained with the Random reduction strategy that eliminates 25, 50 and 75% of the neighbors. The best quality result for this strategy was obtained with a random 50% reduction, producing an average deviation of 11.57%. This result is close to that of the original SRITabu, which produced an average deviation of 10.88%, but it is much faster, with an average of 16 s compared to 29. Columns 5, 6 and 7 contain the results using the Distance reduction strategy. The first test discarded all neighbors that were not closer than 15% of the longest distance on the map (column 5), the second, those that were not closer than 25% of the longest distance (column 6), and the last, those that were not closer than 35% (column 7). Overly restricting the distance produced a significant deterioration in solution quality: the deviation obtained with
Table 3 Best-known solutions. Instance
Best
Instance
Best
Instance
Best
Instance
Best
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
23,113.57 25,391.41 18,268.49 23,585.96 25,173.42 27,788.49 27,581.27 24,156.73 21,783.57 25,171.74 28,200.50 28,937.65 26,087.67 18,151.88 25,250.77 27,917.18 27,729.59 22,731.64 22,629.94 25,401.72 26,403.15 26,359.05 27,477.55 24,901.55 26,618.87
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
24,745.46 27,181.70 25,466.19 23,228.56 29,169.52 24,635.79 25,759.79 22,809.43 23,604.40 26,281.08 24,837.63 26,247.15 25,163.93 19,983.13 20,232.92 24,333.88 24,578.88 24,127.50 24,602.39 20,382.55 20,443.80 22,683.97 26,214.30 24,993.88 24,243.11
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
27,856.98 24,289.10 20,583.28 24,700.93 27,918.29 19,058.70 25,576.13 27,075.20 20,262.76 23,327.91 24,512.32 24,346.14 25,515.31 27,200.36 25,266.38 26,227.86 25,801.10 28,676.45 25,742.31 24,900.09 25,947.90 25,470.19 23,799.70 23,840.44 26,594.19
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
26,828.62 27,672.02 22,771.41 26,817.43 24,737.76 26,776.09 23,887.33 28,674.86 19,024.32 24,307.21 23,319.70 25,379.44 25,831.23 23,950.72 25,659.14 25,024.39 22,336.13 25,274.56 25,586.53 25,879.00 23,210.87 19,520.18 22,771.17 26,181.31 24,832.96
128
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130
Table 4 Results for the initial feasible solution and the improvement phase heuristics. Initial feasible solution
RIP
SRI
Improvement phase
Without
With
Without
With
Average deviation from the best-known solution (%) Average time (s) Number of best solutions
57.23
49.15
64.01
51.28
301 0
301 0
2 0
2 0
Table 5 Comparison of the neighbor reduction strategies for SRITabu with 500 iterations. Original SRITabu
Neighbor reduction strategy Random
% of reduction Average deviation from the best-known solution (%) Average time (s) Number of best solutions
Distance
0 10.88
25 11.68
50 11.57
75 12.16
15 18.48
25 10.36
35 10.57
29 0
23 0
16 0
9 0
18 0
21 0
23 0
distances closer than 15% is 18.48%, while 10.36% was obtained with distances closer than 25%. Both the 25% and 35% distance reduction strategies produced better results than the SRITabu by itself and were also faster, although the distance reduction strategy is, on average, slower than the random strategy. The next subsection presents the results of our tests using the best parameter combination for each neighbor reduction strategy: ‘‘closer than 25%” for the distance strategy and ‘‘50%” for the random reduction strategy. 3.3.3. VRPPDC heuristic: SRITabu results Tables 6 and 7 show the average results obtained by the SRITabu heuristic over the 100 instances, using 500 to 100,000 iterations. The tables are separated horizontally into three main sections: the results obtained with the SRITabu heuristic only, the results obtained with the ‘‘closer than 25%” distance reduction strategy, and then those obtained with the ‘‘50%” random reduction strategy. Each horizontal section shows the worst, best and average deviation from the best-known solution, the average time in
seconds, the number of best solutions found, and the number of solutions found that deviated less than 0.5% from the best solutions. Table 6 uses parameter set A, while Table 7 employs the parameter set B. The best results obtained by SRITabu heuristic with parameter set A are shown in the ‘‘100,000 iterations” column of Table 6. For the SRITabu only (i.e., no neighbor reduction strategy), the average deviation is 2.50%, and the solutions were obtained in an average of 5725 s. This ‘‘no strategy” method produced six best solutions and 24 results with less than 0.5% of deviation. The distance strategy, on the other hand, generated an average deviation of 2.62% and used half the computing time (i.e., an average of 2719 s). This strategy produced four best solutions and 11 results with less than 0.5% of deviation. However, the random strategy generated the best overall results. With only 10,000 iterations, this strategy yielded a 2.51% deviation, a result similar to those obtained by the two others methods with 100,000 iterations. Moreover, this strategy produced these results much faster, within only 300 s on average. With 100,000 iterations, the random strategy obtained the lowest average deviation (1.12%) almost as fast as the distance strategy obtained a much higher average deviation (i.e., an average of 2879 s) and obtained four best solutions and 29 solutions with less than 0.5% of deviation. As these results show, the random reduction strategy halved both the average deviation and the computing time with respect to the original SRITabu application, using parameter set A. From an empirical point of view, these results are significant. Table 7 presents the results for parameter set B. Using this strategy, almost no improvement was obtained over 50,000 iterations. For parameter set B, the best quality results were obtained by the original SRITabu with an average deviation of 1.77%. This strategy generated six best solutions and 14 results with less than 0.5% of deviation. However, this strategy was also the most time consuming of the three, taking 2702 s on average. In fact, both the random and distance reduction strategies obtained results faster, with an average of 1433 and 1420 s, respectively. Unfortunately, these solutions have larger deviation of 2.42% and 2.59%, respectively. With parameter set B, only the original SRITabu performed better than the various scenarios using parameter set A. Examining both tables, strategy by strategy allows parameter set A to be compared to B. In terms of quality, the ‘‘no strategy” results obtained with A at 100,000 iterations were obtained with B at
Table 6 Results of the SRITabu heuristic with parameter set A.
No strategy Deviation
Time Number of solutions Distance 25% Deviation
Time Number of solutions Random 50% Deviation
Time Number of solutions
Iterations
500
1000
2000
5000
10,000
15,000
20,000
50,000
100,000
Worst Best Average Average (s) Best Under 0.5%
18.11% 4.65% 10.88% 29 0 0
17. 03% 2.13% 7.36% 57 0 0
17.03% 0.33% 5.24% 115 0 1
17.03% 0.14% 3.73% 290 0 2
17.03% 0.00% 3.18% 579 1 5
17.03% 0.00% 2.96% 871 1 7
17.03% 0.00% 2.83% 1163 1 11
17.03% 0.00% 2.56% 2892 5 21
17.03% 0.00% 2.50% 5725 6 24
Worst Best Average Average (s) Best Under 0.5%
18.69% 4.36% 10.36% 21 0 0
15.19% 3.53% 7.29% 39 0 0
14.10% 1.86% 5.28% 72 0 0
14.10% 0.80% 3.95% 162 0 0
14.10% 0.05% 3.42% 304 0 1
14.10% 0.00% 3.22% 444 1 2
14.10% 0.00% 3.08% 582 1 3
14.10% 0.00% 2.75% 1.385 4 10
14.10% 0.00% 2.62% 2.719 4 11
Worst Best Average Average (s) Best Under 0.5%
20.34% 5.94% 11.57% 16 0 0
16.09.% 2.87% 7.60% 32 0 0
10.37% 1.29% 5.24% 65 0 0
8.45% 0.13% 3.25% 153 0 2
8.11% 0.13% 2.51% 300 0 3
7.33% 0.13% 2.14% 443 0 4
7.33% 0.13% 1.88% 589 0 4
7.12% 0.00% 1.40% 1.475 2 11
7.12 0.00% 1.12% 2.879 4 29
129
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130 Table 7 Results of the SRITabu heuristic with parameter set B.
No strategy Deviation
Time Number of solutions Distance 25% Deviation
Time Number of solutions Random 50% Deviation
Time Number of solutions
Iterations
500
1000
2000
5000
10,000
15,000
20,000
50,000
100,000
Worst Best Average Average (s) Best Under 0.5%
36.81% 4.15% 11.90% 30 0 0
22.80% 2.92% 8.09% 59 0 0
22.80% 0.97% 5.50% 117 0 0
14.39% 0.00% 3.15% 282 1 3
9.15% 0.00% 2.33% 550 2 6
8.32% 0.00% 2.04% 817 3 9
8.32% 0.00% 1.88% 1095 4 12
8.32% 0.00% 1.77% 2702 6 14
8.32% 0.00% 1.77% 5165 6 14
Worst Best Average Average (s) Best Under 0.5%
30.13% 4.84% 11.16% 18 0 0
16.80% 2.96% 7.82% 34 0 0
10.86% 1.89% 5.32% 64 0 0
7.70% 0.39% 3.49% 145 0 1
7.70% 0.53% 2.84% 274 0 0
7.70% 0.53% 2.70% 409 0 0
7.70% 0.53% 2.63% 550 0 0
7.70% 0.28% 2.59% 1.420 0 1
7.70% 0.28% 2.59% 2.874 0 1
Worst Best Average Average (s) Best Under 0.5%
21.93% 5.42% 11.92% 16 0 0
16.52% 2.87% 8.03% 31 0 0
15.17% 1.37% 5.69% 60 0 0
10.07% 0.48% 3.58% 144 0 1
9.57% 0.48% 2.70% 285 0 1
9.57% 0.48% 2.55% 427 0 1
9.57% 0.48% 2.53% 572 0 1
9.57% 0.48% 2.42% 1.433 0 1
9.57% 0.48% 2.42% 2.800 0 1
10,000 iterations and with only one tenth of the computing time. However, the same improvement is not true for the random and the distance strategies. Overall, the best combination appears to be the random neighbor reduction strategy applied with parameter set A. This is the only combination that produces average results under 2% within less than 600 s, and it leads to the best overall performance, with an average deviation of 1.12%. 4. Conclusion
uikt : upper bound for the load of vehicle k upon leaving customer i in period t; spt : inventory of product p at the DC at the end of period t; qpikt : quantity of product p leaving the DC to be delivered to cusP tomer i with vehicle k in period t. Note that qit ¼ Pp¼1 qpi0t ; t 2 f1; . . . ; Tg. The model is expressed as follows:
Minimize
n X n X m X T X i¼0
We have introduced the split delivery Vehicle Routing Problem with Production and Demand Calendars, a rich variant of the VRP which, to our knowledge, has never previously been studied. The problem was modeled as a mixed integer program and solved by means of a tabu search heuristic with split deliveries. Our results demonstrate that it is possible to obtain an excellent solution within a reasonable amount of time. We have also introduced several neighbor reduction strategies which were proved efficient not only at reducing computing time but also at improving the overall solution quality.
þ
j¼0 k¼1 t¼1 j–i
P X T X
This work was partially supported by the Canadian Natural Sciences and Engineering Research Council under Grants OPG0036509, 39682-05 and OPG0172633. This support is gratefully acknowledged. Thanks are due to three referees for their valuable comments.
n X n X T X i¼0
eit xij0t
j¼0 t¼1 j–i
hp spt
ð2Þ
p¼1 t¼1
subject to spt ¼ sp;t1 þ g pt
n X m X
qpikt
i¼1 k¼0
ðp 2 f1; . . . ; Pg; t 2 f1; . . . ; TgÞ m X t X
qpikh P
k¼0 h¼1
Acknowledgements
cij xijkt þ
t X
dpih
h¼1
ðp 2 f1; . . . ; Pg; i 2 f1; . . . ; ng; t 2 f1; . . . ; TgÞ n X
m X
x0jkt 6 m ðt 2 f1; . . . ; TgÞ
ð4Þ ð5Þ
j¼1 k¼1 n X
xhjkt ¼
j¼0 j–h
n X
xihkt
i¼0 i–h
Appendix A. Problem formulation
ðh 2 f0; . . . ; ng; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
In addition to the notation in introduction, define the following variables:
n X
n X
i¼0
j¼0 j–i
8 > < 1 if vehicle k visits a vertex j immediately xijkt ¼ after vertex i at period t > : 0 otherwise; 1 if vehicle k is used in period t ykt ¼ 0 otherwise;
ð3Þ
P X p¼1
t ij xijkt 6 L ðk 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
qpikt 6 Q k
n X
ð6Þ ð7Þ
xjikt
j¼0 j–i
ði 2 f1; . . . ; ng; k 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ
ð8Þ
130
M.-C. Bolduc et al. / European Journal of Operational Research 202 (2010) 122–130 n X n X i¼0
(10) ensure that the vehicle capacity is never exceeded. Constraints (11) were introduced by Miller et al. (1960) to eliminate subtours. Kara et al. (2004) describe a lifting procedure for (11).
xijkt 6 ðn þ 1Þykt
j¼0 j–i
ðk 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ P X n X
ð9Þ
bp qpikt 6 Q k ykt
p¼1 i¼1
ðk 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
ð10Þ
uikt ujkt þ ðn 1Þxijkt 6 n 2 ði; j 2 f1; . . . ; ng; i – j; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ ð11Þ xijkt 2 f0; 1g ði; j 2 f0; . . . ; ng; i – j; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
ð12Þ
ykt 2 f0; 1g ðk 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ
ð13Þ
uikt P 0 ði 2 f1; . . . ; ng; k 2 f1; . . . ; mg; t 2 f1; . . . ; TgÞ
ð14Þ
qpikt P 0 ðp 2 f1; . . . ; Pg; i 2 f1; . . . ; ng; k 2 f0; . . . ; mg; t 2 f1; . . . ; TgÞ
ð15Þ
spt P 0 ðp 2 f1; . . . ; Pg; t 2 f0; . . . ; TgÞ:
ð16Þ
The objective function minimizes the routing costs of the private fleet, the non-linear costs of the common carrier (see (1)), and the DC’s inventory costs. Constraints (3) maintain the product balance by adjusting the DC’s inventory level according to incoming production and outgoing shipments. Constraints (4) ensure that the contracted customer demand calendars are satisfied. Constraints (5) specify that at most m private fleet vehicles can be used each day, while constraints (6) indicate that the same vehicle k must arrive at and leave from customer h. Constraints (7) require that each private fleet route lasts at most L time units. Constraints (8) and (9) respectively specify that deliveries can only be made if a tour exists already and that a tour exists only if a vehicle is used. Constraints
References Bolduc, M.-C., Renaud, J., Boctor, F.F., 2007. A heuristic for the routing and carrier selection problem. European Journal of Operational Research 183, 926–932. Bolduc, M.-C., Renaud, J., Boctor, F.F., Laporte, G., 2008. A perturbation metaheuristic for the vehicle routing problem with private fleet and common carriers. Journal of Operational Research Society 59, 776–787. Boudia, M., Louly, M.A.O., Prins, C., 2007. A reactive GRASP and path relinking for a combined production–distribution problem. Computers and Operations Research 34, 3402–3419. Campbell, A.M., Clarke, L.W., Savelsbergh, M.W.P., 2002. Inventory routing in practice. In: Toth, P., Vigo, D. (Eds.), The Vehicle Routing Problem, 2002, SIAM Monographs on Discrete Mathematics and Applications, Philadelphia, pp. 309– 330. Cordeau, J.-F., Gendreau, M., Laporte, G., 1997. A tabu search heuristic for periodic and multi-depot vehicle routing problem. Networks 30, 105–119. Fumero, F., Vercellis, C., 1999. Synchronized development of production inventory, and distribution schedules. Transportation Science 33, 330–340. Gendreau, M., Hertz, A., Laporte, G., 1994. A tabu search heuristic for vehicle routing problem. Management Science 40, 1276–1290. Kara, I., Laporte, G., Bektas, T., 2004. A note on the Lifted Miller–Tucker–Zemlin subtour elimination constraints for the capacitated vehicle routing problem. European Journal of Operational Research 158, 793–795. Lin, S., 1965. Computer solutions of the traveling salesman problem. Bell System Technical Journal, 2245–2269. Miller, C.E., Tucker, A.W., Zemlin, R.A., 1960. Integer programming formulation of traveling salesman problems. Journal of the Association for Computing Machinery 7, 326–329. Moin, N.H., Salhi, S., 2007. Inventory routing problems: a logistical overview. Journal of the Operational Research Society 58, 1185–1194. Osman, I.H., 1993. Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem. Annals of Operations Research 41, 421–451. Osman, I.H., Salhi, S., 1996. Local search strategies for the vehicle fleet mix problem. In: Rayward-Smith, V.J., Osman, I.H., Reeves, C.R., Smith, G.D. (Eds.), Modern Heuristic Search Methods. Wiley, Chichester, pp. 131–153. Renaud, J., Boctor, F.F., Laporte, G., 1996. A fast composite heuristic for the symmetric traveling salesman problem. INFORMS Journal on Computing 8, 134–143. Tarantilis, C.D., Kiranoudis, C.T., 2002. BoneRoute: An adaptive memory-based method for effective fleet management. Annals of Operations Research 115, 227–241. Zhao, Q.-H., Chen, S., Zang, C.-X., 2008. Model and algorithm for inventory/routing decision in a three-echelon logistics system. European Journal of Operational Research 191, 623–635.