A Robust Fuzzy Optimization Model for Closed-Loop Supply Chain Networks Considering Sustainability

Supply chain network design (SCND) is an important strategic decision determining the structure of each entity in the supply chain, which has an important impact on the long-term development of a company. An efficient and effective supply chain network is of vital importance for improving customer satisfaction, optimizing the allocation of resources, and increasing profitability. The environmental concerns and social responsibility awareness of the whole society have spurred researchers and managers to design sustainable supply chains (SSCs) integrating the economic, environmental, and social factors. In addition, the innate uncertainty of the SCND problem requires an integrated method to cope. In this regard, this study develops a multi-echelon multi-objective robust fuzzy closed-loop supply chain network (CLSCN) design model under uncertainty including all three dimensions of sustainability. This model considers the total cost minimization, carbon caps, and social impact maximization concurrently to realize supply chain sustainability, and is able to make a balance between the conflicting multiple objectives. Meanwhile, the uncertainty of the parameters is divided into two categories and addressed with two approaches: the first category is missed working days related to social impact, which is solved by the fuzzy membership theory; the second category is the demand and remanufacturing rate, which is settled by a robust optimization method. To validate the ability and applicability of the model and solution approach, a numerical example is conducted and solved using ILOG CPLEX. The result shows that the supply chain network structure and the value of the optimization objectives will change when considering sustainability and different degrees of uncertainty. This will enable supply chain managers to reduce the environmental impact and enhance the social benefits of their supply chain activities, and design a more stable supply chain to better cope with the influence of uncertainty.


Introduction
In the last few decades, with the growing emphasis on resource scarcity and environmental pollution, a closed-loop supply chain (CLSC) that integrates the traditional forward supply chain and reverse logistics has become an important research field in both academia and industry [1][2][3]. In 2002, the Waste Electrical and Electronic Equipment (WEEE) Directive [4,5], which is the first developed by the European Union, obliges original equipment manufacturers (OEMs) to collect scrapped products for subsequent recycling or reuse. Subsequently, similar policies are published in other countries such as the United States and China [6,7]. The incentive from customers also • A mixed integer programming (MIP) model is developed to design a CLSCN in which all three aspects of sustainability (i.e., economic, environmental, and social) are optimized concurrently. In our paper, the economic dimension is measured by calculating the cost of the SCN, while the environmental aspect is evaluated through the limitation of carbon caps. Finally, the social dimension is modeled by minimizing the missed working days due to work damages; • The uncertainty parameters of the proposed model are divided into two categories. The data of missed working days are hard to obtain; therefore, fuzzy membership theory is employed to measure it. Other types of uncertain parameters in our model are the customer demand and the remanufacturing rate, which will obviously influence the tactical plan of the supply chain.
In this context, robust optimization is employed to ensure the stability of the network to some extent [29], and • Lastly, this paper provides a detailed analysis of the computational results to explain the impacts of uncertainty on the Pareto optimal solutions and the variation of the supply chain network structure.
The remainder of our paper is organized as follows. Relevant previous researches are reviewed in Section 2. Section 3 introduces the background of robust optimization and multi-objective methodology. Problem description, assumptions, and the formulation of fuzzy membership degree are described in Section 4. Section 5 presents the developed multi-objective sustainable CLSC optimization model, and then converts it into a crisp formulation. Section 6 provides a numerical example, and the computational results are discussed. Finally, the main conclusions and the directions for further research are presented in Section 7.

Literature Review
In our study, a robust fuzzy mixed integer programming (MIP) mathematical model that considers sustainability is developed to design the CLSC. According to this research content, the literature review section is divided into two parts: (1) CLSC models and (2) a supply chain network design (SCND) problem considering sustainability.

