A fuzzy-stochastic multi-objective model for sustainable planning of a closed-loop supply chain considering mixed uncertainty and network flexibility

Closed-loop supply chain network design (CLSCND) has been increasingly spotlighted over the latest decade. The focus has been given to maximize the economic performance, resource utilization and sustainability through incorporating a holistic decision-making on both forward and reverse logistics. In this paper, a new fuzzy-stochastic multi-objective mathematical model is formulated for sustainable CLSCND. The model aims at balancing the trade-off between cost effectiveness and environmental performance under different types of uncertainty. The environmental performance of CLSCND is measured by carbon emission. Moreover, the network flexibility is modeled and incorporated in the decision-making so that customer demands can be fulfilled by different means. In order to solve the complex optimization problem, the model is first defuzzilized and converted into an equivalent crisp form. Then, a sample average approximation based weighting method (SAAWM) is developed to obtain a set of Pareto optimal solutions between cost and carbon emission under different uncertain environments. The model is validated through a set of numerical experiments. The computational results show, through the incorporation with network flexibility, the proposed mathematical model and solution approach can effectively generate consistent objective values and solutions over different scenario trees and obtain robust strategic decisions on facility locations. Meanwhile, the flexibility and rationality of the decision-making on transportation management, demand allocation and facility operations can be improved as well. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Supply chain management (SCM) refers to the effective and efficient management of the flows of materials, information and capital throughout different actors of a supply chain including supplier, manufacturer, distributor, retailer and customer (Chopra and Meindl, 2016). Conventionally, the objective of SCM is to maximize the supply chain surplus or profit through decision-makings at strategic, tactical and operational levels. In recent years, technological advancement and economic boom have not only resulted in an improvement on people's living standards and lifestyle but also brought significant environmental and social challenges (Darbari et al., 2019;Yu and Solvang, 2016b), e.g., resource depletion, water pollution and waste generation, etc. Due to this reason, tremendous focus on environmental friendliness and sustainable development has been given by worldwide governments, companies as well as the whole society. Defined by the United Nations, sustainable development is "development that meets the present without compromising the ability of future generations to meet their own needs" (Imperatives, 1987), which is supported by economic prosperity, environmental protection and social harmony. To achieve the goal of sustainable development, stringent environmental regulations have been implemented across the globe, with which manufacturers have to be involved in the management of the end-of-life (EOL) returns and to take responsibility of the entire lifecycle of their products. This has dramatically shifted the focus and paradigm of traditional SCM. New concepts and practices, e.g., green logistics (Rodrigue et al., 2008), sustainable SCM (Carter and Liane Easton, 2011), and reverse logistics (Carter and Ellram, 1998), have been proposed and implemented in order to achieve a balance between economic performance and environmental impact. A closed-loop supply chain (CLSC) incorporates the reverse logistics activities into a traditional forward supply chain in order to maximize the utilization of resources and improve both economic and environmental benefits. Compared with the individual design and management of forward supply chain and reverse logistics, which may lead to sub-optimal decision-makings (Choudhary et al., 2015), a CLSC aims at determining the global optimal solution through a holistic analysis.
Supply chain network design (SCND) is one of the most important strategic decisions in SCM (Chopra and Meindl, 2016), which consists of two levels of decisions (Cakravastia et al., 2002). Table 1 gives a comparison between the two levels of decisions in SCND. The first-level decisions determine the optimal network configuration through locating facilities with different functions. The second-level decisions determine the optimal use of the network structure through allocating customer demand to different facilities and formulating the transportation strategy on each itinerary. From decision-making perspective, the first-level decisions are at strategic level and have thus a long-term impact on the performance of a supply chain. Once those decisions are made, it is extremely expensive and time consuming to alter, so the first-level decisions in SCND should be robust in order to withstand the possible changes on external environment for many years to come. On the other hand, the second-level decisions are made based upon the firstlevel decisions. They are at tactical and/or operational levels, which are more flexible and can easily be changed in order to adapt the change of external environment and to optimize the performance of a supply chain. Compared with traditional forward SCND, the close-loop supply chain network design (CLSCND) is more complicated due to the involvement of both forward and reverse logistics activities. In addition, more companies and hybrid facility functions may be included so that the material flows on both channels can be handled in the CLSC.
CLSCND is a complex decision-making problem that aims at simultaneously optimizing several conflicting goals, e.g., economic efficiency, responsiveness to customer demand, and environmental impact. Hence, the proper management of the trade-off among these goals is of significant importance in the planning of a CLSC. Furthermore, CLSCND involves both strategic and operational decisions (Chopra and Meindl, 2016). Once the first-level strategic location decisions are made, it is unlikely to alter them within a short period due to the high expenses. In this regard, a robust first-level decision may achieve long-term benefits in both profitability and sustainability of a CLSC. However, within the planning horizon of a CLSC, the network decisions are inevitably suffered from different types of uncertainty (Talaei et al., 2016;Subulan et al., 2015a). Moreover, some flexible network operating options of capacity adjustment in the second-level, e.g., outsourcing, hire of temporary workers and rent of equipment, etc., may also yield significant influence on the firstlevel location decisions and should thus be considered holistically.
To this end, this papers aims, through advanced mathematical programming approach, optimization and analysis, at answering the following research questions: How to improve the environmental friendliness with a minimum compromise on economic performance in CLSCND? How to manage different types of uncertainty in CLSCND in order to generate robust strategic decisions? How can the first-level location decisions of CLSCND be affected by incorporating the network flexibility at the operational level?
The rest of the paper is structured as follows. Section 2 gives a broad literature review on the recent research works of CLSCND and identifies the contributions of the current research. Section 3 formulates the mathematical model for CLSCND under study. Section 4 develops a solution method in order to effectively solve the complex optimization problem. Computational experiments and analysis are given in Section 5. Finally, section 6 concludes the paper and specifies future improvements.

