Otimização Multiobjetivo da Rede Integrada de Localização-Estoque-Distribuição em VMI via Algoritmo Evolucionário

Expert Systems with Applications 38 2011 67686776 An evolutionary approach for multiobjective optimization of the integrated locationinventory distribution network problem in vendormanaged inventory ShuHsien Liao a ChiaLin Hsieh b PengJen Lai c a Department of Management Sciences and Decision Making Tamkang University Taiwan ROC b Department of Industrial Management and Enterprise Information Aletheia University Taiwan ROC c Department of Mathematics National Kaohsiung Normal University Kaohsiung Taiwan ROC ARTICLE INFO Keywords Integrated locationinventory distribution network system Vendormanaged inventory VMI Supply chain management Multiobjective optimization Genetic algorithm ABSTRACT Vendormanaged inventory VMI is one of the emerging solutions for improving the supply chain efficiency It gives the supplier the responsibility to monitor and decide the inventory replenishments of their customers In this paper an integrated locationinventory distribution network problem which integrates the effects of facility location distribution and inventory issues is formulated under the VMI setup We presented a MultiObjective LocationInventory Problem MOLIP model and investigated the possibility of a multiobjective evolutionary algorithm based on the Nondominated Sorting Genetic Algorithm NSGA2 for solving MOLIP To assess the performance of our approach we conduct computational experiments with certain criteria The potential of the proposed approach is demonstrated by comparing to a wellknown multiobjective evolutionary algorithm Computational results have presented promise solutions for different sizes of problems and proved to be an innovative and efficient approach for many difficulttosolve problems 2010 Elsevier Ltd All rights reserved 1 Introduction Recently two generic strategies for supply chain design have emerged efficiency and responsiveness Efficiency aims to reduce operational costs responsiveness on the other hand is designed to react quickly to satisfy customer demands A crucial question in the supply chain is the design of distribution networks and the identification of facility locations Ballou and Masters 1993 put forward four strategic planning areas in the design of a distribution network system as shown in Fig 1 The first issue deals with customer service levels The second one deals the placement of facilities and demand assignments made to them The third deals with inventory decisions and policies that involve inventory control The fourth deals with transportation decisions of how transport modes are selected utilized and controlled All four of these areas are interrelated and the customer service level is determined by the other three decision areas There are practical challenges for firms when they try to simultaneously reduce operating costs for efficiency and customer service for responsiveness In traditional supply chain network design the optimization focus is often placed on minimizing cost and maximizing profit as a single objective However very few distribution network systems should be considered as intrinsically single objective problems It is not always desirable to reduce costs if this results in a degraded level of customer service Thus it is necessary to set up a multiobjective network design problem Research on integrated locationinventory distribution network systems is relatively new Jayaraman 1998 developed an integrated model which jointly examined the effects of facility location transportation modes and inventoryrelated issues However Jayaramans study did not contain any demand and capacity restrictions Erlebacher and Meller 2000 formulated an analytical joint locationinventory model with a highly nonlinear objective function to maintain acceptable service while minimizing operating inventory and transportation costs Nozick and Turnquist 2001 proposed a joint locationinventory model to consider both cost and service responsiveness tradeoffs based on an uncapacitated facility location problem Miranda and Garrido 2004 studied a MINLP model to incorporate inventory decisions into typical facility location models They solved the distribution network problem by incorporating a stochastic demand and risk pooling phenomenon Sabri and Beamon 2000 presented an integrated multiobjective multiproduct multiechelon model that simultaneously addresses strategic and operational planning decisions by developing an integrated two submodule model which includes cost fill rates and flexibility Gaur and Ravindran 2006 studied a bicriteria optimization model to represent the inventory aggregation problem under risk pooling finding out the tradeoffs in costs and responsiveness Recently Daskin Coullard and Shen 2002 and Shen Coullard and Daskin 2003 introduced a joint locationinventory model with risk pooling LMRP that incorporates inventory costs at dis tribution centres DCs into location problems LMRP solved the problem in two special cases deterministic demand and Poisson demand It assumed direct shipments from DCs to buyers which extended the uncapacitated fixedcharge problems to incorporate inventory decisions at the DCs The uncapacitated assumption at DCs is usually not the case in practice Shu Teo and Shen 2005 solved LMRP with general stochastic demand Shen and Daskin 2005 extended the LMRP model to include the customer service component and proposed a nonlinear multiobjective model including both cost and service objectives In contrast to LMRP and its variants that consider inventory cost only at the DC level Teo and Shu 2004 and Romeijn Shu and Teo 2007 proposed a warehouseretailer network design problem in which both DCs and retailers carried inventory These are actually the two major streams of integrated distribution network design problems Our model builds upon the initial LMRP model but with some differences First a capacitated version of a similar model is estab lished Second to make an original contribution the proposed model incorporates two extra performance metrics corresponding to customer service With these considerations we present a capacitated MultiObjective LocationInventory Problem MOLIP which results in a MixedInteger NonLinear Programming MIN LP formulation Some noteworthy innovative research aspects that are incorporated in our research include i MultiObjective LocationInventory Problem Very few studies have addressed this problem ii multiobjective evolutionary algorithms MOEAs Most previous works have focused on traditional optimization tech niques but few have performed these techniques successfully and efficiently In contrast MOEAs have been successfully devel oped for various optimization problems creating potential for the proposed MOLIP This study is organized as follows Section 2 describes our re search problem and details the model formulation Section 3 pro poses a hybrid evolutionary algorithm with a heuristic procedure for MOLIP Section 4 illustrates our experimental results including i the computational results of a basecase problem ii scenario analysis iii computational evaluation of the proposed algorithm for MOLIP Finally conclusions and suggestions for the direction of future research are provided in Section 5 2 Designing an integrated locationinventory distribution network model In this section we present a mathematical model which pro vides the foundation for our research 21 Problem description 211 VMI coordination mechanism Vendormanaged inventory VMI is one of the most widely discussed coordination mechanisms for improving multifirm supply chain efficiency Evidence has shown that VMI can improve supply chain performance by decreasing inventory costs for the supplier and buyer and improving customer service levels such as reduced order cycle times and higher fill rates Waller Johnson Davis 1999 Fig 2 indicates the system diagram of a VMI system includes its incurred material and information flows Since the sup plier is responsible for managing the inventories at the buyers DC including ordering and inventory holding the supplier ought to re ceive the information about demand directly from the market Since the supplier determines ordering instead of receiving orders from buyers there is no information flow of the buyers orders in the VMI system The main feature of VMI indicates the centralized system with in with which the supplier as a sole decision maker decides the or der quantity based on information available from both buyers and suppliers to minimize the total cost of the whole supply chain sys tem The supplier has full authority over inventory management at the buyers DC to pay all costs associated with the suppliers pro duction cost both the buyers and the suppliers ordering cost the inventory holding cost and distribution cost The supplier mon itors manages and replenishes the inventory of the buyer Thus the decisions on order replenishment quantity and order shipping are given to the supplier in the VMI system rather than to the buyer as in tradition systems Fig 3 presents the operational cost structure between the partners in the VMI system The proposed model is mainly based this cost structure 212 Overview of our research problem In general suppliers and distributors route their products through DCs In practice there are many cases in which each sup plier has its own set of DCs Consider a distribution network config uration problem where a single supplier and DCs are to be established to distribute various products to a set of buyers and both the DCs and buyers are geographically dispersed in a region In this problem each buyer experiences demands for a variety of products which are provided by the supplier A set of DCs must be located in the distribution network from a list of potential sites The DCs act as intermediate facilities between the supplier and the buyers and facilitate the shipment of products between the two echelons The supplier wishes to decide the supply chain distribu tion network for its products such as to determine the subsets of DCs to be opened and to design a distribution network strategy that will satisfy all capacity and demand requirements for the products imposed by the buyers However our problem jointly considers both strategic and tacti cal decisions in the supply chain system The strategic decision in volves the location problem which determines the number and the locations of DCs and assigns buyers to DCs whereas the tactical decision deals with the inventory problem which determines the levels of safety stock inventory at DCs to provide certain service levels to buyers The integrated problem is called a locationinven tory distribution network problem The centralized inventory pol icy is considered under the vender managed inventory VMI mode Waller et al 1999 which refers to the holding safety stocks aggregated at DCs This inclusion acquires especial relevance in the presence of high inventory holding costs and high variability of de mands Fig 4 shows the overall schematic diagram of the hierarchy of the model considered in our study 22 Model assumptions and notations Basic assumptions are used when modeling our problem It is assumed that all the products are produced by a single supplier and one specific product for a buyer should be shipped from a sin gle DC Reverse flows intransit inventory and pipeline inventory are not considered All the buyers demands are uncertain and the Fig 1 Four strategic planning issues in distribution network design SH Liao et al Expert Systems with Applications 38 2011 67686776 6769 storage capacities of the supplier are unlimited but are capacitated at open DCs More assumptions will be stated when we illustrate the mathematical model Here the mathematical notations used in the model are described as follows Indices i is an index set for buyers ði 2 IÞ j is an index set of po tential DCs ðj 2 JÞ k is an index set for product classifications ðk 2 KÞ Decision variables Q k wj is the aggregate economic order quantity for DC j for product k shipped from the supplier Yj 1 if DC j is open 0 otherwise Xk ji ¼ 1 if DC j serves buyer i for shipping product k 0 otherwise Model parameters uj is the capacity volume of DC j dik is the average daily demands for product k of buyer i rik is the standard deviation of daily demands for product k of buyer i fk j is the aver age lead time daily for product k to be shipped to DC j from the supplier w is the number of days per year fj is the fixed annual facility operating cost of opening a DC j h k j is the aggregated inven tory holding unit cost per unit time annually at DC j for product k tck ji is the unit variable transportation cost for shipping product k from DC j to buyer i rck j is the unit variable production and trans portation cost for shipping product k from the supplier to DC j Dk wj is the expected annual demand for product k through DC j Dmax is the maximal covering distance that is buyers within this distance from an open DC are considered well satisfied 23 Mathematical model To begin modeling this problem let us assume for a moment that the assignment of buyers to a DC is known a priori and that all the products are produced by a single supplier We assume that the daily demand for product k at each buyer i is independent and normally distributed ie Ndik r2 ik Furthermore at any site of DC j we assume a continuous review inventory policy Qj rj to meet a stochastic demand pattern Also we consider that the supplier takes an average lead time fk j in days for shipping product k from the supplier to DC j so as to fulfill an order From basic inventory theory Eppen 1979 we know that if the demands at each buyer are uncorrelated then the aggregate demand for product k during lead time fk j at the DC j is normally distributed with a mean of fk j dk wj Fig 2 System diagram of VMI system Fig 3 Cost structure of VMI system 6770 SH Liao et al Expert Systems with Applications 38 2011 67686776 Inputs of Strategic Level Model Number location capacity of existing potential DCswarehouses Location of customersbuyers to be served Aggregated demands of products Strategic Level Supply Chain Network Design Model Outputs of Strategic Level Model Number location capacity of open plantssuppliers Number location capacity of open DCswarehouses Assignment of DCswarehouses to customersbuyers to be served Inputs to Tactical Level Model Customer Service Level Reorder Point for in DCswarehouses and Retailers Customer Demand Variability Replenishment Lead Times at DCswarehouses Tactical Level Inventory Management Distribution System Planning Model Outputs of Tactical Level Model Safety Stock Inventory Levels kept in DCswarehouses and Retailers Cycle Length of periodic review policy Economic Order Quantity Distribution Plans Fig 4 Overview of the integrated strategic and tactical planning model where dwjk sumi dik xjik and a variance of xizetaj sumi sigma2ik xjik Let us consider the centralized supply chain system under the vendor managed inventory VMI mode which refers to aggregating the safety stock pooled at different DCs Then the total amount of safety stock for product k at DC j with risk pooling is z1alpha sqrt sumjk xizetaj sigma2ik Xjik where 1alpha refers to the level of service for the system and z1alpha is the standard normal value with Pz z1alpha 1alpha In the proposed model the total cost is based on the cost structure of the VMI system in Fig 3 and is decomposed into the following items i facility cost which is the cost of setting up DCs ii transportation cost which is the cost of transporting products from the supplier to the buyers via specific DCs iii operating cost which is the cost of running DCs iv cycle stock cost which is the cost of ordering and maintaining inventory at DCs and v safety stock cost which is the cost of holding sufficient inventory at DCs in order to provide a specific service level to their buyers The total cost Z1 is represented as follows Z1 sumj fj Yj psi sumk sumj sumi rck tck dik xjik sumk sumj ok psi dwjk Qwjk sumk sumj hk Qwjk 2 Yjk sumk sumj hk z1alpha sqrt sumjk xizetaj sigma2ik xjik 1 Based on Z1 the optimal order quantity Qwjk for product k at each DC j can be obtained by differentiating Eq 1 in terms of Qwjk each DC j and each product k equal zero to minimize the total supply chain cost We can obtain Qwjk sqrt 2 ok dik hk Yjk for open DC j k hk In this case there is not any capacity constraint for the order quantities Qwjk since we assume the storage capacity at the supplier is unlimited Thus replacing Qwjk in the third and fourth terms of Z1 in Eq 1 we can obtain a nonlinear cost function of Z1 As follows we propose an innovative mathematical model for the MultiObjective LocationInventory Problem MOLIP Min Z1 sumk sumj fj ykj sumk sumk sumj Psijik Xjik sumk sumj rck tck dik xjik sumk sumj hk sqrt sumi dik xjik sumk sumj hk sqrt sumi Deltaj ik Xjik 2 Max Z2 sumk sumi dik Xjik sumk sumi dik 3 Max Z3 sumk sumi dik Xjik sumk sumi dik 4 st sumj Xjik 1 i I k K 5 Xjik Yj i I j J k K 6 sumk sumi dik xjik sumk Lambdajik Xjik uj Yj j J 7 Xjik 01 Yj 01 i I j J k K 8 where Psijik psi rck tck dik Gammajk sqrt 2 ok hk Dik psi dik Lambdajik z1alpha xizetaj sigma2ik Eqs 24 give the objectives While Eq 2 of Z1 is to minimize the total cost Eq 3 of Z2 and Eq 4 of Z3 give the objectives by referring to the maximization of customer service by two performance measurements i volume fill rate VFR defined as the satisfied fraction of total demands without shortage ii responsiveness level RL the percentage of fulfilled demand volume within specified coverage distance Dmax Eq 5 restricts a buyer to be served by a single DC if possible Eq 6 stipulates that buyers can only be assigned to open DCs Eq 7 indicates the maximal capacity restrictions on the opened DCs to enable the capability of holding sufficient inventory for every product that flows through the DC and also part of safety stock so as to maintain the specified service level Eq 8 determines binary constraints The proposed MOLIP model does not only determine the DC locations the assignment of buyers to DCs but also finds out endogenously both the optimal order quantities and safetystock levels at DCs Since two of the three objective functions Z1 and Z3 are nonlinear the formulation results in an intractable multiobjective MINLP model or neither dominates the other One solution x1 is said to dominate the other solution x2 if both the following conditions are true i the solution x1 is not worse than x2 for all objectives ii the solution x1 is strictly better than x2 in at least one objective If any of the above conditions are violated the solution x1 does not dominate x2 NSGAII Deb et al 2002 is one of the best techniques for generating Pareto frontiers in MOEAs First of all for each solution in the population one has to determine how many solutions dominate it and the set of solutions to which it dominates Then it ranks all solutions to form nondominated fronts according to a nondominated sorting process hence classifying the chromosomes into several fronts of nondominated solutions To allow for diversification NSGAII also estimates the solution density surrounding a particular solution in the population by computing a crowding distance operator During selection a crowdedcomparison operator considering both the nondominaton rank of an individual and its crowding distance is used to select the offspring without losing good solutions elitism Whereas crossover and mutation operators remain as usual A NSGAIIbased evolutionary approach for MOLIP is proposed as shown in Table 1 This algorithm starts by generating a random population P1 of size L For each chromosome in P1 the algorithm evaluates its cost and coverage using the encoded solution expressions Then the algorithm applies nondominated sorting of P1 and assigns each chromosome to the front to which it belongs Next the algorithm applies binary tournament selection to form the crossover pool crossover and mutation operators to generate the children population C1 of size L Once initialized the main body of the algorithm repeats for T generations The algorithm applies nondominated sorting to Rt resulting in a population from the union of parents Pt and children Ct The algorithm obtains the next generation population Pt1 after selecting the L chromosomes from the first fronts of Rt Next it applies binary tournament selection crossover and mutation operators to generate the children Ct1 Too make it easy to understand the algorithm is graphically represented in Fig 5 During the selection process of the next generation chromosome fitness depends on the evaluation of the decoded solution in the objective functions and its comparison with other chromosomes The nondomination sorting updates a tentative set of Pareto optimal solutions by ranking a population according to nondomination After that each individual p in the population is given two attributes i nondomination rank in the optimization objectives prank ii local crowding distance in the objectives space directions pdistance If both chromosomes are at the same nondomination rank the one with fewer chromosomes around in the front Table 1 NAGAIIbased evolutionary approach 1 Randomly generate P1 2 Evaluate P1 3 Nondominated sort P1 4 Generate C1 form P1 apply binary tournament selection crossover and mutation 5 Evaluate C1 6 while t T do 7 Rt Pt U Ct 8 Nondominated sort Rt 9 Sort Rt using i see Definition 1 10 Select Pt1 from the first L chromosome of Rt 11 Generate Ct1 from Pt1 apply binary tournament selection crossover and mutation 12 Mutate Ct1 13 Evaluate Ct1 14 t t 1 15 end while Fig 5 Graphical representation of the NSGAII algorithm is preferred Thus a partial order n defined in Definition 1 is used to decide which of the two chromosomes is fitter Definition 1 states that a higher nondomination level chromosome is always preferred If chromosomes are at the same level the one with fewer chromosomes around the front is preferred Definition 1 Let p q Rt be chromosomes in population Rt p is said to be fitter than q p n q either if prank qrank or prank qrank and pdistance qdistance 33 Solution encoding Each solution of the MOLIP is encoded in a binary string of length m J where the jth position bit indicates if DC j is open value of 1 or closed value of 0 This binary encoding only considers if a given DC j is open or closed variables Yj A MOLIP solution also involves the assignment of buyers to open DCs variables Xkji This assignment is performed by a procedure that tries to minimize cost without deteriorating coverage and capacity Since the capacity constraints in MOLIP limit the amount of buyers demands that can be assigned to candidate DCs a greedy heuristic is used to fulfill the buyerDC assignments First of all the buyers are sorted in the descending order of their demand flows After that according to the sorted order they are assigned to a specific DC according to the following rules Rule 1 For each buyer i if the buyer i is covered ie there are DCs within a distance of Dmax it is assigned to the DC with sufficient capacity if one exists which can serve it with the minimal difference between the remaining capacity of an open DC j and the demand flow of the buyer i through DC j That is the DC assignment attempts to be as full as possible Rule 2 If the buyer i cannot be covered or there is no successful assignment from the coverage set it is then assigned to the candidate DC with sufficient capacity that increases the total cost by the least amount regardless of its distance to the DC if possible 4 Model applications and experimental results 41 A baseline problem and its computational results There are no MOLIP instances in the public domain nor are any available in previous studies to serve for benchmarking For this reason a baseline problem was developed by taking the size of a Gammacom supply chain network with 15 DCs and 50 buyers as reference The potential DC locations are randomly generated within a square of 100 distance units of width Other model parameters are given in Table 2 For the sake of simplicity Euclidean distance is used for measuring distribution distances The company intended to determine the number of open DCs needed for order assignments However the capacity limitation of DCs affects the assignments of buyers The managers also need to evaluate tradeoffs among three criteria total cost TC volume fill rate VFR and responsiveness level RL To obtain the approximate Pareto front we attempted to solve the specified problem using the proposed hybrid evolutionary approach Through the GA approach the baseline model of DCs 15 of buyers 50 with product number k 2 resulted in 765 binary variables and 815 constraints In addition defining a reference point is the first step in allowing the MOLIP to obtain tradeoff solutions The reference point is a vector formed by the singleobjective optimal solutions and is the best possible solution that may be obtained for a multiobjective problem With a given reference point the MOLIP problem can then be solved by locating the alternatives or decisions which have the minimum distance to the reference point Thus the problem becomes how to measure the distance to the reference point For the MOLIP problem the decision maker is asked to determine weights by prior knowledge of objectives once all the alternatives in the Pareto front are generated Moreover the reference point can be found simply by optimizing one of the original objectives at a time subjective to all constraints Due to the incommensurability among objectives we measure this distance by using normalized Euclidean distance between two points in kdimensional vector space d tk1wtft ftft212 where f is an alternative solution in the Pareto front f is the reference point and wt is the relative weight for the tth objective Then all alternatives are ranked based on the value of d in descending order The highest ranked alternative with the minimal value of d is then considered as the optimal solution among alternatives for the given MOLIP problem Table 2 Model parameters for the baseline problem Parameters Value Annual cost of operating a DC j U900 1000 Annual holding cost at DC j for product k U2 4 Unit ordering cost at DC j for product k U8 10 Capacity of DC j U500 700 Unit variable transportation cost 1 Unit production and shipping cost for product k from the supplier to DC j U1 3 Maximal covering distance 25 Km Lead time daily U2 4 Working days per year 260 Average daily demand for product k at buyer i U60 80 Standard deviation of daily demands U2 4 Standard normal value service level 095 196 The hybrid evolutionary approach is used for the baseline mod el The input parameters are population size 100 generation number 200 cloning 20 crossover rate 80 mutation rate varies from 5 to 10 The approach program was coded in MAT LAB The algorithm allows the decision maker to rapidly find a set of Pareto solutions that are large in number The decision maker requires the determination of weights by prior knowledge of objec tives to generate the userdefined optimal solution As shown in Fig 6 we illustrate a userdefined optimal solution among the alternatives with equal weights of three objectives ie w1 w2 w3 13 in the baseline model with the minimal total cost of 251112536 the maximal volume fill rate of 7197 and the maximal customer responsive level of 6215 respectively where nine out of 15 candidate DCs are required to open and aggregate It is worth mentioning that most of these aggregated DCs were assigned to the buyers as close as to them as possible within the maximal coverage within 25 kms However about 2903 of buyers were unassigned ðÞ revealing the percentage of the uncovered demands which could possibly result in sales losses Also 3785 of aggregated buyers were assigned to DCs farther than the coverage distance Fig 7 shows the approximate Pareto front of the baseline prob lem obtained from the NSGAIIbased evolutionary algorithm To make it easy to understand the existing tradeoff between the cost and volume fill rate and responsive level respectively we present it as a percentage of the minimal cost instead of using it in an abso lute term As shown in Fig 2 it is possible to increase volume fill rate VFR from 3187 to 6982 and responsiveness level RL from 2539 to 6361 if the percentage over the minimal cost in creases from 204 to 40055 about two times when the number of open DCs is increased from four to nine Thus if the decision ma kers goal is to maintain volume fill rate VFR at a level of about 70 compared to the current status of 3187 extra costs are nec essary increase open DCs up to nine The increase in DCs enhances customers volume fill rate and also increase responsiveness level at the same time 42 Performance evaluation Our goal here is to evaluate the efficiency and the effectiveness of the proposed NSGAIIbased evolutionary approach We establish a set of random instances and try to keep almost all model param eters the same as the basecase problem We generate problem in stances of different sizes of DCs and buyers in the distribution network In addition various capacity and facilitycost scenarios are considered again In this experiment we generated four sets of problem instances SET 1 to SET 4 representing different sizes of problem instances ranging from 15 DCs and 50 buyers to 100 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 Closed DC Unfilled Unassigned buyer lost sales Open DC Assigned Buyer Fig 6 Graphical display of the userdefined optimal solution Fig 7 Approximate Pareto front for the userdefined optimal solution 6774 SH Liao et al Expert Systems with Applications 38 2011 67686776 DCs and 500 buyers problem sizes mn 1550 SET 1 50150 SET 2 75300 SET 3 and 100500 SET 4 Similar to the base case problem all these instances are randomly generated and uni formly distributed to locations within a square of 100 distance units of width for the coordinates of all DCs and buyers However there are two different types of facilitycost structures F1 to F2 Instances labeled F1 and F2 represent different types of facility cost problems There are also problem instances with three different types of DC capacity scenarios C1 to C3 Instances labeled C1 C2 and C3 stand for different DC capacity structures corresponding to tight normal and excess capacity scenarios After combining all the possibilities of problem sizes four types facility cost structure two types and DC capacity structure three types we end up with 24 problem instances Each problem instance is given a name in the following format AmnF1 to F2C1 to C3 For instance the problem instance A50150F2C1 represents a problem in which there are 50 DCs with 150 buyers which are both uniformly distributed within the square area width of 100 distance units The facility cost is expensive and its DC capacity structure is rather tight as compared to others SPEA2 Zitzler Thiele 1999 and NSGAII Deb et al 2002 have been considered two of the most successful and standard evo lutionary approaches among studies on MOEAs To verify the effi cacy of the algorithm we try to make comparisons between NSGAII and SPEA2 in the MOLIP model by using the randomly generated 24 problem instances mentioned above Ten independent runs of each problem instance were conducted for each algorithm The fi nal computational results of each problem instance are obtained by aggregating the approximate Pareto solutions of the 10 indepen dent runs Table 3 summarizes the performance results of the two algorithms considered For each algorithm we report on the identification of the instances the standardized dominatedspace metric w similar to that proposed by Medaglia and Fang 2003 and the number of solutions PF found in the approximate Par eto front respectively Note that the dimensionless metric w ranges between 0 and 1 and as it approaches 1 the closer the approximate and true Pareto front are to each other In addition we represent PF in w a b format where w a and b indicate the worst the average and the best values in the 10 runs of the algorithms Finally the last two columns indicate the incremental percentage of the average execution time T of the 10 independent runs in seconds and the standardized dominatedspace metric w to find out their relative differences between SPEA2 and NSGAII From Table 3 we conclude that for small instances in SET 1 m 15 n 50 NSGAII obtained better results in w Nonetheless SPEA2 runs faster than NSGAII This efficiency in terms of execu tion time is due to the fact that SPEA2 compares the current solu tion with the archive ie one with many as opposed NSGAII which compares many solutions with the current Pareto frontier ie many with many However for larger instances from SET 2 m 50 n 150 to SET 4 m 100 n 500 NSGAII almost always outperforms SPEA2 That is the NSGAII obtained better results than SPEA2 in all w metrics and almost all computing times T ex cept for those instances labeled C1 corresponding to tight capacity scenarios In addition there are significant differences in the qual ity of the solutions between the two approaches in the so called difficulttosolve instances in SET 4 Although NSGAII spends slightly more computation time than SPEA2 for those instances labeled C1 it still outperforms SPEA2 greatly in obtaining better quality of the approximate Pareto fron tiers for larger w For example in the A100500F1C1 instance NSGAII runs slightly slower with a difference 11 compared SPEA2 but favors in solution quality with the value 207 The two columns of PF in Table 3 indicate how many solutions are contained on the approximate Pareto front The results in PF ver ify that NSGAII always provides more Pareto solutions and main tains better diversity properties than SPEA2 Furthermore NSGAII is more stable and robust in computation For all instances SPEA2 is highly variable in PF as compared to NSGAII That is the gaps between the worst value w and the best value b of all experiments obtained by SPEA2 are much larger than for NSGAII Thus we may conclude that NSGAII is a reliable method that provides more ro bust approximate Pareto solutions Fig 8 illustrates the approxi mate Pareto frontiers obtained by the NSGAII and the SPEA2 based algorithms for the problem instance A100500F1U1 For ease of understanding of the existing tradeoffs among the three objectives including total cost TC volume fill rate VFR and responsiveness level RL we normalize total cost instead of using Table 3 Comparisons between NSGAII and SPEAIIbased approaches SH Liao et al Expert Systems with Applications 38 2011 67686776 6775 it in absolute terms Visually the tradeoff curves of these two ap proaches are very similar and partially overlapped However NSGAII results in the solutions covering a larger surface of the approximate Pareto solutions 5 Concluding remarks and research directions This study presented a MOLIP model initially represented as an integrated MOLIP formulation which examines the effects of facil ity location distribution and inventory issues under a vendor managed inventory VMI coordination mechanism The MOLIP model is solved with a proposed hybrid evolutionary algorithm which is preliminarily based on a wellknown NAGAII evolution ary algorithm with an elitism strategy and a nondominated sort ing mechanism We implemented two experiments First we investigated the possibility of a NSGAIIbased evolutionary algorithm solving the MOLIP model Computational results revealed that the hybrid ap proach performed well and presented promising solutions for the MOLIP model in solving practicalsize problems Second we com pared our approach with SPEA2 to understand the efficiency among two approaches The experiment indicates that two algo rithms obtained similar approximations of the Pareto front but our approach outperformed SPEA2 in terms of the diversity quality of the approximation of the Pareto front Moreover SPEA2 was only efficient in terms of execution time in small or tight capacity instances This indicates that the proposed approach could be an efficient approach for providing feasible and satisfactory solutions to largescale difficulttosolve problems In future works we intend to adapt the proposed hybrid evolu tionary algorithm to other location inventory and distribution sys tems that have different characteristics or network structures For instance a network system may have stockpiles or inventories within the suppliers and the customer sites and the shortage pen alty needs to be considered in the overall supply chain operating cost In addition the inclusion of other inventory decisions would be a direction worth pursuing Such inventory decisions could in clude frequency and size of the shipments from plants to the DCs and from DCs to the retailers based on different replenishment pol icies and lead time in addition to safetystock inventory in the model Finding ways to adapt our hybrid evolutionary algorithm into such systems is the task of future research Other possible research directions are to explore more compet itive MOEAs or other existing optimization technologies such as Lagrangian relaxation particle swarm optimization ant colony optimization or other soft intelligent computing techniques Com parative studies of these techniques are worth investigating in the future In addition some possible methods of hybridizations in clude the adaption of new genetic operators for integrated systems and the incorporation of other heuristic search techniques into the evolutionary algorithms such as hillclimbing or local repair procedure References Ballou R H Masters J M 1993 Commercial software for locating warehoused and other facilities Journal of Business Logistics 142 70107 Coello Coello C A 1999 A comprehensive survey of evolutionarybased multiobjective optimization techniques Knowledge and Information Systems 13 269308 Daskin M S Coullard C R Shen Z M 2002 An inventorylocation model formulation solution algorithm and computational results Annals of Operations Research 11014 83106 Deb K Pratap A Agarwal S Meyarivan T 2002 A fast and elitist multi objective genetic algorithm NSGAII IEEE Transactions on Evolutionary Computation 62 181197 Eppen G 1979 Effects of centralization on expected costs in a multilocation newsboy problem Management Science 255 498501 Erlebacher S J Meller R D 2000 The interaction of location and inventory in designing distribution systems IIE Transactions 322 155166 Fonseca CM Fleming PJ 1993 Genetic algorithms for multiobjective optimization formulation discussion and generalization In Proceedings of the Fifth International Conference on Genetic Algorithms pp 416423 San Mateo CA Gaur S Ravindran A R 2006 A bicriteria model for the inventory aggregation problem under risk pooling Computers and Industrial Engineering 513 482501 Jayaraman V 1998 Transportation facility location and inventory issues in distribution network design an investigation International Journal of Operations and Production Management 185 47 Knowles J D Corne D W 2000 Approximating the nondominated front using the Pareto archived evolution strategy Evolutionary Computation 82 149172 Medaglia A L Fang S C 2003 A geneticbased framework for solving multi criteria weighted matching problems European Journal of Operational Research 1491 77101 Michalewicz Z 1996 Genetic algorithms data structures evolution programs Berlin SpringerVerlag Miranda P A Garrido R A 2004 Incorporating inventory control decisions into a strategic distribution network model with stochastic demand Transportation Research Part E Logistics and Transportation Review 403 183207 Nozick L K Turnquist M A 2001 A twoechelon allocation and distribution center location analysis Transportation Research Part E Logistics and Transportation Review 375 425441 Romeijn H E Shu J Teo C P 2007 Designing twoechelon supply networks European Journal of Operational Research 1782 449462 Sabri E H Beamon B M 2000 A multiobjective approach to simultaneous strategic and operational planning in supply chain design Omega 285 581598 Schaffer J D 1985 Multiple objective optimization with vector evaluated genetic algorithms In the 1st International Conference on Genetic Algorithms pp 93 100 Hillsdale NJ Shen Z J Coullard C R Daskin M S 2003 A joint locationinventory model Transportation Science 371 4055 Shen Z M Daskin M S 2005 Tradeoffs between customer service and cost in integrated supply chain design Manufacturing and Service Operations Management 73 188207 Shu J Teo C P Shen Z M 2005 Stochastic transportationinventory network design problem Operations Research 531 4860 Srinivas N Deb K 1994 Multiobjective optimization using nondominated sorting in genetic algorithms Evolutionary Computation 23 221248 Teo C P Shu J 2004 Warehouseretailer network design problem Operations Research 522 396408 Waller M Johnson M E Davis T 1999 Vendormanaged inventory in the retail supply chain Journal of Business Logistics 201 183203 Zitzler E Thiele L 1999 Multiobjective evolutionary algorithms A comparative case study and the strength Pareto approach IEEE Transactions Evolutionary Computing 34 257271 0 02 04 06 08 1 0 05 1 01 02 03 04 05 06 07 08 TC VFR RL SPEA2 NAGA2 Fig 8 Approximate tradeoff curves for problem instance A100500F1U1 6776 SH Liao et al Expert Systems with Applications 38 2011 67686776 Revista Produço On Line Universidade Federal de Santa Catarina wwwproducaoonlineinfbr ISSN 1676 1901 Vol 3 Num 3 Setembro de 2003 O Problema de Roteamento de Estoques um olhar sobre a literatura Fernando Leme Franco Arsenal de Marinha do Rio de Janeiro Doutorando em Engenharia de Produção COPPEUFRJ fernandofrancouolcombr Alberto Gabbay Canen COPPEUniversidade Federal do Rio de Janeiro agcanenpepufrjbr Data de Submissão Mai03 Data de Aprovação Ago03 O Problema de Roteamento de Estoques um olhar sobre a literatura Fernando Leme Franco Arsenal de Marinha do Rio de Janeiro Doutorando em Engenharia de Produção COPPEUFRJ fernandofrancouolcombr Alberto Gabbay Canen COPPEUniversidade Federal do Rio de Janeiro agcanenpepufrjbr Resumo A globalização e o aumento da competitividade nas empresas têm levado pesquisadores a desenvolver dos pontos de vista teórico e prático o problema de roteamento de estoques PRE Este artigo é uma revisão da produção do conhecimento na área do PRE Aborda primeiramente os conceitos envolvidos do PRE assim como as diferenças entre o problema em questão e do roteamento de veículos PRV Os artigos pesquisados são olhados de acordo com o horizonte de planejamento e a modelagem da demanda a partir da abordagem explicitada em cada um Finaliza sugerindo caminhos futuros neste campo de pesquisa Palavras chave Roteamento de Estoques Roteamento de Veículos Controle de Estoques Logística Revisão Bibliográfica The Inventory Routing Problem a glimpse of the literature Abstract Globalization and the increasing competitiveness among the organizations are driving researchers towards developing the Inventory Routing Problem IRP both from theoretical and practical perspectives This article seeks to review knowledge production in the IRP field Firstly outlines the concepts associated with the IRP and its differences as compared to the Vehicle Routing Problem VRP Literature is analyzed according to the planning horizon and demand modeling based on each ones explicated approach It concludes suggesting possible avenues for research in this field Key words Inventory Routing Vehicle Routing Inventory Control Logistics Literature Review Introdução O problema de roteamento de estoques referese ao gerenciamento integrado da distribuição física e do controle de estoques Controle de estoques e distribuição não são problemas novos e têm sido estudados exaustivamente de forma independente Um dos primeiros estudos de controle de estoques remonta o início do século XX Hopp Spearman 1996 O primeiro artigo sobre roteamento de veículos com vistas à distribuição de produtos parece ter sido atribuído a Dantzig e Hamser 1959 e a primazia com relação ao roteamento de estoques é creditada a Federgruen e Zapkin 1984 Além dos aspectos citados roteamento e estoques podese identificar outros a partir da revisão de literatura feita por Baitas et al 1998 sobre o problema dinâmico de roteamento de estoques desta forma incluindo o fator dinâmico como um terceiro aspecto Klewegt et al 1998 argumentam que o roteamento de estoques é o principal problema a ser resolvido visando à implementação da técnica de gerenciamento de estoques de venda Esta é definida como um sistema integrado em que distribuidores e clientes compartilham informações de demanda e estoques visando melhorar o desempenho da cadeia de suprimentos Achabal et al 2000 Fatores como o aumento da competitividade e a globalização têm levado pesquisadores a estudar o PRE e empresas a implementar seus conceitos ainda que parcialmente A implementação prática do gerenciamento de estoques de venda somente se mostrou viável com o grande avanço da tecnologia de informações no final do século XX com conseqüente redução de custos Klewegt et al 1999 O gerenciamento de estoques de venda é apontado por Campbell et al 1999 como um exemplo de valor agregado pela logística incluindo disponibilidade redução de tempos de ciclo entregas consistentes e facilidade de emissão de ordens de compra entre outros elementos reconhecidos como essenciais para a satisfação do cliente As diferenças entre os problemas de roteamento de estoques e de veículos podem ser analisadas sob diversos aspectos sendo que o roteamento de estoques pode ser interpretado como um enriquecimento do roteamento de veículos incluindo nele o fator controle de estoques Ball 1988 O roteamento de veículos referese à distribuição de bens e serviços para diferentes clientes através de um depósito central Canen e Pizzolato 1993 dentre outros A solução deste tipo de problema visa à redução de custos através da definição de rotas em função das ordens de entrega recebidas para um período finito de tempo O intervalo de tempo finito e a necessidade de atender integralmente à quantidade definida em cada ordem de entrega são fatores determinantes do tipo de modelagem utilizada bem como da solução encontrada No problema de roteamento de estoques a quantidade a ser entregue é definida em função da demanda real dos clientes e não da demanda agregada através de ordens de compra A quantidade a ser entregue a cada cliente bem como a sua freqüência passam a ser função desta demanda O objetivo é a redução de custos para um intervalo de tempo infinito A principal restrição é atender integralmente à demanda dos clientes Em alguns casos esta restrição é relaxada com a introdução de custos por demanda não atendida 2 Consideramos dois fatores como fundamentais para definir a modelagem e solução do problema o horizonte de análise e a modelagem da demanda O horizonte de análise se refere ao período de tempo que será considerado para efeito de otimização A modelagem para período de tempo infinito proposta como base do problema é muitas vezes relaxada e o problema é resolvido para períodos finitos de tempo A utilização de um horizonte de tempo finito não descaracteriza o problema além de estar mais próxima de aplicações práticas A modelagem de demanda nos trabalhos considerados vai desde uma demanda constante até uma modelagem por distribuições probabilísticas conhecidas sendo que alguns trabalhos utilizam dados reais para validar a metodologia Campbell et al 1997 Campbell et al 1999 apresentam a solução de um problema prático de roteamento de estoques aplicado a uma distribuidora de gás composta de cerca de 60 fabricas e 10000 consumidores O problema é resolvido em duas fases sendo que a primeira utiliza um modelo de programação linear e a segunda um algoritmo de planejamento de entregas Ao final os referidos autores mostram que o procedimento proposto é melhor do que o procedimento usado pela referida empresa em diversos quesitos de desempenho A partir da importância do roteamento de estoques explicitado o presente artigo é uma revisão da literatura na área objetivando oferecer referências no sentido aplicativo e não em detalhes algorítmicos de modo que pesquisadores interessados possam a vir buscar acesso aos originais das mesmas A Produção do Conhecimento na Área A Figura 1 mostra a distribuição de trabalhos analisados neste artigo considerandose os fatores horizonte de análise e conhecimento da demanda Para um horizonte de análise correspondente a um período o problema de roteamento de estoques se aproxima do problema de roteamento de veículos A diferença básica é que o segundo considera além do custo de transporte os custos de armazenamento e de demanda não satisfeita Federgruen e Zapkin 1984 Quanto maior o horizonte considerado maior a dificuldade de modelagem De forma semelhante quanto menor o conhecimento da demanda maior a dificuldade de solução do problema Horizonte Infinito Vários períodos Um Período Determinística Probabilística Demanda Figura 1 Artigos analisados com relação ao horizonte de análise e modelagem da demanda 3 Para demandas consideradas como determinísticas soluções utilizando a formulação de lote econômico podem ser adaptadas facilitando a solução do problema Anily e Federgruen 1990 Desta forma na figura 1 quanto mais a direita e acima maior será a dificuldade de modelagem e solução do problema Vários trabalhos analisados concentramse na área de demanda probabilística e horizonte de análise infinito no entanto adotam restrições em relação ao problema proposto Algumas destas restrições visam atender características específicas do caso analisado como entregas diretas em BarnesSchuster e Bassok 1997 ou divisão em subproblemas em Kleywegt et al 1999 Artigo Ano Autor Abordagem 2 1990 Anily e Federgruen O algoritmo particiona os clientes por regiao considerando carregamento completo com análise de lote econômico para cada regiao 4 1997 BarnesSchuster e Bassok O algoritmo analisa o impacto da variação de diversos parâmetros nos custos definindo os melhores valores para serem usados 9 1987 Dror e Ball Utiliza algoritmos do caixaeiro viajante incluindo penalidades em função dos estoques de segurança e datas de entregas 10 1996 Dror e Trudeau Algoritmo baseado em fluxo de caixa com ênfase na maximização da eficiência operacional 11 1984 Federgruen e Zapkin O algoritmo mostra primeiramente um particionamento de clientes para resolvêlo como uma extensão do roteamento de veículos 13 1984 Golden Assad e Dahl O algoritmo calcula a urgência de cada cliente roteando em função desta urgência 14 1997 Herer e Levy O algoritmo roteia os clientes minimizando os desvios em relação ao melhor período de entregas calculado previamente 17 1999 Kleywegt et al Soluciona o problema através de cadeias de Markov particionando os clientes Esta partição considera entregas diretas com a restrição de número limitado de veículos 18 1993 Minkoff O algoritmo decompõe o problema gerando uma função de penalidade Estas penalidades são incorporadas ao problema principal que é resolvido através de cadeias de Markov 19 1992 Trudeau e Dror O algoritmo utiliza os dados finais de estoque de cada cliente para efetuar o roteamento do período seguinte sendo então utilizado período a período Tabela 1 Artigos apresentado na Figura 1 de acordo com sua abordagem A Figura 1 e a Tabela 1 acima mostram os artigos cujos algoritmos serão comentados na seqüência deste trabalho A numeração segue a mesma seqüência apresentada nas referências bibliográficas A totalidade dos artigos listados se concentra na solução através de heurísticas Assim estas são utilizadas mesmo por aqueles que usam programação linear em alguma fase do desenvolvimento É importante observar que se trata de trabalhos teóricos sem oferecer aplicações práticas Alguns trabalhos citados na bibliografia não aparecem na Figura 1 pois são revisões bibliográficas Baita et al1998 apresentam métodos de solução utilizando diversos algoritmos Campbell et al 1997 Campbell et al 1999 ou não são específicos do PRE Achabal et al 2000 Canen e Pizzolato 1993 Gallego e SimchiLevi 1990 4 O problema de roteamento de estoque é um problema de difícil solução que na melhor das hipóteses poderá ser resolvido através da adaptação de um problema de roteamento de veículos que é NPHard Baita et al1998 Algumas das soluções propostas resolvem o problema apenas por um curto intervalo de tempo No primeiro trabalho publicado Federgruen e Zipkin 1984 este intervalo de tempo era de apenas um período sendo que posteriormente foi estendido para vários períodos Dror e Ball 1987 Jailet et al 1997 As diferentes abordagens encontradas na literatura em geral referemse à modelagem do efeito das decisões de curto prazo na performance a longo prazo Dror e Ball 1987 A abordagem de curto prazo tende a deixar para períodos posteriores o maior número possível de entregas o que induz um aumento de custos nesses períodos O balanceamento entre a otimização de curto e longo prazo é o fator chave na solução do problema de roteamento de estoques Federgruen e Zipkin 1984 modelam o problema de roteamento de estoques como uma extensão do problema de roteamento de veículos considerando além dos custos de transporte o custo de armazenamento e custo devido à demanda não atendida O problema é modelado para o horizonte de análise de um período sendo que as informações de estoque no final do período para cada cliente são usadas para o cálculo do roteamento de veículos no período seguinte A solução do problema é obtida através de heurística utilizando uma decomposição em subproblemas de alocação e roteamento coordenada apropriadamente Como nem todos os clientes receberão entregas em um determinado período utilizase uma rota fantasma para esses clientes Golden et al 1984 também analisam o problema da minimização de custo para o horizonte de análise de um período no entanto introduzem a restrição de manter um estoque adequado para todos os clientes ao final do período A solução é obtida através de uma heurística que calcula a urgência de cada cliente Desta forma os clientes são primeiramente selecionados em função de sua urgência e o roteamento é feito a seguir Dror e Ball 1987 analisam o problema para um horizonte de análise de vários períodos através de um procedimento no qual reduzem o problema para o horizonte de um período Este procedimento avalia os custos futuros das decisões presentes A política básica estabelece que todos os clientes cujos estoques chegarem a zero no período considerado devem ser reabastecidos Os custos para desvios em relação a esta política são calculados custos de não ser abastecido caso precise e custo de ser reabastecido caso não precise A solução do problema é obtida pela minimização do custo do desvio através de um procedimento heurístico Trudeau e Dror 1992 expandem essas idéias usando uma análise similar consideram o planejamento de alguns períodos No caso específico deste artigo eles utilizam o período de uma semana As informações no final de cada semana são utilizadas como dados de entrada para o planejamento da semana seguinte Herer e Levy 1997 utilizam o conceito de distância temporal ao invés de localização geográfica para definir as rotas No estudo incorporam os custos fixos de armazenamento e relativos à demanda não atendida A solução para um período fixo é obtida usando uma adaptação do algoritmo de Clark e Wright 1964 para roteamento de veículos Nesta 5 modelagem os autores consideram que os custos do capital imobilizado em estoques nos clientes serão de responsabilidade do distribuidor Desta forma o cliente paga somente pelos produtos que efetivamente utilizar Anily e Federgruen 1990 estudam a minimização do custo de transporte e armazenamento para um horizonte longo São estabelecidos padrões para efetuar uma partição dos clientes Após esta partição os clientes são divididos em regiões cuja demanda deve ser igual a uma carga completa do veículo de entrega A estratégia é visitar todos os clientes em uma região caso algum cliente desta região receba uma entrega Usando uma idéia similar Gallego e SimchiLevi 1990 analisam o efeito a longo prazo de entregas diretas A política de entregas diretas considera que a cada viagem o veículo irá visitar um único cliente A conclusão é que entregas diretas são eficientes em 94 dos casos desde que os lotes econômicos dos clientes sejam no mínimo igual a 71 da carga máxima dos veículos de entrega Kleywegt et al 1999 também utilizam o conceito de entregas diretas para formular o problema através de cadeias de Markov O estudo mostra que soluções exatas são obtidas em tempo razoável apenas quando são considerados poucos clientes Para solucionar o problema para instâncias maiores os clientes são separados em grupos e o problema é resolvido para os subproblemas resultantes Soluções utilizando cadeias de Markov foram propostas por Minkoff 1993 que considera a existência de um número ilimitado de veículos Os problemas causados por um número muito grande de estados foram resolvidos via decomposição A heurística incorpora um custo individual de transporte a cada cliente e resolve os diversos subproblemas individualmente Dror e Trudeau 1996 utilizam o conceito de fluxo de caixa para formular um problema para horizonte longo de análise Implicitamente assumese que os clientes pagam o valor relativo à quantidade recebida a cada entrega Desta forma quanto maior for a frequência de entregas mais constante será o fluxo de caixa Esta vantagem deverá ser confrontada com a elevação de custos decorrentes do aumento do número de entregas São considerados tantos os casos de demanda determinística como probabilística BarnesSchuster e Bassok 1997 consideram o problema do ponto de vista de horizonte de análise infinito e demanda estocástica O depósito é considerado apenas um agregador de dados não mantém estoques e são consideradas apenas entregas diretas A solução deste problema estabelece limites inferiores de entrega com metas de otimização Caminhos Futuros Os artigos analisados neste estudo referemse a modelos do problema de roteamento de estoque vistos do ponto de vista matemático Entretanto a implementação prática desses modelos depende de outras considerações Uma análise de cadeias logísticas mostra que a administração das empresas está preparada para enxergar apenas os níveis imediatamente anterior e posterior da cadeia Em outras palavras a empresa tende a enxergar apenas seus 6 fornecedores e seus clientes não tendo vínculos com os fornecedores de seus fornecedores ou com os clientes de seus clientes a não ser que partes da cadeia sejam operadas pela mesma empresa A implementação de conceitos de interligação de atividades logísticas fará com que o relacionamento dos agentes dentro da cadeia seja cada vez mais complexo Desta forma o vínculo entre a empresa e o fornecedor de seus fornecedores tem que ser considerada do ponto de vista administrativo e jurídico Estudos deste tipo são bastante apropriados no momento em que diversas empresas estão em vias de implementar o gerenciamento integrado de elementos da cadeia logística do tipo do roteamento de estoques apresentado neste estudo Conclusões Muito embora existam diversos artigos que apresentam revisões de literatura abordando o problema de roteamento de estoques conforme mencionado neste trabalho este artigo é relevante no sentido de promover uma revisão da produção atual do conhecimento na área focalizando especificamente os aspectos de horizonte de análise e conhecimento da demanda Estes conforme argumentamos são os principais aspectos que diferenciam o problema de roteamento de estoques do roteamento de veículos Mais ainda eles incorporam dificuldades crescentes à modelagem e à solução do problema A literatura pesquisada mostra que o tema roteamento de estoques tem sido estudado por mais de quinze anos porém aplicações práticas não são o foco das mesmas Esta literatura mostra também que existe uma razoável gama de métodos e algoritmos propostos desde adaptações das técnicas de roteamento de veículos e controle de estoques até propostas bastante inovadoras utilizando fluxo de caixa É difícil imaginar modelos abrangentes que possam ser aplicados a uma grande variedade de casos O que parece mais viável é a utilização de modelos que otimizem casos específicos como os apresentados neste artigo Observamos também que a pesquisa com modelos que integram elementos logísticos até então estudados de forma isolada é uma tendência na literatura em função de que hoje já existe a tecnologia necessária para sua implementação Esta implementação está se tornando economicamente viável não somente em função da redução dos custos da tecnologia de informações como também pela publicação de ferramentas de análise mais adequadas O artigo não pretendeu esgotar a questão mas tão somente oferecer subsídios teóricos para a compreensão do tema constituindose em parte de uma agenda de pesquisa em desenvolvimento pelos presentes autores Referências Bibliográficas 01 ACHABAL D D MC INTYRE S H SMITH S A KALYANAM K 2000 A Decision Support System for Vendor Managed Inventory Journal of Retailing Vol76 No4 p430454 7 HERER Y T LEVY R 1997 The Metered Inventory Routing Problem an Integrative Heuristic Algorithm International Journal of Production Economics Vol 51 p 6981 HOOP W J SPEARMAN M L 1996 Factory Physics Irwin JAILLET P HUANG L BARD J F DROR M 1997 A Rolling Horizon Framework for the Inventory Routing Problem Technical Report The University of Texas at Austin TX KLEYWEGT A J NORI V S SAVELSBERGH M W P 1998 A Computational Approach for the Inventory Routing Problem Technical Report Georgia Institute of Technology GA KLEYWEGT A J NORI V S SAVELSBERGH M W P 1999 The Stochastic Inventory Routing Problem with Direct Deliveries Technical Report Georgia Institute of Technology GA MINKOFF A S 1993 A Markov Decision Model and Decomposition Heuristic for Dynamic Vehicle Dispatching Operations Research Vol 41 No1 p 7790 TRUDEAU P DROR M 1992 Stochastic Inventory Routing Route Design with Stockouts and Route Failure Transportation Science Vol26 p 171184 J Sched 2009 12 257280 DOI 101007s1095100800819 The integrated productioninventorydistributionrouting problem Jonathan F Bard Narameth Nananukul Received 14 May 2007 Accepted 4 July 2008 Published online 20 August 2008 Springer ScienceBusiness Media LLC 2008 Abstract The integration of production and distribution decisions presents a challenging problem for manufacturers trying to optimize their supply chain At the planning level the immediate goal is to coordinate production inventory and delivery to meet customer demand so that the corresponding costs are minimized Achieving this goal provides the foundations for streamlining the logistics network and for integrating other operational and financial components of the system In this paper a model is presented that includes a single production facility a set of customers with time varying demand a finite planning horizon and a fleet of vehicles for making the deliveries Demand can be satisfied from either inventory held at the customer sites or from daily product distribution In the most restrictive case a vehicle routing problem must be solved for each time period The decision to visit a customer on a particular day could be to restock inventory meet that days demand or both In a less restrictive case the routing component of the model is replaced with an allocation component only A procedure centering on reactive tabu search is developed for solving the full problem After a solution is found path relinking is applied to improve the results A novel feature of the methodology is the use of an allocation model in the form of a mixed integer program to find good feasible solutions that serve as starting points for the tabu search Lower bounds on the optimum are obtained by solving a modified version of the allocation model Computational testing on a set of 90 benchmark instances with up to 200 customers and 20 time periods demonstrates the effectiveness of the approach In all cases improvements ranging from 1020 were realized when compared to those obtained from an existing greedy randomized adaptive search procedure GRASP This often came at a three to fivefold increase in runtime however Keywords Tabu search Production planning Lotsizing Inventory Vehicle routing problem Allocation model JF Bard N Nananukul Graduate Program in Operations Research Industrial Engineering The University of Texas 1 University Station C2200 Austin TX 787120292 USA email jbardmailutexasedu N Nananukul email naramethmailutexasedu 17 ANILY S FEDERGRUEN A 1990 One Warehouse Multiple Retailer with Vehicle Routing Cost Management Science Vol36 No1 p92114 BAITA F UKOVICH W PESENTI R FAVARETTO D 1998 Dynamic RoutingandInventory Problems A Review Transportation Research Vol33 No8 p585598 BALL MO 1988 Allocationrouting models and algorithms in Herer Y T and Levy R 1997 The Metered Inventory Routing Problem an Integrative Heuristic Algorithm International Journal of Production Economics Vol 51 p 6981 BARNESSCHUSTER D BASSOK Y 1997 Direct Shipping and the Dynamic SingleDepot MultiRetailer Inventory System European Journal of Operational Research Vol101 p 509518 CAMPBELL A CLARKE L KLEYWEGT A SAVELSBERGH M W P 1997 The Inventory Routing Problem Technical Report Georgia Institute of Technology GA CAMPBELL A CLARKE L SAVELSBERGH M W P 1999 An Inventory Routing Problem Technical Report Georgia Institute of Technology GA CANEN A G PIZZOLATO N D 1993 The Vehicle Routing Problem A Managerial Report Logistics Focus Vol1 No2 p 69 DANTZIG G B RAMSER J H 1959 The truck Dispatching Problem Managemente Science Vol6 No1 p8191 DROR M BALL M 1987 InventoryRouting Reduction from an Annual to a shortperiod Problem Naval Research Logistics Vol34 p891905 DROR M TRUDEAU P 1996 Cash Flow Optimization in Delivery Scheduling European Journal of Operational Research Vol88 p504515 FEDERGRUEN A ZIPKIN P 1984 A Combined Vehicle Routing and Inventory Allocation Problem Operations Research Vol32 No5 p10191037 GALLEGO G SIMCHILEVI D 1990 On the Effectiveness of Direct Shipping Strategy for the Onewarehouse Multiretailer Rsystems Management Science Vol36 No2 p 240243 GOLDEN B L ASSAD A A DAHL R 1984 Analysis of Large Scale Vehicle Routing Problem with an Inventory Component Em Campbell A Clarke L Kleywegt A and Savelsbergh M W P 1997 The Inventory Routing Problem Technical Report Georgia Institute of Technology GA 258 J Sched 2009 12 257280 the VMI model a vendor observes and controls the inven tory levels of its customers as opposed to conventional ap proaches where customers monitor their own inventory and decide the time and amount of products to reorder One of the benefits of VMI is that vendors can achieve a more uni form utilization of transportation resources leading to lower distribution costs It also offers them the flexibility to choose the most preferred transportation mode 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 de mand When a single entity is responsible for both planning and scheduling efficiencies are realized at all nodes in the system This has been the underlying motivation for the em phasis on supply chain management that is now widespread in the economy eg see OBrien and Tang 2006 In general the problem of optimally coordinating pro duction and transportation is called the production inventorydistributionrouting problem PIDRP eg see Lei et al 2006 Addressing these components in a single framework offers a holistic view of the logistics network and provides a good starting point for the full integration of the supply chain The PIDRP commonly arises in the retail in dustry where customers or outlets rely on a central supplier or manufacturer to provide them with a given commodity on a regular basis In the version of the problem addressed here a manufacturer must develop minimum cost production and distribution schedules for a single product that are sufficient to meet all customer demand over the planning horizon The purpose of this paper is to outline the full model of the PIDRP along with a composite solution methodol ogy that gives verifiably high quality results within accept able runtimes For planning purposes this means within one or two hours The primary components of the methodology are an allocation model for obtaining initial solutions and lower bounds on the optimum and a tabu search metaheuris tic Glover and Laguna 1997 with path relinking Resende and Ribeiro 2005 for improving the results The tabu search is distinguished by its neighborhood structure short and longterm memory functions and search strategies In the next section the PIDRP is outlined and a por tion of the relevant literature is reviewed with an empha sis on the most recent work In Sect 3 a formal definition of the allocation model is given which takes the form of a mixedinteger program MIP This is followed in Sect 4 by our tabu search algorithm Our lower bounding model is described in Sect 5 and several theoretical statements are made to offer some improvement to the results The compu tational results are highlighted in Sect 6 for benchmark in stances with up to 200 customers and 20 time periods The analysis is discussed in Sect 7 and several ideas are pre sented for extending the work 2 Problem statement and literature review Although the PIDRP can be defined more generally our fo cus is on a single facility that must meet the demand of its customers for a single commodity To ensure timely distribu tion and to avoid shortages excess production can be stored at either the plant or at the customer sites up to some limits however inventory cannot be transferred between sites and stockouts are not permitted It is further assumed that de mand is known for each day of the planning horizon and that all initial inventories are given In the model two lotsizing decisions must be made The first concerns the production inventory tradeoff the second relates to the daily distribu tion decisions over the planning horizon With regard to the latter deliveries are made routinely by a fleet of homogeneous vehicles that begin and end their run at the plant In the most complex scenario investigated a ve hicle routing problem VRP must be solved daily to deter mine in conjunction with the production decisions which customers to visit and how much product to deliver to each Due to either vehicle capacity limits or favorable cost trade offs 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 The PIDRP is different than traditional VRPs because it requires multiple customer visits to satisfy demand spread out over an extended period of time It is most similar to the inventory routing problem IRP Abdelmaguid and Dessouky 2006 Bard et al 1998 Golden et al 1984 and Dror and Ball 1987 and the periodic routing problem PRP Gaudioso and Paletta 1992 Mourgaya and Vanderbeck 2007 and Parthanadee and Logendran 2006 Although there has been much research on these two problems little of it carries over to the PIDRP which has only been studied intermittently The primary reason relates to the formidable complexity of its structure as defined by a combination of a capacitated lotsizing problem Gutiérrez et al 2007 Pochet and Wolsey 2006 and a capacitated multiperiod VRP The full problem has so far proven to be beyond the capability of exact methods By decoupling of the lot sizing and routing decisions though several researchers have had some success in finding good solutions with heuris tics Chandra and Fisher 1994 for example first determine a production schedule without regard to the logistics Next they developed a distribution schedule for each planning pe riod based on the results obtained from the firststage model This approach worked well when there was enough inven tory in the system to buffer production from the distribution operations but consequently led to increased holding costs The IRP and the PRP are relaxations of the PIDRP differ ing in several ways Neither for example takes the produc tion decision and inventory level at the plant into considera tion In addition the PRP assumes that the delivery patterns J Sched 2009 12 257280 259 defined by delivery frequencies or delivery days are given in advance and tries to select the most suitable pattern for each customer Golden et al 1984 were the first to investigate the inter related problem of inventory allocation and vehicle routing The particular application involved an energyproducts com pany 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 Historical data was used to calculate an average consumption rate for each customer and the latest replenishment data was used to cal culate when each customer could be expected to need a re supply The next replenishment was scheduled based on this information The proposed heuristic included a customer se lection algorithm that decides which customers to visit on each day in a costeffective way and a VRP component to construct daily routes Testing was done on instances with up to 3000 customers The results were compared to those from a simulation experiment and showed an improvement of 84 in the number of gallonshour a 50 reduction of stockouts and a 23 reduction in total costs The basic methodology was extended by Bard et al 1998 to better account for demand uncertainty and the use of satellite fa cilities for extending daily tours Parthanadee and Logendran 2006 considered a multi product multidepot periodic distribution problem and for mulated it as a MIP In the model they assumed that the daily demand of each customer was known and that all de liveries could be completed in one day within specified time windows and allowing for multiple vehicle trips Backorder ing of products at the depots was allowed To rationalize de liveries over the planning horizon a set of predefined pat terns was introduced and ranked by the customers Within a tabu search heuristic a penalty scheme was used to direct the search away from those patterns that were deemed unde sirable In their version of the PRP Mourgaya and Vanderbeck 2007 proposed a column generationbased 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 con sidered 1 optimizing the compactness of the geographical regions to which customers were assigned and 2 balanc ing 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 The pricing subproblem decomposed into τ clustering prob lems one for each time period Computational tests were performed using 20 50 and 80 customer data sets that in volved comparisons of various implementation options re lated to initialization column generation strategies number of passes in the rounding heuristic and problem specifica tions 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 longrun average costs were mini mized and all demand was satisfied In their model no in ventory capacity constraints were imposed on the warehouse or on the retailers Testing was done on problem instances with 50 and 75 retailers There has been a vast amount of research in the areas of production planning and inventory management over the last 40 years with the field now including all aspects of the sup ply chain Anily and Federgruen 1993 for example ana lyzed fixed partition policies for the inventory routing prob lem with constant deterministic demand rates and an unlim ited number of vehicles The routing patterns were deter mined by using a modified circular clustering scheme After the customers were grouped those 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 taking into account percentage allocations The routing logic required that all customers in a region be vis ited 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 scheme For extensions to multiple items and warehouses see Feder gruen and Tzur 1999 Lei et al 2006 proposed a twophase solution approach for the PIDRP with multiple plants and a heterogeneous fleet Our methodology is very similar to theirs In phase one a restricted version of the problem that contained all but the routing constraints was solved The results provided a production schedule and the number of items to be deliv ered to each customer in each period The solutions were shown to always be feasible to the original problem how ever they did not allow for the consolidation of lessthan transporter loads LTL which could have reduced overall transportation costs To address this issue a routing heuris tic based on an extended optimal partitioning procedure was used in phase two to consolidate the LTL assignments into more efficient delivery schedules Testing showed that the approach gave good solutions to instances with up to 50 cus tomer sites over 2 to 4 planning periods Boudia et al 2006 2007 proposed both a memetic al gorithm with population management MAPA and a reac tive greedy randomized adaptive search procedure GRASP with pathrelinking Resende and Ribeiro 2005 to solve the PIDRP The model included a single plant and a set of cus tomers located on a grid Holding costs at the customer sites were assumed to be negligible compared to the holding costs at the plant and so were ignored Like ours the objective was to minimize the sum of production holding and transportation costs while ensuring that all demand was satisfied over the planning horizon As an enhancement the authors proposed two strategies to combine path relinking and GRASP The first strategy was to perform path relinking upon the termination of GRASP the second was to perform path relinking every time GRASP improved on the incumbent The computational results showed that on average the first strategy performed better than the second on the 50 and 100 customer instances with 20 time periods but was not as good on the 200 customer instances The algorithms converged within 100 s on a 28 GHz computer in all cases but no lower bounding procedures were provided to gauge the quality of the solutions The reactive GRASP on average provided more savings than the solutions obtained with the GRASP alone by 08 and the solutions obtained with the relinking feature by 042 3 Model formulation We are given a single production facility and a set of n customers geographically dispersed on a grid Each customer i N 1 n has a fixed nonnegative demand dit in time period t of the planning horizon that must be satisfied ie shortages are not permitted If production takes place at the facility in period t then a setup cost ft is incurred t T0 0 1 τ A limited number of items can be produced in each time period and a limited number can store at a unit cost of hP In general it is natural to equate a period with a day which we do but when production is scheduled for more than one shift in a day or when transportation times are measured in weeks a broader interpretation of a period is appropriate In constructing delivery schedules each customer can be visited at most once per day and each of θ homogeneous vehicles can make at most one trip per day The latter restriction implies that all routes overlap in time If cij is the cost of going from customer i to customer j and xijt is a binary variable equal to 1 if customer i is the immediate predecessor of customer j on a route in period t 0 otherwise then the routing costs are given by ijt cij xijt in the full model A limited amount of inventory can be stored at customer is site at a unit cost of hCi but transshipments between customers are not permitted cf Herer et al 2006 Moreover it is assumed 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 on day t is available for delivery only on the following morning this is common in food production and distribution eg see Villegas and Smith 2006 and all inventories are measured at the end of the day Demand on day t can be met from deliveries on day t and from ending inventory on day t1 at the customer site Initial customer inventory on day 0 simply reduces demand on subsequent days while initial inventory at the plant must be routed at the end of the planning horizon all inventory levels are required to be zero The goal is to construct a production plan and delivery schedule that minimizes the sum of all costs while ensuring that each customers demand is met over the planning horizon In so doing four critical decisions have to be made How many items to manufacture each day When to visit each customer How much to deliver to a customer during each visit Which delivery routes to use As part of the solution methodology we do not attempt to solve the full model but instead investigate a relaxation referred to as the allocation model In particular the routing term in the objective function is replaced by a distribution component The following additional notation is used in the developments Parameters Q capacity of each vehicle ImaxP maximum inventory that can be held at the production facility ICmaxi maximum inventory that can be held by customer i C production capacity of the plant fitC fixed cost of making a delivery to customer i on day t eitC variable cost of delivering one item to customer i on day t Decision variables pt production quantity on day t zt 1 if there is production on day t 0 otherwise ItP inventory at the production facility at the end of day t ItC inventory at customer i at the end of day t wit amount delivered to customer i on day t zitC 1 if a delivery is made to customer i on day t 0 otherwise Model ϕIP Minimize tT ft zt tT iN fitC zitC tT iN eitC wit tT0τ hP ItP tTτ iN hiC ItC 1a subject to ItP It1P pt iN wit t T0 1b iN wit It1P t T 1c pt Czt t Tτ 1d p0 iN di1 Ic0 1e ItC It1C wit dit i N t T 1f wit Dmaxit zitC i N t T 1g iN wit 08Qθ t T 1h 0 ItP ImaxP 0 ItC ImaxCi i N t Tτ ItP Icτ 0 i N 1i zt 0 1 zitC 0 1 pt 0 wit 0 i N t T 1j where T T00 and Dmaxit min Q ltτ dil The objective function minimizes the sum of production setup costs a surrogate for the routing costs second and third terms holding costs at the plant and holding costs at the customer sites An explanation of the values used for the coefficients fitC and eitC is given in the next section Because all demand must be met production costs which are assumed to be linear and independent of time can be omitted as can any initial inventory costs Constraints 1b and 1f are inventory flow balance equations in which it is assumed that the initial inventories I0P and I0C are given for all customers i N Constraint 1d limits production on day t to the capacity of the plant A simple way to tighten this constraint is to replace C with Ct minC ImaxP Dmaxt1 where the third term is the demand of all customers for the remainder of the planning horizon The assumption that items produced on day t are only available for delivery on day t 1 implies that pt ImaxP and pt 0 It is possible to strengthen the latter inequality by subtracting from ImaxP the reduction in inventory due to deliveries on day t to get pt ImaxP It1P iN wit but this constraint is dominated by 1b To ensure that demand on day 1 can be met it is necessary to include 1e which allows production on day 0 If I0P Ic0 0 then p0 iN di1 or the problem is infeasible As indicated by constraints 1c the total amount available for delivery on day t is limited by the amount in inventory at the plant on day t 1 The specific amount delivered to customer i is limited by the parameter Dmaxit in 1g which is the smaller of the vehicle capacity Q or the cumulative demand from day t to the end of planning horizon τ Constraints 1h represent an aggregate relaxation of the routing constraints common to capacitated VRPs They simply restrict the total amount that can be delivered on day t to a fixed percentage of the total transportation capacity and provide a hedge against the need for split deliveries Testing showed that using a value of 80 always yielded feasible solutions To conclude the formulation variable bounds are specified in 1i and 1j To obtain the full PIDRP model several modifications to 1a1j are needed First the second and third terms in 1a should be omitted and replaced with ijt cij xijt next constraint 1h should be deleted finally routing constraints that take into account the load on the vehicles must be added Figure 1 depicts an abbreviated network flow diagram of the PIDRP The top portion of the figure represents the production facility denoted by C0 for days 0 through τ Instead of demand driving the decisions on each day t for the plant a series of delivery decisions has to be made for each customer i The amount delivered denoted by wit is combined with customer is inventory It1C at the beginning of day t to satisfy demand dit The corresponding flow is shown in the bottom portion of the figure Because the inventory at the end of the planning horizon at both the plant and the customers sites is required to be zero there is no horizontal flow exiting node τ in either network The size of model 1 is determined largely by constraints 1f and 1g and the number of binary variables zitC all of which grow at a rate proportional to Onτ The majority of the other constraints only grow at a rate proportional to Oτ Problem instances with nτ 30 can be solved with CPLEX 81 in less than 10 min for example with n 5 and τ 2 4 6 corresponding solution times tcplexτ were 4 196 and 494 s However for n 10 and τ 4 CPLEX did not converge within 2 hours on any instances and at the point exhibited an optimality gap of over 10 Although not specified in 1i and 1j our solution methodology requires that the delivery quantities wit be integral The following proposition shows that the allocation model always returns integer values for not only wit but for pt ItP and ItC as well Proposition 1 When the setup variables zt and zitC are fixed at 0 or 1 there exists and optimal solution to model 1a 1j such that wit pt I tP and ItC are integral for all i N and t T0 Proof We show that when zt and zitC are fixed the remaining components of the model are equivalent to a pure net work flow problem whose constraint matrix is known to be totally unimodular The first step is to represent the flows at the plant more accurately by accounting for the fact that production on day t is held in inventory that day and is only available for delivery on day t 1 This can be done by creating three inventory nodes at the plant in each period one for serving customers a second to bind the flow to customers and a third to receive the current days production This division is shown in Fig 2 where primes and double primes are used to distinguish the three nodes At node 1 for example items may be withdrawn and delivered to customers however an upper bound of 08Qθ is placed on the arc from node 1 to node 1 to ensure adherence to the capacity restrictions not shown Items produced on day 1 are channeled to node 1 rather than node 1 At the end of the day whatever inventory remains flows to node 2 as indicated by the variable IPt When zt and zitC are fixed the remaining variables are bounded by integer values Conservation of flow at the nodes with primes is It1P IPt1 iN wit which is equivalent to 1c given the variable bound It1P 0 for all t Conservation of flow at the nodes without primes is essentially 1b with IPt1 replaced with its equivalent The delivery limits in 1h are enforced by the bounds on the arcs between the nodes with single and double primes Finally constraint 1f represents conservation of flow at the customer sites and is already in pure network form Thus we have a bounded pure network flow problem that is guaranteed to have an optimal integral solution Fig1 Network flow representation of productioninventorydistribution problem Production facility C0 Time period p0 p1 p2 pt1 0 1 2 τ2 τ1 τ I0P IP I2 P Iτ2 P Iτ1 P w11 w1n C1 Cn C1 Cn Customer Ci Time period wi1 wi2 wiτ1 wiτ 0 1 2 τ2 τ I0C IC IC IC IC d1i d2i diτ1 diτ 4 Solution methodologytabu search A twophase approach is used in the design of our reactive tabu search algorithm for solving the PIDRP In the first part of phase 1 an initial solution is found by solving the allocation model 1a1j The results provide customer delivery quantities wit for all i 1 n and t 1 τ In the second part these values become the demand for τ independent routing problems An efficient CVRP subroutine also based on tabu search Carlton and Barnes 1996 is called to find solutions It is treated as a black box and will not be discussed here In phase 2 neighborhood search is performed to improve the allocations and routing assignments found in phase 1 41 Initial solution Absent of the routing component model 1 represents a pure lotsizing distribution problem The two terms it fitC zitC and it eitC wit in 1a serve as surrogates for the actual routing costs given by ijt cij xijt In our implementation the specification of the coefficients fitC and eitC depends on the dimensions of the problem For instances where n2τ 500 we set fitC 2ci0 and eitC 0 for all i and t that is the routing costs on any day t are approximated solely by the cost of a round trip between the depot and customer i When n2τ 500 indicating a fairly large instance it is not practical to solve 1a1j with the setup variables zitC included In this case we set zitC fitC 0 and put eitC2ci0dit for all i and t that is the cost of making a delivery to customer i on day t is approximated by the J Sched 2009 12 257280 263 Fig 2 Network flow with detailed inventory nodes round trip cost of visiting customer i directly from the de pot divided by his demand on day t Empirically solutions to the allocation model with these settings were seen to be close to those of the full model when the production and fleet capacities were larger than the total daily demand This is shown in the computations section of the paper Solving model 1 determines the amount of produc tion in each day t pt t 01τ and the quan tity delivered to each customer i on each day t wit t 01τ i 1n Given these values an initial so lution to the PIDRP is found by calling the VRP subroutine to construct delivery routes for each day t 42 Neighborhood definition A neighborhood is a set of points that can be reached from the current solution by performing one or more moves that in general take the form of insertions exchanges or replace ments From any incumbent a new solution is determined by identifying the best solution within its neighborhood For the PIDRP we define the neighborhood as all feasible points that can be reached by two types of moves between periods The first is called a swap and involves a partial exchange of delivery quantities between two customers i1 in period t1 with quantity wi1t1 and i2 in period t2 with quantity wi2t2 where t2 is the first period after t1 in which wi2t2 0 For customer i1 the move considers the maximum portion of wi1t1 that can be reassigned to period t2 without caus ing a shortage in period t1 to be exchanged with the full amount wi2t2 We show below that it is suboptimal to trans fer less than the maximum portion of wi1t1 If customer i1 was not scheduled for a delivery in period t2 then he must be inserted into one of the θ routes In general a swap pro duces a change in holding costs and a change in transporta tions costs in periods t1 and t2 If the net effect is negative then the swap is beneficial Note once again that in the full model the transportation costs are the sum of the routing costs cij and hence do not depend on the delivery quantities Proposition 2 Let wi1t1be the planned delivery quantity to customer i1 in period t1 and let I C i1t1 be the net inven tory position after demand is satisfied ie I C i1t1 I C i1t11 wi1t1 di1t1 For any feasible swap between customers i1 and i2 in periods t1 and t2 the optimal swap amount is equal to min wi1t1 I C i1t1 Proof Consider any swap amount χ for customer i1 in period t1 such that χ min wi1t1 I C i1t1 We will show that there exists a swap amount χ such that χ χ min wi1t1 I C i1t1 that results in a better movevalue The movevalue is equal to the sum of the change in holding costs for customers i1 and i2 the change in trans portation costs in periods t1 and t2 and the change in hold ing cost at the plant For the case where the swap amount is χ the movevalue is hC i1χ ε1 t1 ε1 t2 hC i2 wi2t2 ρ1 where ε1 t1 and ε1 t2 represent the change in transportation costs in period t1 and t2 respectively and ρ1 is the change in holding cost at the plant For the latter we have ρ1 hP χ wi2t2τ 1 t2 τ 1 t1 where τ 1 t1 τ 1 t2 are respectively the periods associated with t1 and t2 in which production must be adjusted to ensure that the swap is feasible Al though the results are stated for a single pair of adjustment periods τ 1 t1 and τ 1 t2 it is easy to show that they hold when production must be adjusted in more than these two periods Similarly for the case where the swap amount is χ the movevalue is equal to hC i1χ ε2 t1 ε2 t2 hC i2 wi2t2 ρ2 The symbols have the same meaning as above but a su perscript 2 is used instead of 1 To compare transportation costs the following two cases must now be considered Case 1 When χ wi1t1 the situation when χ I C i1t1 there is no difference in transportation costs that is ε1 t1 ε2 t1 and ε1 t2 ε2 t2 because the customers who are scheduled for a delivery in period t1 are the same as those scheduled for a delivery in t2 264 J Sched 2009 12 257280 Case 2 When χ wi1t1 it is not necessary to make a delivery to customer i1 in period t1 so ε2 t1 ε1 t1 in period t1 and ε2 t2 ε1 t2 in period t2 Because χ χ more product is moved to a later period which implies a greater reduction in the holding cost at the plant ρ2 ρ1 Finally by noting that the amount of product moved to an earlier period wi2t2 is the same in both cases we have hC i1χ ε1 t1 ε1 t2 hC i2 wi2t2 ρ1 hC i1χ ε2 t1 ε2 t2 hC i2 wi2t2 ρ2 The second move is called a transfer and examines each customer i one at a time and tries to reassign the delivery quantity wit1 scheduled for t1 to the latest period call it t2 preceding t1 in which a delivery is scheduled for at least one customer i that is t2 maxt t t1 wit 0 for some i N In this case we have the following result Proposition 3 For any customer i when considering a transfer move between periods t1 and t2 it is suboptimal to reassign less than the full amount wit1 to period t2 Proof If only a portion of wit1were transferred then the holding costs would increase in t1 for customer i while the transportation costs would remain the same in both peri ods so there would be no benefit in considering intermediate cases Whether a swap or a transfer only moves that result in feasible solutions are permitted so it is necessary to check for violations of the production constraints and the inven tory bounds at the plant and the customer sites Moreover a VRP must be solved in periods t1 and t2 to see whether feasible routes can be found after the move and if so what their optimal costs are The value of a move is determined in part by these costs which must be calculated for each candidate To begin FeasibilityCheckAlgorithm is called to determine whether the move violates any of the produc tion storage or vehicle capacity constraints If not then MoveValueAlgorithm is called to determine the net ben efit see Nananukul 2008 Complexity of neighborhood A swap involves an exchange in delivery quantities between pairs of customers in two different periods t1 t2 so the number of possible can didates in the worst case is On2 in each period Cal culating the portion of wi1t1 that can be exchanged with out causing a shortage in period t1 can be done in O1 The feasibility check is On τ and determining the move value by solving the VRP takes On3 For τ pe riods then the amount of work associated with a swap is On3 nτn2τ On5τ n2τ 2 For transfer moves deliveries to each of the n customers in period t are consid ered for reassignment to period t 1 In the worst case there are On possibilities in each period Taking the feasibility check and the move value computations into account the amount of work associated with a transfer for all τ periods is On3 n τ nτ On4τ nτ 2 Thus the neigh borhood size is On5τ n2τ 2 Example of moves Figure 3 depicts a swap between cus tomers 1 and 3 in periods 2 and 3 respectively The periods are numbered on the far left and the customers on a specific route are represented by circles The depot customer 0 is at the start and end of each route implying that a single vehi cle only is required in each period The parameter values for the five customers in the example are as follows The hold ing cost for customer 1 is hC 1 20 while the holding costs for customers 2 to 5 are hC 2 hC 3 hC 4 hC 5 10 Initial inventories at all customer sites in period 2 are I C 11 I C 21 I C 31 I C 41 I C 51 0 It is also assumed that the storage ca pacity at the customer sites is unlimited and that the vehicle capacity Q 60 Beginning at the depot route costs are calculated by sum ming the individual link costs cij between customers i and j on the route where c01 175c04 10c12 25c02 c20 50c23 75c03 c30 25c45 30 and c51 10 The route cost in period 2 for example is c04 c45 c51 c12 c23 c30 175 Demand for customers 1 to 5 in pe riod 2 is 8 7 7 6 and 8 respectively and in period 3 is 4 4 8 6 and 8 These values are shown in the squares below the circles in Fig 3 only if a delivery is scheduled for the customer in the period of interest Before the swap customer 1 is scheduled for deliveries of 8 and 4 items in periods 2 and 3 respectively and customer 3 is scheduled for a delivery of 15 items in period 2 After the swap the amount to be delivered to customer 1 is in creased to 12 in period 2 and reduced to 0 in period 3 Also the amount to be delivered to customer 3 is decreased to 7 in period 2 in accordance with Proposition 2 and a new de livery of 8 units is scheduled in period 3 Thus for customer 3 the exchange involves the maximum number of items 8 that can be moved in period 2 without causing a shortage customer 3 has a demand of 7 in period 2 and currently 15 items are scheduled for delivery For customer 1 all 4 units that were to be delivered in period 3 are now scheduled for period 2 Note that before and after the swap the total num ber of items scheduled for delivery in either period do not exceed the vehicle capacity Before the swap a total of 58 and 12 items are scheduled to be delivered in periods 2 and 3 respectively while after the swap the delivery quantities are 54 and 16 The results indicate that a big cost reduction is achieved by eliminating the link from the depot to customer 1 c01 175 in period 3 This cannot be realized though by simply rescheduling the 4 items to be delivered to customer 1 in period 3 to period 2 because the vehicle capacity 60 would J Sched 2009 12 257280 265 Fig 3 Example of an exchange or swap move between customers 1 and 3 be exceeded A swap between periods is required Before the swap the transportation cost in period 3 is c01 c12 c20 250 and after the swap it is c02 c23 c30 150 where the solution to the VRP for period 3 indicated that it was optimal to insert customer 3 after customer 2 In period 2 there is no change in route so the transportation cost remains the same as do the holding cots Calculating the difference between the before and after costs gives the movevalue which in this case is 100 Using the same data and starting with the schedule in the bottom portion of Fig 3 Fig 4 gives an example of a sin gle customer transfer move Before the transfer customer 2 is scheduled to receive deliveries of 7 and 4 items in pe riods 2 and 3 respectively where now t1 3 and t2 2 The transfer eliminates the delivery in period 3 by moving all 4 items to period 2 in accordance with Proposition 3 and is feasible because the total load on the vehicle only goes up to 58 which is less than the capacity The move in creases the holding cost for customer 2 from 0 of 40 but the transportation cost in period 3 is reduced by 100 giving a movevalue of 60 More specifically after the transfer the net change in holding cost is 4 10 40 and is asso ciated with customer 2 the total transportation cost in pe riod 3 is c03 c30 50 and the net change in transporta tion cost in period 2 is 0 As a result the movevalue is 50 40 150 60 43 Tabu list and aspiration criterion Tabu search transcends local optimality by forbidding cer tain moves on a temporary basis In general the process is implemented with a shortterm memory structure that pro scribes a subset of the moves in a neighborhood The tabu list stores all forbidden moves Its length normally called the tabu tenure indicates how many iterations a certain move is forbidden The tabu status indicates whether a move on the tabu list is currently forbidden In our implementation we use a reactive tabu list mean ing that the tabu tenure is dynamic We also specify an as piration criterion that allows the tabu status of a move to be overridden Specifically each entry on the tabu list is a com bination of a pair of customers and a pair of time periods represented by the 4dimensional vector customer i1 cus tomer i2 period t1 period t2 Note that for a transfer move where only a single customer is involved a default value of zero is used for customer i2 The length of tabu list is kept constant as long as progress is being made but it is increased when there is no improvement in some fixed number of it erations Because a move associated with a customer in one 266 J Sched 2009 12 257280 Fig 4 Example of a transfer move that assigns a delivery to an earlier period period could create a series of moves in the following peri ods thus affecting up to τ 1 periods we set the tabu tenure proportional to τ Its initial value is set to τ2 and then in creased to τ when there is no improvement for 5 iterations However when a new incumbent is found the tabu tenure is set back to τ2 The process is repeated until the tabu search stopping criteria are met Although some implementations allow infeasible moves as an additional means of overcoming local optimality we do not Preliminary testing showed that when infeasible moves were considered at each iteration runtimes became excessive due to the increase in neighborhood size for even small instances with up to 20 customers After the current neighborhood is searched and a new incumbent is found the move that led to the incumbent is added to the tabu list The tabu status of a move can be overridden though when a certain aspiration condition is met In our case this condition is associated with a move on the tabu list that gives a better solution than the incum bent 44 Search strategy Two important strategies for tabu search are intensification and diversification Intensification focuses on creating so lutions that have good attributes with respect to routing for example solutions that include low cost arcs Diversi fication is the ability of the algorithm to expand the search by generating solutions that have attributes different from those encountered at previous iterations These two strate gies counterbalance and reinforce each other In our algorithm intensification is implemented by using incentives and diversification by using penalties and short term memory The latter is provided by the tabu list Incen tives and penalties are a function of longterm memory and are represented by the n n τ τdimensional ma trices F I and F P respectively For a move that involves customer i1 in period t1 and customer i2 in period t2 an element F I i1i2t1t2 of F I represents the number of times the set i1i2t1t2 has been involved in a move that improves the solution Similarly an element F P i1i2t1t2 of F P repre sents the number of times the set i1i2t1t2 has been involved in a nonimproving move Any move associated J Sched 2009 12 257280 267 with these sets is either rewarded or penalized in accor dance to the values of F I i1i2t1t2 and F P i1i2t1t2 With respect to the above example if a candidate swap involves customers i1 and i2 in periods t1 and t2 then the actual value of the move is movevalue ρIF I i1i2t1t2 ρP F P i1i2t1t2 rather than movevalue alone where ρI and ρP represent incentive and penalty multipliers see next section for more detail Implementing exhaustive search for the neighborhood defined in Sect 42 is possible for small size problems n2τ 500 This was the approach that we took during al gorithmic development For larger instances n2τ 500 including those in the Boudia et al 2007 data sets it is not practical to examine all candidate moves To reduce the computational effort we randomly select a subset of moves and then adaptively decide according to the progress of the algorithm which of two rules to follow The first rule places a 4 to 1 emphasis on transfers over swaps and is used when the most recent tabu iteration resulted in an overall cost re duction The second rule reverses the emphasis and is used when the most recent tabu iteration did not improve upon the incumbent More detail is provided in Nananukul 2008 At this point all moves on the candidate list are processed and the one with the best movevalue is selected The rele vant production and inventory levels are updated the stop ping criteria are checked and if not met the next tabu itera tion is performed 45 Algorithm description Figure 5 summarizes the major steps of the tabu search algo rithm for large instances In the first step an initial feasible solution is created with the procedure described in Sect 41 The result is saved as both the current solution and the best solution found Then the algorithm iterates until either a pre specified maxiterations is reached or there is no improve ment in maxnoimprove iterations Of all the admissible candidates the move selected at each iteration is the one that has the smallest value of movevaluemovepenalty moveincentive A move is considered admissible if it is not on the tabu list or its tabu status has been overridden by the aspiration criterion For the current iteration once the best move is found it is executed and the tabu data structures are updated The parameters and data structures used in the algorithm are as follows numiteration current iteration number maxiterations maximum number of iterations allowed currsoln current solution after executing the best move currcost cost of current solution bestsoln incumbent solution bestcost cost of incumbent solution movevalue difference in the cost before and after the swap move bestmovevalue lowest movevalue among all candi dates freqmatrix longterm memory function that stores the frequency of each possible move movepenalty additional penalty imposed by longterm memory freqmatrix F P when a swap involves customer i in period k and customer j in period lmovepenalty ρP F P ijklρP is set to 5 of the objective function value of the initial feasible solution divided by maxiterations This value gives an approximation of the average amount of penalty per iteration moveincentive additional incentive imposed by long term memory freqmatrix F I when a swap involves customer i in period k and customer j in period l moveincentive ρIF I ijkl where ρI is set to ρP be cause we penalize a move that results in a nonimproved solution the same amount as the incentive of the same move that results in an improved solution bestmovepenalty movepenalty for the best move of all candidates bestmoveincentive moveincentive for the best move of all candidates noimprove number of consecutive iterations during which no better solutions are found maxnoimprove maximum number of consecutive no improvement iterations allowed tabulist shortterm memory function that stores a list of moves that are forbidden tabusize total number of iterations for which a move is held on the tabulist candidaterule candidates that are considered during each iteration in accordance with rules 1 and 2 candidatemove solution that results when the moves as sociated with candidaterule are applied Admissiblemove candidate move that either satisfies the aspiration criterion or is not on the tabulist When evaluating each candidate move it may be pos sible to achieve a further cost reduction by adjusting the production and inventory levels which are not necessar ily optimal for the updated delivery schedule Given the updated values of wit Proposition 1 allows us to solve a linear program to find the optimal production and corre sponding inventory levels however doing so at each it eration is too costly because of the large number of pos sible moves so in the implementation exact solutions are computed once every five iterations At all other times we attempt to improve the solution by calling the heuristic ProductionLevelAdjustmentAlgorithm See Appendix for the details 1 Generate initial feasible solution and save it as currsol and bestsoln 2 Calculate the cost of currsol and set bestcost currcost 3 Initialize tabulist and freqmatrix set noimprove 0 and set candidaterule 1 4 Do bestmovevalue infinity for all candidatemove ifcandidatemove is not on tabulist or if it satisfies aspiration criterion ifFeasibilityCheckAlgorithm return true call MoveValueAlgorithm store result in movevalue ifmovevalue movepenalty moveincentive bestvalue bestmovepenalty bestmoveincentive bestmovevalue movevalue bestmovepenalty movepenalty bestmoveincentive moveincentive bestmove currentmove execute bestmove currcost currcost bestmovevalue update tabulist and freqmatrix and adjust production levels ifcurrcost bestcost bestcost currcost bestsoln currsoln noimprove 0 set candidaterule 1 else noimprove noimprove 1 set candidaterule 2 Whilenumiteration maxiterations or noimprove maxnoimprove 5 Lower bound computations To gauge the quality of the solutions provided by tabu search lower bounds on the true optimum are needed The simplest way to obtain a lower bound is to solve the LP relaxation of the full model Initial testing on small instances showed gaps of roughly 30 to 260not a useful measure As an alternative we developed a procedure based on the allocation model that gives much better results To obtain a lower bound on the true optimum and hence on the feasible solutions provided by tabu search a valid relaxation of the full model must be solved To formulate such a relaxation we begin with 1a1j and make several modifications First the cost coefficients fCit for each customer i are set to the shortest distance to either the depot or any of the other customers that is fCit mincij j in N0i Next the cost coefficients associated with wit are set as follows eCit c0Q where c0 minc0j j in N Finally the righthand side of 1h is increased to Qθ in each period t and the corresponding allocation model solved to get φLB Proposition 4 The solution to the modified allocation model lower bounding model provides a valid lower bound on the true optimum to the PIDRP that is φLB φPIDRP Proof First note that any solution that is feasible to the full model is feasible to the modified allocation model This follows because 1h with its righthand side set to Qθ is a valid relaxation of the routing constraints in the full model We now show that the modified costs sumit fCit zCit sumit eCit wit underestimate the true costs sumijt cij xijt which in any feasible solution contain one arc cost for each customer that receives a delivery in period t and one arc cost for each vehicle used in period t By design the first term underestimates the individual customer arc costs and the second underestimates the vehicle arc costs from the depot With respect to the vehicle arc costs note that the number of vehicles required to deliver sumit wit units in period t is at least ceilingsumit wit Q Removing the integer requirement and multiplying by c0 gives the desired result A slightly improved lower bound can be obtained when the holding costs are small and production setup costs are large with respect to the vehicle cost which is the case here Lemma 1 Assume that the minimum production setup cost fP minft t 0 1 τ 1 is much larger than the cost c0 of using a vehicle in any period and that the customer holding costs are negligibly small ie hCi 0 for all i in N Let Rt ceilingsumi wit Q sumi wit Q and rt sumi wit Q floorsumi wit Q be the fraction of a vehicle corresponding to the difference between the rounded up and rounded down number of vehicles in a solution in period t respectively Under the stated conditions if c0 rt t t hP rtQ tP tP hP rtQ c0 c0 Rt 0 for all t t 2 c0 rt t t hP rtQ tP tP hP rtQ c0 c0 Rt 0 for all t t 3 then the lower bound φLB can be increased by rounding up the number of vehicles used in period t and multiplying the corresponding fraction Rt by c0 Summing over all periods gives the following adjusted lower bound φLB φLB sumt1τ c0 sumi wit Q sumi wit Q 4 Proof The assumption that fP c0 implies that it is never advantageous to set up for production in order to avoid using an additional vehicle in any period Therefore the result follows if this is the only feasible option for eliminating a vehicle in some period t Rounding up number of vehicles used in each period will not change the optimality of the solution if the increase in cost c0 Rt in period t does not exceed the cost of rescheduling the delivery quantity Q rt to another period t Case 1 t t Given φLB the objective function value when rounding up the number of vehicles used without rescheduling in both periods t and t is Cost1 c0 Rt c0 Rt φLB 51 Now let tP and tP be the periods where production of Q rt takes place before and after rescheduling respectively The objective function value after rescheduling is bounded below by Cost2 φLB c0 rt t t hP rt Q tP tP hP rt Q t t maxhCi i in N rt Q c0 sumi wit rtQQ sumi witQ 52 Thus if Cost2 Cost1 for all t t then rounding up the number of vehicles used in period t does not change the optimality of the solution Case 2 t t Using the same logic the objective function value when rounding up the number of vehicles used after rescheduling is bounded below by Cost3 φLB c0 rt t t hP rt Q tP tP hP rt Q t t minhCi i in N rt Q c0 sumi wit rt QQ sumi wit Q 53 The implication of 51 and 53 is that rounding up the number of vehicles used in period t does not change the optimality of the solution if Cost3 Cost1 for all t t Now when hCi 0 for all i in N 51 52 and 53 give Case 1 t t c0 rt t t hP rt Q tP tP hP rt Q c0 c0 Rt 0 Case 2 t t c0 rt t t hP rt Q tP tP hP rt Q c0 c0 Rt 0 A second and third improvement in the lower bound can be obtained by taking into account the fact that exactly ceilingsumi wit Q vehicles must return to the depot in any solution to the full model A fourth improvement can be obtained by recognizing that some minimum number of vehicles must be used in each period in order to meet demand The proofs of Lemmas 2 and 3 similar to the proof of Lemma 1 and can be found in Nananukul 2008 J Sched 2009 12 257280 271 were randomly generated on a 100 100 Euclidean grid For each customer i demand was uniformly distributed be tween 0 and the storage capacity I C maxi and for each period t was positive with a probability that varied from 05 to 1 The value of I C maxi depended on i and the number of cus tomers in the specific data set The vehicle capacities were Q50 8000Q100 8000Q200 12000 and the number of vehicles were θ50 5θ100 9θ200 13 The first set of experiments was designed to gauge CPLEXs performance on the full problem and did not make use of the Boudia et al data sets 61 Preliminary testing To determine the limits of CPLEX to solve the PIDRP di rectly we randomly generated 10 instances on a 100 100 grid The first step was to fix the number of cus tomers n 5101520 and the number of time periods τ 2468 by selecting a range of values from these sets The remaining parameter values and the output can be found in Nananukul 2008 To summarize the results indicated that when either the number of customers or the number of periods is increased the runtime increased dramatically in almost all cases This was to be expected because the PIDRP requires the solu tion of a CVRP in each period When the number of cus tomers was increased the number of binary variables and constraints in each CVRP increased by On2 in the worst case implying a significant increase in both the number of binary variables and the number of constraints in the PIDRP model In all instances the LP solution was obtained in less than a second but CPLEXs MIP solver required anywhere from 4 to 7200 s the runtime limit imposed The instances that could not be solved by CPLEX were the ones that had nτ 40 To test the performance of the VRP subroutine Carl ton and Barnes 1996 we created eight singleperiod prob lem instances in which the numbers of customers ranged from 5 to 150 The same generating procedure used for the PIDRP instances were used here but the inputs were limited to customer locations customer demand vehicle capacity and number of vehicles The results indicated that the VRP subroutine gives high quality solutions denoted by φVRP within a much smaller amount of time compared to CPLEX and is especially ef fective for the larger instances with 30 or more customers For the small size problems n 510 the VRP subroutine provided the same solution as CPLEX within about the same runtime The gap between φVRP and φcplex was less than or equal to zero in all but one case Also the VRP runtimes were substantially less than those of CPLEX with the dif ference increasing dramatically as the number of customers increased We also compared the phase 1 solutions obtained from the allocation model 1 with route optimization in each pe riod with the solutions from CPLEX It was seen that phase 1 reliably provided high quality solutions in a negligible amount of time compared to CPLEX In several instances the optimal solution was found by the allocation model and in the other cases it was at most 6 off however for some instances it provided better results than CPLEX Also the phase 1 runtimes were no more than a few seconds not in creasing much with either n or τ 62 Small instances Table 1 summarizes the results for the 50customer 20time period instances Columns 2 and 3 give the best GRASP solutions φgrasp from Boudia et al and the corresponding runtimes tgrasp for 500 iterations Their computations were performed on a 28 GHz computer with 512 MB of RAM running Windows XP As a word of caution it is always problematic to make comparisons across different platforms and different algorithms There is also some arbitrariness in establishing the termination criteria for metaheuristics like GRASP or tabu search because it is rarely evident when ad ditional iterations will not lead to improvement Therefore all runtimes and other performance measures presented be low should be interpreted with this in mind Columns 4 5 and 6 are phase 1 solutions φPh1 run times tPh1 and percent deviation from the best GRASP solutions respectively Column 7 gives the time to solve the allocation model with CPLEX An optimality gap of 0 was used in all cases Columns 8 9 and 10 report the phase 2 results In those computations tabu search was allowed to run for a maximum of 50 iterations but was terminated ear lier when no improving solution was found in 5 consecutive iterations The 50iteration limit was reached in only two instances Column 9 gives the iteration on which the best phase 2 solution was found From Table 1 we see that the phase 1 solutions are 107 better on average than the best known solutions and were obtained with significantly smaller runtimes In fact tPh1 was 81 less than tgrasp on average whereas φPh1 was better than φgrasp in all cases The gap between φPh1 and φgrasp ranged from 4 to 17 As to be expected φPh2 was smaller than φPh1 in all cases as well provid ing an average improvement of 58 as reported in the last column For the GRASP runtimes averaged 977 s compared to 184 s for phase 1 and 4334 s for phase 2 About 75 of tPh2 resulted from calling CPLEXs LP solver to determine the optimal updated production quanti ties for each candidate move If this option is omitted and only the ProductionLevelAdjustmentAlgorithm is used for this purpose the GRASP is about 1 faster than our twophase approach The 58 improvement between phase 1 and phase 2 however drops to approximately 4 Lemma 2 Let yt ceilingsumi wit Q be the minimum number of vehicles used in period t and let ci mincij j in N be the least cost transition for customer i to another customer Also let Nt i zCit 1 fCit ci0 be the set of customers who receive a delivery in period t and whose corresponding cost is ci0 and let Δci ci ci0 be the opportunity cost of not using the minimum cost arc in a solution Now order the nt Nt customers such that Δci1 Δci2 Δcint and let lt nt yt be the excess number of customers whose delivery cost is ci0 Assume that the minimum production setup cost fP minft t 0 1 τ 1 is much larger than the cost c0 of using a vehicle in any period If i Δci t t wit hP tP tP wit hP t t wit hCi fCit for all t t and ii Δci t t wit hP tP tP wit hP t t wit hCi fCit for all t t then the lower bound φLB can be adjusted as follows φLB φLB sumt1τ sumk1lt Δcik 6 When the holding cost at customer i is negligibly small that is hCi 0 for all i in N the case for our data then conditions i and ii in Lemma 2 become i Δci t t wit hP tP tP wit hP fCit for all t t ii Δci t t wit hP tP tP wit hP fCit for all t t Lemma 3 Expanding on the notation introduced in Lemma 2 let Lt i zCit 1 fCit ci0 be the set of customers who receive a delivery in period t and whose corresponding cost is not ci0 For i in Lt let Δci ci0 ci be the increase in cost that would result if customer i was assigned the arc connected to the depot in a solution instead of ci Assume that yt nt and let lt yt nt be the number of customers whose delivery cost should be set to ci0 Now identify the lt smallest values of Δci for i in Lt that is Δci1 Δci2 Δcilt where i1 ilt in Lt Assuming that the minimum production setup cost fP minft t 0 1 τ 1 is much larger than the cost c0 of using a vehicle in any period if i Δci t t wit hP tP tP wit hP t t wit hCi fCit for all t t and ii Δci t t wit hP tP tP wit hP t t wit hCi fCit for all t t then the lower bound φLB can be adjusted as follows φLB φLB sumt1τ sumk1lt Δcik 7 Lemma 4 Define vt as the minimum number of vehicles needed in period t Let i argminci0 i in N and let Δcit c0i c0 be the incremental cost for customer i in Ni in period t for not using the minimum cost arc with cost c0 from the depot in a solution For period t identify the vt 1 smallest values of Δcit and the corresponding customers The following adjusted lower bound is valid for the modified allocation model φLB φLB sumt1τ sumk1vt 1 Δcikl 8 Proof The number of outbound arcs from the depot to customers in each period must be at least vt in any feasible solution to the full model Because each customer can be visited only once in each period this implies that the outbound arcs from the depot cannot be the same for all vehicles By selecting the vt 1 smallest incremental costs associated with the arcs leaving the depot in period t we increase φLB by sumk1vt1 cikt while simultaneously reducing it by c0 vt 1 Summing over the planning horizon gives a valid lower bound on the total vehicle costs because we are only adjusting the cost for one less than the minimum number of vehicles required in each period The bound improvement in 8 can be calculated in a preprocessing step The other improvements given in 4 6 and 7 can be obtained directly by solving the allocation model after several additional modifications are made For 4 it would be necessary to replace ceilingsumi wit Q in the objective function with a nonnegative integer variable yt t in T and set the righthand side of 1h to Q yt For 6 it would be necessary to introduce n τ binary variables ξit i in N t in T to identify whether cost coefficient ci0 or ci is to be used for customer i in period t and n τ constraints of the form xi ξit 1 to ensure that at most one coefficient is selected A similar modification would be required for 7 Testing has shown that these additions can measurably extend computation times so we have relied on Lemmas 14 to improve the initial bound 6 Computational results All computations were performed on a 253 GHz processor with 512 MB of RAM The optimization models and the tabu search code were implemented in Java Netbean 41 and linked to the CPLEX 81 libraries CPU times were obtained through both CPLEX and the time function in Java For testing purposes we used the three data sets provided by Boudia et al 2007 containing 30 instances of 50 100 and 200 customer problems all with a 20period planning horizon and holding costs hP 1 hCi 0 for all i in N These instances 272 J Sched 2009 12 257280 Table 1 Comparison of solutions for problem instances with 50 customers and 20 periods Prob no φgrasp tgrasp φPh1 tPh1 φPh1 tALLOC φPh2 Iteration tPh2 φPh2 Best Runtime Phase 1 Runtime φgrasp Alloc Phase 2 φPh2 Runtime s φPh1 GRASP s solution s φgrasp 100 time s solution found φPh1100 solution 1 440505 10236 418570 918 5 56 399125 25 22645 5 2 448695 13817 391366 1665 13 130 373581 29 38555 5 3 419730 9535 385897 2018 8 1625 353058 22 23416 9 4 456398 6842 401851 1817 12 170 361309 19 21641 10 5 434466 9937 396977 2029 8 1389 365035 16 72557 8 6 452564 9807 382417 1668 15 146 368082 50 46714 4 7 436812 9890 388935 2021 11 1497 369963 30 40473 5 8 420935 8732 383705 1807 9 1591 370822 30 32934 3 9 434789 14254 391442 1385 10 1525 379401 26 39290 3 10 436221 15840 388957 1530 11 1537 370805 25 40805 5 11 433890 8177 384722 2084 11 1783 357107 30 31656 7 12 452705 8583 382746 1990 15 1678 355199 45 58590 7 13 440771 9985 381645 1774 13 1458 366547 20 34097 4 14 419412 8445 399040 2014 5 1723 364115 35 33542 9 15 453875 8667 403862 1841 11 1443 367659 30 55509 9 16 457310 9170 377530 1916 17 1568 360534 25 49380 5 17 455663 9108 405292 1684 11 1333 398442 20 15388 2 18 441685 8853 404254 1925 8 160 368600 39 5746 9 19 418896 10427 394187 1752 6 1407 377073 32 48867 4 20 452183 9412 403547 1637 11 1343 372141 39 63594 8 21 409677 7827 393013 1675 4 1361 374743 16 16018 5 22 429116 10826 380357 2061 11 1761 347449 50 103196 9 23 443184 10699 387351 1652 13 1276 362619 31 79467 6 24 426113 10140 388221 2034 9 1637 375022 23 23229 3 25 462245 8699 386524 1945 16 1601 374926 19 27334 3 26 442029 8215 397620 2029 10 1665 366733 36 81139 8 27 444695 8569 385085 2030 13 1670 375261 12 15674 3 28 449894 18746 388354 2284 14 1931 373155 26 23068 4 29 461555 9393 400043 1909 13 1562 379320 25 54383 5 30 434006 9373 397217 2000 8 1655 369223 33 49571 7 The final steps of our methodology involved path relink ing and solving the modified version of the allocation model 1 to obtain a lower bound φLB The results for the 30 customer instances are reported in Table 2 Path relinking is a procedure used to explore the opportunity to improve the solution obtained from any methodology that provides mul tiple candidate solutions It is based on the fact that paths be tween solutions give rise to neighborhoods that contain new solutions with attributes similar to those of the endpoints Our algorithm explores paths between all pairs of elite solu tions that were uncovered during tabu search An elite solu tion is defined as the first improved solution following any unimproved solution The second column of Table 2 lists the number of elite solutions found which averaged 44 Columns 35 give the path relinking solution the corresponding runtime and the gap with the phase 2 solution respectively In most cases the latter is either negative or negligibly small calling into question the effort required to apply this procedure The lower bound φLB obtained from the modified version of the allocation model is reported in column 6 The adjusted lower bound φALB is presented in column 7 and was ob tained by applying 4 68 On average it showed an improvement over the lower bound φLB of 05 with a standard deviation of 123 The gap between the best so lution and the adjusted lower bound averaged 128 and is J Sched 2009 12 257280 273 Table 2 Path relinking and lower bound results for instances with 50 customers and 20 periods Prob no No φPR τPR φPR φLB φALB φBest φLP tLP φALB elites Solution Time φPh2 from LB Adjusted φALB LP soln LP time φLP solns from path s φPh2 100 modela LB φALB100 s φLP100 relinking 1 8 398795 180 008 324463 326630 22 218243 472 50 2 2 373374 71 006 325469 326324 14 219206 469 49 3 6 353058 460 0 326365 329607 7 217697 556 51 4 6 361176 300 004 328090 328467 10 220366 584 49 5 5 364819 290 006 324451 325044 12 220178 38 48 6 3 368082 222 0 329875b 332449 11 221852 53 50 7 9 369963 720 0 326712b 328176 13 218911 408 50 8 1 370822 0 324958 325893 14 216969 692 50 9 2 379379 36 001 325737 326277 16 218951 738 49 10 5 370655 260 004 325846b 326483 14 217392 492 50 11 7 354025 540 086 329208 330539 7 217596 1045 52 12 5 354981 180 006 320151b 326316 9 215297 633 52 13 3 365432 208 030 326512 327017 12 217946 638 50 14 3 363404 184 020 321141 323498 12 216514 474 49 15 4 367659 56 0 325857 326339 13 219531 569 49 16 1 360534 0 326324 328384 10 216828 533 51 17 5 398442 120 0 328513 329996 21 221842 594 49 18 5 368533 189 002 323263 324394 14 219581 408 48 19 6 377255 480 005 326077 327712 15 216825 566 51 20 9 372361 670 006 325181 326112 14 218650 431 49 21 2 375228 30 013 326764 327367 14 218483 936 50 22 7 347329 810 003 322818 324015 7 214260 661 51 23 4 362619 220 0 326345b 327824 11 219942 544 49 24 2 375609 100 016 320055b 321520 17 214406 498 50 25 3 374682 163 007 330338b 331818 13 220062 566 51 26 4 366167 129 015 331722 334062 10 220034 516 52 27 1 375261 0 321615 322816 16 216087 514 49 28 3 373464 120 008 323233 324047 15 215604 412 50 29 5 379320 255 0 334784b 336159 13 223175 511 51 30 6 370012 420 021 331837 335212 10 218849 633 53 Only one elite solution aStop at 500 s bLower bound model includes integer variables yt Lemma 1 not satisfied shown in column 8 With respect to size the lower bounding MIP contained 2083 constraints and 3183 variables of which 1021 were binary In all cases a 500 s limit was placed on CPLEX for these runs which terminated with an optimality gap that averaged 22 The last three columns of Table 2 give the LP lower bound φLP for the full model with the routing constraints the corresponding solution times and the gaps between the LP solution and the adjusted lower bound φALB obtained from the modified allocation problem in 500 s The size of this gap which was about 50 implies that the solution to the LP relaxation by itself is not very useful As discussed by Laporte 1992 and others poor performance like this is due primarily to the weak relaxation associated with the subtour elimination constraints of the form yjt yit wit D t max 1 xijt i N j N0 t T where yit is the load on the vehicle just after departing from customer i in period t and D t max minQ iN t l t dil is an upper bound on the load on any vehicle in period t These constraints require the load on a vehicle to be monotoni Springer 274 J Sched 2009 12 257280 Table 3 Comparison of solutions for problem instances with 100 customers and 20 periods Prob no φgrasp tgrasp φPh1 tPh1 φPh1 tALLOC φPh2 Iteration tPh2 φPh2 Best Runtime Phase 1 Runtime φgrasp Alloc Phase 2 φPh2 Runtime s φPh1 GRASP s solution s φgrasp100 time s solution found φPh1100 solution 1 790972 41342 737241 11844 7 9650 711671 21 1079 3 2 782906 50619 734043 14327 6 1246 698512 40 1035 5 3 787830 34676 721767 11725 8 9543 683270 40 1299 5 4 779847 48651 741904 11951 5 9531 718252 21 672 3 5 796176 43906 748370 11833 6 9487 731260 14 335 2 6 793216 34902 758969 14154 4 11268 744927 15 1143 2 7 781317 29864 729426 12622 7 10539 695728 36 1250 5 8 780884 49919 729542 14980 7 11944 706058 30 1008 3 9 784646 44230 742733 13016 5 10822 705035 39 1125 5 10 790156 37854 727431 14340 8 12122 696521 31 985 4 11 787596 39390 734881 12643 7 10071 711895 37 835 3 12 785170 36925 730530 13465 7 11385 703162 32 1129 4 13 777705 50597 742061 13613 5 11188 721066 22 599 3 14 789802 46769 730899 11689 7 9530 698548 35 1065 4 15 790132 52702 732911 15743 7 13664 711506 23 1139 3 16 797322 41819 753956 11516 5 9120 714873 37 1226 5 17 799843 52011 729914 12523 9 10611 702314 33 1218 4 18 787371 41909 752665 15207 4 12681 720238 32 720 4 19 806592 35369 773547 11363 4 8587 748734 36 1349 3 20 809340 40309 748000 12380 8 10175 729099 20 1131 3 21 788736 47735 754214 12712 4 10228 738746 14 544 2 22 804538 41297 735752 14492 9 12360 702849 36 998 4 23 781558 42912 741379 12609 5 10182 712717 23 1037 4 24 798428 41688 758063 14242 5 11380 727741 37 1380 4 25 796591 36814 745669 12521 6 10285 725869 30 1240 3 26 791514 40369 724380 12491 8 10494 700719 37 585 3 27 773662 49551 709640 11837 8 9877 686382 36 454 3 28 780492 41689 724556 14313 7 11752 700980 30 1190 3 29 799417 45751 754116 13521 6 11031 725030 32 993 4 30 785906 39811 734725 15412 7 13415 698942 32 1300 5 cally decreasing during the delivery sequence and are rarely if ever binding in the LP solution The corresponding MIP for the full model contained 54063 constraints and 54183 variables of which 51021 were binary 63 Medium and large instances The results for the 100customer instances with 20 time pe riods each are presented in Tables 3 and 4 The column head ings are identical to those of Tables 1 and 2 respectively On average the GRASP took 421 s while tabu search took a to tal of 11338 s or 174 longer The gap between the GRASP solution and the phase 1 solution obtained from the alloca tion model after running the VRP subroutine is identified in column 5 and represents roughly a 64 improvement An additional 36 is realized in phase 2 giving a total im provement of 10 As seen in Table 4 path relinking takes about 465 s and in only for problem no 2 was an improvement realized On average the path relinking solutions were 12 worse than the phase 2 solutions The modified allocation model had 4083 constraints and 6283 variables of which 2021 were bi nary A total of 850 s was allotted for the computations At termination CPLEX exhibited an optimality gap of roughly 35 The gap between the best feasible solution and the adjusted lower bound is shown in column 8 and averaged 229 The data in the last three columns are as expected J Sched 2009 12 257280 275 Table 4 Path relinking and lower bound results for instances with 100 customers and 20 periods Prob No φPR tPR φPR φLB φALB φBest φLP tLP φALB no elites Solution Time s φPh2 from LB Adjusted φALB LP soln LP time φLP solns from path φPh2100 modela LB φALB100 s φLP100 relinking 1 1 711671 000 577137b 582090 22 281684 2418 107 2 4 694694 240 055 581028b 584847 19 284882 2202 105 3 2 693600 80 151 578440 581317 18 283960 3605 105 4 1 718252 000 571459b 574748 25 283138 2875 103 5 1 731260 000 579452 583117 25 287534 2146 103 6 1 744927 000 581760 586735 27 286161 4387 105 7 3 708958 540 190 571178 573184 21 283575 2728 102 8 3 713267 720 102 561926b 564981 25 278847 3479 103 9 2 713779 360 124 580835 583295 21 286051 2454 104 10 2 710427 270 200 573489b 577928 21 284019 3282 103 11 3 729324 700 245 574951b 580831 23 284048 2095 104 12 6 712832 550 138 577810 581643 21 285148 2062 104 13 5 728004 650 096 568636 572655 26 281696 2324 103 14 2 705304 500 097 578958b 586741 19 282791 2175 107 15 3 720050 560 120 577760 580764 23 286629 2273 103 16 3 719974 627 071 576178 581498 23 283683 3038 105 17 2 725757 160 334 587829b 591956 19 289069 1140 105 18 2 731189 140 152 577616b 580835 24 285723 2217 103 19 3 759505 780 144 573342 576318 30 283849 2352 103 20 3 737988 610 122 589667b 593649 23 290279 3305 105 21 1 738746 000 584901 590280 25 289479 9478 104 22 2 719634 400 239 574965 577327 22 285958 2514 102 23 1 712717 000 567679 569334 25 281909 2382 102 24 3 738178 720 143 574446 577506 26 285820 2057 102 25 6 739029 590 181 573169 576299 26 284631 2170 102 26 3 707568 520 098 567682b 571351 23 280753 2337 104 27 2 695961 200 140 569319 572684 20 283016 2537 102 28 2 710866 400 141 564894 569282 23 280900 4195 103 29 2 740172 520 209 585442 588796 23 288490 6383 104 30 2 712385 320 192 577862b 582331 20 286316 4736 103 Only one elite solution aStop at 850 s bLower bound model includes integer variables yt Lemma 1 not satisfied The LP relaxation of the full model solves quickly with CPLEX but the gap between the solution φLP and the ad justed lower bound φALB averaged 104 The results for the 200customer instance with 20 time periods each exhibited the same patterns as the 100customer instances Table 5 presents the comparisons with GRASP and Ta ble 6 presents the path relinking and allocation model re sults The average improvement of tabu search over GRASP was 3 but average runtimes increased from 1903 to 2879 s or 51 However as the number of customers increases so do the objective function values so small percentage reduc tions in cost often translate into be large reductions in ab solute terms From Table 6 we see that path relinking solutions are on average 05 worse than the phase 2 solutions so it is doubt ful whether the procedure is worthwhile especially with runtimes averaging 2230 s To compute the lower bound from the modified allocation model CPLEX was allowed 1450 s At the time the average optimality gap was 36 276 J Sched 2009 12 257280 Table 5 Comparison of solutions for problem instances with 200 customers and 20 periods Prob φgrasp tgrasp φPh1 tPh1 φPh1 tALLOC φPh2 Iteration tPh2 φPh2 no Best Runtime Phase 1 Runtime φgrasp Alloc Phase 2 φPh2 Runtime s φPh1 GRASP s solution s φgrasp100 time s solution found φPh1100 solution 1 1075528 207839 1051861 280 2 179 1030684 20 2965 2 2 1070340 178575 1038017 320 3 233 1024558 11 1737 1 3 1070505 168852 1036408 265 3 131 1036282 4 1528 0 4 1068959 219475 1067202 289 0 174 1057654 14 2190 1 5 1060220 167819 1046100 266 1 189 1031422 18 2573 1 6 1065700 194855 1052150 360 1 204 1033233 19 1016 2 7 1091538 179672 1061457 393 3 206 1043536 35 1942 2 8 1060164 223727 1072610 349 1 222 1066068 10 2978 1 9 1055447 155442 1047509 321 1 190 1036179 16 2454 1 10 1069590 209038 1057926 377 1 214 1038559 38 1855 2 11 1069280 246475 1054845 388 1 197 1037705 21 3428 2 12 1057631 190958 1052501 397 1 221 1040220 29 2288 1 13 1074180 208702 1075537 483 0 233 1063024 13 2719 1 14 1076460 226477 1054986 397 2 216 1041786 21 3544 1 15 1065340 177681 1045066 376 2 155 1029908 25 2554 2 16 1067550 202239 1052812 410 1 231 1033656 27 3200 2 17 1067007 182645 1053600 382 1 187 1027433 31 2800 3 18 1095350 171627 1089858 377 1 172 1063306 29 4740 2 19 1063445 192098 1079858 372 2 166 1065705 23 2797 1 20 1049854 191066 1057882 360 1 148 1034195 37 2600 2 21 1055436 250695 1065159 477 1 205 1044771 32 3900 2 22 1066185 178139 1059796 369 1 152 1045790 13 2822 1 23 1073265 220173 1033344 459 4 214 1027042 15 1797 1 24 1063585 199914 1078993 381 1 151 1051610 37 3510 3 25 1054230 190287 1055464 410 0 202 1027772 36 3000 3 26 1057443 168531 1059502 335 0 129 1044315 26 2720 1 27 1076798 1484 1060084 409 2 213 1047267 12 2498 1 28 1054225 150817 1052961 417 0 184 1042891 16 4200 1 29 1088853 146353 1040105 409 5 239 1030156 24 1947 1 30 1051195 161392 1051980 325 0 155 1035703 21 2351 2 indicating the increased difficulty in solving the correspond ing IP The lower bound lemmas improved the results by 061 on average not shown in table The gap between the best solution found and the adjusted lower bound was ap proximately 205 slightly better than for the 100customer problems For problem instances of this magnitude these gaps are within an acceptable range 7 Discussion and future directions Good solutions to the PIDRP can yield substantial benefits throughout the supply chain as manufacturers and customers become more integrated A decade ago the primary chal lenges facing manufacturers were establishing online com munications and delivery channels Keeping too much in ventory was discouraged and warehouses equaled waste Today there is a shift in attitude that allows for increased inventory levels where necessary Lean inventory is a lux ury many companies can no longer afford because timely deliveries depend on a far wider range of factors than can be predicted or controlled If your supply chain is global those factors include international transportation providers as well as infrastructure with varying degrees of quality and importexport regulations If your supply chain is domestic you are still affected because certain transportation lanes are J Sched 2009 12 257280 277 Table 6 Path relinking and lower bound results for instances with 200 customers and 20 periods Prob No φPR tPR φPR φLB φALB φBest φLP tLP φALB no elites Solution Time s φPh2 from LB Adjusted φALB LP soln LP time φLP solns from path φPh2100 modela LB φALB100 s φLP100 relinking 1 1 1030684 0 853521 859151 20 442266 12061 94 2 2 1010158 1250 141 851442b 856489 18 440625 12181 94 3 2 1016681 1200 189 843605 847677 20 438910 12388 93 4 4 1042854 2141 140 852993 858620 21 440440 13246 95 5 3 1023680 2423 075 853415 858986 19 440789 7834 95 6 2 1025262 2205 077 859458 863514 19 443041 75 95 7 3 1038746 2610 046 855652 859755 21 442400 11127 94 8 1 1066068 0 856121b 862333 24 442292 7840 95 9 2 1018420 1250 171 838995 845013 21 434188 8301 95 10 8 1035240 2625 032 851118 854605 21 441476 7688 94 11 2 1037767 2730 001 861136 866478 20 444524 10401 95 12 5 1035350 2242 047 856766 861683 20 443012 10569 95 13 1 1063024 0 851622 855760 24 442760 7843 93 14 2 1024491 1920 166 856951 860535 19 443699 10396 94 15 6 1026787 2408 030 848998 855899 20 440104 7598 94 16 2 1043917 1957 099 855689 859910 20 442252 10119 94 17 5 1022250 2018 050 853635 859039 19 440699 8985 95 18 5 1065250 2242 018 854938b 860527 24 441318 9058 95 19 1 1065705 0 859526 867188 23 442439 8206 96 20 4 1027134 2425 068 847312 852918 20 437861 7392 95 21 3 1049028 3003 041 849209 853126 22 440479 8013 94 22 1 1045790 0 848907 854260 22 439190 6717 95 23 4 1034198 2403 070 855591 862561 19 441020 7351 96 24 6 1045014 2908 063 854029 858675 22 439128 9447 96 25 3 1024239 2601 034 859551 864960 18 442202 10296 96 26 4 1043128 1673 011 848613b 855543 22 438277 9342 95 27 2 1030753 2037 158 863295b 869693 19 445644 7968 95 28 3 1032478 2909 100 854208 859167 20 440210 9327 95 29 5 1019371 2179 105 862445 866925 18 444642 9906 95 30 6 1027915 2400 075 859578 865404 19 441228 8831 96 Only one elite solution aStop at 1450 s bLower bound model includes integer variables yt Lemma 1 not satisfied more crowded than ever and the competition for transporta tion services has become even more intense In this paper we have provided an efficient reactive tabu search algorithm for finding high quality solutions to the PIDRP as measured by the optimality gap for small in stances and by the solution to our lower bounding model otherwise Although we included a path relinking feature and several theoretical ways of improving the lower bound none was very effective If there is room for improvement in the methodology the place to begin is with the alloca tion model 1 and its modified version that was used to ob tain lower bounds Adapting the cut generation procedures of Pochet and Wolsey 2006 for the capacitated lotsizing problem offers several opportunities for tightening the LP relaxation of 1a1j A second option that is now under way is to apply branch and price to the full PIDRP using a time period decomposition Because a VRP must be solved for each period there are some inherent limitations to this approach A third option is to cluster the customers first and then solve the PIDP and routing components separately 278 J Sched 2009 12 257280 Each of these ideas offers computational challenges that the research community has yet to address Appendix Adjustment of production levels As mentioned in Sect 45 it may be possible to improve the tabu search results when evaluating candidate moves by adjusting the production levels In the implementation a lin ear program is solved every 5 iterations to find the optimal values of pt In the remaining cases we adjust production levels by calling ProductionLevelAdjustmentAlgorithm The inputs to this algorithm are the production levels pt as sociated with the most recent tabu search solution and the output generated by FindAdjustmentPeriodAlgorithm which is called when checking the feasibility of a solution For a move that involves customer i1 in period t1 and customer i2 in period t2t1 t2 let ηi1t1be the number of items that were originally scheduled to be delivered to cus tomer i1 in period t1 that have been rescheduled for delivery in period t2 and let wi2t2be the number of items to be de livered to customer i2 in period t1 after the move The logic included in the primary algorithms called in the adjustment of production levels is given below The first step is to check feasibility of a candidate move Only the inputs and outputs are stated for this algorithm more detail can be found in Nananukul 2008 FeasibilityCheckAlgorithm complexity Oτ n Input Tabu search solution pt at current iteration periods t1 and t2t1 t2 customer i1 and i2 demand of customer i1 for all t t1 demand of customer i2 from all t t2 delivery amounts wi1tt t1 and wi2tt t2 and swap amounts ηi1t1 and wi2t2 Output Output flag true feasible or false not fea sible If output flag true then return number of de crease adjustment periods decrease adjustment pe riods decrease adjustment amounts number of in crease adjustment periods increase adjustment pe riods and increase adjustment amounts comment The return values are the output from FindAdjustmentPeriodAlgorithm FindAdjustmentPeriodAlgorithm complexity Oτ Input Period t search type 1 search for periods t t such that production level needs to be de creased 2 search for periods t t such that the production level needs to be increased adjust ment quantity AQ and tabu search solution output pt at current tabu search iteration comment t t1 or t2 AQ ηi1t1 wi2t2 Only positive value of adjustment quantity AQ is used in the algorithm Output Feasible flag if feasible periods is found then re turn 1 otherwise return 1 If search type 1 then return number of de crease adjustment periods decrease adjustment pe riods and decrease adjustment amounts otherwise return number of increase adjustment periods in crease adjustment periods and increase adjustment amounts NumAdjPeriods 1 For t t 10 If pt 0 comment Only periods t with pt 0 are considered in the adjustment If Search type 1 TAQ AQ pt else TAQ AQ I P max pt TAQ stores remaining adjustment quantity after considering adjustment in period t end if If TAQ 0 In this case it means period t is the last period that needs to be considered AdjPeriodsNumAdjPeriods t AdjAmountNumAdjPeriods AQ NumAdjPeriods NumAdjPeriods 1 Set feasible flag 1 If search type 1 Store results in NumDecAdjPeriods DecAdjPeriodsn n 1 NumDecAdjPeriods DecAdjAmountn n 1 NumDecAdjPeriods else Store results in NumIncAdjPeriods IncAdjPeriodsn n 1 NumIncAdjPeriods IncAdjAmountn n 1 NumIncAdjPeriods end if return else AdjPeriodsNumAdjPeriods t If Search type 1 AdjAmountNumAdjPeriods pt else AdjAmountNumAdjPeriods I P max pt end if NumAdjPeriods NumAdjPeriods 1 end if end if J Sched 2009 12 257280 279 AQ TAQ If TAQ 0 In this case it means that the adjustment is not feasible Set feasible flag 1 and return end if After the feasibility check is done the movevalue for each feasible candidate is calculated by calling MoveValueAlgorithm Because the logic is straightfor ward we only give the inputs and outputs MoveValueAlgorithm complexity On3 Input Periods t1 and t2t1 t2 costs of current VRP so lutions in period t1 and t2 CVRP t1 CVRP t2 customer i1 and i2 and swap amount ηi1t1 ie number of items rescheduled from period t1 to period t2 for customer i1 number of items that were to be de livered to customer i2 in period t2 but have been rescheduled for delivery in period t1 wi2t2 adjust ment quantity AQ number of decrease adjustment periods decrease adjustment periods decrease ad justment amounts number of increase adjustment periods increase adjustment periods and increase adjustment amounts Output movevalue Once all candidate moves are evaluated the one with the best movevalue is selected Finally the ProductionLevel AdjustmentAlgorithm is called ProductionLevelAdjustmentAlgorithm Input Tabu search solution output pt at current tabu search iteration Number of decrease adjustment periods decrease adjustment periods decrease ad justment amounts number of increase adjustment periods increase adjustment periods and increase adjustment amounts Output Updated values of pt Step 1 Update production levels pt in each period t where the new tabu search solution indicates a decrease in production level For n 1 NumDecAdjPeriods pDecAdjPeriodsn pDecAdjPeriodsn DecAdjAmountn Step 2 Update production levels pt in each period t where the new solution indicates an increase in production level For n 1 NumIncAdjPeriods pIncAdjPeriodsn pIncAdjPeriodsn IncAdjAmountn Example Continued Considering the example in Fig 3 assume that the current solution from the allocation model gives p1 66 p2 0 and p3 0 Base on a swap move in Fig 3 i1 3i2 1t1 2 and t2 3 respectively Following Step 1 of the above algorithm NumDecAdjPe riods 1 DecAdjPeriods0 1 DecAdjAmount0 4 NumIncAdjPeriods 1 IncAdjPeriods0 1 IncAd jAmount0 4 Continuing p1 is updated to 66 4 62 In Step 2 p1 is updated to 62 4 66 As a result the swap move does not change the production level in period 1 For the example in Fig 4 assume that the current so lution from the allocation model gives p1 64 p2 2 and p3 0 Based on the transfer move in Fig 4 i1 0i2 2t1 2 and t2 3 respectively Following Step 1 of FeasibilityCheckAlgorithm NumDecAdjPeriods 2 DecAdjPeriods0 2 DecAdjAmount0 2 DecAdjPe riods1 1 DecAdjAmount1 2 NumIncAdjPeriods 1 IncAdjPeriods0 1 IncAdjAmount0 4 Fol lowing Step 1 of ProductionLevelAdjustmentAlgorithm p2 is updated to 2 2 0 In Step 2 p1 is updated to 64 2 4 66 References Abdelmaguid T F Dessouky M M 2006 A genetic algorithm approach to the integrated inventorydistribution problem Inter national Journal of Production Research 4421 44454464 Anily S Federgruen A 1993 Twoechelon distribution systems with vehicle routing costs and central inventories Operations Re search 411 3747 Bard J F Huang L Jaillet P Dror M 1998 A decomposition approach to the inventory routing problem with satellite facilities Transportation Science 322 189203 Boudia M Louly M A O Prins C 2006 A memetic algorithm with population management for a productiondistribution prob lem In A Doglui G Morel CE Pereira Eds 12th IFAC sym posium on information control problems in manufacturing Vol 3 pp 541546 SaintEtienne France 1719 May Boudia M Louly M A O Prins C 2007 A reactive grasp and path relinking for a combined productiondistribution problem Computers Operations Research 3411 34023419 Carlton B W Barnes J W 1996 Solving the traveling salesman problem with time windows using tabu search IIE Transactions on Scheduling Logistics 288 617629 Cetinkaya S Mutlu F Lee CY 2006 A comparison of outbound dispatch policies for integrated inventory and trans portation decisions European Journal of Operational Research 1713 10941112 Chandra P Fisher M L 1994 Coordination of production and distribution planning European Journal of Operational Research 723 503517 Dror M Ball M 1987 Inventoryrouting reduction from an an nual to a short period problem Naval Research Logistics Quar terly 344 891905 Federgruen A Tzur M 1999 Timepartitioning heuristics appli cation to one warehouse multiitem multiretailer lotsizing prob lems Naval Research Logistics 465 463486 Gaudioso M Paletta G 1992 A heuristic for the periodic vehicle routing problem Transportation Science 262 8692 280 J Sched 2009 12 257280 Glover F Laguna M 1997 Tabu search Norwell Kluwer Golden B Assad A Dahl R 1984 Analysis of a large scale ve hicle routing problem with an inventory component Large Scale Systems 723 181190 Gutiérrez J SedeñoNoda A Colebrook M Sicilia J 2007 A polynomial algorithm for the productionordering planning problem with limited storage Computers Operations Research 344 934937 Herer Y T Tzur M Yucesan E 2006 The multilocation trans shipment problem IIE Transactions on Scheduling Logistics 38 185200 Laporte G 1992 The vehicle routing problem An overview of ex act and approximate algorithms European Journal of Operational Research 593 345358 Lei L Liu S Ruszczynski A Park S 2006 On the inte grated production inventory and distribution routing problem IIE Transactions on Scheduling Logistics 3811 955970 Mourgaya M Vanderbeck F 2007 Column generation based heuristic for tactical planning in multiperiod vehicle routing Eu ropean Journal of Operational Research 1833 10281041 Nananukul N 2008 Lotsizing and inventory routing for a productiondistribution supply chain PhD dissertation Graduate Program in Operations Research Industrial Engineering The University of Texas Austin OBrien C Tang O Eds 2006 Integrated enterprise and sup ply chain management International Journal of Production Eco nomics 101 special issue Parthanadee P Logendran R 2006 Periodic product distribu tion from multidepots under limited supplies IIE Transactions on Scheduling Logistics 3811 10091026 Pochet Y Wolsey L A 2006 Production planning by mixed integer programming New York Springer Resende M G C Ribeiro C C 2005 GRASP with path relinking recent advances and applications In T Ibaraki K Nonobe M Yagiura Eds Metaheuristics progress as real problem solvers pp 2963 New York Springer Villegas F A Smith N R 2006 Supply chain dynamics analysis of inventory vs order oscillations tradeoff International Journal of Production Research 446 10371054 Zhao QH Wang SY Lai K K 2007 A partition approach to the inventoryrouting problem European Journal of Operational Research 1772 786802 A integração produçãoestoquedistribuiçãoroteirização problema O artigo aborda o desafio de integração das decisões de produção e distribuição para otimização da cadeia de suprimentos O objetivo central é coordenar produção estoque e entrega de forma a atender à demanda do cliente minimizando os custos correspondentes O método utilizado apresenta um modelo que inclui uma unidade de produção clientes com demanda variável um horizonte de planejamento finito e uma frota de veículos para as entregas O problema é resolvido por meio de um procedimento de busca tabu reativa com aplicação de realocação de caminho para melhorar os resultados O uso de um modelo de alocação auxilia na busca de soluções viáveis Testes computacionais demonstraram a eficácia da abordagem com melhorias de 10 a 20 em relação a procedimentos existentes O artigo ressalta que boas soluções para o problema têm o potencial de gerar benefícios substanciais para a cadeia de abastecimento permitindo maior integração entre fabricantes e clientes Em relação ao futuro do trabalho há a sugestão de melhorias na metodologia explorando o modelo de alocação e sua versão modificada para obter limites inferiores A adaptação dos procedimentos de geração de corte e a aplicação de ramificação e preço são opções em andamento para apertar o relaxamento do problema Outra opção é agrupar os clientes antes de resolver os componentes de roteamento e o problema de produção e distribuição separadamente o que pode proporcionar uma abordagem mais eficiente para lidar com o problema em períodos de tempo diferentes Uma abordagem evolutiva para a otimização multiobjetivo do problema de rede de distribuição integrada de localização e estoque em estoque gerenciado pelo fornecedor O artigo aborda o uso do Inventário Gerenciado pelo Fornecedor VMI como uma solução emergente para melhorar a eficiência da cadeia de suprimentos O objetivo central é formular um problema integrado de rede de distribuição de estoque de localização sob a configuração do VMI É apresentado um modelo chamado MOLIP MultiObjective LocationInventory Problem e investigase a aplicação de um algoritmo evolutivo multiobjetivo baseado no NSGA2 NonDominated Sorting Genetic Algorithm II para resolver o MOLIP Os resultados dos experimentos computacionais mostraram que a abordagem proposta é inovadora e eficiente para resolver problemas difíceis de tamanho prático O estudo apresenta o modelo MOLIP que considera os efeitos da localização distribuição e questões de estoque em um sistema de inventário gerenciado pelo fornecedor VMI É proposto um algoritmo evolutivo híbrido baseado no conhecido algoritmo NAGAII com estratégia de elitismo e mecanismo de classificação não dominado Foram realizados dois experimentos O primeiro investigou a aplicação do algoritmo evolutivo NSGAII para resolver o modelo MOLIP obtendo resultados promissores O segundo experimento comparou a abordagem proposta com o SPEA2 indicando que a primeira teve melhor desempenho em termos de qualidade e diversidade das soluções aproximadas O estudo sugere a adaptação do algoritmo híbrido proposto para outros sistemas de localização estoque e distribuição além de explorar outras técnicas de otimização existentes Para pesquisas futuras o artigo sugere a adaptação do algoritmo proposto para sistemas com estruturas de rede diferentes a inclusão de outras decisões de inventário como políticas de reabastecimento e estoque de segurança e a investigação de outras técnicas de otimização existentes Também são mencionadas possíveis direções de pesquisa como explorar MOEAs MultiObjective Evolutionary Algorithms mais competitivos e técnicas de computação inteligente suave além da hibridização com outras técnicas de busca heurística como escalada ou reparo local A integração produçãoestoquedistribuiçãoroteirização problema O artigo aborda o desafio de integração das decisões de produção e distribuição para otimização da cadeia de suprimentos O objetivo central é coordenar produção estoque e entrega de forma a atender à demanda do cliente minimizando os custos correspondentes O método utilizado apresenta um modelo que inclui uma unidade de produção clientes com demanda variável um horizonte de planejamento finito e uma frota de veículos para as entregas O problema é resolvido por meio de um procedimento de busca tabu reativa com aplicação de realocação de caminho para melhorar os resultados O uso de um modelo de alocação auxilia na busca de soluções viáveis Testes computacionais demonstraram a eficácia da abordagem com melhorias de 10 a 20 em relação a procedimentos existentes O artigo ressalta que boas soluções para o problema têm o potencial de gerar benefícios substanciais para a cadeia de abastecimento permitindo maior integração entre fabricantes e clientes Em relação ao futuro do trabalho há a sugestão de melhorias na metodologia explorando o modelo de alocação e sua versão modificada para obter limites inferiores A adaptação dos procedimentos de geração de corte e a aplicação de ramificação e preço são opções em andamento para apertar o relaxamento do problema Outra opção é agrupar os clientes antes de resolver os componentes de roteamento e o problema de produção e distribuição separadamente o que pode proporcionar uma abordagem mais eficiente para lidar com o problema em períodos de tempo diferentes O Problema de Roteamento de Estoques um olhar sobre a literatura O artigo em questão é uma revisão da produção de conhecimento sobre o problema de roteamento de estoques PRE que surge devido à globalização e ao aumento da competitividade nas empresas O objetivo central do artigo é revisar a literatura existente na área do PRE com foco nos aspectos de horizonte de análise e modelagem da demanda O método utilizado é uma análise dos conceitos envolvidos no PRE e das diferenças entre o PRE e o roteamento de veículos PRV Os artigos pesquisados são examinados em relação ao horizonte de planejamento e à modelagem da demanda considerando a abordagem específica de cada um O artigo sugere caminhos futuros de pesquisa nesse campo Dois fatores são considerados fundamentais para definir a modelagem e a solução do problema o horizonte de análise e a modelagem da demanda O horizonte de análise referese ao período de tempo considerado para otimização sendo comum a utilização de um horizonte de tempo finito próximo de aplicações práticas A modelagem da demanda varia desde uma demanda constante até uma modelagem probabilística com alguns estudos utilizando dados reais para validar a metodologia O artigo destaca a solução de um problema prático de roteamento de estoques aplicado a uma distribuidora de gás que utiliza um modelo de programação linear e um algoritmo de planejamento de entregas Os autores mostram que o procedimento proposto é superior ao procedimento usado pela empresa em termos de desempenho O artigo não pretende esgotar o assunto mas sim oferecer subsídios teóricos para a compreensão do tema e contribuir para uma agenda de pesquisa em andamento dos autores A revisão da literatura revela que o tema do roteamento de estoques tem sido estudado há mais de quinze anos mas as aplicações práticas não são o foco principal Existem diversos métodos e algoritmos propostos desde adaptações das técnicas de roteamento de veículos e controle de estoques até propostas inovadoras com uso de fluxo de caixa A pesquisa está evoluindo para modelos que integram elementos logísticos de forma mais abrangente aproveitando a disponibilidade de tecnologia e ferramentas de análise adequadas O Problema de Roteamento de Estoques um olhar sobre a literatura O artigo em questão é uma revisão da produção de conhecimento sobre o problema de roteamento de estoques PRE que surge devido à globalização e ao aumento da competitividade nas empresas O objetivo central do artigo é revisar a literatura existente na área do PRE com foco nos aspectos de horizonte de análise e modelagem da demanda O método utilizado é uma análise dos conceitos envolvidos no PRE e das diferenças entre o PRE e o roteamento de veículos PRV Os artigos pesquisados são examinados em relação ao horizonte de planejamento e à modelagem da demanda considerando a abordagem específica de cada um O artigo sugere caminhos futuros de pesquisa nesse campo Dois fatores são considerados fundamentais para definir a modelagem e a solução do problema o horizonte de análise e a modelagem da demanda O horizonte de análise referese ao período de tempo considerado para otimização sendo comum a utilização de um horizonte de tempo finito próximo de aplicações práticas A modelagem da demanda varia desde uma demanda constante até uma modelagem probabilística com alguns estudos utilizando dados reais para validar a metodologia O artigo destaca a solução de um problema prático de roteamento de estoques aplicado a uma distribuidora de gás que utiliza um modelo de programação linear e um algoritmo de planejamento de entregas Os autores mostram que o procedimento proposto é superior ao procedimento usado pela empresa em termos de desempenho O artigo não pretende esgotar o assunto mas sim oferecer subsídios teóricos para a compreensão do tema e contribuir para uma agenda de pesquisa em andamento dos autores A revisão da literatura revela que o tema do roteamento de estoques tem sido estudado há mais de quinze anos mas as aplicações práticas não são o foco principal Existem diversos métodos e algoritmos propostos desde adaptações das técnicas de roteamento de veículos e controle de estoques até propostas inovadoras com uso de fluxo de caixa A pesquisa está evoluindo para modelos que integram elementos logísticos de forma mais abrangente aproveitando a disponibilidade de tecnologia e ferramentas de análise adequadas Uma abordagem evolutiva para a otimização multiobjetivo do problema de rede de distribuição integrada de localização e estoque em estoque gerenciado pelo fornecedor O artigo aborda o uso do Inventário Gerenciado pelo Fornecedor VMI como uma solução emergente para melhorar a eficiência da cadeia de suprimentos O objetivo central é formular um problema integrado de rede de distribuição de estoque de localização sob a configuração do VMI É apresentado um modelo chamado MOLIP MultiObjective LocationInventory Problem e investigase a aplicação de um algoritmo evolutivo multiobjetivo baseado no NSGA2 NonDominated Sorting Genetic Algorithm II para resolver o MOLIP Os resultados dos experimentos computacionais mostraram que a abordagem proposta é inovadora e eficiente para resolver problemas difíceis de tamanho prático O estudo apresenta o modelo MOLIP que considera os efeitos da localização distribuição e questões de estoque em um sistema de inventário gerenciado pelo fornecedor VMI É proposto um algoritmo evolutivo híbrido baseado no conhecido algoritmo NAGAII com estratégia de elitismo e mecanismo de classificação não dominado Foram realizados dois experimentos O primeiro investigou a aplicação do algoritmo evolutivo NSGAII para resolver o modelo MOLIP obtendo resultados promissores O segundo experimento comparou a abordagem proposta com o SPEA2 indicando que a primeira teve melhor desempenho em termos de qualidade e diversidade das soluções aproximadas O estudo sugere a adaptação do algoritmo híbrido proposto para outros sistemas de localização estoque e distribuição além de explorar outras técnicas de otimização existentes Para pesquisas futuras o artigo sugere a adaptação do algoritmo proposto para sistemas com estruturas de rede diferentes a inclusão de outras decisões de inventário como políticas de reabastecimento e estoque de segurança e a investigação de outras técnicas de otimização existentes Também são mencionadas possíveis direções de pesquisa como explorar MOEAs MultiObjective Evolutionary Algorithms mais competitivos e técnicas de computação inteligente suave além da hibridização com outras técnicas de busca heurística como escalada ou reparo local