CLSC Models
The known forward supply chains include the activities of raw materials purchasing, commodity production, transportation, inventory, and sale [30]. Meanwhile, reverse logistics integrates all the activities of recycling used products involving the collecting, inspecting, sorting, remanufacturing, and disposal of returned products [31]. To avoid sub-optimality, more and more studies integrate the activities of the forward supply chain and reverse logistics to formulate the CLSC [32]. Savaskan et al. [33] presented a CLSC model for collecting and remanufacturing used products; their model considered different options for collecting used products and could obtain the applicable CLSC structures for OEMs. Chouinard et al. [34] proposed a stochastic programming model to optimize the forward supply chain network together with reverse logistics; they utilized supply loops to describe CLSCs and considered uncertain parameters related to recovery, processing, and demand. To solve this problem, they designed several scenarios according to the sample average approximation (SAA) method for the volume of uncertain parameters. Zhang et al. [35] formulated an MIP model for a capacitated production plan problem under CLSC with remanufacturing. In this paper, they considered two kinds of storage levels with different costs and addressed a genetic algorithm-based heuristic approach to solve the large-scale remanufacturing problem. Turki et al. [36] studied the problem of CLSCN design for end-of-life (EOL) products; they developed a system to integrate the manufacturing, transport, inventory, and remanufacturing activities of supply chains simultaneously. The customer demand is supposed to be lost when it cannot be satisfied, and the quantity of returned products will be different when the sale of the past period is a variation in their work. An MIP model was proposed by Fattahi and Govindan [37] to plan and optimize a forward/reverse logistics network. The two-stage stochastic programming method was used to cope with stochastic demand and used products return. Maedeh et al. [38] presented an MIP model to solve the CLSCN design problem of the tire industry. They considered three different objectives: supply chain profit maximization, customer satisfaction maximization, and minimization of the distance between collecting facilities and customers. They defined different scenarios to measure the uncertainty and converted it with a robust optimization method. Ruiz-Torres et al. [39] considered both return volume and the supplier's ability uncertainty when they determined the supplier capacity and returner incentives of CLSCs; they also accounted for the satisfaction of products and loss costs.

SCND Problem Considering Sustainability
In the last decade, increasing attention to environmental and social issues has led to a significant increase to consider sustainability when designing SCNs. In most of the existing studies, only one or two dimensions of sustainability have been optimized in the model presented [40].
Some researchers have considered environmental factors in the field of SCND. The work of Nagurney et al. [41] took the environmental aspects of sustainability into account to model an SCND problem and reformulated it as an elastic demand transportation network equilibrium problem. Chaabane et al. [42] proposed a mixed-integer linear programming (MILP)-based framework in the design of sustainable supply chains considering life cycle assessment (LCA), where the tradeoffs between economic and environmental objectives under various cost and operating strategies are evaluated. Research by Mari et al. [43] emphasized the importance of resilience and the disruption risk of sustainable SCNs. They presented a multi-objective optimization model involving carbon emission and carbon footprints. To handle the model, a goal programming method was developed. Taskhiri et al. [44] developed an MILP model to quantify the tradeoff between costs and the carbon footprint of a wood industry logistics network. They utilized LCA software to investigate the environmental effects of the network. Maditati et al. [45] presented a structured review of GSCM literature and conducted the bibliographic analysis tools and techniques to identify the key journals, influential institutions, and impactful and trending articles; the analysis reveals six underlying research streams related to GSCM. Rohmer et al. [46] addressed the global food system sustainability problem with the objectives of supply chain cost and environmental impact minimization to maintain adequate dietary intake. As mentioned ahead, another way to reflect environmental aspects in supply chain strategy decisions is to recycle or remanufacturing used products through the CLSC network structure. This topic has been studied by numbers of researchers in recent decades (see the papers of Min et al. [47], Wang et al. [48], Özceylan and Paksoy [49], Jindal and Sangwan [50], Ma and Li [51], and Zhen et al. [52]).
The social performance of SSCM has received less concern than economic and environmental dimensions, as it is very complex and has an extensive scope [53,54]. The standard of ISO 26000 [55] categorizes social dimension into seven major groups, which has been widely accepted by the public. Based on this, some researchers tried to evaluate the social dimension of supply chains. Dehghanian and Mansour [56] presented a multiple objective mathematical model for a sustainable recovery network of scrap tires in which the social impact is optimized together with economic and environmental dimensions. In this paper, several social factors were taken into account. Pishvaee et al. [57] explored a mathematical model for a medical needle and syringe supply chain with three pillars of sustainability. They utilized the LCA method to measure the environmental and social effects, and four different aspects were considered to evaluate the social dimension. Govindan et al. [58] developed a generic CLSCN structure to rationalize the effect of product recovery for improving manufacturing sustainability. The social benefits are economic welfare and growth, stakeholder responsibility, extended producer responsibilities, and employment practices. Cao et al. [59] presented a multi-objective 0-1 integer programming model for the disaster supply chain. They considered all three dimensions of sustainability in their forward network, and designed a branch and bound approach to solve it. Taleizadeh et al. [60] introduced an MIP optimization model for designing a Sustainability 2019, 11, 5726 5 of 24 sustainable CLSC. They measure the social and environmental effects of supply chain decisions by Global Reporting Initiative guideline indicators.
The above reviewed literature reveals that most of the previous researches consider the economic and environmental objectives of the sustainable CLSCN design problem. While some of the literature begin to identify the social dimensions, the uncertainties of different kinds of parameters are often ignored. To fill the gap, we develop a robust fuzzy MIP model for the CLSCN problem considering sustainability and uncertainty. The findings and managerial implications of this article may be beneficial for designing and optimizing sustainable CLSCNs.