Literature review
Over the years, extensive research works have been conducted for the development of advanced decision-support methods for both forward supply chain (Saberi et al., 2018) and reverse logistics Kuş akcı et al., 2019). Recently, the incorporation of reverse logistics activities into traditional forward SCND has increasingly been focused in the holistic planning of a CLSC. In this regard, comprehensive literature reviews have been conducted by Govindan et al. (2015) and  in order to summarize the development of sophisticated mathematical models, improved computational methods and real-world case studies for CLSCND. This paper focuses on the latest modeling efforts and practices in CLSCND. In order to make a comparison with the earlier research works in this field, the mathematical models are categorized, based on their characteristics in the development environment (deterministic vs. uncertainty) and model objectives (single vs. multiple), into four types: 1. Deterministic model with single objective 2. Deterministic model with multiple objectives 3. Uncertain model with single objective 4. Uncertain model with multiple objectives Table 1 Comparison between first-level decisions and second-level decisions in SCND.

First-level decisions
Second-level decisions Determining the network configuration through selecting the number and locations of facilities Determining the operations of the network through the selection of transportation mode and allocation of demands and material follows on different facilities and links Strategic level Tactical and/or operational levels Long-term impact on the performance of a supply chain Medium-/short-term impact on the performance of a supply chain Very expensive to be changed Easy and inexpensive to be changed Should be robust to withstand the change of external environment Should be flexible to adapt the change of external environment

Deterministic model with single objective
The primary target for a CLSCND problem is to improve the economic performance through the value recovery from EOL products. € Ozceylan et al. (2017) formulated a linear program to maximize the profit generation of an integrated forward and reverse logistics system in automotive industry. Considering multiple products, alternative plants, retailers and suppliers, Amin et al. (2017) proposed a profit maximization model for designing a tire production CLSC system. Several scenarios were tested for analyzing the performance of the optimal result. From the economic efficiency perspective, the performance of a CLSC system can also be evaluated by the overall cost. In this regard, Yi et al. (2016) investigated a mathematical model for a real-world retailed-oriented CLSC network planning of EOL construction equipment.
The CLSCND is a complex optimization problem that usually requires high computational efforts. Hence, several research works have been conducted in order to provide enhanced solution approaches for improving the computational efficiency. Chen et al. (2015) formulated an integer program for CLSCND considering the recycling decisions, and an improved two-stage genetic algorithm (GA) was developed and validated through the comparison with the result obtained from LINGO solver in small instances. Combining with two metaheuristics: particle swarm optimization (PSO) and GA, Soleimani and Kannan (2015) investigated an enhanced algorithm to solve a multi-product multi-period CLSCND problem. The result was compared with that obtained from MATLAB and CPLEX in small sized problems.

Deterministic model with multiple objectives
Nowadays, the increasing focus on environmental protection, sustainable development and customer satisfaction has shifted the traditional CLSCND towards finding the optimal balance among economic benefits, environmental impact and other influencing factors. In order to balance the cost and carbon emission, Taleizadeh et al. (2019) proposed a bi-objective mixed integer model for CLSCND considering pricing decisions and discount of returned products. The objective functions were combined with the fuzzy Torabi-Hassini (TH) method. Hasanov et al. (2019) formulated an optimization model for balancing the trade-off between cost and greenhouse gas (GHG) emission of a four-stage CLSC. The model considers the remanufacturing in reverse logistics and incorporates with the inventory decisions. Garg et al. (2015) developed a biobjective nonlinear optimization model for CLSCND considering both profit and GHG emission. The GHG emission was reduced by minimizing the total use of trucks in the forward channel and an interactive fuzzy approach was used to combine different objectives.
Taking into account of the balance between profit and energy consumption, Kadambala et al. (2017) investigated a bi-objective model for CLSCND, and both multi-objective particle swarm optimization (MOPSO) and non-dominant sorted genetic algorithm (NSGA-II) were applied to solve the proposed model. Zarbakhshnia et al. (2019) proposed a multi-objective model for sustainable CLSCND. The model aims at simultaneously minimizing cost, carbon emission and number of machines in line, and the three objective functions are combined with ε-constraint method. Paksoy et al. (2019) investigated a multi-objective optimization problem for balancing different cost components, carbon emission and percentage of late delivery of raw materials through the decisionmaking on facility location, demand allocation and tour of vehicles in CLSCND. Zohal and Soleimani (2016) formulated a multiobjective integer program for simultaneously maximizing the revenue, minimizing the cost and carbon emission in the planning of a CLSC in gold industry. In addition, an enhanced ant colony algorithm was investigated for improving the computational efficiency of the proposed model.

Uncertain model with single objective
CLSCND is a long-term strategic decision, so exact values of the model input within the planning horizon are difficult to forecast. Therefore, modeling efforts have been given to improve the decision-making of CLSCND under uncertainty; among which stochastic programming is the most extensively used technique. Considering the uncertainty of product quality, Jeihoonian et al. (2017) formulated a two-stage stochastic model for CLSCND and a scenario reduction approach was adopted in order to reduce the size of large optimization problems. Zhen et al. (2019) investigated a stochastic nonlinear optimization problem for minimizing the total cost of CLSCND considering uncertainty from both demand and return. The model was solved by an enhanced Tabu search algorithm. Baptista et al. (2019) developed a multi-period twostage stochastic program for CLSCND aiming at maximizing the profit expectation under uncertainty. The hybrid chance-constraint and second order stochastic dominance risk averse measures were investigated for risk management, and a fixed-and-relax decomposition method was proposed to solve the model. Hajipour et al. (2019) proposed a stochastic optimization model for traceable CLSCND and two meta-heuristics: greedy randomized adaptive search procedure (GRASP) and PSO, were employed to solve the problem.
Almaraj and Trafalis (2019) formulated a robust optimization model for minimizing the total cost of CLSCND. An integrated stochastic and robust model was investigated by Keyvanshokooh et al. (2016), where a Latin Hypercube Sampling method was adopted for scenario reduction and a Benders decomposition method was developed for accelerating the computation. Considering the uncertainty from demand, return as well as carbon tax rate, Haddadsisakht and Ryan (2018) proposed a hybrid stochastic and robust optimization model for minimizing the total cost in CLSCND. The model formulated multiple alternative transportation modes and a Benders decomposition method was applied to solve the problem.
Different from the stochastic perspective that focuses on random uncertainty, Torabi et al. (2016)

