Fuzzy Bi-Objective Closed-Loop Supply Chain Network Design Problem with Multiple Recovery Options

A network design of a closed-loop supply chain (CLSC) with multiple recovery modes under fuzzy environments is studied in this article, in which all the cost coefficients (e.g., for facility establishment, transportation, manufacturing and recovery), customer demands, delivery time, recovery rates and some other factors that cannot be precisely estimated while designing are modeled as triangular fuzzy numbers. To handle these uncertain factors and achieve a compromise between the two conflicting objectives of maximizing company profit and improving customer satisfaction, a fuzzy bi-objective programming model and a corresponding two-stage fuzzy interactive solution method are presented. Applying the fuzzy expected value operator and fuzzy ranking method, the fuzzy model is transformed into a deterministic counterpart. Subsequently, Pareto optimal solutions are determined by employing the fuzzy interactive solution method to deal with the conflicting objectives. Numerical experiments address the efficiency of the proposed model and its solution approach. Furthermore, by comparing these results with the CLSC network design in deterministic environments, the benefits of modeling the CLSC network design problem with fuzzy information are highlighted.


Introduction
In recent years, the closed-loop supply chain (CLSC) has attracted attention for its economic value, environmental impacts and governmental attention in the area of supply chain management. As a comprehensive chain, it involves not only the traditional forward services provided by the suppliers, manufacturers and distributors, but also the reverse logistics related to collecting, recycling, remanufacturing, repairing and disposing of end-of-life products. This new branch of supply chain management has gained traction both in theory and practice, and the key research directions can be listed as follows. First, there has been extensive research into the concept and solution methods laying the foundation of the CLSC network and focusing on the network hierarchy, facilities classification, processing methods, optimization methods and so forth (e.g., [1,2]). Second, research has explored the design of a multi-objective CLSC network taking into account its economic, efficient and sustainable effects, most of which has aimed to build an economical, environmentally friendly and efficient supply chain (e.g., [3,4]). Third, the design of the CLSC network under uncertain environments considering the high complexity, supply-demand imbalance, and other indeterministic factors in the structure has been explored to simulate real environments (e.g., [5,6]).
The design of a remarkable supply chain should incorporate reverse logistics with forward logistics, plan the nodes and routes so that the resources can flow through and make decisions The tabulated summary presented in Table 1 shows different study prospectives for CLSC network design under the condition of uncertainty. Drawing upon Pishvaee and Torabi [6], we measure the customer satisfaction level by the delayed delivery time. Our paper expands upon the listed articles by presenting a CLSC network which integrates multiple recovery options, including remanufacturing, recycling, repair, selling to secondary markets and disposal, making this the most comprehensive study compared with other papers as shown in Table 1. Moreover, we incorporate the uncertainty of cost, recovery ratio, delivery time, product price and quality in the model. To the best of our knowledge, very few research works have considered all of these complex uncertain factors together in one model.
In summarize, this paper differentiates itself from the existing related research and contributes to the study of the CLSC network design from the following aspects: (1) considering the most comprehensive recovery options simultaneously, as compared to the other works in Table 1; (2) establishing a fuzzy bi-objective MILP model which maximizes both company profit and customer satisfaction level for the CLSC network design under fuzzy environments; (3) presenting a two-stage fuzzy interactive solution framework to solve the formulated fuzzy optimization model by integrating the fuzzy expected value operator, fuzzy ranking method and fuzzy interactive solution method; (4) comparing the performance of the fuzzy CLSC model with the deterministic one in [3] by several numerical experiments to demonstrate the importance and necessity of designing the CLSC network under fuzzy environments.
The remainder of this paper is organized as follows. Section 2 presents the problem description and mathematical formulation of the CLSC network design. Section 3 discusses how to process the fuzzy parameters in the formulated fuzzy bi-objective MILP model. Following that, a two-stage fuzzy interactive solution framework for solving the MILP model is constructed in Section 4. In Section 5, some numerical experiments and comparisons are conducted to verify the performance of the proposed model and solution framework. Finally, some pivotal conclusions are addressed in Section 6.

