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

A multi-objective two stage stochastic programming model is proposed to deal with a multiperiodmulti-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.


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 2 Mathematical Problems in Engineering 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.
The a 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 Section 2, 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.

Problem Description
The proposed multiobjective multiproduct multisite production-distribution problem can be described as follows.
There are J 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 product i manufactured at factory j to fulfill stochastic demand of customer zone c by k-level labor, ii the number of k-level labors would be employed, laid off or trained at each factory, iii the quantity of product i 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.

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 .TC n : Total cost of supply chain under scenario n.First objective function 3.1 aims to minimize the expected total cost of supply chain, where TC n is total cost of supply chain under the realization of scenario n and defined as follows:

Multi-Objective Stochastic Production-Distribution Model
Tr n kk j • XU kk jt i,j,t 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 kskill 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.

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.

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.

Mathematical Problems in Engineering
Step 1. Select one of the objective functions Z j 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 where θ is an adequate small number usually between 10 −6 and 10 −3 , X is the feasible region of the original problem s k is a slack variable.Also, r k , g k , Nd k 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 twostage 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.

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 subproblem 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 .
Select one of the objective functions as the main one and calculate ranges for ε k by using the payoff table

Record all feasible solutions obtained in previous step
Set number of grid points g k for the k − 1 objective functions' ranges Generate N scenarios using extended Monte Carlo sampling method Run L-shaped method to construct payoff table: Min Z j : (x), k = 1, . . ., K Convert other objective functions into constraint and run L-shaped method to solve modified ε-constraint model for each vector of ε k

Numerical Experiments
Consider the supply chain network problem depicted in Figure 2. 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 F 1 , F 2 , F 3 , and F 4 which are spread geographically, and three customer centers located in three different cities C 1 , C 2 , and C 3 .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.
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 Table 2.
Figure 3 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 .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 Table 3 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.

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  function and the best value obtained from the proposed method in worst case is not more than 2.78 percent.

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.

Figure 1 :
Figure 1: Flowchart of the proposed method.

Figure 4 :
Figure 4: Pareto curve for the expected and variance of total cost.

Table 1 :
Cost items distribution functions.

Table 2 :
The payoff table for modified ε-constraint method.

Table 3 :
State of the staff upgrading versus the average workers productivity.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

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 The behavior of Z 1 against Z 2 .

Table 4 :
Comparison of the performance of the proposed algorithm with different scenario numbers.