Simulation-Optimization Approach for the Production and Distribution Planning Problem in the Green Closed-Loop Supply Chain

With the growth of multinational companies, increasing international and domestic competition between companies, upgrading information technology, and increasing customer expectations, accurate supply chain (SC) planning is essential. In such an environment, pollution has become more severe in recent decades, and with the weakening of the environment and global warming, green SC management (GSCM) strategies have become more attention in recent decades. In this research, we consider the integrated production and distribution (PD) planning problem of a multi-level green closed-loop SC (GCLSC) system, which includes multiple recycling, manufacturing/ remanufacturing, and distribution centers. We present a three-level bi-objective programming model to maximize profit and minimize the amount of greenhouse gas emissions. A hierarchical iterative approach utilizing the LP-metric method and the non-dominated sorting genetic algorithm (NSGA-II) is introduced to solve the proposed model. Also, the Taguchi approach is applied to find optimum control parameters of NSGA-II. Moreover, Monte Carlo (MC) simulation is applied to tackle uncertainty in demand, and the NSGA-II algorithm is fusioned with MC simulation (MCNSGA-II). The results obtained show that the simulation-optimization approach presented better results than the deterministic approach.

The analysis and design of PD systems have been a vital area of research over the years. Geoffrion and Graves (1974) studied a single-period PD problem and suggested a solution approach based on the Benders decomposition algorithm. This is probably the first article to suggest a general mixed-integer programming (MIP) model for the mentioned problem. Haq et al. (1991) utilized an integrated production-inventory-distribution planning formulation in a large fertilizer industry in North India that incorporates many realistic conditions such as lead times, set-up cost, and recycling of production losses as well as backlogging. Chen and Wang (1997) suggested an integrated framework for steel PD planning at a major steel producer company in Canada. Ozdamar and Yazgac (1999) introduced a hierarchical PD planning method for a multinational plant with multiple warehouses. Lee and Kim (2000) presented a general PD model in the literature and offered a solution method using a hybrid approach combining analytical and simulation techniques. Varthanan et al. (2012) investigated PD problem with stochastic demand to minimize regular, overtime and outsourced production costs. They used a simulation based heuristic discrete particle swarm optimization (DPSO) algorithm to solve the problem. Nasiri et al. (2014) suggested a PD model for a three-level program with uncertain demands. They presented a solution framework based on the Lagrangian Relaxation algorithm, which is developed by a heuristic to solve sub-problems. Niknamfar et al. (2015) considered a three-level SC with multiple production and distribution centers, and different customer areas. They presented a robust counterpart model in PD planning to minimize the total cost. Sarrafha et al. (2015) suggested a bi-objective mixed-integer non-linear programming (MINLP) formulation to design an SC network that includes suppliers, production and distribution centers, and retailers. They used the multi-objective biogeography based optimization (MOBBO) algorithm for solving the problem. Seyedhosseini and Ghoreyshi (2015) introduced a mathematical formulation for integrated PD planning for perishable products. An innovative framework is used to solve the formulation in which the suggested method first solves the production problem and afterward, the distribution problem. Devapriya et al. (2017) considered a PD planning problem of a perishable product and provided a MIP model to minimize cost. They presented a solution approach using evolutionary algorithms to solve the model. Zamarripa et al. (2016) introduced a rolling horizon approach for coordinating the PD of industrial gas SC. Zheng et al. (2016) presented a penalty function-based method to solve a riskaverse PD planning problem. The method changes the formulation into some optimization problems which can be solved by traditional optimization software. Ma et al. (2016) proposed a PD planning model applying bi-level programming for SCM and developed a genetic algorithm to solve the model. Rezaeian et al. (2016) suggested an MINLP model for the integrated PD and inventory planning for perishable products with a fixed lifetime all over a two-echelon SC by integrating production centers and distributors. Moon et al. (2016) introduced a bi-objective MIP model to design a four-stage distribution system under a carbon emission constraint. They proposed a two-phase method to solve the model and find a non-dominated solution. Osorio et al. (2017) suggested a simulation-optimization formulation to make strategic and operational decisions in production planning. They used discrete event simulation to demonstrate the flows across the SC, and MIP model to assist daily decisions. Wei et al. ( 2017) studied the integrated PD planning problem and presented a model with a two-stage production structure. They applied relax-and-fix and fix-and-optimize approaches to solve the problem. Ensafian and Yaghoubi (2017) considered an SC that consists of procurement, production, and distribution of platelets. They presented a bi-objective mathematical model to maximize the freshness of the platelets and minimize the total cost. Moreover, a robust optimization approach is used to tackle uncertain demand. Farahani and Rahmani (2017) proposed a MIP model that includes production planning, allocation-location facilitation, and distribution planning to maximize the profit of a crude oil network. Nourifar et al. (2018) introduced a multi-period decentralized SC network model with uncertainty. Uncertainty parameters such as demand and final product prices were defined by stochastic and fuzzy numbers. They presented a solution framework based on the Kth-best algorithm, chance constraint approach, and fuzzy approach. Rafiei et al. (2018) investigated a PD planning problem within a four-echelon SC. The problem is modeled in two non-competitive and competitive markets to minimize total chain cost and maximize service level. Casas-Ramírez et al. (2018) studied an SC, including factories and depots. They proposed a MIP model to balance the total workload and minimize the total cost of the SC. To solve this model, they used an adapted bi-objective GRASP to find non-dominated solutions. Jing and Li (2018) considered a multi-echelon closed-loop SC planning problem involving a joint recycling center, multiple remanufacturing/manufacturing centers, and multiple distribution centers decentralized to various areas. The solution framework was designed by a hierarchical iterative approach based on the self-adaptive genetic algorithm. Goodarzian et al. (2021) studied a novel multi-objective formulation is devised for the PD problem of a supply chain that consists of several suppliers, manufacturers, distributors, and different customers. Due to the NP-hardness of problems, NSGA-II and Fast PGA algorithm are applied. Pant et al. (2021) developed a bi-objective CLSC model for the paper industry under uncertain environment. The first objective of the proposed model is to maximize SC surplus, and the second objective is to incorporate sustainability through minimizing carbon content by reducing the number of trucks between various echelons of CLSC network. They considered uncertainty at demand points and applied MC simulation to handle it.
The reviewed articles on integrated production-distribution planning are summarized in Table 1. The last row of Table 1 is for the present study. Specifically, the followings are the significant contributions of this paper:  A three-level bi-objective programming model is presented to maximize profit and minimize the amount of greenhouse gas emissions.  NSGA-II algorithm is developed to solve the bi-objective model at each level.  A hierarchical iterative approach is applied to solved the three-level model.  NSGA-II algorithm is combined with MC simulation (MCNSGA-II) to tackle uncertainty in demand. Geoffrion and Graves (1974) 1974 * * * * -Bender method - Haq et al. (1991) 1991 * * * * -Exact - Chen and Wang (1997) 1997 * * * * -Exact - Ozdamar and Yazgac (1999) 1999 * * * * -Exact GAMS Lee and Kim (2000) 2000

