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

Authors: S. M. Seyedhosseini; S. M. Ghoreyshi
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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.


---

## Body

## 1. 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].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.


## 2. 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 production—distribution 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.


## 3. 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 intoT equal discrete time interval that each one is a planning period. The production capacity in each planning period is limited to p
max
⁡. If we have production at the facility in planning period t, then a considerable setup cost f
p
c
t  (
t
=
1
,
…
,
T
) is incurred. A limited amount of inventory can be temporarily stored in the production facility with unit holding cost of h
0. There is a set of n distribution centers geographically dispersed around the production facility that the products are to be delivered to them. Each distribution center i  (
i
=
1,2
,
…
,
n
) has a nonnegative and deterministic demand d
i
,
t in planning period t 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 i with unit holding cost of h
i. A fleet of K 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. Consideringsl as the shelf life of the perishable product, the inventory upper bound for production facility in planning period t 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 i in planning period t 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.Integrated Production and Distribution Model. Consider
(1)
Min
⁡
∑
t
f
p
c
t
·
Z
t
+
∑
t
v
p
c
t
·
P
t
+
∑
t
h
0
·
I
0
,
t
+
∑
i
,
t
h
i
·
I
i
,
t
+
∑
k
,
t
f
t
c
k
·
Y
k
,
t
(2)
I
0
,
t
=
I
0
,
t
-
1
+
P
t
-
∑
i
∑
k
W
i
,
k
,
t
∀
t
∈
T
(3)
I
i
,
t
=
I
i
,
t
-
1
+
∑
k
W
i
,
k
,
t
-
d
i
,
t
∀
i
∈
N
,
t
∈
T
(4)
I
0
,
t
≤
u
0
,
t
∀
t
∈
T
(5)
I
i
,
t
≤
u
i
,
t
∀
i
∈
N
,
t
∈
T
(6)
P
t
≤
p
max
⁡
·
Z
t
∀
t
∈
T
(7)
W
i
,
k
,
t
≤
q
·
X
i
,
k
,
t
∀
i
∈
N
,
k
∈
K
,
t
∈
T
(8)
∑
i
W
i
,
k
,
t
≤
q
∀
k
∈
K
,
t
∈
T
(9)
∑
k
X
i
,
k
,
t
≤
1
∀
i
∈
N
,
t
∈
T
(10)
∑
i
X
i
,
k
,
t
≤
1
∀
k
∈
K
,
t
∈
T
(11)
∑
i
X
i
,
k
,
t
≤
q
·
Y
k
,
t
∀
k
∈
K
,
t
∈
T
(12)
Z
t
,
X
i
,
k
,
t
,
Y
k
,
t
∈
{
0,1
}
∀
i
∈
N
,
k
∈
K
,
t
∈
T
(13)
P
t
,
I
0
,
t
,
I
i
,
t
,
W
i
,
k
,
t
≥
0
∀
i
∈
N
,
k
∈
K
,
t
∈
T
.In the above model, objective function (1) minimizes the sum of all costs. Note that the term ∑
k
,
t
f
t
c
k
·
Y
k
,
t 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 t 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 X
i
,
k
,
t, which both grow at a rate proportional to O
(
N
K
T
). A typical problem with 20 vehicles, 50 distribution centers, and a 10-day 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.


