An Integrated Model for Production and Distribution Planning of Perishable Products with Inventory and Routing Considerations

In many conventional supply chains, production planning and distribution planning are treated separately. However, it is now demonstrated that they are mutually related problems that must be tackled in an integrated way. Hence, in this paper a new integrated production and distribution planning model for perishable products is formulated. The proposed model considers a supply chain network consisting of a production facility and multiple distribution centers. The facility produces a single perishable product that is storable only for predetermined periods. A homogenous fleet of vehicles is responsible for delivering the product from facility to distribution centers. The decisions to be made are the production quantities, the distribution centers that must be visited, and the quantities to be delivered to them. The objective is to minimize the total cost, where the trip minimization is considered simultaneously. As the proposed formulation is computationally complex, a heuristic method is developed to tackle the problem. In the developed method, the problem is divided into production submodel and distribution submodel. The production submodel is solved using LINGO, and a particle swarm heuristic is developed to tackle distribution submodel. Efficiency of the algorithm is proved through a number of randomly generated test problems.


Introduction
During the past decades, mangers have faced big global changes in the business environment as results of advances in technology, globalization of markets, and new situations in economy and politics.With the increase of the number of world-class competitors, organizations have been under pressure to improve their interorganizational processes.In such situations, managers understood that such changes not only are not enough in long term, but also they must enter in management of their supply, distribution, and after sales companies.With such an approach, the terms "supply chain" and "supply chain management" were created [1].
A supply chain includes all facilities, works, and activities that are engaged in production and delivering a product or service from suppliers (and their suppliers) to customers (and their customers).It includes but not limited to planning and management of material requisition, manufacturing, warehousing, distribution, and after sales services.The supply chain management should coordinate such activities so that customers can have high quality products with minimum possible cost [2].
In any supply chain, two main areas that are very important to the managers are production and distribution planning.Production planning tackles decision of how to transform raw materials into finished products respecting to meet demands on time with minimum cost.Determining lot sizes, that is, calculating the quantity to be producecd for each item at each time, is an important decision in tactical production planning.Considerations such as planning horizon, number of levels, number of products, capacity or resource constraints, deterioration of items, demand type, setup structure, and inventory shortage may complicate the lot sizing problem [3].For detailed information regarding classifications and characteristics of lot sizing problem, see, for instance, Ben-Daya et al. [4], Robinson et al. [5], and Buschkühl et al. [6].

Mathematical Problems in Engineering
On the other hand, distribution planning tackles decision of how to deliver the finished products to the customers respecting to meet their demands on time with minimum cost.Routing problem is among one of the model types in the literature that are applied to the distribution planning.In such problems, distribution of a product(s) from a central location to multiple geographically dispersed customers using definite/indefinite number of capacitated vehicles is considered.The decisions to be made are (1) when to deliver to each customer, (2) how much to deliver to each customer at each time, and (3) how to serve customers using the vehicles [7].
In a conventional supply chain, independent manufacturers, wholesalers, and retailers were separate business entities seeking to maximize their own profits.However, it is now demonstrated that the production and distribution decisions are mutually related problems that need to be dealt with in an integrated way.Integrated production and distribution planning in the context of supply chain management has been under attention of many researchers over the past few years.There might be two primary reasons behind this trend: (1) positive effect of integrated production and distribution planning on profitability of supply chains and (2) positive effect of integrated production and distribution planning on reducing lead times and offering quicker responses to market changes [8].In this paper, we focus on the integrated production and distribution planning of perishable products.The organization of this paper is as follows.In the next section, related literature is reviewed and then in Section 3 the studied problem is thoroughly defined and formulated.As the proposed model is hard to solve specially in large problem sizes, a heuristic method is developed which is efficient in light of both quality of the solutions and the running time.Finally, concluding remarks and future research remarks are given in Section 6.