Robust Optimization
The general framework of robust optimization summarized by Li and Floudas [61,62] is introduced in this part. Normally, an optimization problem under uncertainty can be defined as follows: where x j is a continuous or an integer variable. The objective and right-hand side (RHS) uncertainty can be transformed into left-hand side (LHS) uncertainty as follows: So, consider only the LHS constraint under uncertainty: and a ij are subject to uncertainty, which can be defined below: where a ij are the nominal value of parameters,â ij are the positive constant perturbations, ξ ij are the independent random variables, and J i are the index subset that conclude the uncertainty variables. Then the formulation of constraint (3) can be reformulated as follows: In the set induced robust optimization method, after determining the uncertainty set U of ξ, constraint (5) will be as follows: The establishment of the robust optimization model depends on uncertain sets, among which the Box Uncertainty Set, Ellipsoidal Uncertainty Set, Polyhedral Uncertainty Set, and combinations of them have been proved to be feasible. In this paper, the commonly used "Box + Polyhedral" Uncertainty Set with strict equality is selected to deal with uncertain sets. For expression simplicity, the constraint index i is eliminated from in the random vector ξ, so the uncertainty set U is: where Γ and Ψ controls the size of the uncertainty sets, and can be represented as: When the formulation of uncertainty set U is as shown in constraint (7), the robust counterpart, constraint (6), will be converted into: where z i and w ij are the dual auxiliary variables. Note that when Ψ = 1, the set U is the "interval + polyhedral" uncertainty set, and the robust counterpart can be simplified into the formulation proposed by Bertsimas and Sim [63]: This approach has proved that if the numbers of uncertain parameters that violate their nominal values are less than Γ i , the robust solution will remain feasible; if the violation number is more than Γ i , the robust solution will remain feasible with the probability of: where: x * j is the optimal solution of model (27), and ϕ(θ) is the cumulative distribution function of a standard normal random variable.

Multi-Objective Technique
As mentioned before, this paper considers two conflicting optimization objectives, which require a multi-objective technique to get the tradeoff between the objectives. In this way, decision makers will obtain a set of Pareto optimal solutions. Several popular methods, such as the ε-constraint approach [64], the weight-sum method [65], and fuzzy multi-objective programming [66] are frequently applied to solve the multi-objective model. Mavrotas [67] has proved that the ε-constraint approach has several advantages over the weight-sum method, and the procedure of fuzzy multi-objective programming is more complex. Hence, in our paper, the ε-constraint method is employed to solve the proposed model. Consider the following multi-objective programming model [68]: are the p objective functions, and S is the feasible region. In the ε-constraint method, the main objective of the model is optimized, and other objectives are taken as constraints: The Pareto optimal solutions are obtained by the parametric variation of the RHS of the constraint objective function (ε i ).

Problem Description
The CLSCN studied in our study includes forward supply chain and reverse logistics (as illustrated in Figure 1). In the forward part, plants produce new products and then deliver them to distribution centers, and then the products are delivered to customers. After customers use them, they recycle the waste products and transport them to collection centers; after testing and evaluation, collection centers ship the recyclable products to remanufacturing centers for remanufacturing. The remanufactured products will be transported to the distribution center for unified distribution and sent to the new customer area for reuse. Meanwhile, transportation between facilities adopts two different kinds of transport tools: transport vehicle type A is a large truck, which is used for the transportation process from factory and remanufacturing center to distribution center and collection center to remanufacturing center; transport vehicle type B is a medium-sized truck, which is used for the transportation process from distribution centers to customers and from customers to collection centers.

Model Assumptions
To formulate the model, consider the following assumptions: • The model is a single period, single product, but multi-echelon.

•
The potential locations of each echelon are predetermined.

•
All returned products should be collected, and all the customer demands must be met.

•
The capacities of all facilities are limited to an interval if they are established.

•
The quality of the remanufacturing products is the same and sold to the same market.

