# A Multiobjective Stochastic Production-Distribution Planning Problem in an Uncertain Environment Considering Risk and Workers Productivity

Authors: S. M. J. Mirzapour Al-e-Hashem; A. Baboli; S. J.  Sadjadi; M. B. Aryanezhad
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

A multi-objective two stage stochastic programming model is proposed to deal with a multi-period multi-product multi-site production-distribution planning problem for a midterm planning horizon. The presented model involves majority of supply chain cost parameters such as transportation cost, inventory holding cost, shortage cost, production cost. Moreover some respects as lead time, outsourcing, employment, dismissal, workers productivity and training are considered. Due to the uncertain nature of the supply chain, it is assumed that cost parameters and demand fluctuations are random variables and follow from a pre-defined probability distribution. To develop a robust stochastic model, an additional objective functions is added to the traditional production-distribution-planning problem. So, our multi-objective model includes (i) the minimization of the expected total cost of supply chain, (ii) the minimization of the variance of the total cost of supply chain and (iii) the maximization of the workers productivity through training courses that could be held during the planning horizon. Then, the proposed model is solved applying a hybrid algorithm that is a combination of Monte Carlo sampling method, modifiedε-constraint method and L-shaped method. Finally, a numerical example is solved to demonstrate the validity of the model as well as the efficiency of the hybrid algorithm.


---

## Body

## 1. Introduction

One of the problems that could be addressed in the scope of supply chain management is production-distribution planning which is an operational activity that does a plan for the production process, to give an idea to management as to what quantity of materials and other resources are to be procured and when, so that the total cost of operations of the organization is kept to the minimum over that period. Production-distribution planning has attracted the attention of many researchers from several years ago [1]. Numerous production-distribution models with varying degrees of difficulties have been proposed in the last decades. Since Holt et al. [2] proposed the approach for the first time; scholars have developed numerous models to help solving the production planning problems, each with their own supporters and detractors. A rough classification of modeling approaches for production planning is presented by Soyster [3]; techniques that find the exact and mathematical solutions and the techniques search for numerical solutions. As a comprehensive remark, Nam and Logendran [4] reviewed production planning models from 140 journal articles and 14 books and categorized the models into optimal and near-optimal classifications. Hanssmann and Hess [5] developed a model based on the linear programming approach using a linear cost structure of the decision variables. Haehling [6] extended the Hanssmann and Hess [5] model for multiproduct, multistage production systems in which optimal disaggregation decisions can be made under capacity constraints. Masud and Hwang [7] presented three multicriteria decision making (MCDM) methods, which were goal programming (GP), the step method and sequential multiobjective problem. These methods apply to solve production planning problem with maximizing profit, minimizing changes in workforce level, minimizing inventory investment and minimizing backorders. A set of data consisting of two products, a single production plant and eight planning periods was generated to compare the results. Goodman [8] developed a GP model which approximates the original nonlinear cost terms of the Holt’s model by linear terms and solves it using a variant of the simplex method. Baykasoglu [9] developed Masud and Hwang’s model with more constraints such as subcontractor selection and setup issues. A tabu search algorithm was designed to solve the pre-emptive GP model. The integration of production planning problems with other planning problems were considered, for instance scheduling problems [10, 11], manpower planning problems [12] and long set up time problems [13]. Production planning in many production environments is based on some parameters with uncertain values.In spite of the fact that the concepts of variance has been considered in other areas, but to the best of our knowledge, it is the first time workers productivity is considered in a multiobjective scheme to model robust production-distribution planning under uncertainty. Moreover, the idea of involving the human-related issues such as workers' skill level and workers' training is also incorporated into the model. Using this idea, we have the option of training the workers instead of firing them and then hiring new full-skilled ones. Since the expected total cost, the variance of the total cost and the workers productivity are in conflict with each other, it is proposed to model a multiobjective production-distribution problem whose solution will be a set of Pareto-optimal possible plan alternatives representing the trade-off among different objectives rather than a unique solution. Some approaches to deal with solving a multiobjective production-distribution planning under uncertainty are developed such as Possibilistic linear programming method [14] and fuzzy goal programming approach [15, 16].According to Masud and Hwang [7], the methods for solving MOMP problems can be classified into three categories, based on the phase in which the decision maker engages in the decision making process expressing his preferences.Thea priori methods, the interactive methods and the a posteriori or generation methods. In a priori methods the decision maker expresses his preferences before the solution process (e.g., setting goals or weights to the objective functions). The criticism about the a priori methods is that it is very difficult for the decision maker to know beforehand and to be able to accurately quantify (either by means of goals or weights) his preferences. In the interactive methods phases of dialogue with the decision maker are interchanged with phases of calculation and the process usually converges, after a few iterations, to the most preferred solution. The decision maker progressively drives the search with his answers towards the most preferred solution. The drawback is that he never sees the whole picture (the Pareto set) or an approximation of it. Hence, the most preferred solution is ‘‘most preferred” in relation to what he has seen and compare so far. In the a posteriori methods the efficient solutions of the problem (all of them or a sufficient representation) are generated and then the decision maker engages, in order to select among them, the most preferred one.We formulate the proposed model as a multiobjective robust stochastic mixed-integer nonlinear programming problem, after linearization, it is solved by using a hybrid algorithm that is a combination of the extended Monte Carlo sampling method, modifiedε-constraint technique (type of the a posteriori methods) which is a new version of the traditional famous multicriteria decision making method (ε-constraint) for solving multiobjective problems with conflicting objectives simultaneously and the L-shaped method which is on of the efficient heuristic method to solve two-stage stochastic optimization problems. This formulation takes into account not only the expected total cost of supply chain, but also the risk reflected by the variability of the total cost. The result of the model suggest a set of Pareto-optimal solutions and give this chance to the decision maker in order to find the best production-distribution configuration according to his viewpoint.The rest of the paper is organized as follows: in Section2, the problem is described. Section 3 presents the mathematical formulation considering workers productivity. Then the solution procedure is presented in Section 4. Next, the robustness and effectiveness of the proposed model are demonstrated by the computational experiments in Sections 5 and 6. Finally, conclusions are presented in Section 7.