Problem Description
In Jiang et al. [3], the CLSC with multiple recovery options considering enterprise profit and customer satisfaction under deterministic environments have already studied in a detailed manner, and in this paper the problem is extended to a fuzzy one by inquiring into the impact of uncertainty on the network design of the CLSC. Specifically, the discussed CLSC network consists of two customer markets and various facilities serving these customers, such as plants, disposed material processing centers and raw material recycling points, wherein the information of the customers (i.e., demands, expected delivery time and degree of satisfaction), the cost coefficients (e.g., for start-up cost, manufacturing cost or processing cost and unit transportation cost) and capacity of these facilities and various recovery rates of recycled products are all fuzzy numbers. The structure of the CLSC and its corresponding resource flow are depicted in Figure 1.
As we can see, in the forward logistics, the plants process raw materials into new products and transfer them to the primary customers through the distribution centers. In the reverse logistics, then, the disassembly centers collect waste products from the primary customers, test or classify them and determine their recovery modes within the following four options. The products with less damage are directly repaired and then delivered to the secondary customers through the redistribution centers, which are responsible for the delivery of the reprocessed products. The partially usable products are disassembled at the disassembly centers and sent to the plants for remanufacturing and delivered to the secondary customers through the redistribution centers afterwards. For the waste products that cannot be repaired or used for remanufacturing, some useful raw materials are disassembled and recycled at the raw material recycling points, and the worthless materials are sent to the processing centers for disposal. The above-mentioned network provides a variety of recycling options, and thus can maximize the recycling value of waste products. Meanwhile, the design separates the demand markets into two parts, making the network more practical for the operation of many enterprises. Furthermore, as a general network configuration with manufacturing-repairing-remanufacturing-recycling mode, the proposed CLSC network can be diffusely used in the household electrical appliance manufacturing, automobile industry, and other manufacturing industries.
The purpose of the CLSC network design is to make decisions on the quantity and locations of the facilities/the resource flow in the network, and the desirable configuration we aim to design is the one that can bring economic benefits and increase operating efficiency, so as to help the enterprises improve supply chain management, reduce energy consumption, waste emission and reinforce enterprise image. In a network consisting of numerous entities, a series of activities involved will bring revenues to the enterprises, but simultaneously increase operating costs. Consequently, the enterprises had should more attention to the actual economic benefits brought by implementing the CLSC management. Meanwhile, improving the efficiency of the network is conductive to reduce the delivery time and gain customer satisfaction level, and ultimately enhance the economic value of the enterprises. Therefore, this paper studies the design of a bi-objective CLSC network considering company profit and customer satisfaction measured by the delayed delivery time.

Mathematical Formulation
This section simplifies and abstracts the problem of the CLSC network design into a quantitative fuzzy bi-objective optimization model. All of the parameters, decision variables, objective functions and constraints involved are briefly introduced as follows (all notations used in the paper are collected in Appendix A).
Referring to the CLSC network structure displayed in Figure 1, the sets of the potential locations of the plants, distribution centers, redistribution centers, disassembly centers and disposal centers, and the fixed locations of the primary and secondary customers are denoted by I, J, M, L, P, K, and N, respectively. Accordingly, the elements in these sets are respectively indexed by i, j, m, l, p, k, and n.
Constructing and operating this network suffers three types of costs, i.e., the fixed facility establishing cost, the variable processing or handling cost at these facilities and the transportation cost between these facilities (and customers). The fixed start-up costs for plant i, distribution center j, disassembly center l, redistribution center m and disposal center p are, respectively, represented by f i , f j , f l , f m and f p , with the maximum manufacturing or processing capacities of p i , p j , p l , p m , and p p , respectively. Correspondingly, the unit variable manufacturing or processing costs at these facilities include the unit manufacturing cost mc i , unit remanufacturing cost vc i , unit handling costs hc j , hc l , hc m , unit repairing cost rc l and unit disposing cost dc p . Furthermore, the unit transportation costs between the facilities and customers are expressed as tc ij , tc jk , tc li , tc lm , tc l p , tc im and tc mn . Especially, the unit collecting and transporting costs from primary customer k to disassembly center l are denoted by cc kl . The maximum fraction of collecting waste products is ω, and the fractions of disposing and repairing are θ 1 and θ 2 .
Selling products and recycled raw materials are the main revenue sources of the CLSC. s 1 , s 2 and s 3 denote the unit revenues of new products, remanufactured products and recycled raw materials, respectively. The demands of primary customer k (purchases new products) and secondary customer n (purchases remanufactured ones) are denoted by d k and d n . In addition, a primary customer's satisfaction level is dependent on whether the products are delivered on time, which can be judged by comparing the delivery time from distribution center j to customer k, dt jk , with the customer's expected delivery time, et k .
To find out the optimal locations of the facilities in the network and the flow distribution among them are the main tasks of the CLSC network design. Therefore, two kinds of decision variables are included in the model. One is the continuous and positive variables used to represent the flows between the facilities (and customers) in the network, and the other is the 0-1 variables used to indicate whether a candidate site is selected for establishing the facility. These two kinds of decision variables are expressed as X ij , X jk , X kl , X li , X lm , X l p , X im , X mn , X l and Y i , Y j , Y l , Y m , Y p , respectively.
Considering the objectives of maximizing company profit ( Z 1 ) and improving customer satisfaction level by decreasing the delayed delivery time ( Z 2 ), the model of the fuzzy bi-objective CLSC network design is expressed explicitly as follows, where TI, FC, TC, MC, HC, CC, RMC, RC and DC are the total revenue, fixed establishment cost, transportation cost, manufacturing cost, handling cost, collection cost of new products, remanufacturing cost, repair cost and disposal cost of waste products, respectively. Detailed formulas for the calculation of these terms as well as the constraint functions (A10)-(A26) are presented in Appendix B.
Note that the problem of the fuzzy CLSC network design expounded in our paper is an extension of the work presented in [3]. The settings of the two problems are identical except for the assumption that the above parameters, whose notations have tilde symbols on their heads, are uncertain and modeled as triangular fuzzy numbers here. Therefore, we only briefly introduce the fuzzy version of the mathematical model and put the detailed formulas of the model in the appendix to keep this paper self-contained. For more detailed descriptions and discussions on the CLSC network design problem under deterministic environments, readers are recommended to [3].
It should be also noted that although the above model (1) is similar to the work of [3] in form, they are quite different from each other in essence, since the model under fuzzy environments is of higher complexity. Specifically, owing to the existence of the fuzzy parameters, the model cannot be directly solved by the classic optimization techniques for deterministic programming models. Consequently, in what follows, we discuss how to handle these fuzzy parameters and present a two-stage solution framework to solve this fuzzy bi-objective programming model.