## 4. 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(
Z
t
,
P
t
,
I
0
,
t
), and the variables pertaining to distribution decisions (
X
i
,
k
,
t
,
Y
k
,
t
,
W
i
,
k
,
t
,
I
i
,
t
). 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.Algorithm 1:Detailed steps of proposed method.
!Decomposition phase
Solve production sub model using LINGO
Solve distribution sub model usingParticle Swarm Heuristic
Save the solution asCURRENT SOLUTION and BEST SOLUTION
CalculateCURRENT COST and set BEST COST = CURRENT COST
!End of decomposition phase
!Integration phase
While the stopping criterion is not met
UpdateARTIFICIAL DEMAND
RunPerturbation Mechanism
Solve production sub model using LINGO consideringARTIFICIAL SETUP COST and
ARTIFICIAL HOLDING COST
Solve distribution sub model usingParticle Swarm Heuristic
Save the solution asCURRENT SOLUTION
Calculate theCURRENT COST
If (CURRENT COST < BEST COST)
SetBEST SOLUTION =CURRENT SOLUTION
SetBEST COST=CURRENT COST
End If
End While
!End of integration phase
ReturnBEST SOLUTIONandBEST COSTThe 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.Additional Parameters.
D
t is the sum of distribution centers’ demands in planning period t
(
D
t
=
∑
i
d
i
,
t
;
∀
t
∈
T
).Production Submodel. Consider
(14)
Min
⁡
∑
t
f
p
c
t
·
Z
t
+
∑
t
v
p
c
t
·
P
t
+
∑
t
h
0
·
I
0
,
t
(15)
I
0
,
t
=
I
0
,
t
-
1
+
P
t
-
D
t
∀
t
∈
T
(16)
I
0
,
t
≤
u
0
,
t
∀
t
∈
T
(17)
P
t
≤
p
max
⁡
·
Z
t
∀
t
∈
T
(18)
Z
t
∈
{
0,1
}
∀
t
∈
T
(19)
P
t
,
I
0
,
t
≥
∀
t
∈
T
.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 (P
t) from production submodel work as input parameter to this submodel.Distribution Submodel. Consider
(20)
Min
⁡
∑
i
,
t
h
i
·
I
i
,
t
+
∑
k
,
t
f
t
c
k
·
Y
k
,
t
(21)
I
0
,
t
=
I
0
,
t
-
1
+
P
t
-
∑
i
∑
k
W
i
,
k
,
t
∀
t
∈
T
(22)
I
i
,
t
=
I
i
,
t
-
1
+
∑
k
W
i
,
k
,
t
-
d
i
,
t
∀
i
∈
N
,
t
∈
T
(23)
I
i
,
t
≤
u
i
,
t
∀
i
∈
N
,
t
∈
T
(24)
W
i
,
k
,
t
≤
q
·
X
i
,
k
,
t
∀
i
∈
N
,
k
∈
K
,
t
∈
T
(25)
∑
i
W
i
,
k
,
t
≤
q
∀
k
∈
K
,
t
∈
T
(26)
∑
k
X
i
,
k
,
t
≤
1
∀
i
∈
N
,
t
∈
T
(27)
∑
i
X
i
,
k
,
t
≤
1
∀
k
∈
K
,
t
∈
T
(28)
∑
i
X
i
,
k
,
t
≤
q
·
Y
k
,
t
∀
k
∈
K
(29)
X
i
,
k
,
t
,
Y
k
,
t
∈
{
0,1
}
∀
i
∈
N
,
k
∈
K
,
t
∈
T
(30)
I
0
,
t
,
I
i
,
t
,
W
i
,
k
,
t
≥
0
∀
i
∈
N
,
k
∈
K
,
t
∈
T
.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 Algorithm2.Algorithm 2:Details of particle swarm heuristic.
Run LINGO to find an initial feasible solution
Set the initial feasible solution asLocalbest 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 particleX, calculate V
=
V
+
α
·
(
X
l
o
c
a
l
b
e
s
t
-
X
)
+
β
·
(
X
g
l
o
b
a
l
b
e
s
t
-
X
).
For each bitx
j in particle X, If (rand < Sigmoid(v
j)) then x
j = 1, else x
j = 0
Check feasibility for each particle and repair the particle
Calculate fitness function for each particle
UpdateLocalbest Solution and Globalbest Solution
End while
Run LINGO usingGlobalbest Solution to obtain real valued variablesIn particle swarm heuristic, how to encode the problem to set of particles is of great importance. Consider the problem in whichN distribution centers are to be served by K vehicles in T planning periods; we have to map a 2-dimensional array of (
N
×
K
,
T
) for each particle. The first dimension includes K sections, where each section has N binary points. The second dimension includes T 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.Table 1
An example of encoding structure.
Served by vehicle 1
Served by vehicle 2
DC-1
DC-2
DC-3
DC-4
DC-5
DC-1
DC-2
DC-3
DC-4
DC-5
Period 1
1
1
0
0
0
0
0
1
1
1
Period 2
1
0
1
0
0
0
1
0
1
0
Period 3
0
0
0
0
1
0
1
1
0
0
Period 4
1
0
0
1
0
0
0
0
0
0
Period 5
0
0
0
0
0
1
0
0
0
0Based 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 periodsa and b are 1 in a particle. Without loss of generality, let us assume that a
<
a
+
sl
<
b. We set the corresponding value to planning period b
=
0 and set the value corresponding to planning period a
+
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∑
i
,
k
,
t
f
t
c
k
·
X
i
,
k
,
t is used as fitness function. Finally, when the binary variables X
i
,
k
,
t and Y
k
,
t 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 demandsARTIFICIAL 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 Algorithm3.Algorithm 3:Details of perturbation mechanism.
If (NUMBER OF SETUPS >T
1)
ARTIFICIAL SETUP COST = 2 * SETUP COST
ARTIFICIAL HOLDING COST = HOLDING COST/2
End If
If (NUMBER OF SETUPS <T
2)
ARTIFICIAL SETUP COST = SETUP COST/2
ARTIFICIAL HOLDING COST = 2 * HOLDING COST
End If


