A reactive GRASP and path relinking for a combined production–distribution problem

Computers & Operations Research 34 (2007) 3402 – 3419 www.elsevier.com/locate/cor

A reactive GRASP and path relinking for a combined production–distribution problem M. Boudia∗ , M.A.O. Louly, C. Prins ISTIT, University of Technology of Troyes, BP 2060, 10010 Troyes, France Available online 15 March 2006

Abstract An NP-hard production–distribution problem for one product over a multi-period horizon is investigated. The aim is to minimize total cost taking production setups, inventory levels and distribution into account. An integer linear model is proposed as a compact problem specification but it cannot be solved to optimality for large instances. Instead of using a classical two-phase approach (production planning and then route construction for each day), metaheuristics that simultaneously tackle production and routing decisions are developed: a GRASP (greedy randomized adaptive search procedure) and two improved versions using either a reactive mechanism or a path-relinking process. These algorithms are evaluated on 90 randomly generated instances with 50, 100 and 200 customers and 20 periods. The results confirm the interest of integrating production and distribution decisions, compared to classical two-phase methods. Moreover, reaction and path-relinking give better results than the GRASP alone. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Greedy randomized adaptive search procedure; Path relinking; Production planning; Vehicle routing

1. Introduction and literature review Distribution costs often represent a significant fraction of the actual cost of a finished product. In some sectors like mineral waters, they even dominate production costs. This is why, ideally, any cost reduction effort should encompass production planning and distribution planning. However, in industry, decisions concerning production and distribution to customers or retailers are often made separately, and most companies handle these two levels sequentially. Even when distribution costs are considered at production time, they are merely assumed to be proportional to the amount produced and to the distance to customers. Such approximation may lead to suboptimal decisions in less-than-truckload (LTL) transportation, in which routing costs critically depend on the order of visits to customers in the trips. The same dichotomy can be observed in research. Production planning and distribution planning have been widely but separately investigated. On one hand, many vehicle routing problems have been analysed (see [1] for a survey), even for stochastic cases [2], additional inventory constraints [3,4] or multi-period horizons [5], but production decisions are usually ignored or only partially handled. On the other hand, production planning has been studied in detail but without taking distribution into account [6,7]. It is true that supply chain management research tries to design more ∗ Corresponding author.

E-mail addresses: [email protected] (M. Boudia), [email protected] (M.A.O. Louly), [email protected] (C. Prins). 0305-0548/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2006.02.005

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3403

global approaches that encompass both suppliers and end-customers, but it mainly concentrates on qualitative or strategic issues like supplier and market selection or outsourcing [8]. In comparison, there are few quantitative studies concerning the operational and tactical decision levels in supply chains and they are relatively recent [9,10]. The number of such studies dealing with hard combinatorial optimization problems is even more limited [11,12]. This paper belongs to the last category. A production plan and distribution trips are simultaneously determined for a plant that manufactures one single product, over several periods, to minimize total cost (production, inventory and transportation). This approach is desirable and profitable in sectors with relatively stable or predictable demands and when production costs are of the same order as transportation costs, e.g., livestock feed, fertilizers, or milk collection for a dairy product factory. Note that such collection problems are easily tackled in our approach, by inverting the flow of product between the plant and customers. The literature dealing with coordination problems between production and distribution is still scarce. Some papers present general but qualitative methods to choose the suppliers and the kind of relationship in a supply chain [13]. Other papers consider the case of truckload deliveries, thus avoiding the difficulties raised by the vehicle routing aspects [14,15]. Inventory routing also deal with the distribution of a product to a set of customers, over a multi-period horizon, to avoid shortages. For instance, Bard et al. [16] studied a problem inspired by propane distribution, in which vehicles may reload at the plant or at satellite depots. However, such problems do not include production activities. Chandra and Fisher in [11] are very close to our problem defined in Section 2. The main differences in their paper are the production of several products, an unlimited fleet size, no inventory holding cost at the plant, and split deliveries (i.e., multiple deliveries can be made by different vehicles to the same customer). These authors propose an uncoupled (two-phase) approach that computes one production plan and then one distribution plan, and a coupled (integrated) approach. In fact, the latter consists of searching cost-reducing changes in the two plans computed by the first approach. Savings between 3% and 20% are reported on instances with up to 10 products, 50 customers and 10 periods. It should be noted that the instances they use are weakly constrained: one third of them have a total demand per day limited to 85% of production capacity, another third 60%, while the last third considers an unlimited production capacity. Another paper by Chandra [17] concerns the preparation of orders in a regional warehouse to satisfy the demands of customers in the same region. If an order cannot be satisfied, the warehouse may transmit it to a higher echelon (e.g., a factory) but this induces a fixed cost. The tests conducted on 33 test-problems show a cost reduction ranging from 5% to 14% when distribution and order preparation are coordinated. Fumero and Vercellis [12] deal with a problem quite close to the one studied by Chandra and Fisher. Their solution method, based on Lagrangean relaxation, can be applied only to smaller instances, with up to 12 customers, 10 products and 8 periods. Here again, the algorithm was compared with an uncoupled approach and significant savings were obtained. Ertogral et al. [18] also handle the same kind of problem. Like Chandra and Fisher, they start with a decomposition of the global problem into a production planning problem and a distribution problem. Afterwards, they progressively reintroduce coordination constraints to make the results of the two subproblems compatible. Other problems of coordination between production and distribution in a broader sense can be found in postal activities, e.g., Metters [19] investigated the coordination between a sorting centre and mail distribution. The objective not only includes the total cost, but also the reduction of routing delays. Bramel et al. [20] treated a problem of cooperation between several factories that make components or subassemblies for each other. However, the routing aspects are discarded, since truckload transportation is assumed between factories. The remainder of this paper starts with the problem statement and an integer programming model in Section 2. A preliminary greedy randomized adaptive search procedure (GRASP) is introduced in Section 3. This GRASP is improved using a reactive mechanism in Section 4 and a path relinking method in Section 5. Computational results are presented in Section 6. A conclusion with some research perspectives closes the paper.

2. Problem statement and integer programming model The problem is defined on a complete weighted and undirected network G = (N, E, C). N is a set of n + 1 nodes indexed from 0 onwards. Node 0 is a plant with a limited fleet F of identical vehicles of capacity W. The set N  = N \{0} contains n customers. E is the edge set and the weight Cij = Cj i on each edge (i, j ) is the travelling cost between nodes i and j. See Fig. 1 for an example with 11 customers and 3 days.

3404

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

6

3

6

3

5 4

2

5

4

2

4

2 7

7 1

7

1

11

6

3

5

1

11

8

11

8

10

10

10 9

9

8 9

Day 3 Day 2 Day 1 1 2 3 4 5 6 7 8 9 10 11 Demands produced and delivery trips of day 1

1

2 3 4 5 6 7 8 9 10 11 Demands produced and delivery trips of day 2

1

2 3 4 5 6 7 8 9 10 11 Demands produced and delivery trips of day 3

Fig. 1. Graphical problem representation.