Uncertain model with multiple objectives
Several research works have been done for simultaneously managing the conflicting objectives under an uncertain environment. Considering the optimal trade-off between cost and additional demand facility of service, Saedinia et al. (2019) proposed a robust bi-objective mathematical model for CLSCND in oil and gas industry. The fuzzy TH was employed to generate Pareto optimal solutions, and the model was validated through a case study in Iran. Ahmadi and Amin (2019) investigated a chance-constraint stochastic program for designing a mobile phone CLSC. The model aims at balancing the profit generation and the weights of eligible suppliers, which are estimated through a fuzzy method. Vahdani and Mohammadi (2015) formulated an uncertain and responsive CLSCND problem with a bi-objective interval-stochastic robust programming approach. The model aims at simultaneously minimizing the total cost and the waiting time. Das and Posinasetti (2015) proposed a multi-objective stochastic model for simultaneously maximizing profit generation, minimizing energy consumption and minimizing harmful gas emission. Mohammed et al. (2019) developed a stochastic optimization model for planning a multi-period green CLSC. Four different carbon policies: carbon cap, carbon tax, carbon tax-and-trade and carbon offset were incorporated with the cost minimization objective in order to test their effectiveness in the reduction of carbon emission.
Taking into account of the balance among the total profit, the fulfillment of customer demand and the missing working days due to occupational accidents under epistemic uncertainty, Soleimani et al. (2017) formulated a fuzzy multi-objective model for sustainable CLSCND. Asim et al. (2019) proposed a fuzzy multiobjective goal programming approach for the planning of a CLSC. The model aims at determining the optimal trade-off among total cost, total delivery time and defect rate. Combining with a AHP-TOPSIS method and a fuzzy optimization model, Darbari et al. (2019) investigated the sustainable design of a CLSC for an Indian laptop manufacturer. The model aims at simultaneously balance the economic, environmental and social sustainability through the optimal decision-making on facility locations and transportation under fuzzy demand and capacity. Zhalechian et al. (2016) formulated a possibilistic-stochastic model for the decision support on facility locations, rout planning and inventory management of a CLSC. The model aims at minimizing both the cost and the carbon emission, while at the same time, maximizing the job creation. Talaei et al. (2016) proposed a fuzzy-robust optimization model for green CLSCND. The model balances the cost and the environmental impact.
Considering both fuzziness and randomness related to the input parameters, Subulan et al. (2015b) developed a fuzzystochastic multi-objective mixed integer program for balancing the trade-off among total cost, total collection coverage and risk in the design of a CLSC for lead battery. Tosarkani and Amin (2019) investigated a fuzzy-stochastic bi-objective optimization model for sustainable CLSCND. The model balances the trade-off between profit generation and environmental compliance of different actors. Table 2 summarizes the recent publications in CLSCND. Compared with the earlier models (Govindan et al., 2015), the recent research works have increasingly focusing on the proper treatment of the conflictions among multiple objectives (e.g., cost or profit, carbon emission, risk, job creation, etc.) and the uncertainties related to the model input in CLSCND. Considering the different features, e.g., randomness vs. epistemic uncertainty, and statistical dependent vs. non-statistical dependent, etc., fuzzy programming and stochastic programming/robust programming have been applied to formulate different types of uncertainty.

Summary and contribution of the research
Even though extensive research works have been done to tackle the deficiencies of the earlier models, the literature review illustrates several gaps in CLSCND.
CLSCND is a strategic decision that involves multiple conflicting objectives and different types of uncertainty (randomness and epistemic uncertainty). However, due to the increasing model complexity, most research works deal with these problems separately and may thus result in sub-optimal solutions. There is still a lack of optimization models for simultaneously tackling all the aforementioned modeling challenges. To our knowledge, only two research works have been done to address this issue (Subulan et al., 2015a;Tosarkani and Amin, 2019). Even though modelling efforts have been given by Subulan et al. (2015a) and Tosarkani and Amin (2019) to tackle the aforementioned challenges, these as well as many other stochastic models formulated for CLSCND only consider one fixed scenario tree and neglect the inherent uncertainty related to the generation of scenario trees. Without a rigorous validation over several scenario trees generated from the same probability distribution, the quality of solution may be compromised and the use of the stochastic models becomes thus limited. Besides, both models were solved by exact solution methods that are incapable to deal with large problems. Network flexibility of a supply chain is usually considered at the operational planning level (Fiorotto et al., 2018). However, recent research works have revealed the network flexibility may yield significant impact on the location decisions and should thus be considered holistically in the network design (Yu and Solvang, 2018). In this respect, there is a lack of research efforts for analyzing the impact of network flexibility in both model development and practical decision-making of CLSCND.
To overcome the aforementioned shortcomings and fill the literature gap, a new fuzzy-stochastic multi-objective model for CLSCND is first proposed in this paper, which simultaneously considers the balance between different objectives and the treatment of different types of uncertainty. In order to effectively solve the complex fuzzy-stochastic multi-objective optimization problem and to test the model's behavior over a number of different scenario trees, the fuzzy objectives and constraints are first converted into their equivalent crisp forms. Then, a sample average approximation based weighting method (SAAWM) is developed to reduce the size of the original problem and to heuristically obtain a set of Pareto optimal solutions with a high level of confidence. Finally, the network flexibility is modeled and its impact on the first-level location decisions is analyzed via numerical experiments.  (Sahyouni et al., 2007). The multi-layer CLSC includes both forward and reverse channels. In the forward direction, the products are assembled at the production plants and then transported to customers via distribution centers. In the reverse channel, the EOL products generated at customer locations are first collected, sorted and inspected at the collections centers. Based on the quality, the EOL products will be sent to respective facilities for repair/refurbishing, remanufacturing/recycling, and waste disposal. The repaired/refurbished products and components will be re-entered and sold through the forward distribution channels. The remanufactured/recycled parts and components will be sent to and reused at the production/assembly plants, through which the amount of raw materials and components purchased from the suppliers can be reduced. The repair, refurbishing, remanufacturing, recycling and re-entering of the forward supply are value recovery processes.  On the other hand, the non-recoverable EOL products and components will be sent to waste disposal centers for proper treatment.