## 2. Problem Description

The proposed multiobjective multiproduct multisite production-distribution problem can be described as follows.There areJ factories and C customers. Each factory produces several products. Each factory characterized by its own available time for production and warehouse capacities. The available time is limited to the number of k-level labors beside the allowed amount of regular time and overtime. Every factory could subcontract an allowed proportion of its product to subcontractors. The transportation cost between factories and customers’ zones as well as the production cost of a certain item at different factories can be different.The present work formulates the production-distribution problem as a robust multiobjective nonlinear programming and tries to minimize the expected total cost of supply chain, the variability of the total cost of supply chain and the workers productivity, simultaneously, and take decisions for each period as follows:(i)
the quantity of producti manufactured at factory j to fulfill stochastic demand of customer zone c by k-level labor,(ii)
the number ofk-level labors would be employed, laid off or trained at each factory,(iii)
the quantity of producti stored at factory j,(iv)
the amount of demand in each customer’s zone is not met.In our proposed model the scenario-based approach is used to represent the uncertain parameters. Due to the multiperiod multisite multiproduct nature of the model, the problem includes a large number of uncertain parameters; a resulting challenge is that a large number of scenarios are required. To reduce the model size and the number of scenarios, we use an extended Monte Carlo sampling method to generate the scenarios. Each scenario is then associated with the same probability with the summation of the probabilities for all the scenarios equal to 1. The extended Monte Carlo sampling method is an extension of the conventional Monte Carlo sampling method in which interaction between uncertainties is analyzed. Therefore value assignment for dependent uncertain parameters is controlled regarding the type and the level of possible dependencies.


## 3. Mathematical Formulation

In this paper, a novel multiobjective stochastic robust optimization approach is presented in which uncertainty is represented by a set of discrete scenarios (n).


