Computers & Operations Research 36 (2009) 2311--2319 Contents lists available at ScienceDirect Computers&OperationsResearch journal homepage: www.elsevier.com/locate/cor Productionschedulingandvehicleroutingwithtimewindowsforperishable foodproducts Huey-KuoChena,Che-FuHsuehb,∗,Mei-ShiangChangc aDepartment of Civil Engineering, National Central University, Taoyuan 32001, Taiwan bDepartment of Marketing and Distribution Management, Ching Yun University, Taoyuan 32097, Taiwan cDepartment of Civil Engineering, Chung Yuan Christian University, Taoyuan 32023, Taiwan A R T I C L E I N F O A B S T R A C T Available online 7 October 2008 Keywords: Perishable Production scheduling Vehicle routing Time windows We propose a nonlinear mathematical model to consider production scheduling and vehicle routing with time windows for perishable food products in the same framework. The demands at retailers are assumed stochastic and perishable goods will deteriorate once they were produced. Thus the revenue of the supplier is uncertain and depends on the value and the transaction quantity of perishable products when they are carried to retailers. The objective of this model is to maximize the expected total profit of the supplier. The optimal production quantities, the time to start producing and the vehicle routes can be determined in the model simultaneously. Furthermore, we elaborate a solution algorithm composed of the constrained Nelder–Mead method and a heuristic for the vehicle routing with time windows to solve the complex problem. Computational results indicate our algorithm is effective and efficient. © 2008 Elsevier Ltd. All rights reserved. 1. Introduction Perishable goods, such as food products, vegetables, flowers, living animals and ready-mix concrete, often deteriorate during the production and delivery processes. More and more suppliers adopt just-in-time production and delivery strategy to fulfill their orders from retailers because they can reduce the loss of their own profit due to deterioration of perishable goods. Perishable goods involve two different concepts of deterioration. First, all items become simultaneously obsolete at the end of the planninghorizon.Second,itemsdeterioratethroughouttheplanning horizon. This category can be further divided into two classes: (1) items with a fixed shelf life such as blood or ready-mix concrete; (2) items with continuous decay such as food, vegetables, flowers, living animals and so on. Among these perishable goods, fresh food products usually deteriorate rapidly. The value or quality of perishable food products will decrease rapidly once they are produced and will keep decaying when being delivered. The revenue of the food supplier will depend ontheconditionoftheproductswhentheyarereceived.Thus,timely production and delivery of perishable foods significantly affect the supplier's revenue. The rapid and continuous decaying characteristic ∗ Corresponding author. Fax: +88634683994. E-mail address: [email protected] (C.-F. Hsueh). 0305-0548/$-see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2008.09.010 of food products makes the suppliers postpone the time to produce and deliver the products as fast as possible. Furthermore, it is important nowadays that perishable foods must be delivered within allowable delivery times, or time-windows. If a vehicle arrives late, a penalty cost may be incurred. Duetotheaforementionedfacts,anintegratedandwell-designed production schedule and delivery routes must be made so that the supplier can ensure the provision of freshest foods and satisfy customers' requirements in a cost-effective way. However, little has been done to investigate the coordination of production scheduling and delivery planning for perishable goods. Traditional production scheduling focuses on the determination of schedules for production such that some performance measures are optimized without considering the delivery plan. On the other hand, the distribution planning usually focuses on minimizing the transportation cost of finished goods and will not affect the production schedule. Therefore, the coordination of production scheduling and delivery planning becomes an important issue in perishable food industry and urgently needs further studies. In this study, we assumed the value or quality of perishable food products will decrease throughout their lifetimes. The rate of deterioration is given and fixed. A supplier produces perishable food productsandsellstoretailers,whosedemandsareassumedstochastic and follow some known probability distributions. We consider both production scheduling and vehicle routing with time windows for perishable food products, under stochastic demands, in the same framework. The supplier decides production quantities, time to start 2312 H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 producing and the delivery routes to maximize the expected total profit. The remainder of the paper is organized as follows. Section 2 features a brief literature review, after which, in Section 3, the production scheduling and vehicle routing problem with time windows for perishable goods (PS-VRPTW-P) is formulated as an integer nonlinear programming model. In Section 4, an elaborated solution algorithm is proposed, and in Section 5, numerical results are provided. The paper is concluded in Section 6. 2. Literature review Most of the literature regarding perishable goods is in the fields of pricing, return policies and inventory control to a retailer (see [1–4]). Papers discussing production scheduling and/or distribution of perishable goods are relatively rare. In the following paragraphs, we will review some literature in three classes: production scheduling of perishable goods, production scheduling with delivery and vehicle routing for food or perishable goods. First we review literature about production scheduling of perishable goods. Arbib et al. [5] consider a three-dimensional matching model for perishable production scheduling, which is studied under two independent aspects: (i) the relative perishability of products and (ii) the feasibility of launching/completion time. Garcia et al. [6] consider a ready-mix concrete production and distribution problem, in which the selection of orders to be processed by a readymix concrete manufacturing plant has to be made and orders have a fixed due date and are delivered directly to the customer site. Silva etal.[7]alsostudytheproductionanddistributionproblemofreadymix concrete. Sana et al. [8] consider a volume flexible manufacturing system for a deteriorating item with a time-dependent demand rate, allowing shortages in inventory. Entrup et al. [9] develop three mixed-integer linear programming models that integrate shelf-life issues into production planning and scheduling. Literatureaboutproductionschedulingwithdeliveryisdiscussed as follows. Zdrzalka [10] deals with the single-machine scheduling problem in which each job has a processing time and a delivery time; the objective is to find a sequence of jobs which minimizes the time by which all jobs are delivered. Cheng et al. [11] study the single-machine scheduling problem with batch deliveries, in which jobs in a batch are delivered to the customer together upon the completion time of the last job in the batch; the objective is to minimize the sum of the batch delivery and job earliness penalties, which derive from the difference between the delivery time of the batch and the job completion time. However, the travel time is not considered. Chang and Lee [12] also study the machine scheduling and product delivery, in which jobs require different amounts of storage space during delivery; the objective is to find a schedule for processing jobs such that the time required for all jobs to be processed and delivered to the respective customers is minimized. Most of above research assume that a job is delivered directly to the customerandthedeliverytimeisfixed.Lietal.[13]developasinglemachine scheduling model that incorporates the routing decisions of a delivery vehicle which serves customers at different locations. The objective is to minimize the sum of job arrival times. Finally, we focus on the research about vehicle routing for food or perishable goods. The well-known vehicle routing problem with/without time windows (VRPTW/VRP) has been discussed in many forms. However, not many papers consider VRPTW/VRP for perishable goods until recent years. Hwang [14] develops an effective distribution model for determining optimal patterns of food supply and inventory allocation for famine relief area. The proposed vehicle routing problem incorporating inventory allocation and optimal distribution aims at minimizing the amount of pains and starving instead of travel distance or time. Prindezis et al. [15] present an application service provider that would offer services of distribution logistics for central food market enterprises by appropriately solving the fleet management problem. Tarantilis and Kiranoudis [16] propose a threshold-accepting-based algorithm to solve the fresh milk distribution problem. The problem they formulated is essentially a heterogeneous fixed fleet vehicle routing problem and no information about perishable goods is considered. Hsu et al. [17] consider the randomness of the perishable food delivery process and construct a stochastic VRPTW model to obtain optimal delivery routes, loads, fleet dispatching and departure times for delivering perishable food from a distribution center. Naso et al. [18] consider the problem of scheduling the production and distribution activities of a network of plants supplying rapidly perishable materials. They propose a strategy that combines genetic algorithms and schedule construction heuristics for job scheduling and truck routing. Osvald and Stirn [19] develop a heuristic algorithm for the distribution of fresh vegetables in which the perishability represents a critical factor. The problem was formulated as a VRPTW with time-dependent travel times. The model considers the impact of the perishability as part of the overall distribution costs. To the best of our knowledge, it seems that no research simultaneously consider both production scheduling and vehicle routing with time windows for perishable goods under stochastic demands. 3. Model formulation In this section we propose a mathematical model of PS-VRPTWP, in which a supplier has to decide how many and when he/she should start producing perishable food products as well as the optimal vehicle route to deliver the products to retailers. The objective of the supplier is to maximize its expected total profit. The supplier has only one production line but produces several kinds of perishable products. In other words, it produces one kind of products at a time. Given only one production line available, the batch production containing several products for a to-be-dispatched vehicle is arranged one product by one product in order and priority is always given to the product with lower decay rate. After all the products for the to-be-dispatched vehicle are produced, the production of products for the next vehicle can then begin. The value of perishable food products will decay once they are produced,anddifferentkindsofperishableproductsmayhavedifferent decay rates. These perishable products need to be sent to several retailers, which may have different time window requirements (i.e., VRPTW for perishable products). In this paper, soft time windows are considered. Any vehicle that arrives late will incur a penalty. Any vehicle that arrives early has to wait until the beginning of the time window. In addition, each vehicle has a capacity constraint. The demand quantity of a retailer is assumed to be a random variable with a probability density function and is realized when a vehicle arrives at the retailer. If the supply quantity is less than the actual demand quantity, there will be a cost of goodwill loss. The probability density function can be any function with a positive domain and an integral equal to one. Other assumptions in formulating the PS-VRPTW-P are given as follows: 1. The price of products (without deterioration) is known and given. 2. All products on the same vehicle are produced continuously as the same batch and the commodity with lower decay rate in that batch is produced earlier. 3. Thevehicledepartsatthetimewhenallproductsassignedtothat vehicle have been produced and loaded. 4. The setup cost for producing different kinds of products is not considered. H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 2313 commodity 1 commodity 2 commodity 1 commodity 2 t1,1 t1,2 ˆ t1,1 ˆ t2,1 t2,1 t2,2 ˆ t1,2 time manufacturer ˆ t2,2 arrival time a3 departure time of vehicle 1 departure time of vehicle 2 41 3 2 5 0 0 0 0 vehicle 1 vehicle 2 Fig. 1. Illustration of production schedule and vehicle routing plan. Tofacilitatetheunderstandingofourmodel,consideranexample with two kinds of perishable goods, five retailers and two vehicles, wherevehicle1servesretailers1,3,4andvehicle2servesretailers2 and 5. Fig. 1 depicts the production schedule and the vehicle routing plan of a feasible solution. Note that the vehicle's departure time will be affected by (a) production quantities and (b) which retailers the vehicle serves. For easy reference, notations that will be used in this paper are listed as follows. Decision variables: Qin production quantity of commodity n for retailer i, tk,1 time to start producing first commodity for vehicle k and xijk binary decision variable, taking a value of 1 if the arc (i, j) is part of a vehicle route k and 0 otherwise. Non-decision variables and parameters: ai time arriving at retailer i, cij travel time between retailers i and j, cn production cost of unit commodity n, ei starting time of the time window at retailer i, fin probability density function of demand n at retailer i, g1 goodwill loss for shortage of supply, g2 goodwill loss for violation of time windows, K total number of vehicles available, li ending time of the time window at retailer i, N number of kinds of commodities, pin original price of commodity n at retailer i, Qveh capacity for all vehicles, R number of retailers, si service time at retailer i, tk,n time to start producing commodity n for vehicle k, ˆ tk,n ending time of producing commodity n for vehicle k, tn production time of unit commodity n, in value decayed per unit quantity of commodity n when retailer i receives the commodity, n value decayed per unit time per unit commodity n and time value of travel time. The PS-VRPTW-P is now formulated as follows: max {Qin,xijk,tk,1}∈ z= in Qin 0 (pin −in)fin()d + ∞ Qin [(pin −in)Qin −g1(−Qin)]fin()d − in cnQin− ijk cijxijk−g2 i max{ai −li,0} (1) where is the feasible region represented by the following constraints. Flow conservation constraints: j k xijk =1 i=1,...,R (2) i k xijk =1 j=1,...,R (3) i xihk − j xhjk =0 h=1,...,R; k=1,...,K (4) j x0jk1 k=1,...,K (5) xijk ∈{0,1} i,j=0,...,R, ij; k=1,...,K (6) Production scheduling constraints: ˆ tk,N tk+1,1 k=1,...,K−1 (7) ˆ tK,N l0 (8) Vehicle capacity constraints: ijn QjnxijkQveh k=1,...,K (9) Definitional constraints: in =[max{ai,ei}−1 2(tk,n +ˆ tk,n)]n if i is on route k ∀i,n (10) aj = max{ai,ei}+si +cij if xijk =1 i,j=1,...,R, ij; k=1,...,K (11) aj =ˆ tk,N +c0j if x0jk =1, j=1,...,R; k=1,...,K (12) ˆ tk,n =tk,n +tn ij Qjnxijk n=1,...,N; k=1,...,K (13) tk,n+1 =ˆ tk,n n=1,...,N−1; k=1,...,K (14) Nonnegative constraints: Qin0 i=1,...,R; n=1,...,N (15) t1,10 (16) 2314 H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 tk,n tk,n ˆ max{ai, ei} Qin time quantity Fig. 2. Quantity variation of commodity n in time for retailer i. Objective (1) maximizes the expected total profit of the supplier. Thefirsttermintheobjectivefunctionistheexpectedrevenuewhen the retailer's demand is less than the supplier's supply quantity; the second term is the expected revenue minus goodwill loss when theretailer'sdemandismorethanthesupplier'ssupplyquantity;the third and fourth terms are production cost and transportation cost, respectively; the last term is a penalty, which accrues if the arrival time of the serving vehicle exceeds a retailer's time window limit. Constraints (2)–(6) are flow conservation constraints of a traditional vehicleroutingproblem.Constraint(2)requiresthatonlyonevehicle can leave from retailer i once. Constraint (3) denotes that only one vehicle can arrive at retailer j once. Constraint (4) states that for each retailer h, the entering vehicle must eventually leave this node. Constraint (5) designates that each vehicle can leave the supplier once at most. Constraint (6) designates xijk as a 0–1 integer variable. Constraints (7) and (8) state that the ending time of producing all commodities for a vehicle must be earlier than the starting time of producing commodities for the next vehicle or the ending time of the supplier's time window. Constraint (9) is a vehicle capacity constraint. Constraints (10)–(14) are definitional constraints. Constraint (10) calculates the average value decayed when retailer i receives commodity n. Fig. 2 shows the quantity variation of commodity n on vehicle k, which has been produced and assigned to retailer i. Since the value of a commodity will decay as soon as it has been produced, the average value decayed per unit commodity n when retailer i receives the commodities at time max{ai,ei}is equal to in= (shadowed area∗n)/Qin, which derives constraint (10). Constraints (11) and (12) define the arrival time at retailer j. NotethatConstraint(11)alsoguaranteestheeliminationofsubtours (i.e., cycles disconnected with the depot) because it implies that the arrival time at a latter node must be greater than the arrival time at the previous node, which is not possible for all nodes in a subtour. Constraint (13) defines the ending time of producing commodity n for vehicle k. Constraint (14) requires that the starting time of producing commodity n+1 is equal to the ending time of producing commodity n for the same vehicle k. Constraints (15) and (16) are nonnegative constraints. Themodelweformulatedaboveisanintegernonlinearprogramming model embedding a vehicle routing problem, which is known to be NP-hard. 4. Solution algorithm Theproposedintegernonlinearprogramming modelis verydifficulttosolve.Thuswedecomposethedecisionvariables{Qin,tk,1,xijk} into two groups: {Qin,tk,1} and {xijk}. The first group is associated with a production scheduling problem and the second group is subject to a VRPTW. With the concept of decomposition, the PS-VRPTW-P model (1) can be rearranged as follows: Upper level: max {Qin,tk,1}∈1 z= in Qin 0 pinfin()d + ∞ Qin [pinQin −g1(−Qin)]fin()d − in cnQin −ZVRPTW (17) where 1 is the feasible region represented by nonnegative constraints (15) and (16) and ZVRPTW is calculated as follows: Lower level: min xijk∈2 ZVRPTW = in Qin 0 infin()d + ∞ Qin inQinfin()d+ ijk cijxijk +g2 i max{ai −li,0} (18) where 2 is the feasible region represented by constraints (2)–(14) with {Qin,tk,1} given. If a given {Qin,tk,1} causes problem (18) to be infeasible, simplylet ZVRPTW equal infinity. The original model (1) is now converted to a nonlinear programming model (17) with nonnegative constraints and a VRPTW in the objective function. Model (17) can be solved using either a sensitivity-analysis based or a direct search algorithm. The former uses sensitivity analysis to obtain the derivative information of the reaction function (either explicitly or implicitly) while the latter employs only function evaluations. Since the interdependence between production variables{Qin,tk,1}and vehicle routes{xijk}are too complicated and the derivative information is not available in this problem, we adopted a direct search algorithm to solve the problem. One of the most widely used direct search methods for solving nonlinear unconstrained optimization problems is the Nelder–Mead simplex algorithm (see [20]). In the next two subsections, the Nelder–Mead method with boundary constraints is adopted to solve the production scheduling problem (17) and a heuristic is proposed to solve the VRPTW (18). 4.1. Solving the production scheduling problem TheNelder–Meadmethodisasimple,intuitiveandrelativelystable method for solving nonlinear optimization problems and can be applied to non-convex problems. If the number of variables{Qin,tk1} is v, the method will construct a simplex of v+1 vertices; each vertex represents a feasible solution. Suppose y0,...,yv are vertices of the current simplex in a decreasing objective value order and s denotes the center of the gravity of {y0,...,yv−1}. There are three construction operations to determine a new point: reflection, expansion and contraction. Reflection of yv at s is determined by the point yr =s+(s−yv) with 0<1. (cf Fig. 3a). Expansion of yr is determined by the point ye=s+(yr−s) with >1. (cf Fig. 3b). Three different contraction steps can be distinguished. Parameter with 0<<1 denotes a contraction constant. Interiorcontractionofyv isdeterminedbythepointyc=s+(yv−s).(cf Fig. 3c). Exteriorcontractionofyr isdeterminedbythepointyc=s+(yr−s).(cf Fig. 3d). H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 2315 yv yr s y0 y1 yv yr s y0 1 y ye yv ry s y0 y1 yc yv yr s y0 1y yc yv y1' y0 y1 yv' Fig. 3. Basic steps of the Nelder–Mead method (v = 2): (a) reflection; (b) expansion; (c) interior contraction; (d) exterior contraction and (e) total contraction. Total contraction is that all vertices of the simplex except y0 are replaced by the point yi =y0 +(yi −y0). (cf Fig. 3e). Because the Nelder–Mead method is originally applied to an unconstrained problem, an adjustment is necessary that projects its coordinates on the bounds if the new point is out of the domain. We now proceed to formally state the solution algorithm, as follows: Step 1: Initialization Step 1.1: Find an initial solution{Qin,tk1}(designated as y0) of (17) as follows, and solve the corresponding vehicle routing problem (18). Qin =average demand of commodity n at retailer i tk1 =e0 +(k−1)(l0 −e0)/K The initial production quantity Qin is usually set as the mean value of the stochastic demand quantity. The initial time to start producing first commodity for all to-be-dispatched vehicles is scheduled sequentially with equal time increment in the range of working hours. Calculate the value of objective function (17). Step 1.2: Determine other vertices y1,...,yv of the initial simplex by disturbing y0 as follows: yi =y0 + y0Ii ∀i=1...v where is a turbulence factor and Ii is a unit base vector. Project its coordinates on the bounds, if yi is out of the domain. Solvethecorrespondingvehicleroutingproblem(18)andcalculate the value of objective function (17), respectively. Step 2: Identify the vertices with the highest function value as yh, the vertices with the lowest function value as yl, the vertices with the second lowest function value as ys, the center of gravity of the simplex (without yl) ass, and the corresponding objective function values as z(yh),z(yl),z(ys). Step 3: Apply a reflection with respect to yl: yr =s+(s−yl) 2316 H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 Generate initial solution y0 Generate other solutions y1 ~ yv Solve corresponding VRPTW for each yi to calculate objective value z (yi) Apply reflection, expansion and contraction operations Convergent? End Y N Solve corresponding VRPTW for new yi to calculate objective value z (yi) Fig. 4. Flowchart of solution algorithm. Project its coordinates on the bounds, if yr is out of the domain. Step 4: Update the simplex. We distinguish between three cases: (a)Ifz(yr)>z(yh),itmeansthatthereflectioncreatedabetter solution. We attempt to get an even better point through expansion of yr: ye=s+(yr−s).Projectitscoordinatesonthebounds,ifnecessary. Replace yl with ye if z(ye)>z(yr); otherwise, replace ylwith yr. (b) If z(yr)z(ys), replace yl with yr. (c) If z(yr)z(yl), it was probably wrong to do the reflection along that direction. An interior contraction from yl in direction s−yl will be applied: yc=s+(yl−s),projectitscoordinatesonthebounds,ifnecessary. Else, if z(yl)z(yl), replace yl with yc. Otherwise, a total contraction is performed since all attempts to get improvement failed. yi =yh +(yi −yh) ∀ih Step 5: Check convergence. If the distance between yh and any other vertices is smaller than a certain tolerance , then stop; yh and its corresponding vehicle route is the best solution. Otherwise, go to Step 2. Note that in Step 5 the difference of z(yh)−z(yl) less than a preset tolerance could be another choice of stopping criteria. The flow chart for the solution algorithm is depicted in Fig. 4. 4.2. Solving the VRPTW When {Qin,tk1} are temporarily fixed, the problem (18) reduces to a special VRPTW. There are two significant features which differ from the traditional VRPTW, as follows: a. The vehicle routing schedule will affect the departure time of the vehicles because the more loads a vehicle carries, the more production time the supplier needs and as a result the later it departs. b. In addition to transportation cost, the cost of the value decayed and the goodwill loss due to time window violation are also considered in the objective function (18). As in most VRPTW algorithm, we first construct initial routes and then improve them. A modified insertion method is adopted to construct initial routes. When a retailer is inserted into an existing route, the insertion will postpone the departure time (as mentioned above). As a result, the arrival times and the value decayed at all retailers along this route may change. Suppose O is a set of retailers that compose a route and ZVRPTW(O), which only contains retailers in route O, represents the partial cost of objective function (18). ZVRPTW(O)= i∈O,n Qin 0 infin()d+ ∞ Qin inQinfin()d + i∈O,j∈O,k cijxijk +g2 i∈O max{ai −li,0} (19) Denote the routes before and after insertion as, respectively, O and Onew. The insertion cost can then be represented by ZVRPTW(Onew)−ZVRPTW(O). We now propose the solution algorithm as follows. Step 1: Input data Input the production quantities and the time to production {Qin,tk1} calculated from the production scheduling problem, as well as the transportation cost. Step 2: Insert the retailer with the minimum insertion cost into an existing route or perform a new route Calculate the insertion cost for each retailer. Insert the retailer with the minimum insertion cost into an existing route at the optimal position. If the vehicle capacity constraints cannot be fulfilledortheinsertioncostishigherthanthecostofdispatching a new vehicle, then perform a new route. Repeat Step 2 until all retailers have been assigned to vehicles. Step 3: Route improvement Many existing route improvement methods are available, such as 2-opt and Or-opt [21], which may combine with other metaheuristics. Note that in this vehicle routing problem, the possible objective improvement is mainly due to node (or arc) exchange operations, i.e., insertion or removal of a node, which will inevitably affect the arrival times and the value decayed at all retailers along the route. 4.3. Enhancing the algorithm Note that only the nonnegative constraints are considered in the upperlevelmodel.Thereasonthatproductionschedulingconstraints and vehicle capacity constraints are not considered in model (17) is that they will be affected by the vehicle routes, which will be determined in the lower level model (18). Therefore when{Qin,tk,1} are temporarily fixed in the upper level, it may render the lower level model infeasible and results in inefficiency of the algorithm of the upper level model. To avoid this circumstance and increase the efficiency of the algorithm, we modified both the upper level and lower level problems as follows. In the modified upper level problem, the original objective function(17)ismodifiedsuchthatthestartingproductiontimesbetween two successive batches are in order. H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 2317 Upper level—modified: max {Qin,tk,1}∈1 z= in Qin 0 pinfin()d + ∞ Qin (pinQin −g1(−Qin))fin()d − in cnQin −M K−1 k=1 max{tk,1 −tk+1,1,0} −ZVRPTW (20) where M is a very big number and 1 is the feasible region represented by the nonnegative constraints (15) and (16) and the following constraints: n QinQveh i=1,...,R (21) tK,1l0 (22) (21)isalooservehiclecapacityconstraintand(22)isalooserproduction schedule constraint; both of which contain no vehicle routing variables, {xijk}. These two additional looser constraints along with the modified objective function may result in a better production schedule {Qin,tk,1}, which has less opportunity to render the lower level model infeasible as compared with that of original upper level model (17). The modified lower level problem is simply obtained by moving the production schedule constraints (7) and (8), which are side constraints to a traditional VRPTW, to the original lower level objective function as penalty terms. Lower level—modified: min xijk∈2 ZVRPTW = in Qin 0 infin()d+ ∞ Qin inQinfin()d + ijk cijxijk +g2 i max{ai −li,0} +M⎛ ⎝ K−1 k=1 max{ˆ tk,N −tk+1,1,0}+max{ˆ tK,N −l0}⎞ ⎠ (23) where 2 is the feasible region represented by constraints (2)–(6) and (9)–(14) with {Qin,tk,1} given. 5. Computational results 5.1. Numerical examples The proposed algorithm was coded with Visual C++ 6.0. Retailers' information (including locations, demands, time windows and service time) is created based on Solomon's problem set with some modification [22]. The deterministic demands are extended to a range and the probability density functions of demands are assumed tobeuniformdistributions.Therearethreekindsofperishableproducts, which is associated with different production times, production costs and decay rates. All examples were tested on a personal computer with an Intel Core2 CPU and 1GB RAM. Table 1 shows the Table 1 Test results for different problem sizes. Number of retailers 5 10 15 20 25 50 75 100 Objective value 546.7 1046.5 1607.0 2118.8 2622.1 5219.0 8117.0 10245.2 CPU time (s) 0 1 5 11 21 149 578 3153 test results of different problem sizes, in which the tested retailers are selected from the 100 retailers in the problem set. The results show that the proposed algorithm can solve the PS-VRPTW-P with the problem size up to 75 retailers within ten minutes. In the case of 100 retailers, it takes about 53min. For comparison, we also solve the PS-VRPTW-P with mathematical software LINGO 10.0. Since the complexity of PS-VRPTW-P is NP-hard, only small problems can be solved. Ten numerical examples are created for testing. Case 1 to Case 5 have five retailers and Case 6 to Case 10 have six retailers. The test results are summarized in Table 2. In each case, LINGO 10.0 returns a local optimal solution. Table 2 shows that the solutions of proposed algorithm are better than the local optimal solutions found by LINGO 10.0 in eight cases. In other two cases, the gap from the local optimal solution is less than 3.74%. Moreover, the CPU times of the proposed algorithm are less than 1s in all cases except Case 10, which takes 10s. However, LINGO10.0takeshourstofindalocaloptimalsolutioninmostcases. In Case 8, LINGO 10.0 returns a very poor local optimal solution in less than 12min. Fromtheabovetestresults,wefoundthattheproposedalgorithm can solvePS-VRPTW-Pefficientlyandreturnsareliablesolution.The efficiency of the algorithm makes it suitable for solving a real case, which is normally large in scale. 5.2. Sensitivity analyses Twenty-five retailer problems are adopted for conducting sensitivity analyses in terms of the rate of decay and the fleet size. First, wevariedtherateofdecayofProduct3,3,from0.02to0.14.Theobjective value and the change of supply quantity ratio,iQi3/inQin, against decay rate are drawn in Fig. 5. The result shows that the objective value decreases when the rate of decay increases, while the supply quantity ratio doesn't have an apparent change. Second, we varied the number of vehicles from 3 to 14. For the same fleet size, we also run a series of testing which has the different percentage of retailers with time window requirements, including TW(0%), TW(20%), TW(40%), TW(60%), TW(80%) and TW(100%). For example, TW(0%) refers to no retailers have time window requirements and TW(40%) refers to 40% of retailers have time window requirements. The impact of the fleet size on objective value is summarized in Fig. 6. Fig. 6 shows that the objective value improves when the fleet size increases in all series of testing with different percentages of time window requirements. Wecanfurtherderivetherelationsbetweenaverageloadingratio and the number of vehicles from the above testing and illustrate in Fig. 7. For a given number of vehicles, we calculate the average loading ratio of these vehicles in all cases with different percentages of time-window requirements. The loading ratio is defined as the production quantity assigned to a vehicle divided by the capacity of that vehicle. Fig. 7 shows that the average loading ratio declines when the number of vehicle increases. From Figs. 6 and 7, we found that there is a tradeoff between the objective value and the average loading ratio; this may be due to the fact that the vehicle operating cost is not considered in the objective function. Therefore, using more vehicles can usually decrease the deterioration and increase the objective value but results in a lower average loading ratio. For the problem of perishable goods, a supplier usually cares more about the decay of perishable goods 2318 H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 Table 2 Comparisons between proposed algorithm and LINGO 10.0. Case no. 1 2 3 4 5 This study 596.81(00:00:01)a 790.67(00:00:00) 546.10(00:00:00) 662.34(00:00:00) 569.09(00:00:00) LINGO 10.0 357.01(01:29:59) 403.74(03:11:39) 311.62(01:19:43) 688.04(02:25:20) 572.73(02:07:18) Gap (%) 0 0 0 3.74 0.64 Case no. 6 7 8 9 10 This study 593.90(00:00:00) 629.87(00:00:00) 714.76(00:00:00) 614.52(00:00:00) 938.32(00:00:10) LINGO 10.0 254.11(06:51:51) 499.61(12:02:06) −718781.7(00:11:35) 550.47(01:17:28) 350.39(06:55:42) Gap (%) 0 0 0 0 0 aObjective value (CPU time: hh:mm:ss). 0 1000 2000 3000 decay rate objective Supply quantity ratio objective supply quantity ratio 0.252 0.250 0.248 0.246 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Fig. 5. Objective value and change of supply quantity ratio versus decay rate. TW (0%) TW (20%) TW (40%) TW (60%) TW (80%) TW (100%) 0 2000 2500 3000 3500 3 4 5 6 7 8 9 10 11 12 13 14 vehicles objectives Fig. 6. Objective values under different number of vehicles and time windows. than the average loading ratio and is willing to use more vehicles to prevent the deterioration. The more vehicles are used, the less time of production and delivery for each vehicle is required and therefore it results in the less deterioration. 6. Conclusions and suggestions Inthispaper,wehavesuccessfullyformulatedthePS-VRPTW-Pas a mixed-integer nonlinear programming model and proposed a solution algorithm associated with it. Computational results indicated that our algorithm is effective and efficient. Other contributions and some concludsions are summarized as follows: Production scheduling and vehicle routing for perishable goods are integrated into a unified framework which is also applicable to the fields like food, vegetables, flowers, living animals and so on. The PS-VRPTW-P is tactically decomposed into two mutually dependent “simpler” subproblems and each subproblem can be readily solved by some existing algorithms. 0 0.2 0.4 0.6 0.8 3 4 5 6 7 8 9 10 11 12 13 14 vehicles average loading ratio Fig. 7. Relations between average loading ratio and the number of vehicles. Even with a small number of retailers, solving the PS-VRPTWP by means of optimization software may take a very long time. Apparentlysuch optimization softwareisnotsuitable forsolvingPSVRPTW-P. However, the proposed algorithm in this paper is much more efficient and effective in solving PS-VRPTW-P. The number of vehicles is an important factor in PS-VRPTW-P. To reducing deterioration, one may use more vehicles to transport perishable goods. However, the increased number of vehicles may lead to a higher transportation cost. In the future study, the proposed PS-VRPTW-P model can be further improved by including the setup and vehicle operating costs into the model. In addition, different kinds of deterioration, such as nonlinear or fuzzy nature, may be considered in the PS-VRPTW-P. References [1] Pasternack BA. Optimal pricing and return policies for perishable commodities. Marketing Science 1985;4:166–76. [2] Raafat F. Survey of literature on continuously deterioration inventory models. Journal of Operational Research Society 1991;42:27–37. [3] Abad PL. Optimal pricing and lot-sizing under conditions of perishability and partial backordering. Management Science 1996;42:1093–104. [4] Teng JT, Ouyang LY. An EOQ model for deteriorating items with power-form stock-dependent demand. Information and Management Sciences 2005;16: 1–16. [5] Arbib C, Pacciarelli D, Smriglio S. A three-dimensional matching model for perishable production scheduling. Discrete Applied Mathematics 1999;92:1–15. [6] García JM, Sánchez LS, Guerrero F, Calle M, Smith K. Production and vehicle scheduling for ready-mix operations. Proceedings of the 29th international conference on computers and industrial engineering. Montreal, Canada: 2001. p. 70–6 Http://Www.umoncton.ca/CIE/Conferences/Index.htm. [7] Silva CA, Faria JM, Abrantes P, Sousa JMC, Surico M, Naso D. Concrete delivery using a combination of GA and ACO. Proceedings of the 44th IEEE conference on decision and control, and the European control conference. 2005. p. 7633–38. [8] Sana S, Chaudhuri KS. On a volume flexible production policy for a deteriorating item with time-dependent demand and shortages. Advanced Modeling and Optimization 2004;6:57–74. [9] Entrup ML, G¨unther H-O, Van Beek P, Grunow M, Seiler T. Mixedinteger linear programming approaches to shelf-life-integrated planning and scheduling in yoghurt production. International Journal of Production Research 2005;43(23):5071–100. H.-K. Chen et al. / Computers & Operations Research 36 (2009) 2311--2319 2319 [10] Zdrzalka S. Analysis of approximation algorithms for single-machine scheduling with delivery times and sequence independent batch setup times. European Journal of Operational Research 1995;80:371–80. [11] Cheng TCE, Gordon VS, Kovalyov MY. Single machine scheduling with batch deliveries. European Journal of Operational Research 1996;94:277–83. [12] Chang YC, Lee CY. Machine scheduling with job delivery coordination. European Journal of Operational Research 2004;158:470–87. [13] Li CL, Vairaktarakis G, Lee CY. Machine scheduling with deliveries to multiple customer locations. European Journal of Operational Research 2005;164:39–51. [14] Hwang HS. A food distribution model for famine relief. Computers & Industrial Engineering 1999;37:335–8. [15] Prindezis N, Kiranoudis CT, Kouris DM. A business-to-business fleet management service provider for central food market enterprises. Journal of Food Engineering 2003;60:203–10. [16] Tarantilis CD, Kiranoudis CT. A meta-heuristic algorithm for the efficient distribution of perishable foods. Journal of Food Engineering 2001;50(1):1–9. [17] Hsu CI, Hung SF, Li HC. Vehicle routing problem with time-windows for perishable food delivery. Journal of Food Engineering 2007;80:465–75. [18] Naso D, Surico M, Turchiano B. Scheduling production and distribution of rapidly perishable materials with hybrid GAs. Studies in Computational Intelligence 2007;49:465–83. [19] Osvald A, Stirn LZ. A vehicle routing algorithm for the distribution of fresh vegetables and similar perishable food. Journal of Food Engineering 2008;85: 285–95. [20] Nelder JA, Mead R. A Simplex for Function Minimization. Computer Journal 1965;7:308–13. [21] Or I. Traveling salesman—type combinatorial problems and their relation to the logistics of blood banking. PhD thesis, Northwestern University, Evanston, IL, USA, 1976. [22] Solomon MM. Algorithms for the vehicle routing and scheduling problems with time windows constraints. Operations Research 1987;35:254–65.