The plant manufactures one single product to supply the customers over a planning horizon T subdivided into periods or “days”. If production is decided on day t, a setup cost s is incurred. The daily production capacity is limited to a value PPC (plant production capacity). The product may be either stored at the plant, with a storage cost h per unit of product and per period, or shipped to customers. The plant cannot store more than PSC (plant storage capacity) units of products. Like in vehicle routing literature, it is important to note that in reality there can be several products, but it is assumed that they do not create resource conflicts in production (they have dedicated production lines for instance) and that they can be mixed in the same vehicles, like bags or pallets: hence, they can be aggregated into one single equivalent product in the model. Each customer i has a known demand qit for each day t and a limited storage capacity CSCi (customer storage capacity). The local storage cost is paid by the customer and not considered in the problem. Shortages are not allowed. Each customer can be serviced at most once per day and each vehicle may perform at most one trip per day. By convention, without loss of generality, deliveries take place in the morning, before a production run at the plant and before the customers’ daily consumption. Each day’s production is available the next morning. This kind of convention is common in production planning models, to write inventory balance equations and to know when storages costs must be counted. In our case, the goal of the convention is also to impose an order between production and distribution activities in each period, to avoid possible hard scheduling problems if vehicles wait for the product or vice versa. The convention can easily be changed, with minor indexing changes. The periods are indexed from 0 to . The first possible production day is in fact day 1, day 0 being a dummy period for the initial inventories. In particular, the eventual deliveries in day 1 come from the initial stock of the plant. In the sequel, T  = T \{0} denotes the genuine periods, i.e., without the dummy one. The new inventory levels (plant and customers) are measured at the end of each day. For safety reasons, at the beginning of day t each customer must have enough stock to ensure an integer number of consumption days, the minimum amount being the demand qit . This implies that each demand can be delivered early but not late. To avoid shortages, for each customer i demand for day t must be delivered before the demand for day t + 1. This FIFO rule is not mandatory but recommended in many cases, e.g., perishable products. The objective is to determine, for each period, the amount to be produced, the quantities to be shipped to each customer and the associated delivery trips, in order to minimize the total cost (setups, inventories and transportation). Note that the problem is NP-hard, since it reduces to the vehicle routing problem (VRP), a well-known NP-hard problem, in the single period case. Our problem can be modelled as the integer linear program (1)–(16), with the following variables: PIt , inventory level at the plant on day t; CIit , inventory level of customer i on day t; pt , total production of day t; Xij vt , 0–1 variables equal to 1 if vehicle v visits node i just before node j on day t, 0 otherwise; Zitvu , 0–1 variables

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3405

equal to 1 if demand qit is brought by vehicle v on day ut, 0 otherwise; yt : 0–1 variables equal to 1 if there is a production setup on day t, 0 otherwise.    Cij Xij vt + hPIt + sy t (1) min i∈N j ∈N v∈F t∈T 

s.t.



Xij vt 1

j ∈N



Xj ivt =

j ∈N



t∈T 

t∈T 

∀i ∈ N  , ∀v ∈ F, ∀t ∈ T  ,



(2)

∀i ∈ N, ∀v ∈ F, ∀t ∈ T  ,

(3)

∀v ∈ F, ∀u ∈ T  ,

(4)

Xij vt

j ∈N

Zitvu qit W

i∈N  t  u



Xij vt |S| − 1

∀v ∈ F, ∀t ∈ T  , ∀S ⊆ N  : |S| 2,

(5)

i∈S j ∈S

 

Zitvu = 1 ∀i ∈ N  , ∀t ∈ T  ,

(6)

v∈F 1  u  t



Xij vu Zitvu

j ∈N

CIit = CIi,t−1 +

∀i ∈ N  , ∀t ∈ T  , ∀v ∈ F, ∀u ∈ T  : u t,



Ziuvt qiu − qit

∀i ∈ N  , ∀t ∈ T  ,

(7) (8)

u  t v∈F

PIt = PIt−1 + pt −  

 

Ziuvt qiu

∀t ∈ T  ,

(9)

i∈N  u  t v∈F

Zi,t+1,vu 

v∈F 1  u
0 pt PPCyt 0 CIit CSCi 0 PIt PSC Xij vt ∈ {0, 1} Zitvu ∈ {0, 1} yt ∈ {0, 1} ∀t

 

Zitvu

∀i ∈ N  , ∀t ∈ T  : t < ,

v∈F 1  u  t 

∀t ∈ T , ∀i ∈ N  , ∀t ∈ T , ∀t ∈ T , ∀i ∈ N, ∀j ∈ N : i  = j, ∀v ∈ F, ∀t ∈ T  , ∀i ∈ N  , ∀t ∈ T  , ∀v ∈ F, ∀u ∈ T  : u t, ∈ T .

(10) (11) (12) (13) (14) (15) (16)

The objective function (1) represents the total cost to be minimized. Constraints (2) to (5) come from the VRP. Constraints (2) ensure that each customer is visited at most once per day. Trip connectivity is guaranteed thanks to constraints (3). Constraints (4) forbid exceeding vehicle capacity W: the sum cumulates all demands delivered in day u, which may include future demands. The exponential number of constraints (5) prevents all subtours defined over the set N  (customers only). They can be replaced by a polynomial number of weaker constraints (17), derived from the improved version of Miller–Tucker–Zemlin constraints proposed by Desrochers and Laporte for the travelling salesman problem [21]. These alternative constraints require n integer variables (18). eit − ej t + nX ij vt + (n − 2)Xj ivt n − 1 ∀i ∈ N  , ∀j ∈ N  , i  = j, ∀v ∈ F, ∀t ∈ T  ,

(17)

eit ∈ N 

(18)

∀i ∈ N  , ∀t ∈ T  .

In constraints (6), the demand qit of each customer i for day t must be delivered entirely by one single vehicle in one day u t. If vehicle v does not visit customer i on day u, it cannot bring this customer a demand for day u or beyond, which is expressed by (7). Constraints (8) and (9) are the inventory balance constraints for each customer and for the plant, respectively. For each customer i, the demand qi,t+1 cannot be delivered before the demand qit . To enforce this FIFO rule, constraints (10) must be satisfied for each customer and any two consecutive days t and t + 1. In fact, if demand of customer i for day t + 1 is delivered before its demand for day t, the left sum is greater than the right sum and the

3406

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

constraint is violated. Each day that a production setup occurs, the amount produced must not exceed plant production capacity, constraints (11). Constraints (12) and (13) define the range of plant and customer inventories. Constraints (14)–(16) are integrality constraints for variables Xij vt , Zitvu and yt . Only very small instances can be solved using commercial MIP solvers. The model using subtour elimination constraints (17)–(18) instead of (5) has been tested with LINGO on a Pentium 4 PC clocked at 2.8 GHz. Instances with five customers and five periods are solved in 90 min. In 1 day of computation, it is impossible to go beyond n =  = 6. Concerning the relaxed linear program, we succeeded in solving instances with n = 50,  = 8 and |F | = 5 (100,000 variables Xij vt !), but the resulting deviation of our best metaheuristics to the resulting lower bound is important (61%). This situation is not surprising because even a simplified version like the periodic VRP or PVRP (our problem without the production planning part), still lacks exact methods based on integer linear programs. Moreover, the first non-trivial lower bound for the PVRP was presented orally and recently by Mingozzi at CO’2002 [22] and it is not yet published. Since our goal is to solve large instances (e.g., 200 customers and 20 periods), metaheuristics are tools of choice to compute high-quality solutions in a reasonable time. 3. Preliminary GRASP Greedy randomized adaptive search procedures were introduced in 1989 by Feo and Bard [23], see also [24] for a general presentation of the method. It is based on an iterative randomized sampling of the search space, in which each iteration computes one feasible solution to the problem. Each GRASP iteration consists of two main phases, a construction phase based on a randomized heuristic and a local search phase that improves the trial solution provided by the first phase. After one given number of iterations the best solution is kept as the final result. GRASP was preferred here to other metaheuristics because of its relatively simple structure and its small number of parameters. However, like any metaheuristic, it simply gives a general structure: the components are problem-specific and their design is not trivial in our case. For instances, several attempts were required to find a greedy production planning heuristic that can be randomized without degrading solution quality too much. Moreover, the local search is complex: to be effective, it must be able to change production and delivery days for some demands and this affects both the production plan and vehicle routing. In the next section, a third component (either a reactive mechanism or path relinking) is even added to reinforce this combination. 3.1. Construction phase The proposed heuristic can be specified by Algorithm 1, using the notation introduced for the linear model of Section 2, plus dt that denotes the total amount distributed to customers in day t. It comprises two phases. Phase 1 determines day by day a preliminary production plan and associated trips, but without storage at the plant. Phase 2 tries to shift some production days to achieve the best compromise between setup and storage costs at the plant. Phase 1 considers (line 2) each day t with unsatisfied demands to determine the subset of customers to be visited, the amount to be delivered to each of them (this amount may cover several consumption days), and the associated trips. The total amount shipped is produced in day t − 1. The computations for one given day t comprise three steps. Since shortages are not allowed, the goal of step 1 is to fulfil all unsatisfied demands qit for day t. This is always possible because it is assumed that the total demand of a day never exceeds production and fleet capacities. One empty trip is initialized for each vehicle and unsatisfied demands are gathered in a candidate list CL (lines 5–7). Each iteration of loop 8–14 uses the insertion costs of the customers into the trips to satisfy one demand. The insertion cost of one customer i between two consecutive nodes j, k in a trip is defined as Cj i + Cik − Cj k . The best insertion position and associated cost are computed for each customer i in CL. The  demands with cheapest insertion costs are transferred in the restricted candidate list (RCL), one demand is randomly drawn in this RCL and the corresponding insertion is executed. In the loop 8–14, note that the insertions costs and the RCL must be recomputed: indeed, if one customer i has just been inserted between two nodes j and k in one trip, the trip residual capacity decreases and some moves may become impossible. Moreover, new insertions become possible before and after i. Step 2 also tries to satisfy in day t some demands for future days. Contrary to the demands qit , which can be completely satisfied, only the earliest eligible demand qiu for each customer i is now included in CL (lines 16–19): the two storage capacities (at the plant and at customer i) and the fleet capacity must be respected. The loop at lines 20–28 is similar to the one for day t, except that each demand treated is replaced in CL by the next eligible demand for the same customer.

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3407