### 3.1. Notations

ParametersDictn: demand for product i(1,2,…,I) in demand point c(1,2,…,C) in period t(1,2,…,T) in scenario n(1,2,…,N);Cqjn: production cost per hour in regular time (q=1),overtime (q=2), and subcontracting (q=3) at factory j(1,2,…,J) in scenario n;Lkjtn: manpower cost of k-level (k=1,2,…,K) labors at factory j in period t in scenario n;aij: production time of product i at factory j;Fkjtn: firing cost of k-level worker at factory j in period t in scenario n;Hkjtn: hiring cost of k-level worker at factory j in period t in scenario n;Trkk′jtn: training cost for k-level worker trained to level k′ at factory j in period t in scenario n;I1ijtn: inventory holding cost for product i at factory j in period t in scenario n;I2ictn: inventory holding cost for finished product i in customer's zone c in period t in scenario n;Tictξ: transportation cost from factory j to demand point c at period t in scenario n;αt: fraction of the workforce variation allowed in period t;υk: productivity of k-level labors (0≤υk≤1);TIqjt: available regular time (q=1),overtime (q=2) and capacity of subcontracting (q=3) in terms of time unit at factory j in period t;P1j: product storage capacity at factory j;P2c: product storage capacity in customer's zone c;LTjc: lead time required for shipping end product from factory j to demand point c;UPkk′: 1 if training from skill level k to skill level k′ is possible; 0 otherwise;πictn: shortage cost of product i in customer’s zone c in period t in scenario n;ρn: Occurrence probability of scenario n.VariablesXijgt: number of product i produced at factory j using method g in period t;XLkjt: number of k-level workers at factory j in period t;XFkjt: number of k-level workers at factory j fired in period t;XHkjt: number of k-level workers at factory j hired in period t;XUkk′jt: number of k-level workers at factory j trained to level k′ in period t;XPijt: inventory level of product i at factory j in period t;XIictn: inventory level of product i in customer's zone c in period t in scenario n;YSijctn: number of units of product i provided by factory j for demand point c in period t in scenario n;Bictn: shortage of product i in demand point c in period t in scenario n;TCn: Total cost of supply chain under scenario n.


### 3.2. Multi-Objective Stochastic Production-Distribution Model