Modeling the uncertainty
CLSCND is a long-term decision and inherent uncertainty may exist in all the input parameters. Neglecting the impact from uncertainty may significantly compromise the decision-making of the network configuration and limits the use of decision-support models. Considering the nature of uncertainty, two techniques have been extensively used in modeling different types of uncertainty: fuzzy programming and stochastic programming (Govindan et al., 2015). Fuzzy programming is used to formulate the epistemic uncertainty that depicts the imprecision of the information due to inaccurate, incomplete and insufficient historical data (Pishvaee and Torabi, 2010;Pishvaee and Razmi, 2012) and is thus nonstatistical dependent (Zhu, 2014;Subulan et al., 2015a). On the other hand, stochastic programming is applied to model the uncertainty featured with randomness, where statistical distribution based on historical data is used to estimate future conditions. Table 3 presents a comparison between the two techniques in CLSCND.
As shown, the meaning and the conversion method of fuzzy parameters and stochastic parameters are by no means identical. A fuzzy programming model can be converted into a crisp model that is at the same size and level of complexity (Jim enez, 1996).
However, on the other hand, a stochastic program is converted to a deterministic form by adding scenario-dependent components and constraints, which will significantly increase the size and complexity of the problem (King and Wallace, 2012). In this paper, considering the nature of uncertainty and the computational complexity, the customer demands, the rate of EOL returns and the quality of EOL products are modeled as stochastic parameters. The uncertainty related to the other parameters in CLSCND are formulated as fuzzy parameters.
The benefits of including both fuzzy and stochastic parameters in CLSCND are summarized as follows: Model's flexibility: improving the model's flexibility to better depict different types of uncertainty and the model's capability to deal with the incompleteness and inaccuracy of historical data. Optimality of the objective value: enabling the flexibility in the objective value of the model by formulating fuzzy parameters related to the cost and the carbon emission of CLSCND. Feasibility of the constraints: enhancing the robustness of the first-level location decisions under a dynamic environment by modeling the uncertainties featured with high randomness in the constraints. Computational efficiency: the problem size and the computational efficiency can be maintained at a manageable level without excessive use of stochastic parameters.  Table 3 Comparison between fuzzy programming and stochastic programming in CLSCND.

Fuzzy programming Stochastic programming
Epistemic uncertainty Randomness Non-statistical dependent on historical data Statistical dependent on historical data Measures the degree to what extent, for example, A is larger than B Measures the probability of, for example, A is larger than B Static uncertainty and can be converted to a crisp model Dynamic uncertainty and can be converted by adding scenario-dependent objectives and constraints Same problem size and complexity compared to its deterministic counterpart Increased problem size and complexity comparted to its deterministic counterpart 3.3. Network flexibility CLSCND is a two-level decision-making that deals with: 1) the location of facilities at the first level; 2) the operations of the network at the second level. After the facility locations and transportation links have been established, the network flexibility is considered separately at the operational stage, i.e., manufacturing and distribution planning (Kopanos et al., 2012;Fiorotto et al., 2018). This has led to extensive modeling efforts of capacitated location/location-allocation problems for CLSCND (Talaei et al., 2016;Subulan et al., 2015a;Tosarkani and Amin, 2019). However, formulating rigid capacity constraints and demand satisfaction requirements for CLSCND suffers from two main shortcomings when the models are tested in a complex stochastic environment.
Loss of stability: the objective value may vary drastically over different scenario trees due to the changing number of facilities opened in order to fulfill all the customer demands in any conditions (Yu and Solvang, 2017). From the modeling perspective, the poor stability and consistency of a stochastic program significantly limits the use of the decision-support models (King and Wallace, 2012). Improper first-level decisions: satisfying all the customer demands in any conditions requires the setup of a large enough network capacity, which consequently leads to a redundant network structure and a waste of resources in low-demand scenarios (Yu and Solvang, 2017). For instance, with rigid capacity and demand fulfillment constraints, a new facility has to be opened to deal with a tiny increase on customer demands even if the probability of occurrence could be extremely low. However, maintaining an expensive facility for a small and highly fluctuate demand is clearly an improper decision. In practice, instead of building and operating a new facility, supply chain managers prefer to setup temporary or flexible capacity to deal with such demand.
To overcome the aforementioned shortcomings, the network flexibility is formulated in the proposed mathematical model, which enables a higher flexibility of network capacity and improves the robustness of the first-level location decisions in CLSCND. Under a variety of practical conditions, the meaning of network flexibility can be interpreted in several ways, e.g., outsourcing, hire of temporary workers or rent of equipment in demanding seasons, or even the lost sales due to capacity limits and inventory shortage.