•
Products that cannot be recycled and reused are buried nearby, and special landfill sites are not In this article, we consider the following three aspects to achieve sustainability in the SCN: the first is to minimize the total cost of the supply chain; the second is the carbon cap: consider the carbon emissions in the operation of logistics facilities and transportation process, and set the upper cap of carbon emissions. Third, the social impact: select the number of days of missed work as a representative factor for analysis.
In addition, the influences of the uncertainties considered in our model are the number of missed working days, demand, and remanufacturing rate. Then, the fuzzy membership degree theory is used to deal with the uncertainty of missed working days, and the robust optimization method is utilized to handle the uncertainties of the demand and remanufacturing rate.

Model Assumptions
To formulate the model, consider the following assumptions: • The model is a single period, single product, but multi-echelon.

•
The potential locations of each echelon are predetermined.

•
All returned products should be collected, and all the customer demands must be met.

•
The capacities of all facilities are limited to an interval if they are established.

•
The quality of the remanufacturing products is the same and sold to the same market.

•
Products that cannot be recycled and reused are buried nearby, and special landfill sites are not considered.

•
The price of recycled products from customers is the same.

•
A worker should only build one type of facility.

•
Transportation costs from one facility to another one are per product unit.

•
There are two types of transport vehicles that existed in the supply chain to execute different transportation tasks.

Representation of Social Impact Uncertainty
The social objective is to maximize the social sustainability of SCNs considering different social categories [69]. In our study, one of the most important social impact indicators (i.e., the number of missed working days) related to work's damage is considered. This indicator is closely related to the location of the facility, work safety, and human rights issues [70]. We assume that the maximum number of days is permitted. As for the influence of various factors such as the production environment, staff quality, and pre-job training, it is almost impossible for decision makers to obtain more accurate data records. Therefore, the concept of the fuzzy member function is utilized to minimize the total missed days of the CLSC model [71].
As shown in Figure 2, LD ideal is the ideal number of missed working days, LD max is the maximum allowable number of missed working days, and the social impact membership formula of the factory, distribution center, collection center, and remanufacturing center is as follows:   Figure 2. Number of lost work days per worker due to occupational accidents.

Mathematical Model
According to the above description and assumption of the problem, the mathematical model for a CLSCN considering sustainability will be presented in detail. These formulations are developed by extending previous works (mainly [72]).

Mathematical Model
According to the above description and assumption of the problem, the mathematical model for a CLSCN considering sustainability will be presented in detail. These formulations are developed by extending previous works (mainly [72]).

I
Set of candidate plants, i ∈ I. j Set of candidate distribution centers, j ∈ J. K Set of candidate collection center, k ∈ K. L Set of candidate remanufacturing center, l ∈ L. M Set of fixed customers, m ∈ M.

Parameters
Return rate of used products from customer m. θ Remanufacturing rate of used products at remanufacturing centers. σ Carbon cap of supply chain set by the government. cp i Cost of producing one unit product at plant i. cd j Cost of distributing one unit product at distribution center j. cc k Cost of collecting one unit product at collection center k. cr l Cost of remanufacturing one unit product at remanufacturing center l. α 1 Cost of transporting one unit of new product per mile by transport vehicle type A. β 1 Cost of transporting one unit of used product per mile by transport vehicle type A. α 2 Cost of transporting one unit of new product per mile by transport vehicle type B. β 2 Cost of transporting one unit of used product per mile by transport vehicle type B. dpd ij Distance between plant i and distribution center j. ddc jm Distance between distribution center j and customer m. dcc mk Distance between customer m and collection center k. dcr kl Distance between collection center k and remanufacturing center l. drd lj Distance between remanufacturing center l and distribution center j. Carbon emission for producing one unit of product at plant i. cod j Carbon emission for distributing one unit of product at distribution center j. coc k Carbon emission for collecting one unit of product at collection center k. cor l Carbon emission for remanufacturing one unit of product at remanufacturing center l. λ 1 Carbon emission for transporting one unit of new product per mile by transport vehicle type A. γ 1 Carbon emission for transporting one unit of used product per mile by transport vehicle type A. λ 2 Carbon emission for transporting one unit of new product per mile by transport vehicle type B. γ 2 Carbon emission for transporting one unit of used product per mile by transport vehicle type B. s1 i Numbers of missed working days due to work's damages during establishing plant i. s1 max i Maximum averages of missed days due to work's damages when establishing plant i. s2 j Numbers of missed working days due to work's damages when establishing distribution center j. s2 max j Maximum averages of missed days due to work's damages when establishing distribution center j. s3 k Numbers of missed working days due to work's damages when establishing collection center k. s3 max k Maximum averages of missed days due to work's damages when establishing collection center k. s4 l Numbers of missed working days due to work's damages when establishing remanufacturing center l. s4 max l Maximum averages of missed days due to work's damages when establishing remanufacturing center l.