(3.1)Min⁡Z1=∑nρnTCn,(3.2)Min⁡Z2=∑nρn|TCn-∑nρnTCn|,(3.3)Max⁡Z3=∑t∑j∑kukXLkjt∑t∑j∑kXLkjt,(3.4)subjecttoXPijt=XPij(t-1)+∑qXijqt-∑cYSijctn,∀i,j,t,n,(3.5)XLkjt=XLkj(t-1)+XHkjt-XFkjt+∑k′XUk′kjt-∑k′XUkk′jt,∀k,j,t,(3.6)∑kXLkjtυk(TI1jt+TI2jt)≥∑i,q∈{1,2}xijqt⋅aij,∀j,t,(3.7)∑ixij3t⋅aij≤TI3jt,∀j,t,(3.8)XIictn=XIict-1n+∑jYSijc[t-LTjc]n-Dictn-Bic(t-1)n,∀i,c,t,n,(3.9)∑iXPijt≤P1j,∀j,t,(3.10)∑iXIictn≤P2c,∀c,t,(3.11)∑k(XFkjt+XHkjt)≤α(t-1)∑kXLkj(t-1),∀j,t,(3.12)XFkjt+∑k′XUkk′jt≤XLkj(t-1),∀k,j,t,(3.13)∑k′XUk′kjt⋅XFkjt=0,∀k,j,t,(3.14)XUk′kjt≤M⋅UPkk′,∀k,k′,j,t,(3.15)Xijqt,XPijt,YSijct,Bictn≥0,XLkjt,XFkjt,XHkjt,XUkk′jt≥0,andinteger,∀i,j,c,n,k,s,m,t.
First objective function (3.1) aims to minimize the expected total cost of supply chain, where TCn is total cost of supply chain under the realization of scenario n and defined as follows: (3.16)TCn=[∑i,j,q,taijCqjn⋅Xijqt+∑k,j,tLkjtn⋅XLkjt+∑k,j,tFkjtn⋅XFkjt+∑k,j,tHkjtn⋅XHkjt+∑k,k′,j,tTr⁡kk′jn⋅XUkk′jt+∑i,j,tI1ijtn⋅XPijt+∑i,c,tI2ictn⋅XIictn+∑i,j,c,tTictn⋅YSijct+∑i,c,tπictn⋅Bictn]
and including production cost, labor cost, firing cost, hiring cost, training cost, inventory holding cost, transportation cost and shortage cost, respectively. Second objective function (3.2) aims to minimize the variability of total cost of supply chain. Third objective function (3.3) aims to maximize the average manpower productivity along the factories over the planning horizon. Constraint (3.4) is a balance equation for the product inventory at factory j. Constraint (3.5) is also a balance equation for workforce level and ensures that the available k-skill level labors equals the workforce with the same skill level in previous period in addition to the change of workforce level in current period. Constraint (3.6) limits the available production time to available workforce regular and overtime, considering their productivity. Constraint (3.7) restricts the amount of products manufactured by subcontractor. Constraint (3.8) is an inventory balance equation for demand point c. Constraints (3.9) limit the raw material and product inventory levels of factories to their related inventory storage capacities. Constraint (3.10) restricts the product inventory levels at each customer's zones to their related inventory storage capacities. Constraint (3.11) guarantees that the change in workforce level cannot exceed the proportion of workers in previous period. Constraint (3.12) ensures that the number of k-level workers who are fired or trained for upper skill levels in current period cannot exceed the available k-level workforce in previous period. Constraint (3.13) denotes that the labors that are trained for skill level k should not be fired in the same period. Constraint (3.14) guarantees that training from skill level k to level k′ is possible, once this training program exists. Constraint (3.15) denotes the variable types.


## 3.1. Notations

ParametersDictn: demand for product i(1,2,…,I) in demand point c(1,2,…,C) in period t(1,2,…,T) in scenario n(1,2,…,N);Cqjn: production cost per hour in regular time (q=1),overtime (q=2), and subcontracting (q=3) at factory j(1,2,…,J) in scenario n;Lkjtn: manpower cost of k-level (k=1,2,…,K) labors at factory j in period t in scenario n;aij: production time of product i at factory j;Fkjtn: firing cost of k-level worker at factory j in period t in scenario n;Hkjtn: hiring cost of k-level worker at factory j in period t in scenario n;Trkk′jtn: training cost for k-level worker trained to level k′ at factory j in period t in scenario n;I1ijtn: inventory holding cost for product i at factory j in period t in scenario n;I2ictn: inventory holding cost for finished product i in customer's zone c in period t in scenario n;Tictξ: transportation cost from factory j to demand point c at period t in scenario n;αt: fraction of the workforce variation allowed in period t;υk: productivity of k-level labors (0≤υk≤1);TIqjt: available regular time (q=1),overtime (q=2) and capacity of subcontracting (q=3) in terms of time unit at factory j in period t;P1j: product storage capacity at factory j;P2c: product storage capacity in customer's zone c;LTjc: lead time required for shipping end product from factory j to demand point c;UPkk′: 1 if training from skill level k to skill level k′ is possible; 0 otherwise;πictn: shortage cost of product i in customer’s zone c in period t in scenario n;ρn: Occurrence probability of scenario n.VariablesXijgt: number of product i produced at factory j using method g in period t;XLkjt: number of k-level workers at factory j in period t;XFkjt: number of k-level workers at factory j fired in period t;XHkjt: number of k-level workers at factory j hired in period t;XUkk′jt: number of k-level workers at factory j trained to level k′ in period t;XPijt: inventory level of product i at factory j in period t;XIictn: inventory level of product i in customer's zone c in period t in scenario n;YSijctn: number of units of product i provided by factory j for demand point c in period t in scenario n;Bictn: shortage of product i in demand point c in period t in scenario n;TCn: Total cost of supply chain under scenario n.