In Step 3, the Clarke and Wright savings heuristic [25] is used to rebuild in a better way the trips built by successive insertions for day t. This classical method starts from a trivial solution with one trip per customer. Then, it merges the pairs of trips that bring the best saving, until further mergers would violate vehicle capacity or increase total cost. In our problem, the fleet is limited and we can have more trips than vehicles available. In that case, we continue with the least cost-increasing mergers, if any, because each merger saves one vehicle. However, even this strategy can fail. This is why fleet capacity is artificially reduced by 10% when selecting eligible candidates in lines 17 and 20. This restriction can be removed in the local search. The preliminary production plan is improved in phase 2 by a method inspired by Wagner and Whitin’s algorithm for production planning [26]. The aim is to reduce total cost by finding the best compromise between setups and holding costs, while respecting production and storage capacity constraints at the plant. This method proceeds day by day. The setup in day 1 cannot be eliminated. For each day t with a production setup, the algorithm computes the earliest production day u > t (line 34). It then checks whether the inventory level at the plant and the remaining production capacity on day t allows us to bring forward (and to store at the plant) the production of u to day t (line 35). Such a move is performed if it results in a negative cost variation. The computational complexity of the whole greedy heuristic is O(n2 ). Algorithm 1. Randomized heuristic algorithm for the construction phase 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: 35: 36: 37: 38: 39:

//Phase 1: build a preliminary plan without storage at the plant for each day t with unsatisfied demands do //Step 1 : satisfaction of demands of day t (always possible) initialize pt−1 and dt to 0 and prepare one empty trip per vehicle for each customer i with qit unsatisfied do store the demand as a triple (i, t, qit ) in CL end for while |CL| > 0 do compute the best insertion cost of each demand of CL in the trips of day t RCL := the  cheapest demands of CL in terms of insertion cost remove randomly one demand (i, t, qit ) from RCL perform the best insertion for (i, qit ) in the trips add qit to pt−1 , dt and CIit end while //Step 2: also try to satisfy future demands in day t for each customer i do find the earliest day u with qiu unsatisfied, pt−1 + qiu PPC, CIit − qit + qiu CSCi , dt + qiu 0.9|F |W if u exists then add the demand (i, u, qiu ) to CL end if end for while |CL| > 0 do compute the best insertion cost of each demand of CL in the trips of day t RCL := the  cheapest demands of CL remove randomly one demand (u, i, qiu ) from RCL perform the best insertion for (i, qiu ) in the trips add qiu to pt−1 , dt and CIit find the earliest day v with qiv unsatisfied, pt−1 + qiv PPC, CIit − qit + qiv CSCi , dt + qiv 0.9|F |W if v exists then store the demand as a triple (i, v, qiv ) in CL end if end while //Step 3: improve the trips obtained by step 1 and 2 using the Clarke and Wright algorithm end for //Phase 2: find the best balance between setups and storage at the plant for each day t with pt > 0 do for u :=t + 1 to  with pu > 0 do if (pt + pu PPC) and (PIt + pu PSC) and (hp u (u − t) − s < 0) then add pu to pt and PIt pu := 0 end if end for end for

3408

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3.2. Local search phase The local search proposed here is original and rather involved. Indeed, compared to improvement procedures for classical vehicle routing problems, in which one or two trips in one period can be modified at a time, our moves can change, for each customer in a day, the quantity delivered, the day of production, the day of delivery, the delivery trip and the position in this trip. For instance, suppose one delivery for a customer i is moved to another day t. As a customer cannot be visited more than once, the amount moved must be added to the one already delivered in day t, if any. But this may be impossible if the receiving trip has not enough residual capacity. In such cases, the total amount can be moved to another trip, or a new trip may be created if some vehicles are not yet used. Of course, the production plan must be modified accordingly. This example illustrates the kinds of moves which must be enumerated, tested for feasibility and evaluated in terms of total cost variation. The proposed method is a first-improvement local search applied to each day t of delivery and to each customer i visited that day. Two groups of moves are considered: in the same day t and from day t to another day. In the sequel, s(i) is the successor of client i in its trip. In the first group, all pairs (i, j ) of clients visited are considered, even if they belong to distinct trips, and the following moves are evaluated if vehicle capacity is respected: • 2-OPT moves. If i and j are in the same trip, the sub-sequence i → j is inverted. Otherwise, the shortest paths (i, s(i)) and (j, s(j )) are cut and replaced by either (i, j ) and (s(i), s(j )) or (i, s(j )) and (j, s(i)). • Move i, i.e., remove i and reinsert it after j (j may be the depot for this move). • Swap customers i and j. Concerning the second group, the constraint of visiting each client at most once per day must be respected. Then, we check if the quantity delivered to client i on day t can be moved to another day while respecting all constraints (storage capacities of customer and plant, transportation and production capacities, no shortage). It can be shown that the possible destination days define two intervals of consecutive days before and after day t. For each destination day u, there are two cases when moving to day u the amount delivered to customer i in day t: either i is already visited in day u, say by one trip R, or i is not visited in day u. In the first case, we first try to preserve the trip sequence, i.e. the amount for day t is added to the amount already delivered by R. If this violates the capacity of R, two kinds of moves are considered: • exchange the total quantity of client i with another client j in another trip of day u; • insert this amount into another existing trip or one new trip, if fleet size is not exceeded. There is an exception to the first-improvement rule here: the neighbourhood defined by these two last moves is explored completely to give the best improvement. In the second case (i not already visited in day u), we simply try to insert it into an existing trip or into a new trip in day u, without exceeding the number of available vehicles. Concerning complexity, the instances used in the computational testing (Section 6) use real numbers for the routing costs (Euclidean distances). In such conditions, the number of successive improved solutions constructed by the local search is not polynomial in the worst case. However, each iteration (enumerating possible moves for one given solution) is in O(n2 2 ). In practice, the local search is much faster because many tests are included in the implementation to avoid testing some moves. For instance, the possible insertions of customer i from day t into the trips of another day u are not tested if the residual capacity of the fleet in day u is not sufficient.