Literature Review
An integrated production and distribution system often includes facility(s) producing the product(s) and number of distribution centers warehousing products.An integrated production and distribution planning problem is the problem of simultaneously finding the decision variables from different functions that have traditionally been optimized independently.Some of these variables are as follows: (i) quantity of product(s) produced in facility(s) at each planning period, (ii) inventory amount of product(s) temporarily stored in facility(s) at the end of each planning period, (iii) quantity of product(s) shipped from facility(s) to distribution centers in each planning period, (iv) inventory of finished products stored in distribution centers at the end of each planning period.
Due to the number of decision variables to be determined, the integrated production and distribution planning problem is so complex that optimal values are very hard to obtain.
In addition, considerations such as complex structure of the network, geographical span of the supply chain, and involvement of different entities with conflicting objectives can further complicate the problem [9].It is worthy to note that integrated production and distribution planning problems arise in environments where vendor managed inventory (VMI) is implemented.In vendor managed inventory (VMI) environments, the supplier/manufacturer manages inventories at customers to ensure that they do not face any shortages.The integrated production and distribution planning problem has received a lot of attention among researchers particularly in the recent years.Sarmiento and Nagi [10], Chen [11], and Fahimnia et al. [8] provided comprehensive reviews on the general subject.Below we mention some more recent researches.Bard and Nananukul [12] formulated an integrated lot sizing and inventory routing problem as a mixed integer program with the objective of maximizing the net.They developed a two-step procedure that first estimated daily delivery quantities and then solved a vehicle routing problem for each day of the planning horizon.Boudia and Prins [13] studied an NP-hard multiperiod productiondistribution problem to minimize sum of three costs: production setups, inventories, and distribution.This problem then was solved by a memetic algorithm with population management (MAPM).Bard and Nananukul [14] presented a model that included a single production facility and a set of customers with time varying demand.They developed a procedure centering on reactive tabu search for solving the problem.Bilgen and Günther [15] presented an integrated production scheduling and truck routing model for a supply chain of fruit juice.They considered different transportation modes.Bard and Nananukul [16] investigated their previously developed model aimed at minimizing costs.They also developed a hybrid methodology that combined exact and heuristic procedures within a branch-and-price framework.Bolduc et al. [17] proposed a tabu search heuristic for the split delivery vehicle routing problem with production and demand calendars.Jolai et al. [18] considered a supply chain network consisting of a manufacturer, with multiple plants, products, distribution centers, retailers, and customers.They developed a multiobjective linear programming model for integrating production-distribution decisions.They also proposed three metaheuristics to tackle the problem.Toptal et al. [19] considered a joint production and transportation planning problem where two vehicle types were available for outbound shipments.They also presented formulations for three different solution approaches.Nananukul [20] introduced an enhanced clustering model for the latter problem.An algorithm based on a reactive tabu search for solving the clustering problem was also proposed.Liu and Papageorgiou [21] addressed production, distribution, and capacity planning of global supply chains considering cost, responsiveness, and customer service level simultaneously.They developed a multiobjective mixed-integer linear program (MILP) with total cost, total flow time, and total lost sales as key objectives.
Considering perishable products, production/inventory and distribution decisions are often studied separately.Many researchers focused on extending economic production quantity (EPQ) and economic order quantity (EOQ) models for such products.The reader may refer to Nahimias [22], Raafat [23], Goyal and Giri [24], and Bakker et al. [25] for comprehensive reviews.
On the other hand, some researchers focused on extending distribution planning models with concentration on vehicle routing problem for perishable products.Tarantilis and Kiranoudis [26] proposed a fast and robust algorithm to find effective delivery schedules for one of the biggest diary companies in Greece.They formulated the milk delivery problem as a heterogeneous fixed fleet vehicle routing problem (VRP).Chen and Vairaktarakis [27] studied an integrated scheduling model of production and distribution operations in food catering service industries.In their model, a set of customer orders was first processed in a facility and then delivered to the customers directly.The problem was to find a joint schedule of production and distribution such that customer service level and total distribution cost objectives were optimized.Hsu et al. [28] developed a stochastic VRP with time windows formulation for delivery of perishable food products.The goal of their model was to find delivery paths, workloads, and departure times.Similarly, researches like Osvald and Stirn [29], Chen et al. [30], and Gong and Fu [31] applied VRP with time windows to distribution planning of perishable products.Osvald and Stirn [29] formulated a vehicle routing problem with time windows and time-dependent travel times (VRPTWTD) where the travel times between two locations depend on both the distance and on the time of the day.Their model considered the impact of the perishability as part of the overall distribution costs and a tabu search heuristic was used to solve the problem.Chen et al. [30] proposed a nonlinear mathematical model for production scheduling and vehicle routing with time windows for perishable food products.They assumed stochastic demands at retailers.Gong and Fu [31] proposed a multiobjective vehicle routing problem with time window model that includes fixed vehicle cost, operation cost, shelf life loss, and default cost.They applied an ant colony optimization with ABC customer classification (ABC-ACO) to solve the problem.
To the best of our knowledge, no research has ever tackled the integrated problem of production and distribution planning of perishable products with inventory and routing considerations.Therefore, in this paper we deal with an integrated model for production and distribution planning of perishable products.