## 3.2. Multi-Objective Stochastic Production-Distribution Model

(3.1)Min⁡Z1=∑nρnTCn,(3.2)Min⁡Z2=∑nρn|TCn-∑nρnTCn|,(3.3)Max⁡Z3=∑t∑j∑kukXLkjt∑t∑j∑kXLkjt,(3.4)subjecttoXPijt=XPij(t-1)+∑qXijqt-∑cYSijctn,∀i,j,t,n,(3.5)XLkjt=XLkj(t-1)+XHkjt-XFkjt+∑k′XUk′kjt-∑k′XUkk′jt,∀k,j,t,(3.6)∑kXLkjtυk(TI1jt+TI2jt)≥∑i,q∈{1,2}xijqt⋅aij,∀j,t,(3.7)∑ixij3t⋅aij≤TI3jt,∀j,t,(3.8)XIictn=XIict-1n+∑jYSijc[t-LTjc]n-Dictn-Bic(t-1)n,∀i,c,t,n,(3.9)∑iXPijt≤P1j,∀j,t,(3.10)∑iXIictn≤P2c,∀c,t,(3.11)∑k(XFkjt+XHkjt)≤α(t-1)∑kXLkj(t-1),∀j,t,(3.12)XFkjt+∑k′XUkk′jt≤XLkj(t-1),∀k,j,t,(3.13)∑k′XUk′kjt⋅XFkjt=0,∀k,j,t,(3.14)XUk′kjt≤M⋅UPkk′,∀k,k′,j,t,(3.15)Xijqt,XPijt,YSijct,Bictn≥0,XLkjt,XFkjt,XHkjt,XUkk′jt≥0,andinteger,∀i,j,c,n,k,s,m,t.
First objective function (3.1) aims to minimize the expected total cost of supply chain, where TCn is total cost of supply chain under the realization of scenario n and defined as follows: (3.16)TCn=[∑i,j,q,taijCqjn⋅Xijqt+∑k,j,tLkjtn⋅XLkjt+∑k,j,tFkjtn⋅XFkjt+∑k,j,tHkjtn⋅XHkjt+∑k,k′,j,tTr⁡kk′jn⋅XUkk′jt+∑i,j,tI1ijtn⋅XPijt+∑i,c,tI2ictn⋅XIictn+∑i,j,c,tTictn⋅YSijct+∑i,c,tπictn⋅Bictn]
and including production cost, labor cost, firing cost, hiring cost, training cost, inventory holding cost, transportation cost and shortage cost, respectively. Second objective function (3.2) aims to minimize the variability of total cost of supply chain. Third objective function (3.3) aims to maximize the average manpower productivity along the factories over the planning horizon. Constraint (3.4) is a balance equation for the product inventory at factory j. Constraint (3.5) is also a balance equation for workforce level and ensures that the available k-skill level labors equals the workforce with the same skill level in previous period in addition to the change of workforce level in current period. Constraint (3.6) limits the available production time to available workforce regular and overtime, considering their productivity. Constraint (3.7) restricts the amount of products manufactured by subcontractor. Constraint (3.8) is an inventory balance equation for demand point c. Constraints (3.9) limit the raw material and product inventory levels of factories to their related inventory storage capacities. Constraint (3.10) restricts the product inventory levels at each customer's zones to their related inventory storage capacities. Constraint (3.11) guarantees that the change in workforce level cannot exceed the proportion of workers in previous period. Constraint (3.12) ensures that the number of k-level workers who are fired or trained for upper skill levels in current period cannot exceed the available k-level workforce in previous period. Constraint (3.13) denotes that the labors that are trained for skill level k should not be fired in the same period. Constraint (3.14) guarantees that training from skill level k to level k′ is possible, once this training program exists. Constraint (3.15) denotes the variable types.


