Heuristics for a multiperiod inventory routing problem with production decisions

Computers & Industrial Engineering 57 (2009) 713–723

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Heuristics for a multiperiod inventory routing problem with production decisions Jonathan F. Bard a,*, Narameth Nananukul b a b

Graduate Program in Operations Research & Industrial Engineering, 1 University Station C2200, The University of Texas, Austin, TX 78712-0292, United States TidalTV, 3701 N. Lamar Blvd., Suite 300, Austin, TX 78705, United States

a r t i c l e

i n f o

Article history: Received 24 April 2008 Received in revised form 11 November 2008 Accepted 27 January 2009 Available online 5 February 2009 Keywords: Inventory routing Vehicle routing Heuristics Branch-and-price

a b s t r a c t Manufacturers who resupply a large number of retailers on a periodic basis continually struggle with the question of how to formulate a replenishment strategy. This paper presents a comparative analysis of a series of heuristics for an inventory routing problem (IRP) that arises in a manufacturing supply chain. The IRP is formulated as a mixed integer program with the objective of maximizing the net benefits associated with making deliveries in a specific time period to a widely dispersed set of customers. It is assumed that inventory can accumulate at the customer sites, but that all demand must be met without backlogging. Because optimal solutions were not within reach of exact methods, a two-step procedure was developed that first estimates daily delivery quantities and then solves a vehicle routing problem for each day of the planning horizon. As part of the methodology, a linear program is used to determine which days it is necessary to make at least some deliveries to avoid stockouts. The IRP is investigated in the context of an integrated production–inventory–distribution–routing problem (PIDRP). The full model takes into account production decisions and inventory flow balance in each period. For the computations, a previously developed branch-and-price algorithm is used that requires the solution of multiple IRPs (one in each period) to generate columns for the master problem. Testing showed that PIDRP instances with up to eight time periods and 50 customers can be solved within 1 h. This level of performance could not be matched by either CPLEX or an exact version of the branch-and-price algorithm.  2009 Elsevier Ltd. All rights reserved.

1. Introduction In vendor managed inventory (VMI) replenishment systems (e.g., see Cetinkaya, Mutlu, & Lee, 2006; Zhao, Lai, & Wang, 2007) the manufacturer observes and controls the inventory levels of its customers, as opposed to conventional approaches where customers monitor their own inventory and decide the time and amount of products to reorder. One of the benefits of VMI is that it permits a more uniform utilization of transportation resources leading to lower distribution costs. Customers benefit from higher service levels and greater product availability due to the fact that vendors can use existing inventory data at their customer sites to more accurately predict future demand. When a single company is responsible for both planning and scheduling, efficiencies are realized at all nodes in the supply chain (Lee, Kang, & Lee, 2008; O’Brien & Tang, 2006). To achieve optimal performance of a VMI or related system, three key decisions must be made on a periodic basis. The first involves the selection of customers to visit each day. Depending on their demand, their location, storage options, costs, and fleet capac* Corresponding author. Tel.: +1 512 471 3076; fax: +1 512 232 1489. E-mail addresses: [email protected] (J.F. Bard), [email protected] (N. Nananukul). 0360-8352/$ - see front matter  2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.cie.2009.01.020

ity, it may be best to schedule a delivery to restock inventory rather than just to meet that day’s demand. The second involves the quantity to deliver to each customer, and the third concerns the construction of routes. The purpose of this paper is to investigate a series of heuristics that can be used by planners in support of these decisions. In the model, we assume that demand is deterministic, that the planning horizon is finite, and that the fleet consists of a given number of homogeneous vehicles. The resultant problem arises in many inventory routing contexts and has been studied widely (Abdelmaguid & Dessouky, 2006; Bard, Huang, Jaillet, & Dror, 1998; Dror & Ball, 1987; Golden, Assad, & Dahl, 1984). When the goal is to coordinate the primary components of the supply chain, we get what is referred to as the production–inventory–distribution–routing problem (PIDRP). Lei, Liu, Ruszczynski, and Park (2006) were the first to formulate the PIDRP as a mixed integer program (MIP) and proposed a two-phase solution approach that avoided the needed to address lot-sizing and routing simultaneously. Boudia, Louly, and Prins (2006, 2007) developed a similar MIP and proposed both a memetic algorithm with population management (MAPA) and a reactive greedy randomized adaptive search procedure (GRASP) with path-relinking (Resende and Ribeiro 2005) as solution methodologies. Following their lead, Bard and Nananukul (2009) developed a tabu search procedure that provided slightly improved results. Their approach involved

714

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

the solution of an allocation model patterned after the phase one model of Lei et al. coupled with a vehicle routing heuristic to find an initial production and delivery schedule. Path-relinking was also used in a post-processor phase to seek out marginal improvements. In the next section, the problem is formally stated and the related literature is reviewed. In Section 3, we present the mathematical model for our version of the inventory routing problem (IRP) and then discuss its extension to the PIRDP, which is the problem that we ultimately wish to solve. In Section 4, we outline our column generation scheme for the PIRDP and introduce three heuristics for solving the IRP component, the main contribution of this paper. Although the IRP can be solved optimally with a commercial code for relatively small instances, runtimes turned out to be much too long for practical purposes. As an alternative, we developed a two-step approach in which the delivery quantities are first estimated and then a VRP heuristic is run to find routes. Test results are reported in Section 5 for a wide range of instances based on standard data sets. We close with an assessment of the methodology in Section 6.

2. Problem statement and literature review Although the IRP can be defined more generally, our focus is on meeting the periodic demand of a set of customers for a single commodity. Examples might include fuel, beverages and retail goods. Initially, it is assumed that an unlimited amount of inventory can be stored at the customer sites as well as the warehouse from which it is delivered; however, transfers between sites and stockouts are not permitted (cf. Herer, Tzur, & Yucesan, 2006). Demand is considered to be deterministic and known. Deliveries are made by a fleet of homogeneous vehicles that begin and end their run at the warehouse. In the scenarios investigated, a VRP must be solved daily to determine which customers to visit and how much product to deliver to each (for an example, see Kant, Jacks, & Aantjes, 2008). Due to either vehicle capacity limits or favorable cost tradeoffs, it may be desirable to visit a customer on a day in which ample stocks exist at his site in order to build up inventory. In fact, this might be the only feasible option if all demand is to be met. In a more structured version of the problem, customers would be clustered or grouped in a preprocessing step and serviced by the same vehicle each day a delivery is deemed necessary (e.g., see Gaur & Fisher, 2004). Golden et al. (1984) were the first to investigate the interrelated problem of inventory allocation and vehicle routing. Their particular application involved an energy-products company that distributed liquid propane to its customers. They developed a heuristic for designing an integrated delivery planning system aimed at comparing the distribution rule used by the company with their approach. In Dror and Ball’s (1987) version of the IPR, the objective was to minimize the long-term average delivery costs subject to no customers running out of stock at any time. Because the full problem was too large to solve directly, a decomposition heuristic was proposed that had a look-ahead feature to account for costs beyond the current day. This allowed the long-term effects to be transposed to the short-term, thus yielding a much smaller problem to solve. One component of the heuristic included a route improvement scheme that involved node interchange operations for a given set of routes. Anily and Federgruen (1990) analyzed fixed partition policies for an IRP with constant deterministic demand rates and an unlimited number of vehicles. The routing patterns were determined by using a modified circular partition scheme. Based on the results, customers within a partition were assigned to one or more regions. Demand was calculated by summing the demand of the individual