Problem Statement
In this section, the studied problem is defined with details and the proposed model is formulated.We are given a single production facility that produces a single perishable product.The perishable product has a fixed lifetime that is measured by a number of planning periods that can be stored either in the production facility or in distribution centers.The planning horizon is divided into  equal discrete time interval that each one is a planning period.The production capacity in each planning period is limited to  max .If we have production at the facility in planning period , then a considerable setup cost   ( = 1, . . ., ) is incurred.A limited amount of inventory can be temporarily stored in the production facility with unit holding cost of ℎ 0 .There is a set of  distribution centers geographically dispersed around the production facility that the products are to be delivered to them.Each distribution center  ( = 1, 2, . . ., ) has a nonnegative and deterministic demand  , in planning period  of the planning horizon that must be fully met; that is, shortages are not allowed.A limited amount of inventory can be stored in distribution center  with unit holding cost of ℎ  .A fleet of  capacitated homogeneous vehicles is responsible for shipping of the product from the facility to the distribution centers.The company is not owner of fleet and incurred costs are calculated based of number of trips vehicles make, not the distance they travel.Therefore, we aim to minimize the number of vehicles' trips.In addition, two other rules must be applied to deliveries; each vehicle can make at most one delivery per planning period and each distribution center can be visited at most once per planning period.
Our objective is to construct a production-distribution plan by integrating production, inventory, and routing decisions that minimizes total costs while ensuring that all demands are met and there is no shortage.The model must determine, for each planning period, whether there must be production or not, the amount to be produced, the distribution centers that must be visited, and the quantities to be delivered to them.Considering the setup cost, variable production cost, and inventory holding costs over the planning horizon, the model must decide when to overproduce and when to hold items as inventory.
In standard models, the inventory levels are only restricted by the physical storage capacities, while in our model, the perishability dominates the physical storage capacity.We defined an upper bound for inventory levels at production facility and distribution centers.Considering sl as the shelf life of the perishable product, the inventory upper bound for production facility in planning period  is equal to the sum of all distribution centers' demands in that period and next sl planning periods.Similarly, the inventory upper bound for distribution center  in planning period  is equal to the sum of its demands in that period and next sl planning periods.These upper bounds are our main contribution of the formulation.
In the development of the model, see Nomenclature section.
The formulated model is as follows.
In the above model, objective function (1) minimizes the sum of all costs.Note that the term ∑ ,   ⋅  , aims at minimizing the total number of vehicles trips in planning horizon.Constraints (2) are inventory balance equations at production facility that relate its inventory levels to the production quantities and deliveries to distribution centers.Similarly, constraints (3) are inventory balance equations at distribution centers.Constraints (4) and ( 5) guarantee that inventory levels at production facility or distribution centers are never greater than the total demands in next sl following periods.These constraints guarantee that perishable products will never be thrown away.Production quantities in planning period  are limited to the production capacity using constraints (6).Constraints (7) demonstrate that if there is a delivery to a distribution center in a planning period, then the pertaining binary variable must be 1.Vehicle capacity is respected by constraints (8).Constraints (9) guarantee that each distribution center is visited mostly once per planning period.Similarly, constraints (10) guarantee that each vehicle is used mostly once per planning period.Constraints (11) demonstrate number of vehicles' trips during planning horizon.Finally, relations (12) and (13) present decision variables types.
The size of proposed model is determined largely by constraints (7) and the number of binary variables  ,, , which both grow at a rate proportional to ().A typical problem with 20 vehicles, 50 distribution centers, and a 10day planning horizon contains roughly 12000 constraints, 10000 binary variables, and 10000 continuous variables.Initial attempts to solve instances of this size with LINGO commercial optimizer were not encouraging.This led to development of a more efficient algorithm presented in the next section.