## 4. Solution Procedure

In order to overcome the complexity of multiobjective stochastic programming problems, we combine two techniques; the modifiedε-constraint method and the L-shaped method. The modified ε-constraint method offers an overall framework to obtain the optimal Pareto solutions for the multiobjective optimization problems. Within this framework, the L-shaped method is called to solve two-stage stochastic programming model.


### 4.1. Modifiedε-Constraint Method

In this paper, we applied a modified version ofε-constraint method [17]. In this method, one of the objective functions with some changes is selected as the main objective function to be optimized, and all other objective functions are transformed into constraint by considering an upper bound for each of them. The proposed modified ε-constraint method consists of the following stages.Step 1.
Select one of the objective functions (Zj) as the main objective function to be optimized and convert other objective functions into constraint. Then form the payoff table by the individual optimization of eachobjective functions separately. The interval between the ideal value and the worst value over the Pareto set for each objective function is named here as the range of that objective function.Step 2.
Determine the grid points(εk). Then we divide the range of each objective function to m equal intervals using m-1 intermediate equidistant grid points; that are used to vary parametrically the RHS (εk) of that objective function.Step 3.
The modifiedε-constraint model is solved for each value of ε which is obtained in the previous step(4.1)Min⁡Zj(x)-θ∑k≠jskrks.t.Zk(x)+sk=εk,∀k≠j,x∈X,sk∈R+,
where θ is an adequate small number usually between 10-6 and 10-3, X is the feasible region of the original problem sk is a slack variable. Also, rk, gk, Ndk are the range, the number of grid points and the nadir value for objective function k, respectively.Note that, at each iteration of the internal loop of modifiedε-constraint method, a two-stage stochastic programming model must be solved. To achieve the optimal solution for this model, another algorithm is embedded in modified ε-constraint method which is L-shaped method and we describe it the next subsection.


### 4.2.L-Shaped Method

As mentioned before our proposed model is a multiobjective stochastic robust optimization model in which uncertainty is represented by a set of discrete scenarios, we use an extended Monte Carlo sampling approach to descretize the continuous distribution functions and generate the scenarios [18]. This method is characterized by an exponential increase in the problem size with the number of scenarios as well as the number of uncertain parameters due to the nested structure of the two-stage formulation. In this paper, we use L-shaped method, a popular method for solving stochastic programming models, to take advantage of the special decomposable structure of the two-stage stochastic optimization model.The idea behind of the L-shaped method is to first solve the master problem (the model with those constraints that do not include the second stage variables) to obtain a lower bound of the objective value. We then fix all the first stage decisions and solve each scenario sub-problem (inner-model that include second stage decisions) to get an upper bound. If the lower bound and the upper bound fall into a pre specified tolerance, then the algorithm stops. Otherwise, we add a cut by using of the duals of the scenario sub-problems and return to the master problem. We use this method whenever needed in the inner loops of modifiedε-constraint method (see Figure 1).Figure 1
Flowchart of the proposed method.


## 4.1. Modifiedε-Constraint Method

In this paper, we applied a modified version ofε-constraint method [17]. In this method, one of the objective functions with some changes is selected as the main objective function to be optimized, and all other objective functions are transformed into constraint by considering an upper bound for each of them. The proposed modified ε-constraint method consists of the following stages.Step 1.
Select one of the objective functions (Zj) as the main objective function to be optimized and convert other objective functions into constraint. Then form the payoff table by the individual optimization of eachobjective functions separately. The interval between the ideal value and the worst value over the Pareto set for each objective function is named here as the range of that objective function.Step 2.
Determine the grid points(εk). Then we divide the range of each objective function to m equal intervals using m-1 intermediate equidistant grid points; that are used to vary parametrically the RHS (εk) of that objective function.Step 3.
The modifiedε-constraint model is solved for each value of ε which is obtained in the previous step(4.1)Min⁡Zj(x)-θ∑k≠jskrks.t.Zk(x)+sk=εk,∀k≠j,x∈X,sk∈R+,
where θ is an adequate small number usually between 10-6 and 10-3, X is the feasible region of the original problem sk is a slack variable. Also, rk, gk, Ndk are the range, the number of grid points and the nadir value for objective function k, respectively.Note that, at each iteration of the internal loop of modifiedε-constraint method, a two-stage stochastic programming model must be solved. To achieve the optimal solution for this model, another algorithm is embedded in modified ε-constraint method which is L-shaped method and we describe it the next subsection.