3-Mathematical modeling
In this paper, a three-level bi-objective programming model is presented for a GCLSC, which includes multiple recycling centers, multiple manufacturing/remanufacturing factories, and multiple distributors.
Multi-level programming is an optimization approach that has a multi-layer hierarchical form. In this form, decision-makers of different levels have various decision-making authority and objectives. The set of strategies and the goal attainment of lower levels can be affected by the decisions from the upper levels. Nevertheless, lower levels have considerable autonomy, so the upper levels cannot totally control them (Jing and Li 2018). Regarding the basic concept of multi-level programming, each subordinate level in this model is situated in different region and has its policies and strategies. The first-level decision maker (recycling centers) sets own goals and/or decision, and then asks each subordinate level of the organization for their optima, that is calculated in isolation. Because of it, Multi-level programming is used.
In the considered GCLSC, returned products are disassembled to components in the multiple recycling centers. Testing activities to know the defective components perform in the centers. The components can be saved in the inventory of eligible components if they have a remanufacturing standard; otherwise they will be excreted. Eligible components are carried to various remanufacturing/ manufacturing centers.
The remanufacturing/manufacturing centers test the quality and efficiency of Eligible components, or raw materials are converted to new components. Remanufactured and new components are used to produce remanufactured and new products, respectively. The products are added to the inventory of remanufactured and new products and are carried to various distribution centers according to orders.
Each distribution center can select one or more remanufacturing/manufacturing centers for satisfying demand, and provide remanufactured or new products to retailers for sales. Moreover, the returned products from retailers are gathered and recycled by distribution centers and carried to multiple recycling centers. Note that there are several kinds of vehicles to transport parts and products between the GCLSC components.
Furthermore, the assumptions for the considered problem are as follows: 1. The SC is multi-product, multi-period, and multi-level. 2. Demand for remanufactured and new products is uncertain and can only be met by themselves, and remanufactured products have lower selling prices than new products. 3. The production cost of each factory is different. 4. Vehicles have variable capacity in each level. 5. The cost of transfer between the nodes is different. 6. Processing and reprocessing parts in the manufacturing/remanufacturing factories lead to greenhouse gas emissions into the environment. 7. The transportation of vehicles at each level in each period has a greenhouse gas emissions limit. The notations are introduced in Table 13 to Table 19 in the Appendix. The model of each level is as below.