Solution Method
In this section, a two-phase algorithm is presented to solve the proposed model.We face two types of decision variables in the model, the variables pertaining to production decisions (  ,   ,  0, ), and the variables pertaining to distribution decisions ( ,, ,  , ,  ,, ,  , ).A natural way to solve an integrated model is to decompose the problem and consider each part dependently.In this way, we can separate the integrated model into two dependent submodels, the production submodel and the distribution submodel.In proposed method, the production submodel is solved first and the pertaining variables are determined.Then, the results are fed back into the distribution submodel and the pertaining variables are determined.In the next step, the production submodel is solved again considering the previous step.The algorithm continues this iterative solve-feedback-solve procedure until a stopping criterion is met.The details of the proposed method are given in Algorithm 1.
The algorithm consists of two phases, decomposition phase and integration phase.In the decomposition phase, the first iteration of the algorithm, the optimal lot sizes are found considering the sum of distribution centers' demands in each planning period as demand of that period.The production submodel is demonstrated in relations ( 14) to (19).It is a variation of lot sizing model and it is easy to solve by any commercial optimizer such as LINGO.In this submodel, the following additional parameter is used.
Then, the algorithm proceeds with finding the quantity of products to be delivered to each distribution center and the distribution centers that each vehicle must visit in each planning period.This is done through the distribution submodel.Relations (20) to (30) demonstrate the distribution submodel.The production quantities (  ) from production submodel work as input parameter to this submodel.
During initial test, we found that LINGO is not able to find optimal or good quality solution for the above submodel in a reasonable time.Therefore, we used LINGO to obtain initial feasible solution and then a particle swarm based heuristic to improve it.
Particle swarm optimization (PSO) simulates the social behavior of natural organisms such as bird flocking and fish schooling to find a place with enough food.Indeed, in those swarms, a coordinated behavior using local movements emerges without any central control.In PSO, a swarm consists of number of particles.Each particle is a candidate solution to the problem.A particle has its own position and velocity.Optimization takes advantage of the cooperation between the particles.The success of some particles will influence the behavior of the others.Each particle successively adjusts its position according to the following two factors: the best position visited by itself and the best position visited by the whole swarm [32].Originally, PSO has been successfully designed for continuous optimization problems in [33,34].However, Kennedy and Eberhart [35] firstly introduced discrete version of PSO.
The details of proposed heuristic are given in Algorithm 2.
In particle swarm heuristic, how to encode the problem to set of particles is of great importance.Consider the problem in which  distribution centers are to be served by  vehicles in  planning periods; we have to map a 2-dimensional array of (×, ) for each particle.The first dimension includes  sections, where each section has  binary points.The second dimension includes  binary points.If a value is equal to 1, it represents that the corresponding distribution center is served by the relevant vehicle in relevant planning period.An example of encoding structure for problem with five distribution centers, two vehicles, and five planning periods is given in Table 1.
Run LINGO to find an initial feasible solution Set the initial feasible solution as Localbest Solution and Globalbest Solution Generate particles and set all equal to initial feasible solution Set all the particles equal to initial feasible solution While the stopping criterion is not met For each particle , calculate  =  +  ⋅ (  − ) +  ⋅ (  − ).
For each bit   in particle , If (rand < Sigmoid(V  )) then   = 1, else   = 0 Check feasibility for each particle and repair the particle Calculate fitness function for each particle Update Localbest Solution and Globalbest Solution End while Run LINGO using Globalbest Solution to obtain real valued variables Algorithm 2: Details of particle swarm heuristic.Table 1: An example of encoding structure.

Served by vehicle 1
Served by vehicle 2 DC-1 Based on the formulation, the following rules must be respected in each solution: (1) each distribution center must be served at most once per planning period, (2) each vehicle can make at most one delivery per planning period, and (3) the perishability of the product must be respected; that is, time between two consecutive delivery to a distribution center must be at most equal to the shelf life of the product.For instance, if the shelf life of the product is 2 days, a delivery on Monday and another on Thursday is not allowed, because it results in shortage.During the algorithm run, if any of the above rules was violated, the solution becomes infeasible and must be repaired.For rules 1 and 2, if the value of more than one position in the corresponding positions in any particle is 1, we randomly select one position and set its value to 1 and the others to 0. For rule 3, consider that values corresponding to planning periods  and  are 1 in a particle.Without loss of generality, let us assume that  <  + sl < .We set the corresponding value to planning period  = 0 and set the value corresponding to planning period  + sl = 1.
During algorithm run, each particle must be measured according to a fitness function.Because we aim at minimizing the total number of vehicles trips, the term ∑ ,,   ⋅  ,, is used as fitness function.Finally, when the binary variables  ,, and  , are determined via above procedure, we use them as input to distribution subproblem in LINGO to find real variables.
It is clear that approach of decomposing an integrated problem into two submodels and then solving each submodel separately does not necessarily lead to a good quality solution for the integrated problem.So during integration phase that starts from iteration two, new lot sizes and new delivery quantities are found considering the quantities determined in the previous iteration.In a more specific word, new lot sizes are found considering the amount of products delivered to distribution centers in the last distribution plan as demand, instead of the sum of demands of all distribution centers.We call these new demands ARTIFICIAL DEMANDS.Similarly, new delivery quantities are found using the later lot sizes.This iterative process continues until stopping criterion is met, that is, when the outputs of submodels become unchanged.
During the preliminary tests, we often observed a fast convergence of the algorithm in a few iterations.This shows that the algorithm is trapped in the local optimal point.To escape from such point, we implemented a perturbation mechanism in the production submodel.After the computation of such a perturbed lot sizes, the algorithm proceeds with the distribution submodel as usual.Details of perturbation mechanism are given in Algorithm 3.