Decision Variables
A binary variable, it is 1: if collection center k is established; otherwise, 0. R l A binary variable, it is 1: if remanufacturing center l is established; otherwise, 0. Qpd ij Volume of products shipped from plant i to distribution center j. Qdc jm Volume of products shipped from distribution center j to customer m. Qcc mk Volume of used products shipped from customer m and collection center k. Qcr kl Volume of used products shipped from collection center k and remanufacturing center l. Qrd lj Volume of products shipped from remanufacturing center l and distribution center j. µ s1 i Membership degree for the number of missed working days due to work's damages for each worker during establishing plants. Membership degree for the number of missed working days due to work's damages for each worker during establishing distribution centers. µ s3 k Membership degree for the number of missed working days due to work's damages for each worker during establishing collection centers. µ s4 l Membership degree for the number of missed working days due to work's damages for each worker during establishing remanufacturing centers.

Objective Functions
We use the notations described in the previous subsection to develop the objectives step-by-step as follows: The first objective is to minimize the total cost of the whole SCN. The total cost of the supply chain includes the fixed establishment cost, operation cost, and transportation cost. The fixed establishment cost is the total cost of establishing facilities; the operation cost is the processing cost of products in facilities; the transportation cost is the transportation cost of products flowing between facilities. Therefore, the total cost of supply chains is expressed as: Total supply chain cost = fixed establishment cost + operating cost + transportation cost.
The second objective is to minimize the number of missed working days due to work's damages on the CLSC. To realize it, we employ the fuzzy membership introduced in Section 5.3:

Model Constraints
Equations (20)- (39) are constraints. Equation (20) ensures that all customer demands should be satisfied. Equation (20) ensures the collection of used products.  Equations (25)- (28) are the lower capacity constraints. In practice, the facility will not establish if it only needs to handle a few products. Therefore, it is necessary to pre-set a lower bound to avoid the unreasonable solution.
Equation (33) is the carbon caps constraint. The carbon emission considered in our supply chain contains production emission in plants and the remanufacturing center, the handling of emissions in the distribution and collection centers, and the transportation emissions between facilities. i∈I j∈J cop i Qpd ij + l∈L j∈J cor l Qrd lj + j∈J m∈M cod j Qdc jm + k∈K l∈L coc k Qcr kl + i∈I j∈J l∈L m∈M (λ 1 dpd ij Qpd ij + λ 1 drd lj Qrd lj + λ 2 ddc jm Qdc jm ) + l∈L k∈K m∈M (γ 1 dcr kl Qcr kl + γ 2 dcc mk Qcc mk ) ≤ σ Equations (34)-(37) are the fuzzy membership degree constraints related to the missed working days due to work's damage. Equations (38) and (39) are binary and non-negative restrictions.

Robust Counterpart Mathematical Model
To ensure the stability and reliability of supply chain networks, we use the robust optimization method to cope with the uncertainty parameters of customer demand d m and remanufacturing rate θ.
First, according to Equation (4), the uncertain demand d m is redefined as d m with the following value range: where d m and d m represent the nominal demand and the maximum deviation from the nominal demand, respectively, ξ 1 m denotes the independent random variables affected by uncertainty. Based on the above assumptions, Constraints (20) and (21) are rewritten as follows: According to the robust counterpart optimization theory, we introduce parameters Γ 1 m and Ψ 1 m in [0, 1] to control the robustness of customer demands, in which the value ranges are [0, 1]. In this case, constraints (41) and (42) are reformulated in the following robust counterpart: where w 1 m and z 1 m are the dual auxiliary variables for developing the robust counterpart based on formulation (8).
Similarity, according to Equation (4), the uncertain remanufacturing rate θ is redefined as θ with the following value range: where θ and θ represent the nominal remanufacturing rate and the maximum deviation from the nominal remanufacturing rate, respectively, ξ 2 denotes the independent random variables affected by uncertainty. Afterwards, constraint (24) is rewritten as: Then, formulation (49) is converted as follows: where Γ 2 l and Ψ 2 l denote the respected budget of uncertainty.  (24) is converted into expressions (50)- (55).
In our study, the first objective is taken as the objective, and the second objective is dealt with as a constraint of the proposed model. After adopting the robust optimization method and the ε-constraint approach, the mathematics model proposed in our paper is transformed into problem P and represented as follows: minW 1 s.t. Π 2 ≥ ε 2 (56) and constraints (22)