## 4.2.L-Shaped Method

As mentioned before our proposed model is a multiobjective stochastic robust optimization model in which uncertainty is represented by a set of discrete scenarios, we use an extended Monte Carlo sampling approach to descretize the continuous distribution functions and generate the scenarios [18]. This method is characterized by an exponential increase in the problem size with the number of scenarios as well as the number of uncertain parameters due to the nested structure of the two-stage formulation. In this paper, we use L-shaped method, a popular method for solving stochastic programming models, to take advantage of the special decomposable structure of the two-stage stochastic optimization model.The idea behind of the L-shaped method is to first solve the master problem (the model with those constraints that do not include the second stage variables) to obtain a lower bound of the objective value. We then fix all the first stage decisions and solve each scenario sub-problem (inner-model that include second stage decisions) to get an upper bound. If the lower bound and the upper bound fall into a pre specified tolerance, then the algorithm stops. Otherwise, we add a cut by using of the duals of the scenario sub-problems and return to the master problem. We use this method whenever needed in the inner loops of modifiedε-constraint method (see Figure 1).Figure 1
Flowchart of the proposed method.


## 5. Numerical Experiments

Consider the supply chain network problem depicted in Figure2. A typical company is willing to plan its aggregate production planning. The planning horizon of time is assumed to be 12 periods. Also the number of skill levels and products are both assumed to be 5. This company owns four factories F1, F2, F3, and F4 which are spread geographically, and three customer centers located in three different cities C1, C2, and C3. We assume that the demand follows a normal distribution with the expected value and standard deviation equal to 1000 and 100, respectively, (N(μ:1000,σ:100)). Cost items' distribution functions as well as associated parameters are shown in Table 1.Table 1
Cost items distribution functions.
Cost itemProbability distributionProduct  inventory holding cost ($/unit period)Uniform(5,20)Hiring cost (10$/manpower)Normal(μ:6,σ2:3)Firing cost (10$/manpower)Normal(μ:10,σ2:3)Salary cost (10$/manpower)Uniform(20,30)Training cost (10$/manpower)Normal(μ:20,σ2:5)Production cost ($/min)Uniform(0.5,1.5)Transportation cost ($/unit)Uniform(0.015,0.25)Shortage cost ($/period. unit)Normal(μ:2.5,σ2:1.5)Figure 2
Supply chain network.We solve the example with a sampling size of 100 scenarios. The mathematical model has 1920 integer variables, 19680 continuous decision variables and 44148 constraints. Using L-shaped method the ideal and the nadir vales for each objective functions are obtained and reported in Table2.Table 2
The payoff table for modifiedε-constraint method.
The optimal solution forkth single-objective model (k=2,…,K)Z1(x)Z2(x)Z3(x)x1*1,733,21292,705.30.25x2*10,594,46100.5x3*4,414,184281,715.40.95Ideal value (Id)1,733,21200.95Nadir value (Nd)10,594,461281,715.40.25Figure3 shows the convergence of the L-shaped method for the first objective function. We consider the first objective function (total cost of supply chain) as the main objective function and divide the other objectives' ranges by 9 and 3 grid points for variability and workers productivity, respectively. Pareto curve for the expected and variance of total cost is given in Figure 4. As can be seen, there is a significant conflict between expected value and the variability of the total cost of supply chain. This condition arises from this fact that in the case of expected total cost, the model tries to find the solutions that they have a good expected value not regarding to the variability of the total cost under realization of the different scenarios. Conversely, in the case of the variability of total cost, the model tries to find solutions that they have objective values as close as possible to each other under realization of the different scenarios not regarding to the objective values (see Figure 5).Figure 3
Convergence of the L-shaped method.Figure 4
Pareto curve for the expected and variance of total cost.The behavior ofZ1 against Z2.
(a)(b)The size of the resulting stochastic programming problem is very large and increases exponentially as the number of scenarios increases. For the stochastic programming model with 100 scenario case, mathematical programming solvers could not solve the problem in reasonable amount of time due to its huge size, that is, the problem cannot be solved directly, although the deterministic model can be solved to optimality within five minutes. By using the proposed method, we can obtain the optimal solution for the 100 scenario case in around 1 h with 0.01% optimality tolerance. In Table3 we report the state of the staff upgrading versus the average workers productivity, as can be seen, as the average productivity decreases, the diversity of training courses decreases. For instance, in average productivity level 1, 9 training courses are held, but in levels 0.75 and 0.5 the number of courses decreases to 8 and 4, respectively. Finally in the case of level 0.25 training diversity decreases to 3. In other words, workers' training plays a significant role to enhance workers productivity in factories.Table 3
State of the staff upgrading versus the average workers productivity.
Z3Upgraded levelFactoryjPeriodt1234567891011123→5113→52494→52511→33312→43113→5461→51430.751→52282→51101→5180.53→5181→52120.251→51633→514