## 5. 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 MATLAB 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.


### 5.1. 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 Table2. We classified test problems into two groups, small size and large size.Table 2
Number and size of test problems.
Small size
Large size
Problem number
Size (N
×
T
×
K)
Problem number
Size (N
×
T
×
K)
1
4
×
5
×
2
9
30
×
10
×
10
2
4
×
10
×
2
10
30
×
15
×
10
3
6
×
5
×
3
11
40
×
10
×
15
4
6
×
10
×
3
12
40
×
15
×
15
5
8
×
5
×
4
13
50
×
10
×
20
6
8
×
10
×
4
14
50
×
15
×
20
7
10
×
5
×
5
15
60
×
10
×
25
8
10
×
10
×
5
16
60
×
15
×
25We 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.
(vii)
Vehicle capacity was set to(
1.5
∑
i
,
t
d
i
,
t
)
/
N
·
T.
(viii)
Production capacity was set to(
3.5
∑
i
,
t
d
i
,
t
)
/
N
·
T.


### 5.2. Lower Bound Computation

A simple way to validate the quality of solutions provided by proposed method is to compare them to lower bounds on the optimal solution. In order to calculate a lower bound, it is necessary to solve a relaxation of the proposed model. We formulated such a relaxation in relations (31) to (39). In this formulation, the following additional parameter is used.Additional Parameters. v
t
c
k is the variable transportation cost of one item using vehicle k which is equal to f
t
c
k
/
q.The main difference between the original model and the relaxed model is that the binary variablesX
i
,
k
,
t and Y
k
,
t and pertaining constraints are removed, which makes the relaxed model easy to solve.Relaxed Model. Consider
(31)
Min
⁡
∑
t
f
p
c
t
·
Z
t
+
∑
t
v
p
c
t
·
P
t
+
∑
t
h
0
·
I
0
,
t
+
∑
i
,
t
h
i
·
I
i
,
t
+
∑
i
,
k
,
t
v
t
c
k
·
W
i
,
k
,
t
(32)
I
0
,
t
=
I
0
,
t
-
1
+
P
t
-
∑
i
∑
k
W
i
,
k
,
t
∀
t
∈
T
(33)
I
i
,
t
=
I
i
,
t
-
1
+
∑
k
W
i
,
k
,
t
-
d
i
,
t
∀
i
∈
N
,
t
∈
T
(34)
I
0
,
t
≤
u
0
,
t
∀
t
∈
T
(35)
I
i
,
t
≤
u
i
,
t
∀
i
∈
N
,
t
∈
T
(36)
P
t
≤
p
max
⁡
·
Z
t
∀
t
∈
T
(37)
∑
i
W
i
,
k
,
t
≤
q
∀
k
∈
K
,
t
∈
T
(38)
Z
t
∈
{
0,1
}
∀
t
∈
T
(39)
P
t
,
I
0
,
t
,
I
i
,
t
,
W
i
,
k
,
t
≥
0
∀
i
∈
N
,
k
∈
K
,
t
∈
T
.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∑
i
,
k
,
t
v
t
c
k
·
W
i
,
k
,
t that we substituted to the term ∑
k
,
t
f
t
c
k
·
Y
k
,
t in the original model and always underestimate it.


### 5.3. Results