4. Reactive GRASP The GRASP of Section 3 can be made reactive. In a basic GRASP, the size of the RCL  is fixed for all instances and all iterations. However, it is difficult to find the RCL size that gives the best average results. With  = 1, the

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3409

GRASP becomes a deterministic heuristic and all iterations yield the same solution. If  is too large, well-diversified but low-quality solutions are obtained. The principle of reactive GRASP, proposed for the first time by Prais and Ribeiro [27], is to let the algorithm find the best value of  in a small set of allowed values. Let v be the fixed number of allowed values and L = {1 , 2 , . . . v } the list of these values. Both must be determined during some preliminary testing. Let k be the kth value of L, k = 1, 2, . . . , v. The reactive GRASP is very similar to the basic version, except that  is randomly selected in L before building one trial solution and that the probabilities of choosing the different k are progressively modified to favour the RCL sizes that give the best average solution costs. The probabilities Pk of taking  = k are initialized to 1/v (uniform distribution). During the algorithm, the number of solutions built with each value k is recorded in count(k) and the corresponding solution costs (after local search) cumulated in score(k). Let Zmin be the provisional best solution cost. For each block-iteration (sequence of  successive iterations, e.g.,  = 50), the Pk are updated to favour the values yielding the best average costs, using score, count, and one given amplification factor  (e.g.,  = 10): • for each value k of L, compute avg(k) = score(k)/count(k) and then Qk = (Zmin /avg(k)) ; • compute the sum  of the Qk ; • replace each Pk by Qk / to get values in [0, 1]. After some block-iterations, the reactive GRASP tends to use the value of  that has given the best solution costs on average. The whole process is summarized by Algorithm 2, in which ni is the number of GRASP iterations, i an iteration counter, Z(S) the cost of solution S and Smin the best solution. The reactive mechanism adds only a small overhead to the basic GRASP. In particular, the complexities already stated for the randomized heuristic and the local search are preserved. Algorithm 2. General structure of reactive GRASP 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19:

Zmin := +∞ for k := 1 to v do initialize count(k) and score(k) to 0 and Pk to 1/v end for for i :=1 to ni do draw k ∈ [1, 2, . . . , v] with probability Pk compute one trial solution S using Algorithm 1 with  = k improve S using the local search of 3.2 if Z(S) < Zmin then Smin :=S; Zmin :=Z(S) end if count(k) :=count(k) + 1; score(k) := score(k) + Z(S) if iter mod  = 0 then avg(k) = score(k)/count(k) Qk = (Zmin /avg(k)) compute , the sum of all Qk , and replace Pk by Qk / for each k end if end for return Smin

5. Path-relinking Path-relinking (PR) is another possible technique to improve a basic GRASP. It consists of exploring paths between elite solutions found by the GRASP and kept in a small pool. Each path links one initiating solution to one guiding solution. For instance, in Fig. 2, the path from A to B (solid line) does not yield solutions improving A or B, while the path from B to A (dashed line) produces three better solutions.

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

Objective function value

3410

A B

Moves Fig. 2. Principle of path relinking.

Ideally, the elite solutions must be well scattered in the solution space and this can be achieved using some kind of distance measure between solutions. In practice, it may be hard to find a measure satisfying all properties of a distance, especially the triangle inequality. Hence, a dissimilarity measure D(A, B) that grows with the number of differences between two solutions A and B may be used instead. The following subsections present such a measure for our problem, the construction of the pool, the path generation, the general structure of the PR phase and its possible locations in the GRASP. 5.1. Measure definition The proposed measure takes the production and distribution levels into account and uses different weights for the possible kinds of differences between two solutions A and B. The measure is initialized to 0. For the production level, the amounts produced in each period t are compared and • if t is a production day in A or B but not both, 3 is added to D(A, B); • if t is a production day in both A and B, the amounts produced are compared: a weight of 1, 2 or 3 is added if the difference with respect to A is less than 5%, between 5 and 10% or greater than 10%, respectively. For the distribution level, the trips in each period t of a solution are viewed as one single string, using symbol 0 (the depot) as trip delimiter. For instance, if one solution A contains two trips (1, 4, 3) and (5, 2) in day t, the resulting string is At = (0, 1, 4, 3, 0, 5, 2, 0). Here our measure extends a distance between permutations proposed by Campos et al. [28], that counts the number of pairs of consecutive elements on one permutation that are split in the other. In our case, there are several periods and even the strings At and Bt for a given day t in A and B are not permutations if some customers are not visited. We derive a suitable measure by considering each day t in turn: • if customer i is visited in At or Bt but not both, add 2 to D(A, B); • if i is visited in At and Bt but with a different predecessor or successor, add 1. By construction, a pair (i, j ) of consecutive customers in A is not counted as split if (j, i) is found in B, contrary to the original distance proposed by Campos et al. Thanks to this property, D(A, B) = 0 only if the two solutions are identical, even if some trips are inverted: in fact, in an undirected graph, a trip is equivalent to its mirror (inverted trip). Consider the following simple example of distance computation, with five customers and two periods. Solution A p1 = 28 and A1 = (0, 4, 1, 3, 0, 2, 5, 0) p2 = 11 and A2 = (0, 3, 2, 0, )

Solution B p1 = 30 and B1 = (0, 1, 2, 3, 0, 4, 5, 0) p2 = 09 and B2 = (0, 2, 5, 0)

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3411