## 6. Efficiency of the Proposed Method

In this section, five large scaled test problems are generated to evaluate the efficiency of the proposed algorithm. As we described earlier, for these problems, standard mathematical programming solvers cannot solve the problem in reasonable amount of time. Therefore, each test problem is solved four times with 50, 100, 500 and 1000 scenarios and compared with lower bounds of linear programming solvers obtained after one and half an hour.To evaluate the efficiency of the algorithm, the usual relative gap (RG) between the average of best values of first objective function in Pareto set solutions (AB) (obtained from the proposed method) and the average of the lower bounds (AL) of the first objective function in Pareto set solutions (obtained by standard linear programming solver) is used and reported in Table4(6.1)RG=AB-AL⁡AL⁡×100.
Table 4 shows the characteristics of the test problems and compares the performance obtained by the proposed method with different scenario numbers. As can be seen, the proposed method can solve the problem in less than half an hour, even for large scaled ones. This comparison demonstrates that the relative gap between the lower bound of the objective function and the best value obtained from the proposed method in worst case is not more than 2.78 percent.Table 4
Comparison of the performance of the proposed algorithm with different scenario numbers.
Prob. No.Problem info. No. ofj/c50 scenarios100 scenariosCPU Time (min)RG%CPU Time (min)RG%112/20171.253171.358220/25181.668181.681325/30202.446232.588425/35212.477242.356530/40252.302253.01Prob. No.Problem info. No. ofj/c500 scenarios1000 scenariosCPU Time (min)RG%CPU Time (min)RG%112/20181.392201.501220/25201.382211.661325/30242.508272.446425/35272.210282.477530/40292.788322.302


## 7. Conclusion

In this paper a multiobjective two-stage stochastic programming model is developed to deal with production-distribution planning in an uncertain supply chain considering workers productivity. In addition to the traditional production planning problem in which the total cost is considered as the main objective function we added two extra objective functions that are variability and workers productivity. Risk is described in the form of absolute deviation and indicates the variability of the total cost of supply chain and productivity is described in the form of average workers productivity among all factories in all periods. It is assumed that all of the parameters are subject to uncertainty. The proposed model is solved with a novel hybrid algorithm composed of modifiedε-constraint method, extended Monte Carlo sampling method and L-shaped method. Finally a numerical example is generated based on some normal and uniform distributions and the result demonstrates the validity of the model as well as the efficiency of the proposed method.



---

*Source: 406398-2011-04-20.xml*