To run the proposed method, the algorithm parameters were set according to Table3.Table 3
Algorithm parameters.
Parameter
Value
T
1
[T/3]
T
2
[2T/3]
Number of PSO particles
50
α
0.35
β
0.65
Number of PSO iterations
100A summary of the results for small size test problems is reported in Table4. In columns 2–4 of the table, objective function obtained by the proposed method (
Z
Alg
), runtime of the proposed method in seconds, and optimal solution obtained by LINGO (
Z
LINGO
) are reported, respectively. In the final column, the gap between Z
Alg and Z
LINGO that is equal to (
(
Z
Alg
-
Z
LINGO
)
/
Z
Alg
)
×
100 is presented.Table 4
Computational results for small size test problems.
Number
Z
Alg
CP
U
Alg
Z
LINGO
Gap
1
38369
6.3
38369
0.00%
2
76739
14.7
75709
1.34%
3
57554
17.0
56239
2.29%
4
115108
36.9
113258
1.61%
5
76739
32.5
75315
1.86%
6
153478
70.3
149429
2.64%
7
95924
48.9
93993
2.01%
8
191004
118.6
184744
3.28%We observe in Table4 that the proposed method can provide good solutions in short runtime. The gaps between Z
Alg and Z
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 Table5. In columns 2–5 of the table, the average and worst objective functions obtained by the proposed method (Z
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 (Z
LB) and the average and worst gap between Z
Alg and Z
LB that is equal to (
(
Z
Alg
-
Z
LB
)
/
Z
LB
)
×
100 are reported, respectively. Finally in columns 9-10, the best solution obtained by LINGO within 2 hours run (Z
LINGO) and the gap between Z
Alg and Z
LB which is equal to (
(
Z
LINGO
-
Z
LB
)
/
Z
LB
)
×
100 are presented, respectively.Table 5
Computational results (large size test problems).
Number
Z
Alg
CP
U
Alg
Z
LB
Gap 1
Z
LINGO
Gap 2
Avg.
Worst
Avg.
Worst
Avg.
Worst
9
556028
572704
1097
1141
526941
5.52%
8.68%
590796
12.12%
10
834042
860841
1691
1773
769475
8.39%
11.87%
895742
16.41%
11
741370
765830
2175
2244
692744
7.02%
10.55%
799822
15.46%
12
1112056
1147772
3252
3374
1054671
5.44%
8.83%
1180939
11.97%
13
926713
955584
3919
4039
863965
7.26%
10.60%
984028
13.90%
14
1390069
1426843
4702
4924
1293635
7.45%
10.30%
1476745
14.15%
15
1112056
1143883
4892
5088
1030104
7.96%
11.05%
1184692
15.01%
16
1668083
1711109
5791
5965
1549744
7.64%
10.41%
1792531
15.67%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 Table4 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 Section3, 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.Figure 1
Average run time of the proposed method.


## 5.1. 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 Table2. We classified test problems into two groups, small size and large size.Table 2
Number and size of test problems.
Small size
Large size
Problem number
Size (N
×
T
×
K)
Problem number
Size (N
×
T
×
K)
1
4
×
5
×
2
9
30
×
10
×
10
2
4
×
10
×
2
10
30
×
15
×
10
3
6
×
5
×
3
11
40
×
10
×
15
4
6
×
10
×
3
12
40
×
15
×
15
5
8
×
5
×
4
13
50
×
10
×
20
6
8
×
10
×
4
14
50
×
15
×
20
7
10
×
5
×
5
15
60
×
10
×
25
8
10
×
10
×
5
16
60
×
15
×
25We 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.
(vii)
Vehicle capacity was set to(
1.5
∑
i
,
t
d
i
,
t
)
/
N
·
T.
(viii)
Production capacity was set to(
3.5
∑
i
,
t
d
i
,
t
)
/
N
·
T.


## 5.2. Lower Bound Computation

A simple way to validate the quality of solutions provided by proposed method is to compare them to lower bounds on the optimal solution. In order to calculate a lower bound, it is necessary to solve a relaxation of the proposed model. We formulated such a relaxation in relations (31) to (39). In this formulation, the following additional parameter is used.Additional Parameters. v
t
c
k is the variable transportation cost of one item using vehicle k which is equal to f
t
c
k
/
q.The main difference between the original model and the relaxed model is that the binary variablesX
i
,
k
,
t and Y
k
,
t and pertaining constraints are removed, which makes the relaxed model easy to solve.Relaxed Model. Consider
(31)
Min
⁡
∑
t
f
p
c
t
·
Z
t
+
∑
t
v
p
c
t
·
P
t
+
∑
t
h
0
·
I
0
,
t
+
∑
i
,
t
h
i
·
I
i
,
t
+
∑
i
,
k
,
t
v
t
c
k
·
W
i
,
k
,
t
(32)
I
0
,
t
=
I
0
,
t
-
1
+
P
t
-
∑
i
∑
k
W
i
,
k
,
t
∀
t
∈
T
(33)
I
i
,
t
=
I
i
,
t
-
1
+
∑
k
W
i
,
k
,
t
-
d
i
,
t
∀
i
∈
N
,
t
∈
T
(34)
I
0
,
t
≤
u
0
,
t
∀
t
∈
T
(35)
I
i
,
t
≤
u
i
,
t
∀
i
∈
N
,
t
∈
T
(36)
P
t
≤
p
max
⁡
·
Z
t
∀
t
∈
T
(37)
∑
i
W
i
,
k
,
t
≤
q
∀
k
∈
K
,
t
∈
T
(38)
Z
t
∈
{
0,1
}
∀
t
∈
T
(39)
P
t
,
I
0
,
t
,
I
i
,
t
,
W
i
,
k
,
t
≥
0
∀
i
∈
N
,
k
∈
K
,
t
∈
T
.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∑
i
,
k
,
t
v
t
c
k
·
W
i
,
k
,
t that we substituted to the term ∑
k
,
t
f
t
c
k
·
Y
k
,
t in the original model and always underestimate it.


## 5.3. Results

To run the proposed method, the algorithm parameters were set according to Table3.Table 3
Algorithm parameters.
Parameter
Value
T
1
[T/3]
T
2
[2T/3]
Number of PSO particles
50
α
0.35
β
0.65
Number of PSO iterations
100A summary of the results for small size test problems is reported in Table4. In columns 2–4 of the table, objective function obtained by the proposed method (
Z
Alg
), runtime of the proposed method in seconds, and optimal solution obtained by LINGO (
Z
LINGO
) are reported, respectively. In the final column, the gap between Z
Alg and Z
LINGO that is equal to (
(
Z
Alg
-
Z
LINGO
)
/
Z
Alg
)
×
100 is presented.Table 4
Computational results for small size test problems.
Number
Z
Alg
CP
U
Alg
Z
LINGO
Gap
1
38369
6.3
38369
0.00%
2
76739
14.7
75709
1.34%
3
57554
17.0
56239
2.29%
4
115108
36.9
113258
1.61%
5
76739
32.5
75315
1.86%
6
153478
70.3
149429
2.64%
7
95924
48.9
93993
2.01%
8
191004
118.6
184744
3.28%We observe in Table4 that the proposed method can provide good solutions in short runtime. The gaps between Z
Alg and Z
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 Table5. In columns 2–5 of the table, the average and worst objective functions obtained by the proposed method (Z
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 (Z
LB) and the average and worst gap between Z
Alg and Z
LB that is equal to (
(
Z
Alg
-
Z
LB
)
/
Z
LB
)
×
100 are reported, respectively. Finally in columns 9-10, the best solution obtained by LINGO within 2 hours run (Z
LINGO) and the gap between Z
Alg and Z
LB which is equal to (
(
Z
LINGO
-
Z
LB
)
/
Z
LB
)
×
100 are presented, respectively.Table 5
Computational results (large size test problems).
Number
Z
Alg
CP
U
Alg
Z
LB
Gap 1
Z
LINGO
Gap 2
Avg.
Worst
Avg.
Worst
Avg.
Worst
9
556028
572704
1097
1141
526941
5.52%
8.68%
590796
12.12%
10
834042
860841
1691
1773
769475
8.39%
11.87%
895742
16.41%
11
741370
765830
2175
2244
692744
7.02%
10.55%
799822
15.46%
12
1112056
1147772
3252
3374
1054671
5.44%
8.83%
1180939
11.97%
13
926713
955584
3919
4039
863965
7.26%
10.60%
984028
13.90%
14
1390069
1426843
4702
4924
1293635
7.45%
10.30%
1476745
14.15%
15
1112056
1143883
4892
5088
1030104
7.96%
11.05%
1184692
15.01%
16
1668083
1711109
5791
5965
1549744
7.64%
10.41%
1792531
15.67%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 Table4 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 Section3, 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.Figure 1
Average run time of the proposed method.


## 6. 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 making. 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.



---

*Source: 475606-2014-05-12.xml*