Processing Fuzzy Parameters in the CLSC Model
In model (1), both the objective functions and constraint functions contain fuzzy parameters. Thus, the objectives cannot be directly optimized, and whether the constraints are satisfied cannot be judged as well. Considering that, a general fuzzy programming model is formulated below as model (2) for discussing how to manipulate these fuzzy parameters, where x = (x 1 , x 2 , · · · , x p ) is a p-dimensional vector representing the decision variables, A i and C T denote the fuzzy coefficient matrices of the constraints and the objective function, respectively, and B i is the vector representing the fuzzy bound of the constraints. Afterwards, by means of the fuzzy expected value operator and fuzzy ranking method, the model can be transformed into its deterministic counterpart, in which both the fuzzy objective functions and constraint functions in the original model (1) are explicitly expressed.

Fuzzy Objective Functions
For converting the fuzzy objective functions into their crisp counterparts, it is natural to approximate them by their expected values as illuminated by Jimenez et al. [23] and Liu [24]. To this end, the definitions of the triangular fuzzy number (TFN) and the expected value of a TFN as well as its linearity property should be elaborated in advance.
is a fuzzy set on the real line, whose membership function can be expressed as where f ξ (x) and g ξ (x) are the left and right shape functions, and ξ p , ξ m , and ξ o are the lower bound, median value, and upper bound, respectively, of the TFN, as shown in Figure 2.
Definition 2 (Jimenez et al. [23]). The expected interval (EI) and expected value (EV) of a TFN ξ = (ξ p , ξ m , ξ o ) are, respectively, Theorem 1 (Jimenez et al. [23]). Let ξ and ζ be two TFNs. For any real numbers a and b, we have According to the above theorem, we convert the expected value of the fuzzy objective function Z 1 in model (1) into the following form, where EV[ TI] can be further calculated as However, the second objective function Z 2 in model (1) cannot be similarly calculated as presented above by applying the linearity of the EV operator (Theorem 1) owing to the existence of the operator (·) + . In some previous studies, the operator was simply approximated, such that where dt and et are TFNs (see, e.g., [6]). Clearly, this treatment would lead to serious deviation from the real exact value in some cases. Therefore, in this paper, based on the concepts of EI and EV proposed by Jimenez et al. [23], we exactly calculate the value of EV[( dt − et) + ] as follows, where

Fuzzy Constraint Functions
Analogously, with the purpose of determining whether a constraint with fuzzy parameters is valid, the method proposed by Jimenez et al. [23] for comparing two fuzzy numbers is adopted in this paper as follows.
Definition 3 (Jimenez et al. [23]). For TFNs ξ and η, the degree of ξ > η is defined as After applying the above definition into model (2), the constraints A i x ≥ α B i , i = 1, 2, · · · , l, are equivalent to the following group of equations, where α is the confidence level used to adjust the degree of the constraints Similarly, the degree of equality of two fuzzy numbers was defined by Parra et al. [25] below.
Definition 4 (Parra et al. [25]). For any two TFNs ξ and η, There are three kinds of constraints in the fuzzy programming model (2), including the fuzzy inequality, fuzzy equation and crisp constraints. According to Definition 3, the fuzzy inequality constraints can be transformed into the corresponding deterministic inequality ones formulated as Equation (7). Furthermore, the fuzzy equality constraints can be converted into according to Definition 4.