customers assigned to a region. When a customer appeared in multiple regions a percentage of his demand was allocated to each. The routing logic required that all customers in a region be visited as long as there was a need to visit one. A lower bound on the longrun average cost was also determined to provide an understanding of the effectiveness of the routing approach. Anily and Federgruen (1993) extended their earlier work to include the case where the depot can hold inventory. A new algorithm in which the interval between two visits to a region follows the power-of-two (POT) principle was developed. They showed that the difference between the cost from this algorithm and a lower bound on the minimum cost obtained by a similar approach proposed by Roundy (1985) for uncapacitated systems was at most 6% for a sufficiently large number of customers. Herer and Roundy (1997) investigated the one-warehouse multi-retailer distribution problem where the vehicle routing cost was based on the length of the traveling salesman tour beginning and ending at the central warehouse and visiting a subset of retailers. In the model, they assumed that a central warehouse receives deliveries from an outside supplier(s), and that inventories can be held at both the warehouse and at the retailer sites. Demand was specified for the retailers only and was assumed to be deterministic. The objective was to minimize the long-run average costs, which consisted of fixed ordering costs and inventory holding costs. The authors developed a number of polynomial time heuristics, which produced results that were provably close to the cost of an optimal policy. They showed that given a submodular function that is a good approximation of the true ordering cost, there exists a POT policy whose cost is only moderately greater than the cost of an optimal policy. Bertazzi, Paletta, and Speranza (2002) considered a deterministic version of the IRP with a single capacitated vehicle over a longterm period. The objective was to minimize the cost of deliveries to a set of the retailers while ensuring that inventory levels were maintained between some minimum and maximum value. They presented a two-phase heuristic to determine the vehicle route at each discrete time point while using an order-up-to policy to management inventory. In the initialization phase, a feasible solution is built with a sequential procedure that inserts a retailer at each iteration. When a retailer is being considered, a set of delivery time instants is determined by solving a shortest-path problem. Then, for each of the selected time instants, the retailer is inserted into the route at the point that minimizes the incremental cost. In the second phase, the incumbent is improved by performing swaps that removed or inserted retailers at different positions of the route used by the vehicle. Cetinkaya et al. (2006) investigated the impact of alternative outbound dispatch policies on integrated stock replenishment and transportation decisions in a VMI system. In their model, retailers were geographically grouped and it was assumed that replenishment costs were fixed. Demand was generated randomly for each group and deliveries were made in batches of consolidated loads. In the analysis, they compared the performance of timebased and quantity-based policies. Generally speaking, a timebased policy ships accumulated loads every T periods, whereas a quantity based policy ships accumulated loads when a specific dispatch quantity q is reached. Zhao et al. (2007) studied the integration of inventory control and vehicle routing for a distribution system in which a set of retailers with constant rates of demand were resupplied with a single item from a central warehouse. The objective was to determine inventory policies and routing strategies such that the long-run average costs were minimized and all demand was satisfied. In their model, no inventory capacity constraints were imposed on the warehouse or on the retailers. Testing was done on problem instances with 50 and 75 retailers.

715

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

A problem closely related to ours is the periodic vehicle routing problem (PRP), which assumes that the delivery patterns defined by visit frequencies or visit days are given in advance. The aim is to select the most suitable pattern for each customer. Gaudioso and Paletta (1992) proposed several algorithms for the PRP when the objective is to minimize the maximum number of vehicles used simultaneously during any day of the planning horizon. In their model, each customer i requires the delivery of a specified amount of commodity exactly once in each time interval of Ti days during a period of s days. They also required that two successive deliveries be spaced between Ki and Ui days apart. When two or more routes were disjoint, it was permissible to service them with the same vehicle. Gaur and Fisher (2004) developed and implemented a system for a PRP at Albert Heijn, a supermarket chain in the Netherlands with 207 stores. Their model assumes that there is a single homogeneous product with deterministic, time-varying demand. Although the application involved the replenishment of several items, all items were replenished from a single distribution center so they were able to formulate an equivalent single-product problem. This was done by expressing the demand rates of all items in units of weight and volume, translating the result into fractions of containers, and aggregating. To find solutions, a fixed partition policy was applied. For such a policy, customers are partitioned into disjoint subsets or regions and each region is served separately and independently from the others. In their version of the PRP, Mourgaya and Vanderbeck (2007) proposed a column generation-based heuristic to fix dates for customer visits and to assign customers to vehicles. The daily sequencing decisions were left to an operational model. In formulating the problem two objectives were considered: (1) optimizing the compactness of the geographical regions to which customers are assigned, and (2) balancing the workload evenly between vehicles. Using a Dantzig-Wolfe reformulation, they found that the resulting master problem provided substantially stronger lower bounds than the LP relaxation of the original problem and that there were fewer difficulties due to symmetry during branch and bound. Computational tests were performed using 20, 50 and 80 customer data sets that involved comparisons of various implementation options related to master problem initialization, column generation strategies, number of passes in the rounding heuristic, and problem specifications. 3. Model formulation For the basic IRP, we are given a set of n customers geographically dispersed on a grid and a single warehouse where a unique commodity is stored. Each customer i has a fixed nonnegative demand dit that must be satisfied in each time period t of the planning horizon; i.e., shortages are not permitted. A benefit of ait is realized for each unit delivered to customer i in period t, and a cost of cij is realized when the link between customers i and j is traversed. Deliveries are made by a fleet of h homogeneous vehicles but not all need be routed each day. The goal is to construct a delivery schedule that maximizes the total net benefit after the costs are deducted while ensuring that each customer’s demand is met over the planning horizon s. Implicit in the solution will be a daily distribution of additional items to be placed in inventory. For now, production and inventory holding costs are not considered. In the development of the model, we make use of the following notation. Indices and sets i, j indices for customers, where 0 corresponds to the warehouse t index for periods or days

N T