3-1-The recycling centers, the first level
The bi-objective MIP model at the first level (FLM) to maximize profit and minimize adverse environmental effects is as follows: .
Expression (1) maximizes the profit of the recycling centers. The profit is equal to revenue minus total costs, including recycling, disassembly, testing, disposal, inventory, and transportation costs. The objective function (2) minimizes the amount of greenhouse gas emissions of each vehicle per kilometer for a unit of eligible components between the recycling centers and remanufacturing/manufacturing factories. Constraint sets (3) and (4) are the inventory equations for returned products and eligible components. Constraint set (5) shows the quantity constraint for disposed of components. Constraint sets (6) and (7) ensure that the inventory rate of returned products and eligible components do not exceed the maximum level. Constraint set (8) represents that the quantity of returned products to be disassembled and tested in a recycling center does not exceed the maximum amount. Constraint set (9) guarantees that the quantity of eligible components that are transported by a vehicle from the recycling centers to the remanufacturing/manufacturing factories does not exceed the capacity of the vehicle. Constraint set (10) ensures the sum of greenhouse gas emissions for the transportation system between recycling centers and factories does not exceed the maximum allowance. Constraint sets (11) and (12) describe the value ranges of the variables.

3-2-The manufacturing/remanufacturing centers, the second level
The bi-objective MIP model at the second level (SLM) to maximize profit and minimize adverse environmental effects is as follows: 1 1 1 1 1 1 1 1 1 1 1 .. st ,1 11 Expression (13) maximizes the profit of the remanufacturing/manufacturing factories. The profit is equal to sales revenue minus total costs, including manufacturing, remanufacturing, inventory, purchase, and transportation costs. The objective function (14) minimizes the amount of greenhouse gas emissions of each vehicle per kilometer for a unit of new products and remanufactured products between remanufacturing/manufacturing factories and distribution centers, as well as the amount of greenhouse gas emissions during processing and reprocessing. Constraint sets (15)-(19) represent the inventory equations for eligible, new, remanufactured components, and products. Constraint sets (20)-(24) ensure that the inventory quantities of eligible, new, remanufactured components and products do not exceed the maximum quantity. Constraint sets (25)-(28) guaranty that the quantity of processed and reprocessed components, as well as the quantity of new and remanufactured products, do not exceed the capacity of relevant factories. Constraint set (29) considers that the amount of new and remanufactured products that are transported by a vehicle from the manufacturing/remanufacturing factories to distribution centers does not exceed the capacity of the vehicle. Constraint set (30) ensures the sum of greenhouse gas emissions for the transportation system between factories and distributions centers does not exceed the maximum allowance. Constraint set (31) represents that the sum of greenhouse gas emissions for processing and reprocessing operations in factories must be less than the maximum allowance. Constraint sets (32) and (33) describe the value ranges of the variables.