Notations
The sets, parameters and variables used in the mathematical model are first given as follows: Set of manufacturing plants, indexed by p D Set of candidate locations for distribution center, indexed by d F Set of customers, indexed by f C Set of candidate locations for collection center, indexed by c U Set of candidate locations for repair/refurbishing center, indexed by u R Set of candidate locations for remanufacturing/recycling center, indexed by r X Set of products, indexed by x S Set of scenarios, indexed by s Fxr Fixed opening and operating cost of distribution center d2D, collection center c2C, repair/refurbishing center u2U and remanufacturing/recycling center r2R g vbpx, g vb dx , g vbcx, g vbux, g vbrx, g vbwx Variable processing cost for one unit product x2X at production plant p2P, distribution center d2D, collection center c2C, repair/ refurbishing center u2U, remanufacturing/recycling center r2R and waste disposal center w2W g Purpx Purchasing cost of raw materials for producing one unit product x2X at production plant p2P g Tc pdx , g Tc dfx , g Tc fcx , g Tccux, g Tccrx, g Tccwx, g Tc udx , g Tcrpx Transportation cost of one unit product x2X on link pd, p2P;d2D, link df, d2D;f 2F, link fc, f 2F;c2C, link cu, c2C;u2U, link cr, c2C; r2R, link cw, c2C; w2W, link ud, u2U; d2D, link rp, r2R; p2P g Oc fx , g Orc fx Cost of product x2X for customer f 2F under flexible capacity in both forward and reverse channels Carbon emission for processing one unit product x2X at production plant p2P, distribution center d2D, collection center c2C, repair/refurbishing center u2U, remanufacturing/recycling center r2R and waste disposal center w2W g Carbon emission for transporting one unit product x2X on link pd, p2P; d2D, link df, d2D; f 2F, link fc, f 2F; c2C, link cu, c2C; u2U, link cr, c2C; r2R, link cw, c2C; w2W, link ud, u2U; d2D, link rp, r2R; p2P Product-to-component conversion rate (disassembly) of product x2X g OLx, g ROLx Total flexibility limit of product x2X in both forward and reverse channels Objective function (1) minimizes the overall cost of the CLSC. The first four components are facility opening and operating cost of distribution center, collection center, repair/refurbishing center, and remanufacturing/recycling center, which include both fixed cost and scenario-dependent variable cost. The fifth and sixth parts are variable operating cost of existing production plants and waste disposal centers, and the seventh part calculates the purchasing cost of raw materials for production. The eighth component is the scenario-dependent transportation cost on each link between two facilities. The last two components are the cost of demand fulfillment by flexible capacity, e.g., outsourcing cost, etc., in both forward and reverse channels. It is noted that the signs "y" and "~ " are used in fuzzy objectives and constraints in order to represent the imprecision of the objective value and inequality, as argued by Darbari et al. (2019).
The second objective function (2) minimizes the total carbon footprint from the operation of the CLSC. The first six parts inside the parenthesis calculate the carbon emission of the facility operations, and the seventh and eighth parts calculate the carbon emission related to the demand fulfillment by flexible capacity in both forward and reverse logistics. Finally, the last components calculate the carbon emission associated with the transportation within the CLSC.

Constraints
The first two sets are capacity constraints. Constraints (3) and (4) are capacity requirements for the existing production plants and waste disposal centers, which guarantee the amount of products or EOL products processed at respective facilities cannot exceed their capacities at any scenarios.
Constraints (5)e(8) are capacity requirements for distribution centers, collection centers, repair/refurbishing centers, and remanufacturing/recycling centers, respectively. Moreover, when the binary variable equals to 0, this group of constraints also restricts that an unselected candidate location cannot be used as a node in the established network structure of the CLSC.
Constraints (9) and (10) set the upper limits on the overall customer demands fulfilled by the flexible capacity in both forward and reverse channels.
Constraints (11) are the demand satisfaction in the forward channel of the CLSC, which ensure the customer demand should be satisfied in all scenarios by both production and flexible options.
Eqs. (12)-(15) are flow balance requirements at production plants and distribution centers, which specify the relationship between the input and the output at respective facilities in each scenario. It is noteworthy that, in the CLSC, the remanufactured/recycled parts and components are reused in the production/assembly and the repaired/ refurbished products are used to fulfill the customer demand in order to maximize the value recovery of EOL products.
Eqs. (16) estimate the generation of EOL products at each customer location. The amount of EOL products is converted from a portion of the customer demand and is featured with high randomness. Constraints (17) (27) are flow balance constraints in the reverse logistics, which specify the relationships between the input and the output at collection centers, repair/refurbishing centers, remanufacturing/ recycling centers and waste disposal centers, respectively. The repair/ refurbishing fraction and the remanufacturing/recycling fraction are affected by the quality of EOL products. Herein, parameters w xu and w xr are the conversion rate of each product at its best condition. The actual fractions for repair/refurbishing and remanufacturing/recycling are compensated by the level of quality (Qal s x 1), which is a scenario-dependent parameter with a high stochastic nature.

Solution approach
The proposed mathematical model is a complex multi-objective optimization problem with both fuzzy and stochastic parameters. In order to solve the problem, the first step is to convert the fuzzy model into an equivalent auxiliary crisp model. Then, two approaches: sample average approximation (SAA) and weighting method (WM), are combined to manage the stochastic parameters and the multiple objective functions.

Conversation of an equivalent auxiliary crisp model
The fuzzy uncertainty or imprecision related to the parameters can be described by a membership function given in Eq. (31), where c l , c m and c h represent the pessimistic (lower bound) estimation, the most likely estimation and the optimistic (upper bound) estimation of the fuzzy number, respectively. It is noteworthy that the conversion method is established based on the general ranking approach developed by Jim enez (1996) and can be applied to defuzzilize both fuzzy triangular numbers and fuzzy trapezoidal numbers, as shown in Fig. 2. Eq. (31) gives the membership function of a fuzzy triangular number, and it can be converted into the membership function of a fuzzy trapezoidal number through replacing x ¼ c m by an interval x2½c ml ; c mh when m c ðxÞ ¼ 1.
The fundamental concepts of the conversion method are based on expected interval (EI) and expected value (EV), which have been extensively discussed by Yager (1981), Dubois and Prade (1987) and Heilpern (1992). Eqs. (32) and (33) illustrate how of a fuzzy triangular number can be converted to respective EI and EV, and the conversion method of a fuzzy trapezoidal number is given in Eqs. EIðcÞ¼ When converting a fuzzy model into an equivalent auxiliary crisp model, two questions need to be taken into consideration: the optimality of an objective function with fuzzy parameters and the feasibility of the crisp solution under fuzzy constraints (Jim enez et al., 2007). The optimality of the objective function with fuzzy coefficient can be estimated by EV. In order to solve the feasibility issue, the ranking approach (Jim enez,1996) and the aefeasible coefficient (Jim enez et al., 2007) are introduced to convert the fuzzy constraints into the crisp ones with a certain level of feasibility. The basic idea of the method is to evaluate to what extent a fuzzy number is "!" or " ¼ " or " " another fuzzy number.
For instance, given a pair of fuzzy numbersã andb,ã ! ab means the inequalityã is greater than or equal tob is fulfilled in a degree of a. Based on the aforementioned methods, the proposed fuzzy model can be converted into an equivalent auxiliary crisp model by using Eq. (36), whereas ½E ai 1 ; E ai 2 and ½E bi 1 ; E bi 2 are the expected interval of the fuzzy coefficient vectors a i and b i . minimizez ¼ EVðcÞx Eqs. (37) and (38) can be used to convert the fuzzy constraints with "!" and " ¼ " relationships into respective crisp constraints.