set of customers; N0 = N [ {0}; |N| = n set of days in the planning horizon; T0 = T [ {0} and |T| = s

Parameters dit demand of customer i on day t ait benefit of delivering one unit of product to customer i in period t h number of available vehicles Q capacity of each vehicle cij cost of going from customer i to customer j Decision variables xijt 1 if customer i immediately precedes customer j on a delivery route on day t; 0 otherwise yit load on a vehicle immediately before making a delivery to customer i on day t wit amount delivered to customer i on day t

3.1. Model

XX

/IRP ¼ Maximize

t2T

i2N

X

subject to

ait wit 

8i 2 N; t 2 T

xijt  1;

xijt ¼

i2N0 i–j

X

cij xijt

ð1aÞ

i2N 0 j2N 0 t2T i–j j–i

j2N 0 j–i

X

XXX

X

xjit ;

8j 2 N; t 2 T

ð1bÞ ð1cÞ

i2N0 i–j

8t 2 T

x0jt  h;

ð1dÞ

j2N

ð1  xijt Þ; yjt  yit  wit þ Dmax t

8i 2 N; j 2 N0 ; t 2 T wit  Dmax it

X

ð1eÞ

xijt ;

j2N0

8i 2 N; t 2 T

ð1fÞ

s xijt 2 f0; 1g; 0  yit  Q ; wit 2 Zn þ ;

(

8i–j 2 N0 ; t 2 T

¼ min Q; where Dmax it ( Dmax ¼ min Q ; t

s X

)

dil

ð1gÞ

and

l¼t

s XX i2N

)

dil

l¼t

The objective function (1a) maximizes the difference between the benefits associated with making deliveries and the transportation costs. In the application that we have in mind, the coefficient ait is a function of the marginal cost of holding inventory at the warehouse versus holding it at the customer site, as well as the marginal benefit of delivering one more unit to customer i in period t (see Section 4 for its derivation). Constraints (1b)–(1e) represent the routing aspect of the problem. Constraints (1b) ensure that if customer i is visited on day t (xijt = 1) then it can have at most one successor j on its route, which may be the warehouse. This constraint must be written as an inequality because we do not necessarily want to service each customer every day. Route continuity is enforced by (1c); if a vehicle arrives at customer j on day t, then a vehicle must depart customer j on day t. Logic, and the fact that the fleet is homogeneous, requires that it be the same vehicle. The number of vehicles that leave the warehouse on any day t is limited to the number available, h, as indicated by (1d). Constraint (1e) keeps track of the load on the vehicles and guarantees that on

716

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

day t, if customer i is the immediate predecessor of customer j on a route, then the load on the vehicle before visiting customer j must be less than or equal to the load just before visiting customer i minus the amount delivered, which is represented by the variable is specified to be as small as wit. The value of the parameter Dmax t possible while ensuring that (1e) is always feasible. Because the load on each vehicle is monotonically decreasing as customers are visited, (1e) provides the added benefit of eliminating subtours that do not include the factory. The specific amount delivered to customer i is limited by the in (1f), which is the smaller of the vehicle capacparameter Dmax it ity Q or the cumulative demand from day t to the end of planning horizon s. When customer i has no successors on a route, no delivery is possible because the right-hand side of (1f) would be zero. After all deliveries are made, the fleet returns to the warehouse empty so y0t can be set to 0 for all t e T. To conclude the formulation, variables are defined in (1g). It is straightforward to generalize model (1) to include multiple warehouses and multiple products. The size of model (1) is determined largely by constraint (1e) and the number of binary variables, xijt, both of which grow at a rate proportional to O(n2s). The majority of the other constraints only grow at a rate proportional to O(ns). A small problem with three vehicles, 30 customers, and a 5-day planning horizon contains roughly 5000 constraints, 4500 binary variables, and 300 continuous variables. Initial attempts to solve instances of this size with CPLEX 8.1 were not encouraging. When a time limit of 60 min was imposed, optimality gaps between 5% and 10% were the norm. This was not surprising since model (1) is NP-hard in the strong sense due to the fact that when wit is held constant, a VRP results. 3.1.1. Extensions to the PIDRP The integrated version of our problem includes a production facility collocated with the warehouse and inventory considerations. If production takes place at the facility in period t, then a setup cost ft is incurred, t = 0, 1, . . . , s. A limited number of items can be produced in each time period and a limited number can be stored at the warehouse for a unit cost of hP, or at cusC tomer i’s site at a unit cost of hi . An additional assumption is that all deliveries take place at the beginning of the day and arrive in time to satisfy demand for at least that day. All production in period t is available for delivery only on the following morning (this is common in food production and distribution; e.g., see Villegas & Smith, 2006) and all inventories are measured at the end of the day. Demand on day t can be met out of deliveries from the warehouse on day t and ending inventory on day t – 1 at the customer site. Initial customer inventory on day 0 simply reduces demand on subsequent days, while initial inventory at the factory must be routed; at the end of the planning horizon all inventory is required to be zero. The goal for the PIDRP is to construct a production plan and delivery schedule that minimize the sum of all costs while ensuring that each customer’s demand is met each day. Under the common assumption that unit production costs are constant, the decision to overproduce and hold items in inventory will be a function of the setup cost at the factory, system holding costs, production and storage capacities, vehicle routes, and daily demand. Some additional notation is needed to represent the full model. Parameters maximum inventory that can be held at the production I Pmax facility I Cmax :i maximum inventory that can be held by customer i C daily capacity of the factory

ft

setup cost at production facility on day t

Decision variables production quantity on day t pt 1 if there is production on day t; 0 otherwise zt inventory at the production facility at the end of day t IPt inventory at customer i at the end of day t ICit

3.2. Model

XXX

/PIDRP ¼ Minimize þ

X

X

cij xijt þ

t2T i2N0 j2N 0

X X

P

h IPt þ

t2T 0 nfsg

f t zt

t2T C

hi ICit

t2Tnfsg i2N

subject to I Pt ¼ I Pt1 þ pt 

ð2aÞ

X

wit ;

8t 2 T 0

ð2bÞ

i2N C þ wit  dit ; 8i 2 N; t 2 T IitC ¼ Ii;t1 X wit  IPt1 ; 8t 2 T

ð2cÞ ð2dÞ

i2N

pt  Czt ; 8t 2 T 0 =fsg X p0  ðdi1  ICi0 Þ

ð2eÞ ð2fÞ

i2N

0  I Pt  I Pmax ; 0  yit  Q ;

0  I Cit  I iC max ; pt 2 Zsþ ;

IPs ¼ ICis ¼ 0;

8i 2 N;

zt 2 f0; 1g;

s wit 2 Zn þ ;

t2T

ð2gÞ

and (1b)–(1g) The objective (2a) is now to minimize the sum of transportation costs, production setup costs, holding costs at the factory, and holding costs at the customer sites. Because all demand must be met, production costs, which are assumed to be linear and constant, can be omitted, as can any initial inventory costs. Constraints (2b) and (2c) are inventory flow balance equations in which it is assumed that the initial inventories I P0 and I Ci0 are given for all customers i e N. Constraints (2e) limit production on day t to the capacity of the factory. The assumption that items produced on day t are only available for delivery on day t + 1 implies that pt  I Pmax and ps ¼ 0. To ensure that demand on day 1 can be met, it is necessary to include (2f) which allows production on P day 0. If IP0 ¼ ICi0 ¼ 0; then p0  i2N di1 or the problem is infeasible. 4. Branch and price Several heuristics have been developed that provide good solutions to the PIDRP for up to 200 customers and 20 time periods (see Boudia, Louly, & Prins, 2006, 2007; Bard & Nananukul, 2009). Finding exact solutions, though, has been problematic for all but very small instances. Bard and Nananukul (2008) developed a branchand-price (B&P) methodology whose performance was comparable to that of CPLEX 8.1. Table 1 highlights test results for instances with up to 10 customers and 8 time periods using a 30-min time limit. As can be seen, the largest instance (problem no. 2) that could be solved to optimality contained 10 customers and two time periods, where ucplex and /B&P are the best objective function values returned by CPLEX and B&P, respectively, and tcplex and tB&P are the corresponding runtimes. Informally speaking, the idea behind B&P is to represent a linear optimization problem as a set of columns such that each column is associated with a partial solution, in our case, a delivery schedule for one day. In the model, the right-hand-side values are equal to the demand. When all feasible columns are enumerated, forming what is called the master problem (MP), the optimal solution to

717

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723 Table 1 Solutions from CPLEX and B&P algorithm. Problem no.

s

n

(/cplex) CPLEX solution

(tcplex) CPLEX runtime (s)

(/B&P) B&P solution

(tB&P) B&P time (s)

[(/B&P  /cplex)//cplex] 100%

1 2 3 4 5 6 7

2 2 4 4 6 6 8

5 10 5 10 5 10 5

21,961 38,776 24,416 48,287 44,377 71,546 54,558

8 62 >1800 >1800 >1800 >1800 >1800

21,961 38,776 24,416 52,321 44,831 70,804 56,073

34.3 1075 >1800 >1800 >1800 >1800 >1800

0 0 0 8.4 1.02 –0.29 3.77

the original problem can be found by identifying the subset of columns that satisfies all the demand at minimum cost [see Wolsey (1998) for a full description of the methodology]. Rather than trying to enumerate all feasible columns, it is common to take an iterative approach and generate promising columns one or more at a time, and then solve the LP relaxation of the corresponding restricted master problem. When it is no longer possible to generate promising columns, branch and bound (B&B) is used to arrive at an integer solution. At each node of the search tree, column generation is used to find additional candidate schedules. The restricted MP is generally created from those constraints that make the original problem difficult. For the PIDRP, it is natural to separate the production and inventory constraints (2b)– (2g), which cut across two time periods each, from the routing and delivery constraints (1b)–(1f), which can be written separately for each time period. We used the former along with the objective function in (2a) to construct the MP, and the latter to construct s subproblems, SPt, t e T. During column generation, the subproblems act as schedule generators guided by the values of the dual variables associated with the LP solution to a restricted MP. This is done by solving the subproblems given by model (1) to determine whether one or more schedules can be identified that have a negative reduced cost. When added to the restricted MP, such schedules in the form of columns will improve the current LP solution, which serves as a lower bound on the solution to the PIDRP (1)–(2). If all reduced costs are nonnegative but fractional variables appear in the MP solution, B&B is performed. The details of our implementation are given in Bard and Nananukul (2008). It is instructive, though, to see how the benefit coefficient, ait, in (1a) is computed. For subproblem t, let a1t be the dual variable associated with the factory inventory constraint (2b), let a2it be the dual variable associated with customer i’s inventory constraint (2c), and let a3t be the dual variable associated with the delivery constraint (2d). Only a3t is required to be nonnegative. For each day t, the reduced cost of a generic column in MP is

ct ¼ ct  a1t

X i2N

wit 

X i2N

a2it wit þ a3t

X

wit þ a4t

i2N

where ct is the cost of a generic route and a4t is a constant term that arises from the MP formulation. Collecting terms in wit gives ait ¼ a1t þ a2it  a3t . A closer examination of the statistics in Table 1 indicated that the total amount of time devoted to generating the columns was 1275 s, on average, which is approximately 88% of the total runtime, tB&P. This was not surprising in light of the complexity of the subproblem. In an effort to improve on these results, we investigated three heuristics for the IRP (1) along with a model for determining periods in which at least one customer requires a delivery. This model is used in a preprocessing step to fix variables associated with the delivery decisions. Each is discussed in the next section.

5. Heuristic modifications and preprocessing The amount of effort required to solve SPt increases significantly as either the number of customers n or the number of time periods s increases. To ease this burden, we propose a two-step heuristic for solving the subproblems (cf. Choi & Tcha, 2007; Mourgaya & Vanderbeck, 2007; Savelsbergh, 1997). The first step involves developing a model for determining delivery quantities for each customer in each time period. Recall that SPt is not a pure VRP. The second step involves finding actual routes in light of the current set of branching constraints. This is done with a VRP tabu search code that we adjust as necessary. 5.1. Heuristic model (1) for determining delivery quantities VRP algorithms require the delivery quantities to be known, even for cases in which visiting a customer is optional as in the prize collecting version of the problem. One way to estimate these quantities is to solve a linear program whose objective function coefficients are based, in part, on the current values of the MP dual variables. Ideally, we would like to know a priori ci, the cost of visiting customer i. In reality, ci can only be determined after the routing decisions embodied in xijt are made and is equal to Rj cij xijt in period t. Alternatively, we can estimate its value by, say, fitC ¼ 2ci0 for all i and t, which is the cost of a round trip between the warehouse and customer i. For ait > 0, there is some value in making a delivery to customer i in period t. The ratio ait =ðci þ a4t Þ gives the relative value for each customer. We would like this ratio to be greater than or equal to 1 but as noted, it is impossible to know the value of ci a priori. For our purposes, then, it is sufficient to consider all i for which ait > 0 and to remove the constant term a4t from the ratio. The first heuristic model (H1) is as follows.

Maximize subject to

n X s X

maxf0; ait =ci gwit

i¼1 t¼1 n X

wit  hQ;

t ¼ 1; . . . ; s

ð3aÞ ð3bÞ

i¼1 t X

wil 

l¼1

0  wit 

t X

dil ; l¼1 Dmax it ;

i ¼ 1; . . . ; n; t ¼ 1; . . . ; s i ¼ 1; . . . ; n; t ¼ 1; . . . ; s

ð3cÞ ð3dÞ

The objective function (3a) is designed to maximize the value of the items delivered to all customers over the planning horizon. The coefficients in (3a) will change at each MP iteration. Constraints (3b) limit the delivery quantities in each period to the fleet capacity. Constraints (3c) ensure that a sufficient amount is delivered to each customer over the planning horizon to meet his or her demand in each period. Bounds are placed on the delivery variables in (3d). The upper bound limits the amount delivered to any customer i in period t to the minimum of the capacity of a vehicle or the remaining demand for that customer from t through s. It does not restrict the total amount delivered over the planning horizon to be equal to the total demand.

718

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

Model (3) ignores the routing aspect of the PIDRP so there may be some periods in which it is not possible to service all those customers whose wit values are greater than 0. This issue can be addressed in two ways. First, we could reduce the righthand side of (3b) by some fraction to compensate for the routing constraints, or second, the VRP code discussed presently could be modified to accommodate infeasible solutions by heavily penalizing violations of (3b) rather than prohibiting them. If the capacity of a vehicle was exceeded in a solution, then the individual delivery quantities wit could be reduced incrementally, say, starting with the customer whose coefficient ait > 0 was the smallest, until the violation was removed. During the reduction process, it would be necessary to ensure that (3c) remained satisfied and that the reduced cost of the solution remained negative. Proposition 1. The feasible region given by constraints (3b)–(3d) is a relaxation of the feasible region of the original PIDRP. Proof. The load variable in the original problem, yit, is bounded by ^t ¼ maxfyit : i ¼ 1; . . . ; ng and forming Q in each period t. Letting y ^t gives an upper bound on the delivery quantity in the product hy period t, which is equivalent to (3b). Constraints (3c) can be derived by substituting out the inventory variables, ICit , in Eq. (2c) of the PIDRP and noting that ICit  0 for all i and t. Finally, (3d) is a relaxation of (1f) in the original model. h Proposition 2. An optimal solution exists to model (3) in which wit is integral for all i  N and t  T when all the problem data are integral. Proof. We show that if any of the wit are fractional, the solution cannot be optimal. In the simplest case in which there is at most one fractional value of wit in period t, it is always possible to increase its value to dwit e without violating any of the constraints because all the data are assumed to be integral. The objective function would then increase by ðdwit e  wit Þ  ðmaxf0; ait =ci gÞ in period t. For the case in which, say, wi1 t and wi2 t are fractional in period t, let us assume without loss of generality that maxf0; ai1 t =ci1 g  maxf0; ai2 t =ci2 g. Now, increasing and decreasing the delivery quantities to customers i1 and i2 by some arbitrarily small amount e, respectively, produces a marginal improvement in the objective function value without violating either constraints (3b)–(3d). The increase in wi1 t cannot lead to a violation of (3c) nor can the reduction of wi2 t . The latter claim follows from the fact that (3c) must be nonbinding for those constraints in which wi2 t appears due to the fact that the right-hand sides are integral. All other cases are generalizations of these two and so the same arguments apply to show that either the solutions cannot be optimal or that an equivalent integral solution can be constructed. h

t X

wil 

l¼1

Our inability to specify exactly the cost of visiting customer i in period t may distort the selection of delivery quantities wit in the solution to (3a)–(3d). A more accurate model can be formulated by reintroducing the routing variable xijt. In doing so, we limit the reduced feasible region to include only the assignment constraints so that optimal solutions can be obtained quickly. The augmented model, call it H2, is

Maximize

n X s X i¼1

subject to

n X i¼1

ait wit 

t¼1

wit  hQ ;

n X n X s X i¼0

j¼0

cij xijt

ð4aÞ

t¼1

t ¼ 1; . . . ; s

ð4bÞ

i ¼ 1; . . . ; n; t ¼ 1; . . . ; s

dil ;

ð4cÞ

l¼1

0  wit  Dmax it

n X

i ¼ 1; . . . ; n; t ¼ 1; . . . ; s ð4dÞ

xijt ;

j¼0 j–i

routing constraints ð1bÞ—ð1dÞ xijt 2 f0; 1g; 8i; j; t

ð4eÞ

The objective function (4a) is the same as that in (1a). Constraints (4b)–(4d) are the same as their counterparts in (3) except that the upper bound in (4d) has been modified to prevent a delivery to customer i if no routing assignment is made. The routing constraints (1b)–(1d) ensure balance of flow in a tour and that at most h vehicles leave the warehouse in period t. In (4e), binary restrictions are placed on the flow variables. A solution to model (4) is a weak approximation to the collective solutions of the s subproblems, SPt, t e T. The primary differences are that the individual vehicle capacity restrictions are not maintained in (4) and that a solution may have subtours in each time period that do not include the warehouse. Nevertheless, the advantage of H2 is that it turns out to be much easier to solve than the IRP (1), and like H1, it provides reasonable values for the delivery quantities wit, for all i e N and t e T, that can serve as input to the VRP code (Carlton & Barnes, 1996). The third model proposed for determining delivery quantities, call it H3, is similar to H2 (4) except (4c) is removed. The advantage of H3 is that the runtimes are much less than for H2, since (4a), (4b); (4d), (4e) decompose into s separate problems. 5.3. Determining periods with delivery requirements During B&P, subproblem solutions may yield null columns because there is no explicit requirement that a delivery be made in any particular period. This can lead to a form of degeneracy where null columns are contained in a solution, implying that a fraction of a route has been selected. During branch and bound this issue is ultimately resolved, but infeasible solutions may introduce a degree of inefficiency in the computations. One way to strengthen the B&P approach is to identify periods in which at least one customer must be visited. The following sequence of s  1 LPs patterned after model (3) can be used for this purpose. Problem Pt; t = 2, . . . , s

ft ¼ Minimize subject to

n X i¼1 n X

wit

ð5aÞ

c ¼ 2; . . . ; s

wic  hQ ;

ð5bÞ

i¼1

c X l¼1

5.2. Heuristic models (2) and (3) for determining delivery quantities

t X

0 n X

wil 

c X

dil ; l¼1 wic  Dmax ic ;

i ¼ 1; . . . ; n; c ¼ 2; . . . ; s

ð5cÞ

i ¼ 1; . . . ; n; c ¼ 2; . . . ; s

ð5dÞ

wi;c1  fc1 ;

c ¼ 2; . . . t

ð5eÞ

i¼1

where f1 ¼

n X

di1

i¼1

The objective function (5a) determines the minimum quantity that must be delivered in period t subject to the conditions that no shortages occur in any period over the planning horizon and that the load on the h vehicles does not violate the capacity of the fleet. These conditions are ensured by (5b)–(5d) which duplicate (3b)– (3d) except that now the index t has been replaced with c to avoid confusion with sequence index t. The new constraints (5e) require that the aggregate minimum delivery quantities computed for peri-

719

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

ods prior to t be maintained in all subsequent periods; however, the individual quantities wit may change. A solution to Pt with ft > 0 indicates that a delivery must take place in period t so the constraint Riwit P 1 is added to model (1) or the heuristic being used in its place. Nevertheless, the fact that the feasible region (5b)–(5e) is a relaxation of the routing constraints implies that there may exist periods in which ft = 0 but where a delivery is required. Proposition 3. There are no feasible solutions to the PIDRP (1)–(2) which a delivery must take place in period t when the solution of model (5) returns ft = 0. Proof. From Proposition 1, we know that constraints (5b)–(5d) are a relaxation of the original feasible region. Constraint (5e) states that in period t all previously found minimum delivery quantities computed in periods 1, . . . , t  1 cannot be reduced, which is a feasibility condition for the PIDRP. If a feasible solution to the PIDRP P existed in which i wit > 0 for some period t where nt = 0, that solution would also be feasible to (5) and hence contradicts the optimality of nt. h Model (5) is designed to schedule deliveries as early as possible in periods preceding the current period t to achieve a nt as small as possible without regard to production, inventory or routing considerations. When nt = 0, it is often possible to avoid deliveries in subsequent periods by increasing deliveries in periods 1, . . . , t which constraints (5e) permits. This observation has implications for branching when some variables are set to 0. 6. Computational results Two sets of experiments were carried out to evaluate the proposed solution methodology. The first was aimed at determining which of the two-step approaches (H1, H2 or H3 in conjunction with the VRP code) was most suitable for solving the PIDRP (1)– (2), and the second to gain experience with the most effective algorithmic configuration. All computations were performed on a 2.53 GHz processor with 512 MB of RAM. The optimization models and heuristics were implemented in Java Netbean 4.1 and linked to the CPLEX 8.1 libraries. CPU times were obtained through both CPLEX and the time function in Java. The B&P algorithm and its various components were coded with Microsoft Visual Studio .NET 7.0; the MINTO 3.1 libraries were used to manage the search tree.

6.1. Heuristics for the IRP For the initial experiments, we randomly generated seven instances, each on a 100  100 grid. The first step was to fix the number of customers n e {5, 10} and the number of time periods s e {2, 4, 6, 8} by selecting combinations from these sets. The remaining parameter values were generated as follows:  The distance between a pair of customers was measured by the Euclidean norm and the transportation costs were assumed to be proportional to that distance. A proportionality constant of 1 was used.  Setup and holding costs at the factory, and demand and holding costs at the customer sites were generated from a discrete uniform distribution with given upper and lower bounds: – Setup cost at the factory on day t, ft: upper bound = 1000 and lower bound = 500 (values of ft were generated independently for each day t e T) – Holding cost at the factory, hP: upper bound = 200 and lower bound = 100 – Demand on day t, dit: upper bound = 10 and lower bound = 0 (these values imply that for each day t e T, 9.1% of the customers have a chance of having zero demand) C – Holding cost at customer sites, hi : upper bound = 100 and lower bound = 50 for all i e N  Initial inventories at the factory as well as at the customer sites were set to zero; IP0 ¼ ICi0 ¼ 0, "i e N  The number of vehicles h was fixed at 5 and the vehicle capacity Q was determined by calculating the total demand on each day, identifying the largest values, and multiply them by three; P Q ¼ ð3  maxf i2N dit : t ¼ 1; . . . ; sg.  No limit was imposed on the production capacity or the maximum inventory allowed at the factory; C ¼ IPmax ¼ 1  Maximum inventory allowed at a customer site was about three days of demand for that customer The above parameters were chosen to make the model as realistic as possible while avoiding the inefficiency of daily deliveries to a majority of customers. To reduce the implicit delivery frequency, the setup cost at the factory was specified to be higher than the holding costs at both the factory and the customer sites. In addition, the number of vehicles and their capacity were chosen to achieve high vehicle utilization when they are on the road. Low-

Table 2 Results using column generation heuristic H1 to solve the PIDRP at the root node. Problem no.

s

n

(/H1_Root) Root node solution

Strategy 1: Add one column at each iteration 1 2 5 21,223 2 2 10 38,322 3 4 5 31,099 4 4 10 62,896 5 6 5 70,976 6 6 10 126,423 7 8 5 178,832 Problem no.

s

n

ð/H1 Root Þ Root node solution

Strategy 2: Add multiple columns at each iteration 1 2 5 21,223 2 2 10 38,322 3 4 5 22,369 4 4 10 46,771 5 6 5 41,145 6 6 10 67,719 7 8 5 49,267

[(/H1_Root  /Root)//Root] 100%

(jH1_Root) Number of columns

(tH1_Root) Total runtime (s)

[(tH1_Root  tRoot)/tRoot] 100%

0.02 0.05 39.03 36.69 72.62 87.59 263.13

20 31 7 9 11 19 8

2.98 1.95 1.16 4.14 3.19 3 1.78

–8.87 –92.49 –70.78 –98.26 –74.76 –99.56 –91.99

½ð/H1 Root  /Root Þ=/Root  100%

ðjH1 Root Þ Number of columns

ðtH1 Root Þ Total runtime (s)

½ðtH1 Root  tRoot Þ=tRoot Þ 100%

16 31 35 30 40 86 81

2.86 2.03 1.8 1.97 1.95 3.84 3.61

–12.54 –92.18 –54.66 –99.17 –84.57 –99.44 –83.76

0.02 0.05 0.00 1.65 0.07 0.48 0.04

720

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

er vehicle capacities would result in more frequent trips to the customer sites, which translates into less efficient schedules because they would not be able to take advantage of the low holding costs at the customer sites. Preliminary testing showed that when Q was decreased to the point where almost every customer required a visit every day, the computational effort increased noticeably. Increasing either Q or h would result in lower vehicle utilization and an easier problem to solve. Tables 2–4 report the computational results when H1, H2 and H3 are used, respectively, instead of CPLEX to solve the PIDRP at the root node of the B&P search tree. The following two column generation strategies were used: (i) at each iteration begin solving the subproblems and insert only the first column that prices out negatively into MP, and (ii) find as many columns as possible and insert all of them into MP. The results in the top portion of the tables are for the first strategy, and the bottom for the second. Columns 4 and 5 give the root node solutions (/H#_Root) and their deviation from the solutions obtained with exact column generation (/Root) (recall that the results in Table 1 are for the entire problem and not just the root node). The number of columns generated at the root node (jH#_Root) is reported in column 6. The runtime (tH#_Root) and the deviation from the runtime obtained with exact column generation (tRoot) are reported in columns 7 and 8, respectively. For strategy 1, after some additional manipulation of the data in Tables 2–4, we can see that H1 gives the smallest runtimes, tH1_Root, on average, which are 76.67% less than the average values of tRoot; however, /H1_Root is approximately 71.31% larger than /Root. This calls into question the effectiveness of H1 because good heuristics should give solutions relatively close to the optimum, /Root. The deviation of /H3_Root and /H2_Root from /Root are 0.03% and 0.09% on average, implying that H3 gives the best solutions at the root node. Also, since tH3_Root is 48.23% lower than tH2_Root on average, and the number of columns generated (jH3_Root) is roughly the same as jH2_Root, H3 is arguably the best choice of the three. The solutions provided by H3 have the smallest deviations from /Root with 71.17% shorter runtimes when tH3_Root is compared to tRoot, and 23.67% fewer columns when jH3_Root is compared to jRoot. For strategy 2, we found from the data in Tables 2–4 that when H1, H2 and H3 are respectively used to generate multiple columns at the root node (an asterisk  is used to distinguish these results from single column generation), on average, tH1 Root (total runtime), is slightly less than tH1_Root, whereas tH2 Root is 54.67% less than tH2_Root, and t H3 Root is 56.56% less than tH3_Root on average. For H2

and

H3

the

explanation for these results is that decrease with respect to jH2_Root and jH3_Root so the master problem LPs are smaller. For H1, the number of columns generated jH1 Root is greater, on average, than jH1_Root, but the number of LPs solved is less. In addition, we see that the deviation from /Root is reduced from 71.31% for /H1_Root to about 0.33% for /H1 Root , which shows that generating multiple columns provides a large improvement in the performance of H1. However, the deviation of /H1 Root from /Root is still relatively high compared to the deviations of /H3 Root and /H2 Root from /Root, which are only 0.04% and 0.08%, respectively. This analysis indicates that H2 and H3 are superior to H1 with respect to solution quality, and that the deviation of /H3 Root and /H2 Root from /Root is about the same as the deviation of /H3_Root and /H2_Root from /Root. Looking at the number of columns generated, we see that jH1 Root is 203% larger than jH1_Root; however, jH2 Root and jH3 Root are 36.24% and 31.17% less than jH2_Root and jH3_Root, respectively. When we consider only H1 and H2 the results show that generating multiple columns at each iteration is beneficial in term of shortening runtimes, reducing the overall size of MP, and providing good solution quality. The analysis also shows that H2 and H3 with multi-column generation outperform H1 and exact column generation with respect to runtime. In addition, their average statistics for solution quality and the number of columns generated are about the same. A significant difference arises, though, when we consider runtime: tH3 Root is 50.39% less than tH2 Root , which suggests that overall, H3 with multi-column generation is the best choice. On average, H3 gives good solutions, leads to relatively small master problems, and runs quickly. Nevertheless, when H2 and H3 with multi-column generation were incorporated in the B&P algorithm we found that H2 was more efficient as measured by the size of the search tree, runtime, and solution quality. B&P could only solve small instances when H3 was used. Table 5 reports test results comparing H2 with H3 on the first six problem instances. The computations were terminated at 1800 sec. In Table 5, columns 4 and 5 give the respective B&P solutions when H2 and H3 were used. Columns 6 and 7 give the corresponding number of nodes in the search tree at termination, and columns 8 and 9 give the number of columns generated. The last two columns report the total runtime. From these statistics, we see that on average the number of nodes (oH3_B&P) and the number of columns generated (jH3_B&P) are more than 60% larger than oH2_B&P

jH2 Root and jH3

Root

Table 3 Results using column generation heuristic H2 to solve the PIDRP at the root node. Problem no.

s

Strategy 1: Add one column 1 2 2 2 3 4 4 4 5 6 6 6 7 8 Problem no.

s

n

(/H2_Root) Root node solution

at each iteration 5 21,223 10 38,322 5 22,394 10 46,151 5 41,121 10 67,457 5 49,259 n

ð/H2 Root Þ Root node solution

Strategy 2: Add multiple columns at each iteration 1 2 5 21,223 2 2 10 38,322 3 4 5 22,369 4 4 10 46,147 5 6 5 41,124 6 6 10 67,482 7 8 5 49,268

[(/H2_Root  /Root)//Root] 100%

(jH2_Root) Number of columns

(tH2_Root) Total runtime (s)

[(tH2_Root  tRoot)/tRoot] 100%

0.02 0.05 0.11 0.3 0.01 0.09 0.02

16 39 43 172 118 168 98

1.22 2.84 3.83 16.14 7.56 48.06 19.34

–62.7 –89.06 –3.53 –93.22 –40.19 –93.02 –13

½ð/H2 Root  /Root Þ=/Root  100%

ðjH2 Root Þ Number of columns

ðtH2 Root Þ Total runtime (s)

½ðtH2 Root  t Root Þ=tRoot  100%

0.02 0.05 0.00 0.29 0.02 0.13 0.04

13 32 36 86 51 123 76

3.02 3.86 2.38 11.38 3.48 15.83 4.89

–7.65 –85.13 –40.05 –95.22 –72.47 –97.70 –78.00

721

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723 Table 4 Results using column generation heuristic H3 to solve the PIDRP at the root node. Problem no.

s

n

(/H3_Root) Root node solution

Strategy 1: Add one column at each iteration 1 2 5 21,218 2 2 10 38,305 3 4 5 22,369 4 4 10 46,094 5 6 5 41,122 6 6 10 67,419 7 8 5 49,250 Problem no.

s

ð/H3 Root Þ Root node solution

n

Strategy 2: Add multiple columns at each iteration 1 2 5 21,218 2 2 10 38,312 3 4 5 22,369 4 4 10 46,087 5 6 5 41,122 6 6 10 67,441 7 8 5 49,250

[(/H3_Root  /Root)//Root] 100%

(jH3_Root) Number of columns

(tH3_Root) Total runtime (s)

[(tH3_Root  tRoot)/tRoot] 100%

0.00 0.01 0.00 0.18 0.01 0.04 0.00

10 22 34 119 56 273 102

1.11 3.38 2.48 9.55 4.67 19.21 10.81

–66.06 –86.97 –37.53 –95.99 –63.05 –97.21 –51.37

½ð/H3 Root  /Root Þ=/Root  100%

ðjH3 Root Þ Number of columns

ðtH3 Root Þ Total runtime (s)

½ðtH3 Root  t Root Þ=tRoot  100%

0.00 0.02 0.00 0.16 0.01 0.07 0.00

16 52 32 85 53 118 68

2.39 4.76 1.5 3.75 2.22 5.25 2.36

–26.91 –81.66 –62.22 –98.42 –82.44 –99.24 –89.38

Table 5 Performance comparison of H2 and H3 with multi-column generation when incorporated in the B&P heuristic. Problem no.

s

n

/H2_B&P

/H3_B&P

(oH2_B&P) Number of nodes

(oH3_B&P) Number of nodes

(jH2_B&P) Number of columns

(jH3_B&P) Number of columns

(tH2_B&P) Total runtime (s)

(tH3_B&P) Total runtime (s)

1 2 3 4 5 6

2 2 4 4 6 6

5 10 5 10 5 10

21,961 38,776 25,770 56,258 57,509 92,622

24,760 48,489 42,860 84,876 98,720 156,062

7 5 1981 853 1239 769

486 502 2520 1340 1890 1020

18 38 465 1238 443 738

420 445 2200 1180 1682 860

14 23 1315 >1800 >1800 >1800

53 77 >1800 >1800 >1800 >1800

and jH2_B&P. Also, the total runtime tH3_B&P is 43% larger than tH2_B&P, while solution /H3_B&P is 56% higher than /H2_B&P. The relative inefficiency of H3 when compared to H2 is due to the fact that the columns generated by H3 are not necessarily feasible with respect to customer demand over the planning horizon s. Although this is not essential for SPt, it proved to be beneficial. The advantage of H2 is that it takes the consolidation requirement into account with constraints (4c) so the columns provided by H2 at each iteration include all customers who must be serviced to ensure feasibility. Without these constraints the size of the B&B tree increases by an order of magnitude. Generally speaking, it takes a significant amount of branching effort to consolidate deliveries into the same column. In summary, H2 demonstrated the best performance when implemented with the B&P algorithm and will be used in the experiments discussed in the following section. 6.2. Solving the PIDRP A fine-tuning of the B&P algorithm determined the best rules to use for branching and the frequency with which to call our feasibility heuristic [see Bard and Nananukul (2009) for the details]. Although the preprocessing model (5) did not significantly improve the results from our initial testing, we felt that it could still have some impact on performance since its effectiveness depends on the actual data. In addition, its computational costs are low, at most 40 seconds for the largest instances, so it was included in the methodology. With this in mind, we now provide test results for a series of instances derived from the three data sets developed by Boudia et al. (2007). The original data sets contain 30 instances of 50, 100 and 200 customer problems, all with a 20-period planning horizon

C

and holding costs hP = 1, hi ¼ 0 for all i e N. Customers were randomly located on a 100  100 Euclidean grid. As expected, none of these instances could be solved optimally, at least within 24 h, with either CPLEX or our exact B&P algorithm. Therefore, we selectively created five data sets from the original 150 problems to test the limits of our methodology. The first data set was created from the 50 customer problems, the second and third from the 100 customer problems, and the fourth and fifth from the 200 customer problems. Each contains 20 instances. The first step was to fix the number of customers n e {10, 20, 30, 40, 50} and the number of time periods s e {2, 4, 6, 8}. By selecting a range of values from these sets, an instance had anywhere from 2 to 8 planning periods and 10–50 customers. For each customer i, demand was uniformly distributed between 0 and the storage capacity ICmax :i , and for each period t the probability of that demand being realized varied uniformly in the interval [0.5, 1]. The value of ICmax :i was a function of i and the number of customers in the specific data set. Because these new data sets have fewer customers than the original, the number of vehicles was proportionally adjusted downward as follows: h10 = 1, h20 = 2, h30 = 3, h40 = 4 and h50 = 5 for corresponding n e {10, 20, 30, 40, 50}. Vehicle capacities were Q50 = 8000, Q100 = 8000, Q200 = 8000 for the data sets with 50, 100 and 200 customers, respectively. The maximum storage capacity for customer i, ICmax :i , for data sets 1, 2 and 3 follow the default values from the original data sets. For data sets 4 and 5, ICmax :i is set to 1000 for all customers. In the computations, a maximum of 3600 s was allowed for a run. In a few cases, termination took place prematurely when the out-of-memory error appeared. The H2 heuristic was used for all but the largest instances; for data sets with either n = 40 and s P 6 or n P 50, model (4) was replaced with model (3) to deter-

722

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

Table 6 Results using CPLEX and tabu search to solve the PIDRP. Problem no.

s

n

(/cplex) CPLEX solution

(tcplex) CPLEX runtime (s)

(/TS) Tabu search solution

(tTS) Tabu search runtime (s)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8

10 10 10 10 20 20 20 20 30 30 30 30 40 40 40 40 50 50 50 50

80,333 95,815 132,545  88,548 142,847 207,587 255,352 97,814 159,946 282,244  102,549    138,234   

0.64 >3600 >2821a >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600 >3600

80,333 95,086 127,183 170,281 88,542 122,377 178,096 227,518 96,452 131,453 221,178 278,383 95,491 142,480 233,825 288,903 102,894 144,871 239,250 325,537

0.48 6.24 16.71 1.87 0.52 10.52 5.25 5.13 0.6 24.19 18.24 1.98 0.76 26.91 20.08 6.78 0.96 21.45 13.94 8.01

* Cplex cannot find integer solution. Out of memory.

a

Table 7 Results using the B&P algorithm to solve the PIDRP. Problem no.

s

n

(/B&P) B&P solution

(tB&P) Total runtime (s)

[(/B&P  /TS)//TS] 100%

[(/B&P  /cplex)//cplex] 100%

(oB&P) Number of nodes processed

(jB&P) Number of columns

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8

10 10 10 10 20 20 20 20 30 30 30 30 40 40 40 40 50 50 50 50

80,333 85,196 127,183 132,265 88,542 122,377 178,096 186,581 96,452 131,453 214,716 268,375 95,491 142,480 222,990a 288,903a 102,894a 144,871a 236,148a 296,062a

3.7 83.7 100.89 638.21 8.7 54.84 879.55 1017.38 8.17 288.24 >3600 2890 8.77 173.55 750.38 2934.28 9.94 189.27 3100 226.05

0 –10.40 0 –22.33 0 0 0 –17.99 0 0 –2.92 –3.62 0 0 –4.63 0 0 0 –1.30 –9.05

0 –11.08 –4.05  –0.01 –14.33 –14.21 –26.93 –1.39 –17.81 –23.93  –6.88    –25.57   

1 1 37 105 1 1 3 10 1 4 6 11 1 1 15 63 1 1 81 1

0 53 146 245 0 53 189 210 0 235 374 198 0 77 1214 3232 0 710 3173 257

* Cplex cannot find integer solution. a Heuristic H1 used in subproblem.

mine delivery quantities prior to calling the VRP code. This modification was necessary because the amount of time required to solve model (4) alone at the root node did not permit column generation to be completed within the allotted 3600 s. Table 6 reports the results for one of the five data sets for which CPLEX and our tabu search algorithm were used to solve the PIDRP (Bard & Nananukul, 2009). These statistics are provided for comparative purposes. The column headings should be self-explanatory. Table 7 reports the results for the B&P algorithm. The fourth and fifth columns give the B&P solution (/B&P) and corresponding runtime (tB&P) at termination. Columns 6 and 7 give the percentage gap between B&P and tabu search, and B&P and CPLEX, respectively. The eighth column identifies the total number of nodes in the search tree (oB&P), and the final column reports the total number of columns generated (jB&P). In all but one case, B&P terminated within the 1-h limit, in all but two cases, CPLEX reached the 1-h limit, and in all 20 cases tabu

search required at most 26.91 s. On average, B&P improved upon tabu search by 3.6% and by 12.2% when compared to CPLEX. The statistics for the remaining four data sets are similar and can be found in Bard and Nananukul (2008). 7. Discussion and future directions The use of heuristics to find solutions to large-scale combinatorial optimization problems is common, but little work has been done on coupling them with decomposition methods. In this paper, we have shown how the inventory routing component of the PIDRP can be solved efficiently within a branch-and-price framework. The results indicate that high quality solutions can be obtained for instances of the PIDRP with up to eight time periods and 50 costumers. One idea for reducing runtimes would be to replace the VRP component of the IRP with hTSPs. This could be achieved by first

J.F. Bard, N. Nananukul / Computers & Industrial Engineering 57 (2009) 713–723

clustering the customers into groups whose aggregate delivery requirements could be satisfied by a single vehicle over the planning horizon. This does not mean that the daily aggregate demand of each group would have to be within the capacity limits of the vehicle as long as the option of building inventory at the customer sites was taken into account. Another area of possible research involves the expansion of the PIDRP model and the development of new algorithms for accommodating uncertainty in demand. A simple case might consider a probability distribution over each customer’s demand with an objective aimed at minimizing expected long-run costs. A more complex case might be to postulate individual scenarios for daily demand and formulate the problem as a two-stage stochastic integer program, where production decisions were made in the first stage and distribution decisions were made in the second. These ideas are now under investigation. References Abdelmaguid, T. F., & Dessouky, M. M. (2006). A genetic algorithm approach to the integrated inventory–distribution problem. International Journal of Production Research, 44(21), 4445–4464. Anily, S., & Federgruen, A. (1990). One warehouse multiple retailer system with vehicle routing costs. Management Science, 36(1), 92–114. Anily, S., & Federgruen, A. (1993). Two-echelon distribution systems with vehicle routing costs and central inventories. Operations Research, 41(1), 37–47. Bard, J. F., Nananukul, N. (2008), A two-stage supply chain planning problem with inventory routing. Working paper, Graduate Program in Operations Research & Industrial Engineering, The University of Texas, Austin, TX, USA. Bard, J. F., Nananukul, N. (2009). The integrated production–inventory– distribution–routing problem for a single commodity. Journal of Scheduling. Bard, J. F., Huang, L., Jaillet, P., & Dror, M. (1998). A decomposition approach to the inventory routing problem with satellite facilities. Transportation Science, 32(2), 189–203. Bertazzi, L., Paletta, G., & Speranza, M. G. (2002). Deterministic order-up-to level policies in an inventory routing problem. Transportation Science, 36(1), 119–132. Boudia, M., Louly, M. A. O., Prins, C. (2006). A memetic algorithm with population management for a production–distribution problem. In A. Doglui, G. Morel, C. E. Pereira (Eds.), 12th IFAC symposium on information control problems in manufacturing, May 17–19, Saint-Etienne, France (Vol. 3, pp. 541–546). Boudia, M., Louly, M. A. O., & Prins, C. (2007). A reactive grasp and path relinking for a combined production–distribution problem. Computers & Operations Research, 34(11), 3402–3419.

723

Carlton, B. W., & Barnes, J. W. (1996). Solving the traveling salesman problem with time windows using tabu search. IIE Transactions on Scheduling & Logistics, 28(8), 617–629. Cetinkaya, S., Mutlu, F., & Lee, C.-Y. (2006). A comparison of outbound dispatch policies for integrated inventory and transportation decisions. European Journal of Operational Research, 171(3), 1094–1112. Choi, E., & Tcha, D.-W. (2007). A column generation approach to the heterogeneous fleet vehicle routing problem. Computers & Operations Research, 34(7), 2080–2095. Dror, M., & Ball, M. (1987). Inventory/routing: Reduction from an annual to a short period problem. Naval Research Logistics Quarterly, 34(4), 891–905. Gaudioso, M., & Paletta, G. (1992). A heuristic for the periodic vehicle routing problem. Transportation Science, 26(2), 86–92. Gaur, V., & Fisher, M. L. (2004). A periodic inventory routing problem at a supermarket chain. Operations Research, 52(6), 813–822. Golden, B., Assad, A., & Dahl, R. (1984). Analysis of a large scale vehicle routing problem with an inventory component. Large Scale Systems, 7(2–3), 181–190. Herer, Y. T., & Roundy, R. (1997). Heuristics for a one-warehouse multiretailer distribution problem with performance bounds. Operations Research, 45(1), 102–115. Herer, Y. T., Tzur, M., & Yucesan, E. (2006). The multilocation transshipment problem. IIE Transactions on Scheduling & Logistics, 38(3), 185–200. Kant, G., Jacks, M., & Aantjes, C. (2008). Coca-Cola enterprises optimizes vehicle routes for efficient product delivery. Interfaces, 38(1), 40–50. Lee, B. K., Kang, K. H., & Lee, Y. H. (2008). Decomposition heuristic to minimize total cost in a multi-level supply chain network. Computers & Industrial Engineering, 54(4), 945–959. Lei, L., Liu, S., Ruszczynski, A., & Park, S. (2006). On the integrated production, inventory, and distribution routing problem. IIE Transactions on Scheduling & Logistics, 38(11), 955–970. Mourgaya, M., & Vanderbeck, F. (2007). Column generation based heuristic for tactical planning in multi-period vehicle routing. European Journal of Operational Research, 183(3), 1028–1041. O’Brien, C., & Tang, O. (Eds.). (2006). Integrated enterprise and supply chain management. International Journal of Production Economics, 101(1) (special issue). Roundy, R. (1985). 98%-Effective integer-ratio lot-sizing for one-warehouse multiretailers systems. Management Science, 31(11), 1416–1430. Savelsbergh, M. (1997). A branch-and-price algorithm for the generalized assignment problem. Operations Research, 45(6), 831–841. Villegas, F. A., & Smith, N. R. (2006). Supply chain dynamics: Analysis of inventory vs. order oscillations trade-off. International Journal of Production Research, 44(6), 1037–1054. Wolsey, L. A. (1998). Integer programming. New York: Wiley. Zhao, Q.-H., Lai, S.-Y., & Wang, K. K. (2007). A partition approach to the inventory/routing problem. European Journal of Operational Research, 177(2), 786–802.