3-3-The Distribution centers, the third level
The bi-objective MIP model at the third level (TLM) to maximize profit and minimize adverse environmental effects is as bellows: 1 1 1 1 1 1 1 Expression (34) maximizes the profit of the distribution centers. The profit is equal to sales revenue minus total costs, including the purchase of new and returned products, shortage, inventory, and transportation costs. The objective function (35) minimizes the amount of greenhouse gas emissions of each vehicle per kilometer for a unit of returned product between distribution and recycling centers. Constraint sets (36)-(38) show the inventory equations for new, remanufactured, and returned products. Constraint set (39) ensures that the number of collected returned products does not exceed the number of returned products available in the market. Constraint sets (40) and (41) represent the quantity shortage for new and remanufactured products. Constraint sets (42)-(44) guaranty that the inventory rate of remanufactured, new, and returned products do not exceed the maximum level. Constraint set (45) ensures that the quantity of returned products, which are transported by a vehicle from distribution centers to recycling centers, does not overstep the capacity of the vehicle. Constraint set (46) guarantees that the sum of greenhouse gas emissions for the transportation system between distribution and recycling centers does not exceed the maximum allowance. Constraint set (47) defines the value ranges of the variables.

4-Solution approach
The multi-level linear programming is NP-hard (Hansen et al. 1992), so finding its solutions are usually hard. In this paper, the proposed three-level bi-objective programming formulation is solved based on a hierarchical iterative approach using the LP-metric method for small-size instances and the NSGA-II algorithm for large-size instances, respectively. Moreover, the NSGA-II algorithm is combined with MC simulation (MCNSGA-II) to tackle uncertainty in demand.

4-1-LP-metric method
LP-metric is a classical approach which is applied to solve multi-objective models. The method attempts to find the best solution, which has the shortest distance from the ideal solution. Therefore, the LP-metric function (48) is applied to compute the distance of an available solution to the ideal solution .
fx and j w represent the th j objective function and its importance degree, respectively.
p demonstrates the emphasis degree on available deviations. To calculate the ideal value of the th j objective function ( * () j j fx ), the model with the th j objective function and the existing constraints is solved. Moreover, to find the anti-ideal value of the th j objective function ( () j j fx ), the objective function is reversed. Finally, we minimize the LP-metric function (48) subject to the constraints to obtain the solution.

4-2-The NSGA-II algorithm
In this paper, the NSGA-II algorithm (Deb et al. 2002) is applied for solving the considered bi-objective models at each level. Note that the algorithm is adopted to solve medium and largescale instances.

4-2-2-Crossover
The crossover operator is used to transfer the same characteristics from parents to nextgeneration children. The simplest types of crossovers are single-point, two-point, and uniform crossover. In this paper, the two-point crossover is applied. First, two parents are selected, and then two numbers are randomly selected to determine the cut-off area. The first child inherits the cut region from the second parent, and the second child inherits it from the first parent. Figure 2 shows the crossover operator in this study.

4-2-3-Mutation
The mutation operator is usually performed with a very low probability. In the present study, two chromosomes from the parents are randomly selected, and the position of the two chromosomes is changed, as shown in Figure 3.

4-3-The hierarchical iterative approach
The proposed three-level model is solved based on a hierarchical iterative algorithm. The description of the algorithm is as below: Step 1: Determine the initial value of variable jkpvt da . It is randomly generated from formula (45) and (46) at iteration 0 ( 0 Iter  ).
Step 2: Solve the FLM using the values of  will be considered for the FLM as a parameter in a conceivable new iteration.
Step 5: Check that the stopping condition has been met, that is: where * w is an iteration accuracy given in advanced. If so, stop iteration; if not, start Step 2. Note that the NSGA-II algorithm is implemented in the context of the above procedure. The NSGA-II is a population-based algorithm, so the Mean Ideal Distance (MID) criterion is applied to calculate w ̅ 1 , w ̅ 2 , and w ̅ 3 . By considering the equation (51) to calculate MID, the equation (52) is used to calculated w ̅ 1 , w ̅ 2 , and w ̅ 3 in Step 5.