Sample average approximation
In the proposed mathematical model, the uncertainty featured with randomness is formulated by scenario-based stochastic parameters, which may involve a large number of scenarios and thus may lead to significant computational challenges. Due to this reason, an exact method may not be able to determine the optimal solution within reasonable time, so, in this paper, the sample average approximation (SAA) is used to solve the stochastic optimization problem. The SAA is a Monte-Carlo simulation based sampling approach and is capable to test the quality of solution of a stochastic program over a variety of scenario trees. SAA aims at approximating the optimal solution of the original large stochastic optimization problem through repeatedly solve several smaller samples in order to have a high level of confidence (Kleywegt et al., 2002).
Considering the single cost minimization objective model for CLSCND obtained in previous section, a general stochastic form can be written in Eq. (49), whereas x are the first-level location decision vectors and y are second level allocation decision vectors with respect to the feasible solution space Q . The expected value of the resource function E P ½Fðx; xðyÞÞ is dependent on the first-level decisions, the parameters with the realization of each scenario, as well as the respective probability. With the SAA, instead of solving gðx; yÞ, the expected value of the original problem is approximated byg K ðx; yÞ with the sample size K, as shown in Eq. (50). min x; y2Q According to Kleywegt et al. (2002) and Verweij et al. (2003), the SAA method can be implemented with the following steps. First, given a sample size K, J individual problems are randomly generated from a given probability distribution P . Then, a reference sample with size H ( H[K) is generated based on the same P , which is considered as the original problem g. Repeatedly solving the SAA problems formulated in Eq. (50) for J times and obtaining the non-scenario-dependent first-level decisions x 1 , …, x J and the values of the objective functions g 1 K , …, g J K . Based on which the statistical lower bound estimators of the original problem can be calculated using the average g K;J and variance s 2 g K;J (Mak et al., 1999), as shown in Eqs. (51) and (52).
Next, the upper bound estimators are calculated by Eqs. (53) and (54), where a feasible first-level decision b x (e.g. one obtained from the previous step) is used in the reference sample H. Even if the sample size H is much larger than K, the current optimization problem only involves the second-level continuous variables and becomes thus a linear program which is computationally manageable.
Finally, the optimality gap gap K;J;H ðb xÞ and variance s 2 gap are calculated by Eqs. (55) and (56) in order to assess the quality of the optimal solution obtained from the SAA. If the estimators fulfill the pre-determined convergence criteria, the result of the SAA g K can be used to approximate the optimal result of the original stochastic optimization problem g. Otherwise, increased sample size and repetition need to be tested until all the pre-determined convergence criteria is satisfied. gap K;J;H ðb xÞ ¼ b g H ðb xÞ À g K;J (55) In addition, one of the most important purposes of using a stochastic optimization model is to make a robust first-level decision. With the SAA method, the optimal first-level decision is determined by the comparison of the performance among the candidates x 1 , …, x J in the reference sample H.

Weighting method
The proposed mathematical model is a multi-objective optimization problem that aims at determining the optimal trade-off between cost and carbon emission in CLSCND. Eq. (57) gives a general form of a multi-objective optimization problem that includes n objective functions. The optimal solution of such a problem is called a Pareto optimum at which the performance of one objective z i¼l ðxÞ; l2½1; n cannot be improved without a compromise on the others z isl ðxÞ; l2½1; n (Sakawa et al., 2013).
Min zðxÞ ¼ ðz 1 ðxÞ; z 2 ðxÞ; …; z n ðxÞÞ T S:t: x2X (57) Weighting method (WM) is a scalarization method and is put forward by Zadeh (1963). It has been used to solve multi-objective optimization problems in a large variety of industries and businesses (Sheu, 2007;Yu and Solvang, 2016c). As an a priori method, the weight of each objective function should be given in advance, which determines the relative importance in decision-making. The weighted objective functions can then be combined in a weighted sum (WS), as shown in Eq. (58). For a minimization problem, the Pareto optimum is the one that has the smallest WS.

Min wzðxÞ
In the proposed mathematical model, the cost and the carbon emission are measured by different units. Hence, as shown in Eq. (59), each objective function is first normalized by its optimal value that is obtained from solving the individual single objective optimization problem. Similar to the principle of goal programming, the normalized objective ziðxÞ z min i measures the deviation of the objective value from its individual optimal performance and finds out the one with the least deviations in the WS. The advantage of the WM is that it is very simple to implement and the computational complexity is at the same level of the individual single objective optimization problem.

Algorithmic procedures
The proposed fuzzy-stochastic multi-objective model is first converted into an equivalent auxiliary crisp model and is then solved with a hybrid SAA based WM (SAAWM) approach. The algorithmic procedures are given as follows.

Computational experiments
Computational experiments are conducted in order to validate the proposed fuzzy-stochastic multi-objective model and the solution approach.

Parameter generation
In the experiment, a CLSCDN problem consists of two production plants, seven potential locations for distribution center, ten customer zones, five potential locations for collection center, three potential locations for repair/refurbishing plant, four potential locations for remanufacturing/recycling plant as well as one waste disposal center was considered. All the test parameters were generated randomly in accordance with the method used by Pishvaee and Torabi (2010).
In order to generate fuzzy triangular numbers, the most likely value of each fuzzy parameter (c m ) was first generated from a given interval. Table 4 shows the generation interval of the fuzzy parameters of production plants and distribution centers. It is assumed the carbon emission per unit product is inversely related to the variable processing cost due to the increased investment on technological updates for the reduction of carbon emission (Wang et al., 2011). Afterwards, both the pessimistic value (c l ) and the optimistic value (c h ) were calculated from Eq. (60) in order to maintain the generality (Pishvaee and Torabi, 2010). Herein, b is a random number generated from a uniform distribution between 0.1 and 0.3. Tables 5 and 6 present the fuzzy parameters of production plants and distribution centers used in the experiment.
In order to convert the fuzzy constraints into their respective crisp forms, a ¼ 0.7 was used in the experiment. The generation of all the fuzzy parameters as well as the calculation of their crisp numbers are given in Appendix A (supplementary data file). In addition, the uncertainty related to the customer demands, the rate of EOL returns and the quality of EOL products are formulated as stochastic parameters. Table 7 shows the interval from which the scenarios trees of the test samples are generated randomly.
The mathematical model was coded in LINGO 18.0 solver and all the computations were performed on a PC with i5-6400T 2.20 GHz CPU and 8 GB RAM.

Results and discussions
In the experiment, three different sample sizes: 10, 30 and 50 with 10 repetitions each were tested. In the initial step, two single objective problems for cost minimization and minimization of carbon emission were solved individually. Table 8 presents the statistical lower bounds, upper bounds and gap estimators. First, the in-sample stability of the model with respective to different scenario trees were evaluated for both objective functions. Insample stability measures the deviation of the objective value with respect to different scenario trees generated from the same data distribution. According to King and Wallace (2012), it is a very important characteristic to test the internal consistency of a stochastic model and the robustness towards the discretization procedure. Without the in-sample stability, the use of a stochastic model is significantly limited and unreliable. In this paper, the coefficient of variation (CV) is used to evaluate the in-sample stability, which is calculated by CV ¼ s = m . As can be seen in Table 7, with the increase of the sample size from 10 to 50, the CV of cost is reduced by 63.9% and the CV of carbon emission is reduced by 81.5%. Both of which show, based on the same model structure and data distribution for CLSCND, the in-sample stability may be improved with the increase on sample size.
Furthermore, the CVs of cost and carbon emission related to the facility operation, the transportation and the use of flexible capacity with respect to different sample sizes were compared. As shown in Fig. 3, the CVs of all elements decrease gradually with the increase of sample size. However, compared with the CV of the use of flexible capacity related to both cost and carbon emission is much larger than that of the facility operation and the transportation. As discussed in previous literature (King and Wallace, 2012;Yu and Solvang, 2017), a stochastic network design model may lead to bad stability due to the fluctuation on the objective values under Algorithmic Procedures Step 1: Convert the fuzzy objective functions and constraints into an equivalent auxiliary crisp model by Eqs. (32)-(38) Step 2: Generate J individual problems with sample size K and a reference problem with sample size H from a given probability distribution P Step 3: For each problem within J 3.1 Generate two single objective optimization problem for minimizing the cost and carbon emission individually 3.2 Solve the two single objective optimization problem and obtain respective g j K and x j Step 4: Calculate the lower bound estimators g K;J and s 2 g K;J by Eqs. (51) and (52) Step 5: Calculate the upper bound estimators b g H ðb xÞ and s 2 b g H ðb xÞ by Eqs. (53) and (54) Step 6: Evaluate the optimality gap gap K;J;H ðb xÞ and variance s 2 gap by Eqs. (55) and (56). If the convergence criteria is fulfilled, proceed to next step. Otherwise, return to Step 2 and increase K and/or J Step 7: Calculate the individual optimal value of each objective function z min i by solving the reference problem with sample size H using the candidates obtained by the SAA Step 8: Determine the weight of each objective function Step 9: For each problem within J Solve the normalized WS in Eq. (59) and obtain respective g j K and x j and the objective value of the cost and carbon emission Step 10: Calculate the lower bound estimators g K;J and s 2 g K;J for the normalized WS, cost and carbon emission, respectively, by Eqs. (51) and (52) Step 11: Calculate the upper bound estimators b g H ðb xÞ and s 2 b g H ðb xÞ for the normalized WS, cost and carbon emission, respectively, by Eqs. (53) and (54) Step 12: Evaluate the optimality gap gap K;J;H ðb xÞ and variance s 2 gap of the normalized WS, cost and carbon emission, respectively, by Eqs. (55) and (56). If the convergence criteria is fulfilled, proceed to next step. Otherwise, return to Step 2 and increase K and/or J Step 13: Calculate the optimal objective value and determine the first level decisions of the multi-objective problem by solving the reference problem with sample size H using the candidates obtained by the SAA Table 4 Generation interval of fuzzy parameters of production plant and distribution center.    In this respect, the introduction of network flexibility can solve that problem. This is a very important observation that, even if the CV of the use of flexible capacity is instable, it improves the flexibility of the model by relaxing the rigid capacity constraints and is an effective way to improve the in-sample stability of the objective value over different scenario trees. In practice, it means, instead of opening and operating an additional facility, the small and fluctuate increase on customer demands may be satisfied by a variety of flexible options in a CLSC. Thereafter, the gap estimators between the lower bound and the upper bound were tested. They measure the quality of the different SAA problems in comparison with the reference problem. With the increase of sample size from 10 to 50, the optimality gap of cost minimization decreases dramatically from 1.73% to 0.16%. However, the optimality gap of the minimization of carbon emission is not improved, which remains stable at approximately 5%. Considering the quality of the solution in both in-sample stability and optimality gap estimators, the sample size 50 was selected to solve the multi-objective optimization problem in our experiment.

Fuzzy parameters Interval
The individual optimal cost and carbon emission of the reference sample, which are used as the benchmark in the weighting method, were first calculated through comparing the candidate first-level decisions obtained from the SAA. Then, the weight combination of w cost ¼ w emission ¼ 0:5 was tested. The result is presented in Tables 9 and 10. Table 9 compares the statistical lower bounds, upper bounds and gap estimators of the WS, the cost and the carbon emission in the optimal solution. It is noteworthy that not only the WS but also each objective value should comply with the convergence criteria in a multi-objective optimization problem. Table 10 shows the comparison of the cost, the carbon emissions and the first-level decisions in different scenarios. Compared with the two benchmark values, in the optimal solution of w cost ¼ w emission ¼ 0:5, the total cost deviates by 3.4% and the carbon emission deviates by 5.3% from their individual optimums. The result shows the trade-off between cost and carbon emission in CLSCND and, in practice, it can be interpreted as the total carbon emission of the CLSC may be reduced by 8% when 3.4% more investments are budgeted for improving the technological levels as well as other carbon reduction measures. In addition, it is also noted that the network configurations are by no means identical when the focus of system planning varies.

Sensitivity analysis
The model's behavior with respect to the change of weight combinations are of interest in the experiment. Hence, a sensitivity analysis was conducted. The weight of the cost objective increases gradually from 0.1 to 0.9 by the step of 0.1, while at the same time, the weight of the carbon emission objective decreases with the    Fig. 4, was obtained by solving the reference sample in each scenario with the candidates of the first-level decisions obtained from the SAA. Fig. 5 illustrates, in the optimal solution of each scenario, the deviation of cost and carbon emission from their respective individual optimums. In the individual carbon emission optimization scenario, compared with its adjacent point (w cost ¼ 0:1 and w emission ¼ 0:9), the carbon emission can only be reduced by 0.46% with a 68.2% increase on the total cost due to the large cost related to facility operation. Hence, this point is considered as an unrealistic and near dominated solution, which is not included in the comparison in Figs. 4 and 5.
It is observed the slope of the Pareto frontier is much higher when w cost decreases from 1 to 0.6, and this reveals the increased investments> on technological updates as well as other measures are much more effective in the minimization of carbon emission. Besides, when w cost ranges from 0.3 to 0.8, the cost and carbon emission are balanced with each other as can be seen in Fig. 5, so the system behavior of the scenarios within this range may particularly be interested for investigation. Based on the above discussions, for the test instance, the weight combinations with w cost ¼ 0:8, w cost ¼ 0:7 and w cost ¼ 0:6 may be suggested to the decision-makers.

Implications
Through the analysis of the computational results, some general modeling and managerial implications are obtained.
The implications related to the modeling and computations are summarized as follows: From the modelling perspective, the fuzzy parameters enable a higher flexibility in the interpretation of the objective value obtained, while on the other hand, the stochastic parameters in the constraints improve the robustness of the first-level location decisions in CLSCND. The inclusion of network flexibility (penalty) in the objective function can dramatically improve the in-sample stability and the consistency of a stochastic network design problem over a variety of scenario trees. Hence, the quality and the reliability of decision-making is improved. In general, the quality of solution of a SAA problem is improved with the increase of sample size. However, in some conditions, the increased sample size may not yield a significant improvement on the convergence rate.
The managerial implications are given as follows: From the practical perspective, the network flexibility at the operational level, e.g., outsourcing, hire of seasonal workers and rent of equipment, etc., may have significant impact on the location decisions and should thus be taken holistically into account of CLSCND. The carbon emission of a CLSC may be reduced at an increased cost for technological updates as well as other methods. The effectiveness in the reduction of carbon emission may be of significant difference when different investment strategies are used.

Conclusion
The increasing focus on sustainable development accompanied with stringent environmental regulations has caused companies to take responsibility of the whole lifecycle of their products through CLSC management. This paper formulates a new fuzzy-stochastic multi-objective optimization model for CLSCND under both epistemic uncertainty and randomness. The customer demands, the rate of EOL returns and the quality level are featured with high randomness and are thus formulated as stochastic parameters. The other uncertain parameters are formulated as fuzzy triangular numbers in order to deal with the inherent imprecision of forecast. Besides, the network flexibility is considered and is incorporated in the model development. In order to solve the proposed mathematical model, an equivalent auxiliary crisp model is first formulated and a SAAWM is then developed. The model is validated through a set of computational experiments whose results lead to  From the theoretical perspective on model development, this research first provides a thorough elaboration on the significance and reasons why the mixed uncertainty and the network flexibility need to be considered holistically in CLSCND. Based upon which, a new fuzzy-stochastic multi-objective model for CLSCND is the developed to simultaneously deal with multiple objectives, different types of uncertainty and the network flexibility. From the practical implementation perspective, a new solution approach is proposed to effectively and efficiently solve the complex optimization problem formulated and to test the performance of the model over different scenario trees. In addition, general managerial implications are obtained through the analytical study of a set of numerical experiments in order to provide practical insights into CLSCND under multiple objectives, mixed uncertainty and the network flexibility.
For further development of this research, four suggestions are made.
Research efforts may be given to the development of methods to assess, formulate and include social sustainability in CLSCND. The carbon emission objective may be reformulated as different carbon policies in order to better reflect the practical regulations (Mohammed et al., 2019). The model may be expended with more flexible options, for example, a flexible manufacturing system for multi-sourced product flow in CLSCND (Yu and Solvang, 2018). The solution approach may be improved with more advanced methods in order to improve the computational efficiency and the quality of solution, e.g., augmented εeconstraint method for multi-objective optimization (Mavrotas, 2009;Yu and Solvang, 2016a).

Author contributions
Hao Yu:Conceptualization; Literature Review; Modeling; Algorithm; Experiment and Analysis; Writing of the paper. Wei Deng Solvang: Conceptualization; Supervision.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.