D(A, B) is initialized to 0. The contribution of day 1 is 7: all five customers in A1 are also in B1 , but with a different successor or predecessor (contribution 5), and the difference between the amounts produced is 7% (contribution 2). The contribution for day 2 is 8: 2 × 2 = 4 for customers 3 and 5 that are not visited in both solutions, 1 for customer 2 found in both solutions but with different neighbours, and 3 for the production discrepancy of 18%. Hence, D(A, B) = 15. 5.2. Construction of the pool of elite solutions The path-relinking is applied to a small pool  of elite solutions. The pool is empty at the beginning. In a first phase, the GRASP fills the pool: each new best solution is added to , until elite solutions are available. The resulting pool is sorted in increasing order of costs. In a second phase, during the remaining GRASP iterations, the pool is updated. We have compared two techniques. The simplest one is to replace the worst solution in  (the last one) each time a better solution is found. A more sophisticated technique was suggested by Resende and Ribeiro [29] to increase pool diversity using a dissimilarity measure (in our case the measure D described before). Each time a new solution A is better than the worst in , replace the most similar solution (i.e., with the smallest measure) among all solutions with cost worse than A. The second technic was selected because of its better results in preliminary testing. 5.3. Path generation Given one initiating solution A and one guiding solution B, the idea is to modify progressively A to reproduce in each day t the trips present in B. B is scanned day by day, trip by trip in each day and customer per customer in each trip. Let t be the selected day, i the selected customer, ditS the amount brought to i in day t of a solution S, = ditA − ditB , and u the earliest delivery day after t in B (i.e., the days without trips are skipped). At that time, the strings coding the trips are identical in both solutions for days 1, 2, . . . , t − 1 and, in day t, the two strings in A and B are identical until the predecessor of i in B. The first step is to reproduce in A the amount which is delivered to i in B. If = 0, there is nothing to do in that A > 0), the surplus is aggregated with the step. If > 0, this surplus is delayed to day u in A. If i is visited in day u (diu A amount already planned. When diu = 0, if the predecessor  or the successor of i in day u of B i is also visited in A, i is inserted after  or before (the idea is to relink a pair of consecutive clients in B that is split in A). If  and are not visited in day u of A, a best insertion of i is performed in the trips of day u. Due to vehicle capacity, it may be impossible to ship the whole surplus in day u. In that case, the surplus is split in order to bring an integer number of consecutive consumption days to i in day u , to respect our FIFO rule. The process is repeated with the next delivery day u, until the residual surplus becomes zero. If < 0, the shortage is corrected by bringing forward some future deliveries to day t. Let v be the earliest delivery A ) is added to d A and is updated ( ← − R). Of course, i is day to i after day t in A. The amount R = min( , div it A also removed from its trip in day v if R = div . If does not equal zero, the process is repeated with the next delivery day, until the lack is satisfied. The second step consists of correctly placing customer i in day t. The customer is simply moved in order to have the same predecessor in A and B. If the target trip is overloaded, it is split just after customer i. Of course, in parallel, the production plan and the inventory levels (customers and plant) must be updated accordingly. Infeasible solutions (mainly because of the limited fleet capacity) are so rare that they can be discarded instead of being repaired or penalized (this is confirmed in the next section on computational experiments). In general, a path relinking process finds improved solutions only if intermediate solutions generated along a path are improved by local search. In our case, to limit running time, a new feasible solution undergoes the local search only when a trip of B is completely reproduced in A. This option is also justified by the fact that consecutive solutions on a path often lead to the same solution after local search, this is why systematic calls to the local search are useless. To illustrate the method, consider the following example with the same solutions as in Section 5.1, completed by the amounts delivered to each customer. Solution A can be transformed into solution B in five steps. Bold numbers indicate which customer is moved in each step. In the third step, customer 3 is removed from day 2 because its demands for the two days are now delivered in day 1. Customer 5 appears in the fifth step because the surplus of day 1 is inserted in day 2. The other steps move one customer without changing the amount received.

3412

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

Customers Day 1, demands: Day 2, demands:

1 2 5

Instance data 2 4 8

A1 : Amounts: A2 : Amounts:

0 0 0 0

Solution A 4 1 3 7 7 1 3 2 0 3 8 0

A1 : Amounts: A2 : Amounts:

0 0 0 0

1 7 3 3

A1 : Amounts: A2 : Amounts:

0 0 0 0

1 7 2 8

3rd step 2 3 4 4 0 0

0 0 0 0

1 7 2 8

A1 = B1 : Amounts: A2 = B2 Amounts:

1st step 4 3 7 1 2 0 8 0 4 7

3 1 3

4 0 7

5 8 1

0 0

2 4

5 9

0 0

B1 : Amounts: B2 : Amounts:

0 0 0 0

Solution B 1 2 3 7 4 4 2 5 0 8 1 0

0 0

2 4

5 9

0 0

A1 : Amounts: A2 : Amounts:

0 0 0 0

1 7 3 3

0 0

5 9

0 0

5th and last step 2 3 4 4 5 0 1 0

A1 : Amounts: A2 : Amounts: 0 0

0 0 0 0 4 7

1 7 2 8

2nd step 2 4 4 7 2 0 8 0 4th step 2 3 0 4 4 0 0 0 5 8

0 0

4 7

5 8

0 0

3 1

0 0

5 9

0 0

0 0

5 9

0 0

4 7

0 0

5.4. General structure of path-relinking and location One possible first strategy to combine PR and GRASP is to perform the PR after the GRASP, as a post-optimization process. We call this strategy S1. For any pair {A, B} of solutions in the elite pool , two paths A → B and B → A are generated. Let S be the best solution found after having evaluated all pairs. While S improves the best solution in , the following steps are iterated: (a) generate two paths A → S and S → A for each solution A ∈ , (b) record the best solution T obtained, (c) add S to  and (d) assign T to S (S ← T ). Another strategy called S2 consists of interleaving GRASP iterations and PR steps. More precisely, a PR step is applied each time the GRASP discovers a new best solution A, provided the pool  is full. A is added to the pool as usual. If its cost does not exceed the average cost of the pool, two paths A → B and B → A are generated for any pool solution B  = A. Applying the PR when the cost of solution A exceeds the average cost increases the running time without bringing significant improvement of the best solution. It is also useless to execute the PR when the pool is not yet full, because preliminary tests have shown that the average quality of elite solutions may be poor at that stage.

6. Computational evaluation 6.1. Implementation and instances All algorithms designed in this paper were implemented in Delphi 7 [30] and executed on a 2.8 GHz PC under Windows XP. To the best of our knowledge, no instances are publicly available for our problem. This is why the instances randomly generated in our previous work [31] on simpler heuristics are reused here. There are three sets of 30 instances with 50, 100 and 200 customers respectively, all with T = 20 periods.

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3413

Some data in the instances depends on the number of customers: production capacity, customers and plant storage capacities, number of vehicles and setup cost. On average, production capacity corresponds to 3–4 days of customer consumption and the fleet capacity represents 2 days of consumption for the full set of customers. The production setup cost is proportional to production capacity. The storage capacity at the plant varies between one and a half and twice the production capacity. The holding cost at the plant is the same for all instances. One can note that VRP (single-period) instances with 200 customers are already considered as respectable, and adding a 20-period horizon makes the problem much more combinatorial. Note also that our instances are much bigger than the ones used by Chandra and Fisher [11]. 6.2. Algorithms compared and parameters Four metaheuristics are compared with our two earlier heuristics H1 and H2 presented in [31]: the basic GRASP, the reactive version and the two versions with the path relinking (either at the end or interleaved). H1 and H2 are not simple greedy heuristics. Briefly speaking, H1 is a classical two-phase approach: a production plan is elaborated using Wagner and Whitin’s method, then this plan is frozen and the trips are built for each day. H2 implements a weak integration between production and distribution: it corresponds to H1, but with one final feedback that adjusts the production plan after trip construction. In both heuristics, the trips are optimized day by day using a local search that combines 2-opt moves and transfers of customers. Contrary to the metaheuristics, demands cannot be moved to different days because the production plan is definitive. The ideal number of iterations ni depends on the GRASP version. The basic GRASP and the two versions with path relinking work well with ni = 100, while the reactive GRASP requires ni = 500 to become effective. Indeed, recall (Section 4) that this version needs to collect enough statistics about solution costs before finding ad hoc probabilities for each allowed RCL size. Hence, the results in the next subsection are given for ni = 100 and ni = 500. For the preliminary GRASP, the RCL size  was set to 14, 15 and 12 for 50, 100 and 200 customers, respectively. These values give the best average cost for each set of instances. Concerning the reactive GRASP, the additional parameters are the v possible RCL sizes k , the amplification factor  and the block size . The following values were selected after some preliminary experiments. For all instances, v = 5,  = 10 and  = 50. However the values of the v allowed RCL sizes that give the best average results depend on the set of instances and were selected as follows: • n = 50 :  = (4, 9, 12, 14, 16); • n = 100 :  = (15, 16, 17, 19, 20); • n = 200 :  = (4, 7, 8, 11, 12). Concerning the GRASP versions with path relinking, the RCL size  and the size of the elite pool have been set, respectively, to 14 and 8 for n = 50, 15 and 10 for n = 100 and 12 and 7 for n = 200. The percentage of infeasible solutions generated during the PR is negligible, except for n = 200 where it reaches 8%. However, this percentage is always small enough to avoid repair or penalizing procedures: such infeasible solutions are discarded. Due to the disappointing results (mentioned in Section 2) with the integer linear programming model, no lower bound is included in the comparisons. Indeed, finding a tight lower bound seems a very challenging task, even for the PVRP that corresponds to our problem with production constraints removed. For instance, it is difficult to generalize capacity-related cuts because the optimal amount delivered to each customer in each period is not known. Recall also that the first efficient lower bound for the PVRP is quite recent, very involved, and not yet published. 6.3. Results The results are summarized in Tables 1–3 for the instances with 50, 100 and 200 customers, respectively. These tables share the same format. For each instance, the solution value (Cost columns) and the running time in seconds (Time columns) are provided for six algorithms: our two earlier heuristics (columns H1 and H2), the basic GRASP (GRASP column), the reactive version (RGRASP), the GRASP with path relinking at the end or strategy 1 (GRASP+PR/S1) and the one with the PR interleaved or strategy S2 (GRASP+PR/S2). All GRASP versions are executed with ni = 500

0.10

0.16 0.14 0.09 0.09 0.11 0.10 0.10 0.09 0.10 0.10 0.10 0.11 0.09 0.08 0.09 0.10 0.10 0.12 0.09 0.10 0.10 0.10 0.11 0.09 0.09 0.09 0.09 0.09 0.10 0.09

985 355 885 565 855 200 150 835 805 245 970 435 640 030 660 865 020 425 945 480 040 495 615 830 080 155 660 920 170 520

464 294.5 0 9.10%

457 472 464 457 458 465 467 475 463 462 460 461 466 464 461 466 468 468 470 458 457 455 467 459 471 471 458 459 471 463

The best solution value on each line is emphasized in boldface.

511 578.9 0 0

Average cost or time Best solutions found Cost saving from H1

835 636 862 366 549 784 413 887 002 695 978 599 836 893 994 658 234 669 135 517 970 653 725 746 122 779 330 786 026 688

487 517 510 519 516 554 498 506 497 534 510 539 495 494 490 496 531 572 502 511 487 475 524 503 531 500 489 488 532 522

Cost

Cost

Time

Heuristic H2

Heuristic H1

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

Instance

Table 1 Results for instances with 50 customers

0.17

0.14 0.11 0.11 0.17 0.18 0.19 0.18 0.18 0.18 0.18 0.18 0.17 0.18 0.18 0.17 0.19 0.18 0.18 0.18 0.18 0.18 0.11 0.18 0.17 0.18 0.18 0.18 0.19 0.18 0.18

Time 630 005 618 518 686 639 677 583 849 163 055 140 891 018 100 911 878 015 571 798 112 296 882 242 605 321 984 254 111 901

445 848.4 2 12.74%

440 465 455 456 434 452 436 422 434 454 434 453 440 433 454 457 455 442 461 463 410 429 455 430 462 444 444 450 464 433

Cost

GRASP

86.97

51.16 88.42 98.07 56.13 65.65 90.73 110.62 77.34 133.81 121.12 72.79 74.49 97.08 79.90 90.08 90.51 91.08 88.53 99.39 91.25 66.23 83.17 95.90 101.40 76.65 75.03 82.28 104.55 67.52 88.39

Time 565 490 468 398 631 564 872 263 834 488 115 705 771 973 875 385 978 180 896 858 677 281 252 972 245 209 695 194 555 531 443 264 10 13.25%

440 449 455 456 434 452 436 422 434 453 434 452 440 432 453 457 455 442 418 452 409 429 455 429 462 442 444 450 461 434

Cost

RGRASP

93.45

63.93 98.07 104.87 68.42 74.55 98.07 112.74 89.96 123.80 121.25 82.11 85.83 99.85 82.90 86.67 100.39 98.18 94.93 104.27 97.93 78.27 88.76 102.16 81.42 86.99 88.35 85.69 108.92 93.93 100.51

Time 505 005 618 518 019 639 677 153 849 506 255 050 741 133 912 490 738 015 192 963 847 296 275 042 605 194 935 254 320 901 441 735.3 12 13.41%

440 452 455 456 434 452 436 420 434 436 433 453 440 416 453 457 455 442 432 452 409 429 452 429 462 442 444 450 462 433

Cost

Time

87.90

67.88 90.45 95.35 73.74 78.91 89.75 98.80 87.32 107.82 105.00 81.77 88.27 93.99 84.45 81.16 91.70 89.40 86.64 98.03 94.12 72.46 83.25 91.34 78.35 82.20 82.15 81.17 98.19 89.91 93.73

GRASP + PR/S1

580 880 985 386 686 639 677 856 841 076 680 140 771 578 746 610 843 985 481 085 082 116 812 207 440 554 835 113 710 901 442 076.5 10 13.49%

440 451 453 456 434 452 436 421 434 447 433 453 440 415 452 457 455 441 419 445 410 429 455 429 462 442 444 445 462 433

Cost

Time

131.49

102.36 138.17 135.31 77.58 99.37 158.50 150.09 133.45 142.54 158.40 117.42 94.54 163.14 114.13 156.29 80.18 159.98 122.68 257.19 130.55 157.48 108.26 106.99 106.30 58.93 130.33 128.09 187.46 136.77 132.27

GRASP + PR/S2

3414 M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

888 819 994 752 945 350 940 303 867 192 584 656 435 707 347 510 473 923 653 115 144 325 747 608 907 225 204 618 083 149

963 648.8 0 0

954 970 956 954 969 947 941 950 960 963 960 967 961 957 958 961 979 968 1 006 974 959 980 946 974 958 966 959 967 970 959

0.54

0.47 0.52 0.55 0.51 0.55 0.51 0.53 0.53 0.55 0.53 0.53 0.55 0.52 0.53 0.55 0.53 0.53 0.53 0.61 0.53 0.53 0.52 0.62 0.61 0.53 0.55 0.53 0.53 0.55 0.55

445 915 560 000 445 970 520 925 560 170 020 820 560 805 890 875 210 805 785 945 355 650 655 690 690 100 870 230 255 895

841 920.5 0 12.63%

842 843 841 831 848 830 827 839 841 838 841 845 835 839 841 844 852 845 850 849 842 853 833 849 843 845 834 838 846 836

Cost

Cost

Time

Heuristic H2

Heuristic H1

The best solution value on each line is emphasized in boldface.

Average cost or time Best solutions found Cost saving from H1

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

Instance

Table 2 Results for instances with 100 customers

0.67

0.76 0.67 0.65 0.65 0.70 0.64 0.69 0.67 0.66 0.62 0.65 0.64 0.64 0.64 0.67 0.67 0.66 0.66 0.62 0.64 0.64 0.70 0.65 0.68 0.63 0.68 0.69 0.66 0.72 0.70

Time 183 228 768 130 622 960 909 804 180 321 507 520 341 644 775 755 484 922 291 026 644 817 617 762 871 232 542 572 090 952

792 515.6 1 17.75%

796 785 790 782 799 794 781 790 763 792 790 798 778 792 795 799 800 788 807 811 794 806 784 798 796 796 784 782 804 785

Cost

GRASP

415.90

375.09 542.21 395.20 524.45 439.06 429.60 427.33 487.45 429.37 366.27 400.86 316.64 513.70 467.77 454.45 393.74 510.40 381.72 298.10 347.59 359.96 390.68 308.83 431.23 364.84 384.81 464.78 416.89 457.51 396.56

Time 819 038 972 212 181 509 804 884 646 903 961 804 243 802 224 656 159 078 964 867 136 707 508 428 591 611 526 295 725 906 791 805.3 5 17.83%

791 785 787 783 798 793 781 780 784 791 788 790 781 791 795 799 803 788 806 810 790 805 783 798 796 795 778 782 801 785

Cost

RGRASP

415.60

411.47 505.48 391.52 497.42 444.28 423.86 419.23 499.19 442.30 357.20 391.87 347.53 502.98 449.61 432.47 395.84 508.69 400.23 323.80 356.13 362.92 380.39 339.16 416.88 368.14 374.86 461.02 406.33 443.00 398.11

Time 509 048 342 130 190 209 781 804 222 329 964 605 341 900 866 277 301 909 291 517 714 817 736 404 614 449 168 095 608 952 790 236.4 12 17.99%

792 783 786 782 796 794 779 790 758 788 787 796 778 787 793 796 797 785 807 808 792 806 781 795 796 793 781 780 801 785

Cost

436.42

404.30 521.78 407.98 521.90 463.61 456.63 423.53 497.15 455.31 378.54 420.40 369.25 505.97 467.69 451.29 418.19 520.11 419.09 353.69 387.27 392.96 412.97 366.12 444.66 396.02 403.69 495.51 435.10 480.51 421.42

Time

GRASP + PR/S1

609 228 505 122 044 131 845 004 585 214 731 805 341 344 928 746 484 734 291 917 183 133 059 402 137 892 368 657 908 945 789 909.8 14 18.02%

791 785 790 780 799 794 778 789 756 792 789 788 778 792 791 793 800 781 807 808 788 805 782 798 795 793 784 777 795 785

Cost

466.22

413.42 506.19 346.76 486.51 470.19 349.02 298.64 507.42 569.16 391.05 393.90 548.23 557.68 497.19 527.02 398.70 565.42 414.26 513.82 403.09 477.35 450.73 429.12 471.84 428.43 495.64 464.18 599.32 559.06 453.30

Time

GRASP + PR/S2

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419 3415

2.10

2.03 2.06 2.06 2.00 2.03 2.16 2.00 2.36 1.99 2.03 2.14 2.23 2.26 2.33 2.08 2.06 2.00 2.24 2.09 2.03 2.12 2.05 2.05 2.19 2.02 2.02 2.11 2.13 2.17 2.13

115 221 108 100 099 232 242 089 095 099 210 108 113 228 098 102 248 113 103 103 209 100 109 099 106 097 115 101 098 217

430 903 770 835 275 910 915 225 030 065 258 545 135 161 825 035 285 210 355 925 610 160 355 725 085 325 040 075 195 177

1 136 294.6 0 13.40%

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

The best solution value on each line is emphasized in boldface.

1 312 612.3 0 0

Average cost or time Best solutions found Cost saving from H1

885 448 514 302 850 811 975 367 548 798 551 198 792 506 868 485 606 248 273 225 456 313 233 522 588 163 022 872 875 075

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

321 309 311 285 280 302 313 371 286 285 310 319 391 390 304 285 302 387 308 284 308 284 321 301 289 307 318 294 285 313

Cost

Cost

Time

Heuristic H2

Heuristic H1

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

Instance

Table 3 Results for instances with 200 customers

3.42

3.30 3.51 3.21 3.26 3.22 4.11 3.48 3.38 3.33 3.22 3.57 3.33 3.26 3.60 3.27 3.30 3.48 3.22 3.21 3.24 3.60 3.29 3.22 3.35 3.22 4.21 3.33 3.32 3.43 4.06

Time 090 085 095 083 080 077 107 064 079 086 092 076 090 090 080 084 084 109 078 074 079 083 083 081 090 066 111 059 103 080

955 245 374 718 953 213 547 031 408 390 441 313 483 927 625 519 039 761 144 537 881 007 438 914 172 180 924 007 486 461

1 085 069.7 1 17.29%

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Cost

GRASP

1801.80

1802.70 2002.73 1419.91 1908.09 1326.52 2101.69 1631.94 2237.27 1549.83 1804.78 1850.03 1857.58 2087.02 2103.95 1751.50 1609.41 1746.55 1684.84 2143.42 1764.66 2422.88 1681.44 1761.78 1592.80 1739.25 1685.31 1517.36 1508.17 1522.44 2238.44

Time 075 070 073 068 060 066 091 075 055 081 069 057 077 085 066 077 067 095 064 049 055 078 076 063 054 057 076 059 088 056

528 764 938 959 319 636 538 540 447 622 896 631 672 983 161 891 007 350 431 854 436 908 888 772 230 719 798 019 853 982 1 070 025.7 13 18.45%

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Cost

RGRASP

1893.75

2078.39 2282.84 1710.44 2194.75 1365.09 2128.75 1796.72 2245.25 1554.42 1935.06 2151.26 1909.58 2063.05 2068.09 1924.59 1658.39 1826.45 1716.27 2087.14 1910.66 2506.95 1649.62 1803.39 1655.44 1902.87 1782.70 1484.00 1574.81 1463.53 2382.09

Time 079 069 073 068 066 067 074 060 063 070 069 074 090 074 066 067 073 075 071 070 070 065 073 068 072 065 076 059 064 072

895 845 685 930 335 265 450 935 545 010 050 045 483 870 620 685 220 380 510 040 325 825 355 045 965 735 795 007 280 755 1 070 562.8 10 18.40%

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Cost

GRASP + PR/S1

1846.69

1951.86 1785.75 1688.52 1937.23 1678.19 1948.55 1781.17 1893.25 1896.53 1887.64 2075.83 1904.72 1779.33 2052.12 1776.81 2022.39 2006.67 1738.83 1920.98 1816.48 1947.78 1684.47 1878.69 1694.95 1950.41 1683.75 1743.19 1864.83 1795.87 1613.92

Time 079 069 073 071 066 067 107 063 063 069 068 073 090 075 066 067 074 076 064 069 055 064 072 068 074 066 101 059 064 072

820 740 805 750 995 415 547 275 020 005 675 370 483 230 965 850 855 880 331 845 365 400 935 435 375 140 604 007 745 395 1 072 008.5 8 18.29%

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Cost

GRASP + PR/S2

2116.86

2364.11 2106.27 1978.34 2199.55 1924.00 2174.81 2056.58 2170.50 2302.20 2090.38 2464.75 2090.33 2065.41 2264.77 2046.14 2341.55 2157.81 1905.55 2206.37 2216.72 2139.45 1781.39 2201.73 1999.14 2174.92 1998.20 2019.44 2279.80 1921.73 1863.91

Time

3416 M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3417

Table 4 RGRASP (ni = 500) vs other versions with ni = 100 RGRASP

GRASP

PR/S1

PR/S2

446 267.4 13.60 1 12.60%

444 852.2 18.94 7 12.94%

444 941.8 19.80 3 12.91%

793 170.2 112.44 2 17.60%

792 188 155.79 6 17.78%

792 495 160.00 7 17.74%

1 087 107.9 312.78 0 16.92%

1 073 602.7 395.50 4 17.49%

1 073 994.3 587.85 5 17.65%

50 customers Avg. cost Avg. time(s) Nb of best solns Avg. saving from H1

443 264 93.45 19 13.25%

Avg. cost Avg. time(s) Nb of best solns Avg. saving from H1

791 838.6 415.60 15 17.83%

100 customers

200 customers Avg. cost Avg. time(s) Nb of best solns Avg. saving from H1

1 070 025.7 1893.75 21 18.45%

iterations. The three last lines indicate, respectively, for each algorithm the average cost and running time, the number of best solutions found, and the average percentage saved compared to the only uncoupled approach (H1) taken as reference level. Table 1 (50 customers) shows that the only two-phase heuristic (H1) is the fastest but provides the worst results. The other columns demonstrate the positive impact of a coordination between production planning and distribution. Indeed, even the simplest integrated method (H2) saves 9%. The four GRASP do better (around 13%) but at the expense of a larger computational effort: contrary to H1 and H2, they build 500 solutions instead of one and use a more involved local search. The best results are obtained by the GRASP with interleaved path relinking (saving 13.49%), followed by the GRASP with final PR (13.41%) and the reactive GRASP (13.25%). Compared to the basic GRASP (12.74%), the more sophisticated versions save 0.51–0.75% more. The gains are more visible if one considers the number of best solutions found by each method: 10–12 for the improved GRASPs and only 2 for the basic GRASP, while no best solution is obtained by the two simpler heuristics H1 and H2. After a careful analysis, we can explain the longer running times of the three improved GRASP versions as follows. The statistics in the reactive GRASP require a negligible amount of time but, as long as the algorithm has not yet selected the best RCL size, the greedy randomized heuristic produces poor solutions, for which the local search needs more moves before reaching a local optimum. The two other versions take longer because of the paths explored in the path relinking process. The longest time is achieved when the PR is interleaved with the GRASP, because more paths are browsed. Similar conclusions can be drawn from Table 2 (100 customers). The savings brought by coordinated approaches are even larger (13–18%). The four GRASPs are still much more effective than the only other integrated method H2 but reaction or path relinking improve only slightly the basic GRASP (additional saving 0.08–0.27%). This time, the two PR versions outperform the reactive GRASP, especially in terms of best solutions found (12 and 14 versus 5). In Table 3 (200 customers), the deviations to H1 continue increasing. The reactive GRASP is the best algorithm, it obtains 13 best solutions and, respectively, saves 5.05%, 1.16%, 0.05% and 0.16% compared to H2, the basic GRASP and the path relinking variants with strategy 1 and 2. Finally, on all sets of instances, the algorithms that integrate production and routing decisions provide much better results than H1. Among them, all GRASP versions clearly outperform H2, the weakest integration. The savings brought by integration increase with instance size: on the smallest instances, the local search finds less feasible and improving moves. The basic GRASP is outperformed by the improved variants, although the best one (reactive GRASP or GRASP with PR) depends on the set of instances. Concerning running times, they remain reasonable for a tactical problem that does not need to be solved every day: 2 min for n = 50, 7 min for n = 100 and 30–35 min for n = 200, the PR version with strategy 2 being in all cases the lengthiest method.

3418

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

A natural question is to wonder whether combining PR and reaction could improve the results further. We tried this but no significant improvement was obtained, as if each technique reduces the margin for improvement left to the other. Ribeiro [32] told us that such an antagonism between reaction and PR has been observed in several other optimization problems. As already mentioned, all GRASP variants were tested with 500 iterations to let the RGRASP collect enough statistics, but this results in long running times. Table 4 compares the RGRASP with ni = 500 (data copied from previous tables) with the other versions with ni = 100: the two PR versions are still almost as good as the RGRASP but they become 2.5–4.5 times faster. 7. Conclusion and future directions In this paper, different versions of a GRASP improved by a reactive mechanism or by a Path Relinking process have been proposed for solving a combined production–distribution optimization problem. The originality of the approach is to tackle simultaneously production and routing decisions instead of resorting to classical two-phase approaches which are still widely used in practice. The two principal difficulties in the GRASP design were to randomize the construction of a trial solution without altering solution quality too much and to develop a local search able to modify both the production plan and the sets of trips on each period. The results underline the potential benefits that can be achieved by a closer integration between production and distribution in logistic processes. Possible research directions include the study of more complex cases (several products and additional constraints), the design of a non-trivial lower bound and the development of population-based metaheuristics (such as memetic algorithms and scatter search) to try to improve the solutions found by the GRASPs. Moreover, it should be possible to improve the current running times by optimizing the implementations. References [1] Toth P, Vigo D. The vehicle routing problem. Philadelphia: SIAM; 2002. [2] Gendreau M, Laporte G, Séguin R. Stochastic vehicle routing. European Journal of Operational Research 1996;88:3–12. [3] Christiansen M, Nygreen B. A method for solving ship routing problems with inventory constraints. Annals of Operations Research 1998;81: 357–78. [4] Dror M, Ball M. Inventory-routing: reduction from an annual to a short-period problem. Naval Research Logistics 1987;34:891–905. [5] Golden BR, Wasil EA. Computerized vehicle routing in the soft drink industry. Operations Research 1987;35(1):6–17. [6] Graves SC, Rinnooy Kan AHG, Zipkin PH. Handbooks in operations research and management science, vol. 4: logistics of production and inventory. North-Holland: Elsevier; 1993. [7] Nahmias S. Production and operations analysis. 3rd ed., Homewood, IL: Irwin; 1997. [8] Brewer AM, Button KJ, Hensher DA. Handbook of logistics and supply chain management. Oxford: Pergamon; 2001. [9] Tayur S, Ganeshan R, Magazine M. Quantitative models for supply chain management. Massachusetts: Kluwer; 1999. [10] De Kok AG, Graves SC. Handbooks in operation research and management science , vol. 11: supply chain management: design, coordination and operation. North-Holland: Elsevier; 2003. [11] Chandra P, Fisher ML. Coordination of production and distribution planning. European Journal of Operational Research 1994;72(3):503–17. [12] Fumero F, Vercellis C. Synchronized development of production, inventory, and distribution schedules. Transportation Science 1999;33(3): 330–40. [13] Selçuk Erengüc S, Simpson NC, Vakharia AJ. Integrated production/distribution planning in supply chains. An invited review. European Journal of Operational Research 1999;115:219–36. [14] Blumenfeld DE, Burns LD, Daganzo CF. Synchronizing production and transportation schedules. Transportation Research 1991;25B(1): 23–37. [15] Jayaraman V, Pirkul H. Planning and coordination of production and distribution facilities for multiple commodities. European Journal of Operational Research 2001;133:394–408. [16] Bard JF, Huang L, Jaillet P, Dror M. A decomposition approach to the inventory routing problem with satellite facilities. Transportation Science 1998;32(2):189–203. [17] Chandra P. A dynamic distribution model with warehouse and customer replenishment requirements. Journal of the Operational Research Society 1993;44(7):681–92. [18] Ertogral K, Wu SD, Burke LI: Coordination production and transportation scheduling in the supply chain. Technical Report 1998, Lehigh University. [19] Metters RD. Interdependent transportation and production activity at the United States postal service. Journal of the Operational Research Society 1996;47:27–37. [20] Bramel J, Goyal S, Zipkin P. Coordination of production–distribution networks with unbalanced leadtimes. Operations Research 2000;48(4): 570–7.

M. Boudia et al. / Computers & Operations Research 34 (2007) 3402 – 3419

3419

[21] Desrochers M, Laporte G. Improvements and extensions to the Miller–Tucker–Zemlin subtour elimination constraints. Operations Research Letters 1991;10:27–36. [22] Mingozzi A, Valetta A. A lower bound for the periodic vehicle routing problem. Oral communication. International symposium on combinatorial optimisation (CO 2002), Paris, France. [23] Feo TA, Bard j. Flight scheduling and maintenance base planning. Management Science 1989;35(12):1415–1432. [24] Feo TA, Resende MGC. Greedy randomized adaptive search procedures. Journal of Global Optimization 1995;6:109–33. [25] Clarke G, Wright JW. Scheduling of vehicles from a central depot to a number of delivery points. Operations Research 1964;12:568–81. [26] Buffa ES, Sarin RK. Modern production–operations management. 8th ed., New York: Wiley; 1987. [27] Prais M, Ribeiro C. Reactive GRASP: an application to a matrix decomposition problem in TDMA traffic assignment. INFORMS Journal on Computing 2000;12(3):164–76. [28] Campos V, Martí R, Laguna M. Context-independent scatter search and tabu search for permutation problems. INFORMS Journal on Computing 2005;17:111–22. [29] Resende MGC, Ribeiro CC. GRASP with path-relinking: recent advances and applications. In: Ibaraki T, Nonobe K, Yagiura M, editors. Metaheuristics: progress as real problem solvers. Berlin: Springer; 2005. p. 29–63. [30] Borland Software Corporation. Delphi 7, United States of America, 2002. [31] Boudia M, Louly MAO, Prins C. Combined optimization of production and distribution. CD-ROM proceedings of the international conference on industrial engineering and systems management, IESM’05, Marrakech, Morocco, 16–19 May 2005. Mons, Belgium: I4E2. [32] Ribeiro CC. Private communication, August 2005.