Transformation
Referring to the methods presented before, the general form of a fuzzy programming model (2) can be transformed into the following deterministic counterpart, By converting model (2) into model (10), the fuzzy inequality constraints are transformed into a group of crisp inequalities, and the fuzzy equality constraints are converted into two groups of crisp inequality constraints, respectively, with the confidence level of α, whereas the crisp ones are maintained.
Similarly, the fuzzy constraints in the fuzzy CLSC model (1), i.e., (A10), (A11), (A13), (A15)-(A17), (A20)-(A24), are, respectively, transformed into the crisp ones as follows, Consequently, a deterministic counterpart of the fuzzy CLSC model (1) is deduced as follows, where the EVs of (A2)-(A9) are calculated as Equation (4), and EV[( dt jk − et k ) + ] is calculated by Equation (6). So far, the fuzzy bi-objective model (1) for the problem of designing a fuzzy CLSC network has been converted into a crisp bi-objective mixed 0-1 integer linear programming (CB-MILP) model (26). In order to achieve a compromise between the two objectives within the model, a solution framework on the basis of the fuzzy interactive approach is proposed in the following section.

Interactive Approach to Multi-Objective Programming
Multi-objective programming is a model that optimizes multiple objectives simultaneously while the constraints are all satisfied, whose general form is where x = (x 1 , x 2 , · · · , x p ) is a p-dimensional decision vector, f 1 (x), · · · , f m (x) are m objective functions, and the feasible set is denoted by S.
In order to solve this kind of multi-objective problems, Zimmermann [26] proposed a fuzzy interaction method, in which the implementation degree of each exact objective is provided according to the preferences. The method can help to transform a multi-objective programming into a single-objective Min-Max model, thereby balancing multiple goals. While being applied, some expansions have been made by Werners [27], Torabi and Hassini [28] (called the TH method) and Selim and Ozkarahan [29] (called the SO method). idea of the Min-Max model is to construct the membership function of the target value so as to quantify the fuzzy target of decision-makers. Subsequently, based on the membership function of each objective, a synthesis evaluation function of the single objective is constructed to optimize the interaction among different objectives, and then the optimal solution is obtained. Therefore, the key of applying this method lies in the construction of the membership function and synthesis evaluation function.
To do so, we first determine the membership function of objective i(1 ≤ i ≤ m), say µ i (x). As a preliminary step, the positive ideal solution (PIS) and negative ideal solution (NIS) of each objective are calculated. In model (27), there are m objectives of minimization, and thus the minimum value of each objective function is the corresponding PIS, recorded as f PIS i = min f i (x), x ∈ S. After that, fix the objective function and plug the other x values into it, then the maximum objective function value can be deduced as f N IS i . For instance, letting m = 3, we have the PIS of each objective and the corresponding decision vector x, i.e., ( f PIS 1 , , respectively. Then, the NIS of each objective is, Based on the concepts of PIS and NIS, the membership function µ i (x) is constructed by combining the subjective preferences of decision-makers, among which the most commonly used linear form is shown below, Second, combining with the membership function in Equation (29), we employ the synthesis evaluation function λ(x) to convert the bi-objective model (26) into a single one by means of the TH method [28], which adds the compensation factor γ on the basis of Min-Max method as follows, where x is the decision vector consisting of continuous and 0-1 variables, Z g (x) = EV[ Z g (x)], g = 1, 2, and F represent the objective functions and feasible set of model (26), respectively, and w g represents the predetermined importance weights of different objectives. It should be pointed out that the decision variable λ 0 = min g {µ g (Z g (x))} expresses the minimum implementation level of the objective function, and γ ∈ [0, 1] is the compensation factor of the minimum implementation level. Through the TH evaluation function, we find that a trade-off between the weighted and the minimum operators can be achieved, hence the decision-makers can adjust the parameters w g and γ to obtain equilibrium and non-equilibrium solutions between the two objectives according to their preferences.

Two-Stage Fuzzy Interaction Solution Procedure
As elaborated in Sections 3 and 4.1, the solution framework for solving the problem of the fuzzy CLSC network design can be integrated into a two-stage fuzzy interaction solution procedure. In the first stage, transform the original fuzzy bi-objective CLSC model (1) into a crisp bi-objective MILP model (26) with the feasibility level, α, based on the fuzzy expected value theory and fuzzy ranking method introduced by Jimenezi et al. [23], as performed in Section 3. Following that, by means of the fuzzy interaction method proposed by Zimmermann [26] and improved by Torabi and Hassini [28], the crisp bi-objective MILP model is turned into a deterministic single-objective model (30) in the second stage, as presented in Section 4.1. Since the fuzzy interaction method calculates the realization level of each objective function directly, decision-makers can obtain equilibrium solutions between multiple objectives and non-equilibrium solutions with their subjective preferences, and provide a variety of valuable decision-making schemes by adjusting the related parameters including the feasibility level α, the compensation factor γ and the weight coefficient w g . In a word, the specific steps of the two-stage fuzzy interaction solution method can be depicted in Figure 4 and summarized in the following: Step 1: Transform the objective functions with fuzzy parameters, Z 1 and Z 2 in the original CLSC model (1), into the crisp form, Z 1 and Z 2 in model (26), by means of the fuzzy expected value theory elaborated in Definition 2 and Theorem 1.
Step 3: Calculate the PIS Z α-PIS 1 and Z α-PIS 2 for each objective at the feasibility level α and the corresponding decision schemes x 1 α-PIS and x 2 α-PIS , and then figure out Z α-N IS Step 4: Based on the PIS and NIS, compute the realization levels of the objectives in model (26), µ g (x), which are definitely constructed as follows, Step 5: Establish the synthesis evaluation function λ(x) via the multi-objective interactive TH method to convert the CB-MILP model (26) into a deterministic single-objective MILP model (30), and then determine the compensation factor γ and the weight coefficients w g of the evaluation function to solve it.
Step 6: Adjust the parameters α, γ (if necessary, the weight coefficients w g can be changed), and repeat Steps 2-6 until the solution is satisfactory from the perspective of the decision-maker.

Numerical Analysis
To confirm the validity of the transformed CB-MILP model and the proposed solution framework, some numerical experiments are conducted in this section. Afterwards, some comparisons with the deterministic problem solved by Jiang et al. [3] is presented. All the experiments in this paper were run in the computer environment of Intel(R) Core(TM) i5-8265@1.60ghz, and the solutions are derived from CPLEX 12.8.0 called by Matlab 2018 R(b).

Experimental data
The experiments consider a CLSC network with 43 sites, including 4 factory alternative points, 8 distribution center alternative points, 12 primary customer points, 6 disassembly center alternative points, 4 alternatives for redistribution centers, 6 secondary customer points and 3 alternatives for processing centers. More details of the sites are accessible in [3]. The fuzzy parameters in the problem are all assumed as TFNs, the lower bound, median value and upper bound of which are generated following the uniform distribution. The specified process is as follows. First, randomly generate a median value of ξ m for each parameter. Second, in order to avoid loss of generality, generate two random numbers (τ 1 ,τ 2 ) following the uniform distribution between 0.2 and 0.5, respectively, as the left and right spreads. Third, calculate the lower and upper bounds of TFN ξ, ξ p and ξ o , as To facilitate the comparisons between our work with the deterministic problem, the medians of the fuzzy numbers are set as the corresponding deterministic values in [3]. Furthermore, the left and right spreads, (τ 1 , τ 2 ), are randomly generated as shown in Table 2.

Two-Stage Fuzzy Interaction Analysis
For solving the above problem, we plug the parameters into the fuzzy bi-objective programming model (1), carry out the crisp equivalence transformation, and obtain the CB-MILP model (26) for the experiments. Afterwards, by setting different feasibility levels α, the PIS of the profit and the NIS of the delayed delivery time representing the optimal solution results when only one of the two objectives is considered are figured out, as shown in Figure 5. It can be observed that, with the increase of α, the optimal value of the profit marked by the dotted black line decreases, and the optimal value of the delayed delivery time represented by the solid blue line increases. Since the objective functions as well as the constraints of model (1) are all linear combinations of the triangular fuzzy parameters, the two lines are changing in a linear fashion.
Then, substitute the optimal decision scheme of one goal into the other according to Equation (28), we obtain the worst values of the two goals, i.e., the NIS of the profit and the PIS of the delayed delivery time, as presented in Figure 5. Obviously, with the increase of α, the worst cases show a strict monotone tendency and the differences between the PIS and NIS of the objective functions both shrink. As mentioned in Section 4.1, the PIS and NIS prescribe a limit to the range of one objective while another objective is optimized, thus Figure 5 indicates that the range of the two objectives narrows with the increase of α. This trend is consistent with common sense, since the feasibility level, α, is a parameter representing the degree of constraint relaxation, i.e., a higher feasibility level is equivalent to a greater degree of constraints satisfaction and a smaller feasible region of the solution. The experimental results can provide decision-makers with the range of the two objective functions at different feasibility levels and help to avoid the risks caused by uncertain factors.  Based on the above PIS and NIS of the two objectives at different feasibility levels with the weight coefficients of w 1 = w 2 = 0.5, namely, maximizing the profit is as important as minimizing the delayed delivery time, and the compensation factor of γ = 0.5, the realization levels of the two objectives, µ 1 and µ 2 , by using the TH fuzzy interaction method are obtained as shown in Table 3. In general, with the increase of α, the realization levels of the two objectives are relatively high and experience strictly monotone increasing trends. Typically, when α = 0, the most relaxed constraints result in the largest feasible region, and the solution which can immunize against the largest parameter perturbation is achieved, with the minimum realization levels of the two objectives, 80.29%. By contrast, when it comes to α = 1, the tightest constraints lead to the most reliable results with the realization levels of 95.32% and 98.33%. Furthermore, the average realization levels are, respectively, 87.47% and 87.74%, indicating that the TH method is competent to obtain stable and precise solutions. On the other hand, we can see from Table 3 that there is little difference between the realization levels, the reason is that the compensation factor is set to be 0.5, and accordingly the objective function of model (30), λ = 0.75 min{µ 1 , µ 2 } + 0.25 max{µ 1 , µ 2 }, endows a greater weight to the goal with lower realization level, which means the model pays more attention to the short board and is aimed to obtain a compromise and equilibrium solution between the two objectives. For more information, we set the compensation factor as γ = 0, i.e., λ = 0.5µ 1 + 0.5µ 2 , and record the realization levels as well as the gap between them, µ 1 − µ 2 , at different feasibility levels in Table 4. The results explicitly show that the gap tends to decrease with the increase of α and the realization levels of maximizing enterprise profit is constantly higher than that of minimizing customers' delayed delivery time except for the case of α = 1, demonstrating that the most balanced solution can be obtained when α is between 0.9 and 1. Therefore, to further explore the sensitivity of the TH method to γ, we fix the feasibility level at 0.9 after comparing the solutions of different α and figure out the realization levels of the two objectives as well as the gap between them, µ 1 − µ 2 , as illustrated in Table 5. It shows that the TH method tends to seek equilibrium solutions, for the gap narrows down with the increase of γ. Especially, when γ ≥ 0.4, µ 1 is equal to µ 2 and is constantly unchanged. That is, as γ gradually approaches 1, the results pays more attention to optimize the objective with lower realization level, so that a more equilibrium state can be reached. Additionally, we find that µ 1 and µ 2 evince similar trend of monotonic increase, which is consistent with the conclusion of [28], verifying the remarkable monotonicity of multiple realization levels with the increase of γ while introducing the TH method. Besides, it should be pointed out that, at a fixed feasibility level, when the compensation factor γ is small, the TH method can obtain non-equilibrium solutions between multiple objectives, as the realization levels are relatively high. Ultimately, in order to investigate the joint effect of the weight coefficients and the compensation factor (w 1 , w 2 and γ) on the realization levels of the two objectives, we fix α at 0.5 and accordingly depict the results in Figure 6. It can be seen from Figure 6a that the weight coefficients of the two objectives have no effect on the realization level of maximizing the enterprise's profit when γ > 0.6 (with µ 1 remaining steady at 90%), yet have significant influence when it comes to γ ≤ 0.6. Specifically, when w 1 ≥ 0.4 (w 2 ≤ 0.6), µ 1 maintains a high level of more than 90%; and when w 1 < 0.4 (w 2 > 0.6) and γ ≤ 0.6, it decreases rapidly and reaches the bottom of 0 with w 1 = 0 and γ = 0. In Figure 6b, the similar conclusion can be drawn as Figure 6a since the tendency of µ 2 in Figure 6b is close to that of µ 1 in Figure 6a if the x-axis, w 1 , is replaced by w 2 . Hence, to get more equilibrium solutions, it is supposed to set a relatively concentrated value of w 1 and w 2 and a relatively high value of γ. To sum up, the above analyses show that while designing the CLSC network in fuzzy environments, the decision-makers can obtain satisfactory solutions by modifying the parameters, i.e., the feasibility level, the compensation factor, and the weight coefficients, according to their preferences, and the fuzzy settings provide more flexible and diversified options for decision-making. In the next section, we will elaborate this advantage in detail by some comparisons with the deterministic problem presented by [3].

Comparisons with the Deterministic Problem
By fixing the parameters at γ = w 1 = w 2 = 0.5 and solving the proposed fuzzy model and the deterministic one in [3], the comparative results on the optimal solutions at different feasibility levels are depicted in Table 6 and Figure 7. Table 6. Optimal solutions of the fuzzy and the deterministic model at different feasibility levels. As we can see from the four curves representing the optimal solutions of the objective functions obtained by using the TH method under deterministic and fuzzy environments, both the profit and the delayed delivery time in the deterministic problem are less than that in the fuzzy one. Furthermore, it also can be observed that the values of the two objectives under fuzzy environments decrease with the feasibility level α increases, whereas the results under deterministic environment are constantly unchanged.
Explicitly, the delayed delivery time in fuzzy cases is much longer. At the root, under deterministic environments, Jiang et al. [3] calculated the delayed delivery time by adding up the non-negative differences between the actual and the expected delivery time, which is equivalent to the treatment of Equation (5) corresponding to the fuzzy cases. However, in our fuzzy settings, the difference between the fuzzy expected delivery time and the fuzzy actual delayed delivery time are calculated by Equation (6), since we have proved in Section 3 with concrete data that the difference between fuzzy numbers cannot be directly calculated by their EVs and Equation (5) may result in a huge deviation. Thus, a big gap between the delayed delivery time under different environments arises, as shown in Figure 7.
On the other hand, minimizing the delayed delivery time and total profit are two conflicting objectives, that is, a shorter delayed delivery time brings a higher service level and higher cost/lower profit simultaneously. Therefore, in Figure 7, the profit under fuzzy environments is also larger than the deterministic case. Specifically, the maximum disparity of the two lines, representing the profits under different environments (4,518,735 and 3,878,633, respectively) at α = 0, occupies about 17% (640,102/3,878,633) of the profit in the deterministic case. Then, with the increase of α, the disparity decreases. When coming to α = 1, it narrows down to the minimum (of 192,911).
Furthermore, it is universally acknowledged that the optimal solution of the deterministic case is a degenerate version of the fuzzy case. Henceforth, with the increase of α, both the profit and delayed delivery time are supposed to decrease to the deterministic values. However, due to the treatment of calculating the EV of the delayed delivery time aforementioned, the downtrend of the disparity between the profits under fuzzy and deterministic environments is remarkably more dramatic than that of the delayed delivery time. The results indicate that under fuzzy environments, while the realization levels of the objectives are relatively high, as the constraints become more stringent, improving customer service level is costly for the enterprise, which is quite different from the deterministic cases. Besides, under fuzzy environments, it is more flexible for managers to make decisions by adjusting the feasibility level.
From the above analyses, it is not difficult to conclude that putting the problem of the CLSC network design under fuzzy environments is more reasonable from the perspective of improving enterprise's profit and flexible management. Since fuzzy numbers are sets on the real line (see Definition 1), it is natural for some researchers to argue that we may set the values of the variables as crisp intervals so that the problem of the CLSC network design can be studied under deterministic environments. Considering that, the following experiment is designed to manifest the essential distinction between the deterministic and fuzzy models.
First, by keeping the other parameters in [3] unchanged and increasing or decreasing 0 to 50% of the total cost, the two lines representing the range of the profits are obtained, as plotted in Figure 8a. Since the deterministic total cost can be regarded as a fuzzy number with the spread of 0, we can realize the equivalent increase and decrease of the total cost by adjusting the spread of the fuzzy number under fuzzy environments, that is, expanding the right spread is equivalent to increasing the total cost under deterministic environments, and expanding the left spread is equivalent to decreasing the deterministic total cost. Henceforth, fix γ at 0.9, adjust the spreads of the cost coefficients to realize the equivalent variations of the total cost, and set α as 0, 0.5, and 0.9, respectively, the three curves standing for the optimal profits under fuzzy environments can be figured out and depicted in Figure 8a. It should be noted that in our solution framework, the fuzzy parameters in the original model (2) are all deffuzified by their EVs. Thus, the profit obtained are a crisp value rather than a fuzzy number and its variation is relatively lower with regard to the variation of the total cost. However, by setting the value of the total cost as a crisp interval and substituting the extreme values into Jiang et al. [3]'s model, we can merely obtain the upper and lower limit of the profit whose range is proportionally changing with respect to the variation of the total cost. As a result, the three curves are in the middle of the two lines, demonstrating that the solutions of the fuzzy settings are less fluctuant than that of the deterministic settings. Additionally, we can see that the three curves representing the solutions at different feasibility levels are quite different from each other. As mentioned in Definition 3, the feasibility level α is the essential parameter of transforming the fuzzy constraints into their crisp counterparts, and a higher level is equivalent to more strict constraints. Therefore, the solutions with higher feasibility levels are less fluctuant, which is consistent with Figure 8a that the dashed curve with α = 0.9 is more stable than the dotted one with α = 0 and the dotted-dashed one with α = 0.5, and the profit of α = 0.9 is much higher even though the variation of the total cost is large.
Analogously, while keeping the other parameters unchanged, a comparative experiment on the influence of total cost over another objective of minimizing the delayed delivery time under deterministic and fuzzy environments is also conducted, and the results are depicted in Figure 8b. We can also conclude that the fuzzy model is superior to the deterministic one, and the curve with α = 0.9 is most stable among the three curves. Especially, the delayed delivery time at α = 0.5 is mostly lowest, showing that 0.5 is more appropriate for the decision-maker from the perspective of improving customer service levels if the variation of the total cost is small.  Above all, the comparative results of the optimal solutions illustrated in Table 6 and Figure 7 confirm the superiority of the fuzzy model proposed in this paper over the deterministic one presented by Jiang et al. [3] in maximizing enterprise's profit and flexible management. In addition, the curves (representing the fluctuation of the profit/delayed delivery time in the proposed fuzzy settings which assumes the uncertain total cost as fuzzy numbers) and lines (standing for the upper and lower limits of the objectives under deterministic environments) in Figure 8 demonstrate the reasonableness and necessity of taking fuzziness into account while designing the CLSC network, since the solutions under fuzzy environments are less fluctuant and can offer more accurate information about the objectives. If the fuzzy approach was not applied and the market's uncertainty was handled by restricting the parameter's range, then we can merely obtain an interval of the profit/delayed delivery time. Such a result could be of high uncertainty and less useful for decision-making. Besides, the comparisons between the curves with different α show that the decision-maker can obtain stable and ideal solutions by adjusting the feasibility level, i.e., the key parameter for transforming the fuzzy constraints into crisp ones.

Discussions
In this section, several experiments are carried out to validate the feasibility and necessity of investigating the CLSC design problem under fuzzy environments and the performance of the proposed solution framework in solving the fuzzy CLSC model. In order to exemplify the superiority of the fuzzy model over the deterministic one, the data set is selected from the work of [3] and modified by replacing the crisp parameters into symmetric triangular fuzzy ones. The results with different feasibility levels α, compensation factor γ, and weights of objectives w indicate that these parameters play crucial rules on the performance of the proposed method. Specifically, a higher feasibility level meaning stricter constraints is conducive to a higher-quality and more stable solution, a lager compensation factor and more evenly weights tend to make the solution more balanced.
Furthermore, the comparative experimental study between the fuzzy model and the model in [3] show that putting the problem of the CLSC network design under fuzzy environments is more reasonable from the perspective of improving enterprise's profit and flexible management since the profit of the fuzzy model is higher than that of the deterministic one, and the differences of the profits and delayed delivery times both decrease with the increase of the feasibility level. In the last two experiments, we also illustrate the difference between handling the market's uncertainty by restricting the parameters' range and assuming the uncertain parameters as fuzzy numbers, that is, the deterministic settings can merely provide a range of the objectives, whereas the fuzzy settings can offer more accurate and stable information.
It should be pointed out that although the proposed solution framework is tested only on a bi-objective case with symmetric fuzzy parameters, it can be extended to solve many other CLSC variants with slight modifications, such as the problems with multiple objectives, more diversified recovery options or complicated risk preferences. Therefore, it is highly flexible and applicable for the decision makers who are in need of solving various types of practical CLSC problems. Moreover, it is easy to understand and implement, and is very successful in returning high-quality solutions since the realization levels of both objectives are demonstrated to be high as long as the above-mentioned parameters are reasonably selected.

Conclusions
This paper studies the network design problem of the CLSC with multiple recovery models by putting it under fuzzy environments and establishing a fuzzy bi-objective mixed linear programming model, wherein the objectives of maximizing enterprise's profit and maximizing service level measured by delayed delivery time are taken into consideration. For dealing with the fuzzy parameters and two objectives in the proposed model, a two-stage solution framework including the equivalent transformation from the original fuzzy model to a crisp one by means of the fuzzy expected value operator and fuzzy ranking method and the switch from the bi-objective programming to a single-objective one via the fuzzy interactive solution method was elaborately introduced.
Then, through the experiments of the CLSC network design under fuzzy environments at different feasibility levels, weights of objective functions and compensation factor, the good performance of the proposed model and the two-stage solution framework was verified. The results also demonstrated the importance of selecting appropriate parameters, which can provide guidance for practical decision-making for obtaining equilibrium and non-equilibrium solutions of the problem. Eventually, the objectives of the fuzzy and deterministic models were compared in two ways, jointly proving the flexibility, accuracy, practicality, and necessity of the proposed fuzzy model.
Practically speaking, by assuming fuzzy parameters, the market's uncertainty can be taken into account while making decisions, therefore the proposed model is more flexible and practical, which can not only offer a broader perspective to scholars, but also provide some insights and inspirations to the enterprises' decision-makers. As a consequence, the model can be easily extended to many kinds of supply chains with diverse properties. In the future work, we will consider more realistic features to design more practical CLSC networks, such that (1) considering the disruptions in the model to enhance the flexibility of the CLSC, e,g., the influence of COVID-19 on most supply chains; (2) considering the inventory management in the reverse flow of the CLSC network; (3) introducing more accurate and fast heuristic algorithms to solve the large-scale problems so that the methods can be more applicable for some practical cases.
Author Contributions: Conceptualization, K.W., W.X., Q.Z. and J.Z.; methodology, W.X., Q.Z., K.W. and J.Z.; software, W.X.; validation, K.W., H.L. and W.X.; formal analysis, W.X., H.L. and J.Z.; investigation, Q.Z.; data curation, W.X. and K.W.; writing-original draft preparation, W.X.; writing-review and editing, H.L. and J.Z.; visualization, W.X.; supervision, K.W. All authors have read and agreed to the published version of the manuscript. There are 4 types of constraints. The demand constraints ensure that the demand of the primary customers must be met, whereas the secondary customers should be satisfied as far as possible. They are formulated as, ∑ m∈M X mn ≤ d n , ∀n ∈ N.
The flow balance constraints keep the amount flowing into the corresponding facilities equaling to the amount flowing out to its downstream, being represented as