Computational Results
In this section, the efficiency of the proposed algorithm is discussed with respect to time performance and quality of results.We implemented the proposed algorithm in MAT-LAB 2010 in which DLL of LINGO 8.0 is contained in the code.All test runs of the algorithm are performed on a Celeron(R) CPU 2.40 GHz PC with 512 MB of RAM.

Test Problems.
To the best of our knowledge, no public test problems are available for our problem.So in order to evaluate the algorithm efficiency, we randomly generated test problems.Sizes of test problems, including number of distribution centers, number of planning periods, and number of available vehicles, are given in Table 2.We classified test problems into two groups, small size and large size.We generated model parameters as follows.
(i) Demand quantities were generated from a uniform distribution with lower bound = 10 and upper bound = 20.
(ii) Variable production costs were generated from a uniform distribution with lower bound = 50 and upper bound = 100.
(iii) Fixed production costs were generated from a uniform distribution with lower bound = 500 and upper bound = 1000.
(iv) Inventory holding costs in production facility and distribution centers were generated from a uniform distribution with lower bound = 5 and upper bound = 10.
(v) Fixed transportation costs were generated from a uniform distribution with lower bound = 200 and upper bound = 300.
(vi) Shelf life of the product was set to 2 or 3 periods randomly.

Lower Bound Computation.
A simple way to validate the quality of solutions provided by proposed method is to It is worthy to note that any feasible solution for the original model is also feasible to the relaxed model.In addition, the objective function of the relaxed model is similar to original model, except for the term ∑ ,, V  ⋅  ,, that we substituted to the term ∑ ,   ⋅  , in the original model and always underestimate it.

Results.
To run the proposed method, the algorithm parameters were set according to Table 3.A summary of the results for small size test problems is reported in Table 4.In columns 2-4 of the table, objective function obtained by the proposed method ( Alg ), runtime of the proposed method in seconds, and optimal solution obtained by LINGO ( LINGO ) are reported, respectively.In the final column, the gap between  Alg and  LINGO that is equal to (( Alg −  LINGO )/ Alg ) × 100 is presented.
We observe in Table 4 that the proposed method can provide good solutions in short runtime.The gaps between  Alg and  LINGO are less than 3% in most cases.The worst performance of the method also can obtain a solution gap that is less than 4%.
The proposed method was run 5 times for each large test problem and a summary of the results is reported in Table 5.In columns 2-5 of the table, the average and worst objective functions obtained by the proposed method ( Alg ) and the average and worst runtimes of the proposed method in seconds are reported.In columns 6-8, the solution obtained by relaxed model ( LB ) and the average and worst gap between  Alg and  LB that is equal to (( Alg − LB )/ LB )×100 are reported, respectively.Finally in columns 9-10, the best solution obtained by LINGO within 2 hours run ( LINGO ) and the gap between  Alg and  LB which is equal to (( LINGO −  LB )/ LB ) × 100 are presented, respectively.
The gaps are computed according to the lower bound not to the developed method.In this way, they represent the maximum deviation from the optimal solution.
It is worthy to note that we run LINGO for the test problem number 9 for a night and did not see any significant improvement in results to when we run it for 2 hours.Therefore, we set the LINGO for other test problems for 2 hours.
We observe in Table 4 that the proposed method can provide good solutions in short runtime.The integrality gaps of the proposed method are less than 8.5% in average and less than 12% in worst.This means that the maximum deviation from optimal solution is very reasonable according to the problems complexity.The deviations between the average and worst gaps are about 3% that means the algorithm is very robust to produce results.Comparing the proposed method and LINGO, we see in the table that gap 1 is about half of gap 2 that means the proposed method can provide better solutions than LINGO.In light of the run time, proposed method can provide solutions in less time than LINGO.As stated in the final paragraph of the Section 3, complexity of the proposed model depends on multiplication value of three parameters; number of distribution centers, vehicle, and planning periods.As shown in Figure 1, with increase of this value, run time of the proposed method also increases.