4-4-Monte Carlo Simulation
MC simulation is premier to a deterministic simulation of a system when the input variables of the system are random. In this simulation approach, a single value is selected for each input random variable based on the best guess by the modeler. The approach is then run, and the output is considered. This output is a single value or a single set of values based on the selected input. MC simulation randomly samples values from every input variable distribution and applies the sample to compute the model's output. This procedure is iterated many times until the modeler achieves a sense of how the output changes given the random input values (Sokolowski and Banks 2010).
In this paper, it is assumed that the demand parameter in the third level (distribution centers) is uncertain. We use MC simulation to deal with the uncertainty in the NSGA-II algorithm. Therefore, MC simulation is applied to appraise the objective functions of each chromosome, considering the probability distribution of the demand that is obtained from historical data. Moreover, the expected value of each objective function related to a chromosome is estimated by simulation. Moreover, we calculate the number of simulation iterations using equation (53).

4-5-Comparison of multi-objective solution methods
Two metrics, set coverage and spacing, are usually applied to evaluate convergence and diversity of a multi-objective solution method. We apply the following metrics to assess the obtained non-dominated solutions the NSGA-II and MCNSGA-II algorithms. Zitzler and Thiele (1999) proposed the set coverage metric   , C A B that computes the ratio of solutions in B are weakly dominated by solution A .

4-5-1-Set coverage metric
C A B  means that all members of B are weakly dominated by A .   ,0 C A B  shows that no member of B is weakly dominated by A . Scott (1995) presented the following metric to calculate the comparative distance between successive solutions in the non-dominated set (Q ). Based on the above-mentioned metrics, an algorithm with a greater C and a smaller S is preferred.

5-Computational results
In this section, numerical experiments are conducted to investigate the performance of the three-level bi-objective MIP model, the NSGA-II algorithm, and the MCNSGA-II algorithm based on random data. First, we present how the test problems are generated. The Taguchi approach is then employed to tune the parameters of the NSGA-II algorithm, and the best level of each parameter is obtained by the signal-to-noise diagram in the Minitab software. Sensitivity analysis is then performed for different weights of LP-metric function, and the best weight for each level is selected. The results obtained by the LP-metric method and the NSGA II algorithm are compared, and the relative distance of each sample is determined. Note that the LP-metric method and the NSGA II algorithm are implemented in GAMS and MATLAB, respectively, and are tested on a computer with a 2.3 GHz CPU and 8.0 GB of RAM. Finally, the MCNSGA-II algorithm is implemented on two instances, and the results are compared with the NSGA-II algorithm.

5-1-Data generation
We need to test the suggested mathematical model and the solution approach using randomly generated test problems. The dimension of the test problems that are categorized into small, medium, and large-scale are introduced in Table 5. Moreover, the values of model parameters are randomly generated based on a uniform distribution that are shown in Table 6.

TE
Uniform 

TE
Uniform 

5-2-Parameter tunning of NSGA-II
The Taguchi experimental design approach is utilized to calibrate the parameters of the NSGA-II algorithm. The approach minimizes the effect of noise and specifies the optimal level of signal factors. To do so, the signal to noise ratio (S/N), that computes the value of variation of the response, is implemented. Then, the method aims to maximize the S/N ratio (Peace 1993). In this paper, MID and maximum scattering (MS) metrics calculate response based on equation (56). Table 7 exhibits different levels of the factors for NSGA-II algorithm parameters at the recycling centers level to run the Taguchi method. where Pc and Pm are the crossover and mutation probability, respectively. Then, applying Minitab Software, the L9 design is applied for the NSGA-II algorithm, and experimental results are illustrated in Table 8. The above approach is applied to parameter tuning of the NSGA-II algorithm at the second and the third levels. Figure 5 indicates the effect plot of the S/N ratio related to each level. According to Figure 5, the best parameters values of the NSGA-II algorithm for each level are indicated in Table 9.