Data Generation
In order to explore the validity and applicability of the proposed robust fuzzy mathematic model and its solution approach for designing and planning CLSCs, a numerical example has been presented in this section. The model will be coded and tested by CPLEX12.8 with Windows 10 (Intel(R) Core (TM) i7-7700HQ @2.80GHz, 8 GB of RAM). The numerical example contains five candidate points for establishing plants, six candidate points for constructing distribution centers, 10 fixed customer points, five potential locations for collection centers, and four remanufacturing centers. We assume two transport vehicle types, and the value of transport vehicle type B parameters is calculated based on type A, which is expressed as α 2 = 1.5α 1 , β 2 = 1.5β 1 , λ 2 = 1.5λ 1 and γ 2 = 1.5γ 1 . The model has 4665 constraints, 31 binary variables, and 8450 other variables. The values of the other parameters u are randomly generated with a uniform distribution, as shown in Table 1.

Multi-Objective Solution Mechanism
In order to solve the robust fuzzy model P converted in the previous section, we need to determine ε 2 (the lower bound of the second objective function). Hence, we calculate the payoff table (Table 2) by individually optimizing each objective function while calculating the value of another objective in the deterministic case. After obtaining the payoff table, we divide the ranges of the objective functions into four equal intervals, and five grid points are utilized (i.e., [6][7][8][9][10] as the values of ε 2 .

Effect of Parameters' Uncertainty
The effect of uncertainty is studied by varying the conservatism degree and data variability in the robust uncertain parameters of demand and the remanufacturing rate. First, we analyze the changes in the values of the Pareto solutions of different levels of conservatism degree. Assume the positive constant perturbations parameterd m is 45,35,40,30,30,40,30,20,35, and 40, respectively;θ is 0.10. As Γ 1 m , Ψ 1 m and Ψ 2 l ∈ [0, 1], Γ 2 l ∈ [0, 5], we consider these parameters, taking 0%, 25%, 50%, 75%, and 100% of its maximum value, respectively, and dividing them into two groups (i.e., Γ 1 m and Γ 2 l , Ψ 1 m , and Ψ 2 l ) to study the behavior of the Pareto optimal solution. Table 3 reports the Pareto optimal solution under different conservatism degree levels for test model P according to the value of ε 2 .  Table 3 indicates the relationship between conservatism degree and Pareto optimal solutions, which have presented good managerial insights. As expected, the Pareto optimal solutions reveal a monotonic increase when robustness is enforced by enlarging the conservatism degree. When the newly added RHS of the constrained objective function ε 2 is less than a certain value, the solution increases slowly, and the change is not obvious; however, when the value of ε 2 increases gradually to a critical value, the Pareto optimal solution changes obviously. The result also shows that the cost of the supply chains will increase when the social impact is taken into account when enterprises formulate strategies. In this example, when the RHS of the constrained objective function ε 2 is below 10, the economic and social benefits of the network are relatively stable. Table 3 also shows that when one of the two conservatism coefficients is 0, the objective function value does not change. When both groups of coefficients are not 0, the overall trend of the Pareto optimal solution value increases gradually with the increase of the conservatism degree, and in most conservatism cases, the increase is obvious. The reason is that with the increase of the conservatism degree, the range of parameter perturbation becomes larger. The target value increases.
Under the condition of different conservatism degrees, the newly added RHS of the constrained objective function ε 2 is also different. When the conservatism degree is very low, the objective function is a definite optimal solution, which has no effect on the value of ε 2 . With the gradual increase of the conservatism degree, the upper bound of ε 2 will become smaller. In this sense, if the value of ε 2 is too small, it does not affect the optimal solution of the objective function; if the value of ε 2 is too large, it will cause a slack solution or no solution. In a certain range (i.e., [6][7][8][9][10], the value of the Pareto optimal solution will increase when the value of ε 2 increases; we called this value range the "sensitive area", which is of great significance for enterprise decision making. Another important experiment is determining the magnitude of impacts of robust uncertain parameters (demand and remanufacturing rate) on Pareto optimal solution. Suppose all the conservatism degree parameters vary by the same level, and for each assigned conservatism degree, the variability of the uncertain parameters is set to 5%, 10%, 15%, 20%, 25%, and 30% relative to their nominal values. We calculate the Pareto solutions based on different values of ε 2 , and the results are reported in Table 4.  Table 4 presents the change of Pareto solutions when the conservatism degree and robust parameters are simultaneously uncertain. It can be seen that the Pareto solutions will increase when either of these two kinds of uncertainty are enhanced. However, this procedure cannot obtain the Pareto solution efficiently when both of the conservatism degrees and parameters variability increases are too large. The reason can be interpreted as follows. In reality, problem P is similar to a multi-objective mathematical programming problem; when we apply the ε-constraint method to deal with this model, one of the objective functions, i.e., here the supply chain cost, is optimized while using the other objective as constraint, i.e., here the fuzzy membership. On this occasion, when the lower bound of the second objective (the value of ε 2 in problem P) and the uncertainty of the model is too large, the solution will be infeasible or out of the solution space. Therefore, our proposed model is helpful for decision makers to identify the bottlenecks of the supply chain, such as the maximum uncertainty levels or the achievable minimum missed working days.
On the other hand, to show the trend of changes in the Pareto optimal solution values and the probability of constraint violation based on different conservatism degrees, we take the change of remanufacturing rate as an example. When the value of Ψ 1 m and Ψ 2 l is 1, the initial robust equivalence becomes the "interval + polyhedral" equivalence proposed by Bertsimas and Sim [63], assumingd m is 45, 35, 40, 30, 30, 40, 30, 20, 35, and 40, respectively,θ is 0.10, and the value of ε 2 is 8. The calculation result is depicted in Figure 3.
is helpful for decision makers to identify the bottlenecks of the supply chain, such as the maximum uncertainty levels or the achievable minimum missed working days.
On the other hand, to show the trend of changes in the Pareto optimal solution values and the probability of constraint violation based on different conservatism degrees, we take the change of remanufacturing rate as an example. When the value of 1 m  and 2 l  is 1, the initial robust equivalence becomes the "interval+polyhedral" equivalence proposed by Bertsimas and Sim [63], assuming ˆm d is 45, 35, 40, 30, 30, 40, 30, 20, 35, and 40, respectively,  is 0.10, and the value of 2  is 8. The calculation result is depicted in Figure 3.  Figure 3 reveals that improving the conservatism degree of the uncertainty remanufacturing rate and the budget of uncertainty will get a worse Pareto optimal solution value, reflecting the tradeoff between robustness and performance. Additionally, the higher conservatism degree will get a lower constraint violation probability, and the network structure will be more stable. In this case, decision makers can select the appropriate conservatism degree value to obtain objective values and optimal assignments of products from different facilities and determine which facilities should be established.

Effect of Network Structure
In this section, we interpret the effect of uncertainty on the location selection and the network structure in the supply chain network by changing the conservatism degree. For consistency, the values of ˆm d ,  , and 2  are the same as in Figure 2, while the conservatism degree parameters are set to 0%, 50%, and 100% of the maximum value. Then, the location selection result is obtained as shown in Table 5, and the network is illustrated in Figures 4-6.   Figure 3 reveals that improving the conservatism degree of the uncertainty remanufacturing rate and the budget of uncertainty will get a worse Pareto optimal solution value, reflecting the tradeoff between robustness and performance. Additionally, the higher conservatism degree will get a lower constraint violation probability, and the network structure will be more stable. In this case, decision makers can select the appropriate conservatism degree value to obtain objective values and optimal assignments of products from different facilities and determine which facilities should be established.

Effect of Network Structure
In this section, we interpret the effect of uncertainty on the location selection and the network structure in the supply chain network by changing the conservatism degree. For consistency, the values ofd m ,θ, and ε 2 are the same as in Figure 2, while the conservatism degree parameters are set to 0%, 50%, and 100% of the maximum value. Then, the location selection result is obtained as shown in Table 5, and the network is illustrated in Figures 4-6.          As shown in Table 5, the location selection varies with different conservatism degree levels. When the conservatism degree changes from 0 to 50%, the factory site selection is changed from 2, 3 to 1, 2; the distribution center site selection is changed from 1, 5, 6 to 1, 4, 6; the collection center site is changed from 2 to 2.5; and the remanufacturing center site selection is changed from 1 to 1,2. When the degree of conservation changes from 50% to 100%, the location of the factory, distribution center, and remanufacturing center remain unchanged, but the collection center changes from 5 to 2, 3. This is because with the increase of conservatism, the robustness of supply chain networks will be gradually enhanced, and the solutions obtained can adapt to more situations. However, it is also necessary to select facilities with better stability and greater processing capacity, or open more facilities of the same category, so as to meet customer requirements and various constraints. Figures 4-6 depict that the network structure varies greatly with different degrees of conservatism. When the conservatism degree changes from 0 to 50%, the location of each echelon and network flows are different. This reveals a weak network structure. When the conservatism degree changes from 50% to 100%, as the location selection is almost the same, the network is capable of remaining most of its structure, but more facilities (collection centers) still need to be established to meet the demand.
In summary, the robust optimization method utilized in this paper is practical. With the increase of conservatism degree, the network structure can still be applied when the degree of uncertainty parameter changes increases. This will need more resources and costs; therefore, the value of Pareto optimal solutions increases continuously. The change reflects the balance between the robustness and optimality of the model. Improving the robustness of the model can make the model more stable, but the optimality of the model may be reduced. On the other hand, the value of ε 2 has a "sensitive area", which expresses the inevitable and maximum value of social factors; decision makers need to choose the affordable value in different degrees of conservatism and choose the most suitable decision-making scheme for enterprise strategy.
Meanwhile, when the range of uncertainty parameters increases, we need to establish a new facility to satisfy demand, which will increase the huge fixed cost and enlarge the objective value. In addition, under different conservatism degrees, the location of facilities in supply chain networks will change, and transportation routes will also be different. Therefore, decision makers can choose the values of conservative and uncertain parameters according to their strategic positioning of supply chain networks. If decision makers pursue the applicability of the model, the more conservative solution results can be obtained by choosing larger parameters, and the robustness pursuit of the decision maker can be satisfied by increasing the value of the objective. If decision makers are more concerned about the realization of the optimization objective of supply chains, then more definite parameters can be selected to obtain a robust solution for satisfying their own requirements.

Conclusions
In this paper, we developed a robust fuzzy MIP optimization model for the CLSCN design problem considering sustainability and uncertainty. We aimed to minimize the costs of CLSCNs, as well as simultaneously minimizing missed working days by maximizing the fuzzy membership degrees. The other kind of uncertainty considered in this article is the customer demand and remanufacturing rate, which are modeled using the robust optimization approach. Then, the ε-constraint method is employed to solve the developed multiple objective model, generating efficient Pareto solutions. The proposed robust fuzzy model along with the solution methods are tested with a numerical example. The results of the example reveal that the developed model and solution approaches can be effectively used in practice to obtain the economic, environmental, and social benefits from the designing of CLSCs. Investors are able to design the most cost-effective and sustainable supply chain to meet the demand of customers and estimate the social impacts generated by constructing the supply chain. Additionally, sensitivity analysis experiments of parameter uncertainty reveal that the changing demands regarding the conservatism degree and uncertainty parameter perturbation will both affect the Pareto optimal solution. Furthermore, there is a "sensitive area" for the solution to remain feasible, which means that decision makers should select the conservatism degree and uncertainty parameter perturbation rate cautiously. Finally, a test is implemented to illustrate the effect of conservatism degree on the location selection and network structure of the supply chain. According to the results, the position of facilities and network structure will be variable when the conservatism degree changes; therefore, decision makers should modify the value of uncertain parameters according to the actual situation in order to adapt to the complex and changeable market environment flexibly.
In the past, many studies about the CLSCN design problem considering sustainability have taken economic or carbon emission as the objective, but few have considered quantifying social factors in their research [41][42][43][44][45][46][47][48][49][50][51][52]73]. Compared with these research studies, our study considers all three dimensions of sustainability, which can provide better theoretical support for supply chains to realize sustainability. Moreover, we also consider two kinds of uncertainty in our paper and handle them with different methods; the numerical result and the discussion show that uncertainty will change the structure of SCNs and influence objectives. This indicates that considering uncertainty can provide more accurate information for the formulation of a supply chain strategy than a deterministic SCND model [36][37][38][39]72].
As for future research, our model and solution method can be applied in real-world instances such as energy, tires, or the food industry. In addition, heuristic algorithms or hybrid metaheuristics can be employed to solve large-scale optimization problems. The proposed model can also be expanded by adding some more real-life constraints (e.g., inventory control or vehicle routing). Moreover, different transportation modes (railway, shipping . . . etc.) and different production technologies can be considered in the model. Lastly, different measurement approaches and indicators can be employed to evaluate and quantify the sustainability of the proposed supply chain model, which could also be promising areas for future research.
Author Contributions: X.Z. developed the concept and drafted the manuscript. Y.Q. revised the manuscript, and G.Z. supervised the overall work. All authors read and approved the final manuscript.