Conclusion and Future Research Remarks
In this paper, a new formulation for integrating production planning and distribution planning of perishable products in a supply chain is presented.Some assumptions about the problem are as follows: the supply chain network includes a single production facility and multiple distribution centers.There is also fleet of homogenous vehicles that are responsible for transporting the product from facility to distribution centers.The company is not owner of the fleet and pertaining incurred costs are fixed and measured based on the number of trip vehicles make not the distance they travel.In each planning period, each vehicle can make one trip at most and each distribution center can be visited at most once.The products are perishable; that is, they are storable only for predetermined periods.The production capacity is limited and shortage in distribution centers is not allowed.The decisions to be made in each planning period are the production quantities, the distribution centers to be visited, and the delivery quantities.Based on the above assumptions, a mixed integer linear program (MILP) is formulated.
Respecting the computational complexity of the model, a heuristic method is developed to tackle it.In the proposed method, LINGO and particle swarm optimization heuristic are combined.In order to validate the performance of proposed method, a number of randomly generated test problems are used.The results are compared with the LINGO and lower bound of the optimal solution.The proposed method is efficient in light of solution quality and time performance.In addition, it is able to provide robust results.In addition, it is shown that time complexity of the proposed method depends on linear multiplication of three parameters, number of distribution centers, vehicle, and planning periods.
Although the studied problem fits into the real world application, it may imposes some limitations into the decision If (NUMBER OF SETUPS > 1) ARTIFICIAL SETUP COST = 2 * SETUP COST ARTIFICIAL HOLDING COST = HOLDING COST/2 End If If (NUMBER OF SETUPS < 2) ARTIFICIAL SETUP COST = SETUP COST/2 ARTIFICIAL HOLDING COST = 2 * HOLDING COST End If Algorithm 3: Details of perturbation mechanism.

Figure 1 :
Figure 1: Average run time of the proposed method.

Table 2 :
Number and size of test problems.

Table 3 :
Algorithm parameters.V  is the variable transportation cost of one item using vehicle  which is equal to   /.The main difference between the original model and the relaxed model is that the binary variables  ,, and  , and pertaining constraints are removed, which makes the relaxed model easy to solve.  ,  0, ,  , ,  ,, ≥ 0 ∀ ∈ ,  ∈ ,  ∈ .

Table 4 :
Computational results for small size test problems.

Table 5 :
Computational results (large size test problems).For example, the problem can be developed when the company owns the vehicle fleet and so the distance each vehicle travels in each trip is important.In this situation, a customized inventory routing model must be developed with respect to problem definition.In addition, other real world assumptions such as inventory transshipments between distribution centers or split and delivery of the inventories can be added to the model.It should be noted that such assumptions complicate the model, which may result in more sophisticated solution methods.Set of planning periods  = 0, 1, ..., : Set of production facility and distribution centers, where 0 corresponds to production facility  = 1, 2, ..., : Set of vehicles.Parameters , : Demand of distribution center  in planning period    : Fixed production cost in planning period  V  : Variable production cost in planning period   max : Production capacity ℎ 0 : Inventory holding cost at production facility ℎ  : Inventory holding cost at distribution center  sl: Shelf life of perishable products which is measured in number of planning periods that the product can be stored  0, : Upper bound of inventory level at production facility in planning period , which is equal to ∑  ∑ ≤≤+sl  ,  , : Upper bound inventory level at distribution center  in planning period , which is equal to ∑ ≤≤+sl  , : Vehiclecapacity   : Fixed transportation cost of using vehicle .Variables   : 1 if there is production on planning period , 0 otherwise   : Production quantity in planning period   0, : Inventory level at production facility in planning period   ,, : 1 if distribution  is served by vehicle  in planning period , 0 otherwise  , : 1 if vehicle  is used in planning period  to serve distribution centers, 0 otherwise  ,, : Amount delivered to distribution center  in planning period  by vehicle   , : Inventory level at distribution center  in planning period .