5-3-The weight coefficients in LP-metric method
It is necessary to determine the suitable weight coefficients of the objective functions for executing the proposed model in GAMS software based on the LP-metric method. Different weight combinations are solved for a small-size determinant problem, and the appropriate weight is generalized to the other problem. Figure 6 illustrates how the value of the objective functions is affected by various weights in different levels of the problem. The most significant variations in the objective functions happen in the weight coefficients between 0.4 and 0.9, and in the other weights, the objective function value becomes zero. So, it is better to weigh the goals between these limits. According to Figure 6

5-4-Validation of the NSGA-II algorithm
The LP-metric method is used for model validation. Thus, six small-scale problems are solved in GAMS software with the LP-metric method. Note that it is not possible for large-scale problems. Moreover, for comparison, the six small problems are solved with the NSGA-II algorithm. Since the NSGA-II is run based on the assumption that the initial population is randomly generated; hence, any test problem is solved ten times. The best non-dominated solution is selected from the ten times of running results, and its LP-metric measure is calculated. The measured value is compared with the value obtained from the exact method that is run in GAMS. In Table 10, the relative gap between the values of the two LP-metric measures, which are achieved from the exact method and the NSGA-II algorithm, is calculated. According to Table 10, the average of the relative gaps is 7%, 7%, and 1% for the three levels, respectively. Therefore, we can declare that the NSGA-II algorithm is suitable to solve the considered problem.

5-5-Analysis of MC Simulation
In this research, it is assumed that the demand for remanufactured and new products in distribution centers at different periods are stochastic. The probability distribution function associated with the value of demand is uniform. The MC simulation approach is used to deal with the uncertainty Equation (53) is used to calculate the number of simulation replications. First, a pilot run of MC simulation is implemented considering an initial sample size of 10 replications to calculate mean and the standard deviation of 1 and 2 for a chromosome. By considering a confidence level of 95% and an error of 0.1, the number of simulation replications is obtained equal to eight based on equation (53).
The medium and large-size problems are solved to compare the NSGA-II and MCNSGA-II algorithms. Figure 7 and Figure Table 11 and Table 12 indicates the value of convergence (C) and spacing (S) metrics of the different optimization methods. As can be seen, the MCNSGA-II algorithm presents better non-dominated solutions than the NSGA-II algorithm.

6-Conclusions
In this paper, we considered an integrated PD planning problem of a multi-level GCLSC system, that includes multiple recycling centers, multiple remanufacturing/manufacturing centers, and multiple distribution centers. A three-level bi-objective MIP model is presented to maximize profit and minimize the amount of greenhouse gas emissions. A hierarchical iterative approach using the LP-metric method and the NSGA-II algorithm is suggested to solve the proposed model. The Taguchi experimental design approach is applied to find optimum control parameters of NSGA-II. Moreover, the NSGA-II algorithm is fusioned with MC simulation (MCNSGA-II) to deal with the uncertainty of the product demand in distribution centers. The results obtained show that the simulation-optimization approach presented better results than the deterministic approach. Furthermore, there are several opportunities in future research. First, to understand the realities of the environment, it is important to consider uncertainty in other influential parameters such as purchase costs, shipping costs, and so on. Second, studying the literature on multi-level production-distribution planning shows that fewer study goals, such as customer satisfaction and reduction of delay time in the model, were considered. Adding the aforementioned goals could be of interest for future research. Third, in this and other studies, multi-level planning has been solved at two or three levels. Adding the fourth level as the supplier level can be an appropriate area for future research.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest regarding the publication of this paper.
Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.      x Inventory quantity of remanufactured product p in distributor j at the end of period t D jpt  Inventory quantity of return product p in distributor j at the end of period t jkpvt da This notation has taken value in the FLM ijpvt fdn This notation has taken value in the SLM ijpvt fdr This notation has taken value in the SLM

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict ofinterest.