Closed-Loop Supply Chain Design with Sustainability Aspects and Network Resilience under Uncertainty: Modelling and Application

In this study, we present a multiobjective mixed-integer nonlinear programming (MINLP) model to design a closed-loop supply chain (CLSC) from production stage to distribution as well as recycling for reproduction. The given network includes production centers, potential points for establishing of distribution centers, retrieval centers, collecting and recycling centers, and the demand points. The presented model seeks to ﬁnd optimal locations for distribution centers, second-hand product collection centers, and recycling centers under the uncertainty situation alongside the factory’s ﬁxed points. The purpose of the presented model is to minimize overall network costs including processing, establishing, and transportation of products and return ﬂows as well as environmental impacts while maximizing social scales and network ﬂexibility according to the presence of uncertainty parameters in the problem. To solve the proposed model with fuzzy uncertainty, ﬁrst, the improved epsilon ( ε )-constraints approach is used to transform a multiobjective to a single-objective problem. Afterward, the Lagrangian relaxation approach is applied to eﬀectively solve the problem. A real-world case study is used to evaluate the performance of the proposed model. Finally, sensitivity analysis is performed to study the eﬀects of important parameters on the optimal solution.


Introduction
Manufacturing supply chains have experienced incredible advancements in knowhow including the adoption of efficient technologies, flexibility, mass production enabled by knowledge management, computerization, and rapid reconfigurations over generations of products. As life cycles are continuously shortening thus increasing competitiveness [1], but also competition, manufacturers have to gain more efficiency in their supply chains (SCs) [2]. ere are two types of flows that need to be managed: the first one is the forward SC including the process from procurement of raw materials until the transfer of products and/or services to customers and the second one is the backward/reverse SC where defective products or those handed back from customers after their use to the manufacturer are disassembled, reworked, and reused for production [3]. Such return flows can be handled in the same facilities sharing material handling equipment and infrastructure [4], so that storage and handling processes for backward and forward flows should be considered in an integrated way. Mathematical models can provide decision support for planning these flows [5].
ey are generally referred to as "close-loop" supply chains. Considering them together assures the avoidance of suboptimal decisions [6].
Environmental pollution and insufficiency of resources have become significant and urgent issues during recent years. Countries and international organizations are emphasizing environmental protection and conservation of energy. Manufacturing as a great consumer of the latter gains a special role of adhering and implementing sustainability measurements. Sustainability can be defined as "using resources to meet the present needs without compromising the ability of future generations to meet their own needs" [7]. It bases on three pillars [8], i.e., (1) environmental, (2) economic, and (3) social, where the environmental pillar encompasses, e.g., measures such as the reduction of environmentally harmful emissions. e economic aspects include profitability objectives of companies, whereas the social pillar encompasses aspects of social prosperity, e.g., the maintenance of jobs or creation of new job opportunities and aspects of work environments including health and safety measures, maternity/paternity leave, flexible work schedules, learning and development opportunities [9], as well as community engagement such as fund raising, sponsorships, or scholarships for local communities [10].
To achieve these goals, the governing tendency in a sustainable and low-carbon economy is the attentive use of resources and, to that end, the motivation for effective and efficient realization of recycling and remanufacturing processes of used products [11] that further create job opportunities for local communities and can be economically viable for companies complying with rules and regulations concerning the environment. In fact, companies implement various approaches associated with one or more product life cycle phases such as production planning integrated with remanufacturing control policies [12], product recovery [4], backward logistics integration building CLSCs [13], and overall carbon emission reduction [14,15]. In that sense, production and transportation activities are crucial targets in implementing sustainable operations.
is is done by reviewing and designing sustainable CLSCs [16]. is leads to implications for other decisions both on the strategic time horizon of 3-10 years and on the tactical encompassing 1-12 months encompassing the selection of remanufacturing/ production technologies, suppliers, facility locations, strategies for product recovery, and transportation networks. ese decisions affect the societal contribution of the company especially in terms of employment affecting local communities [17].
Uncertainty has a great impact on CLSC network design due to inclusion of reverse flows of products. is can be attributed to different factors including quality and quantity of the returning products [18]. Uncertainty in decision-making can be subdivided into two categories, i.e., systemic uncertainty and environmental uncertainty. e systemic uncertainty is associated with some processes, including production and distribution of goods. e environmental uncertainty is associated with the performance of different SC players, such as suppliers and manufacturers [19]. e more players are acting within a CLSC, the more complex and prone to knock-on effects it becomes [1]. Besides, unforeseen events affecting globally operating SCs can lead to disruptions within SCs, e.g., the outbreak of the Icelandic volcano Eyiafjallajökull in 2010 where air traffic in Europe came to an abrupt halt. Production lines had to stop due to shortage of parts from suppliers. e financial crisis in 2008 was another example, or the global pandemic due to the COVID-19 virus where social distancing and the shutdown of places where people could meet and get infected significantly affected whole business areas.
Companies try to minimize negative impacts caused by uncertainty and unforeseen events by implementing measurements that render their SC more resilient [20,21]. is can be done by introducing different kinds of flexibility that take into account overall SC structures [22], which needs to be carefully done [23]. Not surprisingly, research efforts within this topic have significantly increased in recent years. In fact, we find many concepts that aim at introducing actions for rendering SCs less prone to disruptive events. Perspectives and methods of introducing resilience (see also [24]) are discussed within the concepts of disruption of SCs, agility [25,26], fragility [27], antifragility [28], etc. Here, we define a resilient SC by its flexibility of reacting to unforeseen events. Resilience is incorporated in our model by regarding flow complexity and node complexity.
We do not find much work that combines aspects of resilience with sustainability considerations in CLSCs; see, for instance, [26], Mehrjerdi and Shafiee [29], and Mehrjerdi and Lofti [30]. e same is valid for mathematical models integrating these aspects.
We close this gap in providing a multiobjective programming model for design decisions in CLSCs including sustainability expressed by job opportunity creation in rural areas and resilience expressed by network flexibility. We consider uncertainty in input parameters by using a probabilistic programming method, i.e., the epsilon (ε)-constraint method that transforms algorithms for unconstrained into those for constrained problems. Here, we use it to convert a multiobjective into a single-objective model. Furthermore, we use Lagrangian relaxation to solve the resulting high-dimensional model.

Literature Review
CLSC management has received great attention in the last decade; see, e.g., the reviews of Fleischmann et al. [31], Angell and Klassen [32], Dekker et al. [33], Inderfurth et al. [34], Akcali et al. [35], Souza [36], Raza [37], and the references therein. Related research deals with the consideration, design, and integration of return flows and rework of used products including restoring, collecting, and repairing together with forward-related activities within the overall network [2]. is is also motivated by many countries' governmental rules and regulations that enforce manufacturing companies to reduce their wastages [38]. Nevertheless, retrieving valuable materials from recycled products also lead to profit gains as in the case of mobile phones that contain valuable metals such as copper, aluminum, steel, silver, gold, and cobalt, just to mention a few, where some are put on the list of critical raw materials by the EU commission [39].
A first reference that provides a review of quantitative models for reverse logistics is [31]. ey consider three main areas that deal with returns flows, i.e., production planning, distribution planning, and inventory control. Designing CLSCs include strategic decisions on where and how many facilities to open (or close). is significantly depends on decisions on the tactical level regarding the prospect interfacility forward and reverse flows. Akcali et al. [35] provide a review of the literature regarding models accounting for the strategic design including infrastructure decisions of either backward flows or for both forward and backward flows. A great portion of the literature seems to deal with multicriteria/objective formulations. According to Govindan et al. [16], 88% propose biobjective or multicriteria models for CLSCs.
As we are interested in the interaction of resilience and sustainability aspects in a CLSC subject to uncertainty, we review the specific literature that includes at least some of these mentioned features in proposed mathematical models for CLSC design decisions. For studies including only one of these features, we refer the reader to the abovementioned literature reviews.
2.1. Single-Objective Models. As pointed out earlier, uncertainty has a strong impact on the strategic and tactical decisions of the CLSC. Accordingly, we find many stochastic models for designing CLSCs. For instance, Inderfurth [40] presented a CLSC network based on a stochastic programming model. He considers products, demands, and return rates as random to address the uncertainty of products quality.
Ozceylan and Paksoy [41] propose a multiperiod multiproduct mix integer nonlinear program (MINLP) model to optimize strategic decisions for the CLSC specifically accounting for uncertainty in capacities, demand, reverse rates, and objectives, the latter aiming at minimizing total costs for transportation, purchasing, refurbishing, and fixed costs. In fact, reverse rates significantly complicate mathematical models [42]. eir proposed framework accounts for the logistic managers subjective and often imprecise, but desired levels for the uncertain aspects stated above by fuzzy theory. In the same realm, Alamdar et al. [43] examined the decisions uncertainty situations in CLSC. ey consider fuzzy price and sales demand within a game theoretic mathematical model to compare the performance of all facilities in the networks such as production centers composed by a collector and a retailer.
As stated earlier, sustainability addresses three pillars, one being environmental that encompasses all influencing factors of the environment, e.g., CO 2 emissions and overuse of resources. Reducing carbon emissions adds to preserving the environment and thus accounting for sustainability. Accordingly, Zohal and Soleimani [44] develop a multifunctional CLSC model for the gold industry with the objective to minimize carbon emissions besides overall SC costs.
ey apply an ant colony optimization solution method to discover near-optimal solutions. Furthermore, Soleimani et al. [45] propose a green CLSC chain for environmental considerations and lost business days. ey use a genetic algorithm to solve the problem.
Jabbarzadeh et al. [1] develop a robust-stochastic optimization model including a Lagrangian relaxation solution method for addressing the operational and disruption hazards in their CLSC. ey specifically account for supply and demand as well as whole facility disruptions and introduce aspects of fortifying their occurrence through investments. Based on a real case, Monte Carlo simulation is used to evaluate the performance of the proposed model. However, sustainability aspects are not considered in the study.

Biobjective or Multicriteria Models.
Altiparmak et al. [46] propose a multiobjective mixed-integer nonlinear programming (MINLP) model for a CLSC that minimizes total cost of the network and maximizes the rate of fulfilling the customer's demands at the available time as well as the capacity utilization balance for distribution centers. ey use simulated annealing to solve a real case, i.e., a plastic producing company in Turkey.
Amin and Zhang [47] present a three-stage stochastic multiobjective model for CLSC based on qualitative and quantitative criteria. e three stages include the evaluation of suppliers, remanufacturing subcontractors, and rework sites, the network configuration, and the order allocation. Fuzzy set theory is used to account for uncertainty in demand. e objectives are to maximize expected profits while minimizing purchasing costs and costs related to the refurbishing sites as well as remanufacturing subcontractors to determine the amounts of products to be disposed, reworked, and procured from suppliers. e results show that expected profits decrease with the increase in the uncertainty of demand.
Pishvaee et al. [48] develop a mixed-integer linear problem (MILP) fuzzy/possibilistic mathematical model for a CLSC with a biobjective function concerning a single period for a logistic network problem that concurrently minimizes overall costs including opening facilities, variable processing, and transportation in the first objective. In the second objective, the total environmental impact defined by an ecoindicator is minimized. ey further consider uncertainty in demand, rate of returned product, operating cost, shipping cost, and delivery time. e two-stage solution method transforms the model into an equivalent auxiliary crisp model and is then modified using the ε-constraint method to find compromise solutions. e model is tested on a real-life case involving medical equipment, i.e., syringes, that can be either recycled or incinerated in the reverse flow at collection/disassembly centers. e results show that the solution approach can find a reasonable tradeoff regarding the two objectives.
Talaei et al. [49] develop a fuzzy biobjective MILP mathematical model for a (green) CLSC network consisting of manufacturers, remanufacturers, collection/inspection centers, disposal centers, and customers that simultaneously minimizes overall network costs and environmental effects expressed by carbon emissions. ey use a robust fuzzy programming approach to account for uncertainty of variable costs and demand rates and further use the ε-constraint method to solve the model. e results show that the model can well balance network uncertainties and provide a robust price which is the cost allocated to the system to face uncertainty in a robust manner.
Zhu and Hu [50] use an adapted ε-constraint method to solve a multiobjective mathematical model for a three-level SC composed of suppliers, manufacturers, and customers that addresses all three pillars of sustainability in strategic decisions. It employs life cycle assessment to measure the environmental factor and sustainability reporting guidelines to quantify the social component. e model supports the decision-maker in selecting facility locations, technology, product manufacturers, and demand allocation by providing trade-off solutions under different environmental and social requirements.
Sreekumar and Rajmohan [51] deal with complexity of the evaluation and selection process of sustainability strategies in manufacturing SCs in a multicriteria setting. For that purpose, they combine analytical hierarchical processing (AHP), the subjective and objective weight integrated approach (SOWIA), and the "Technique for Order of Preference by Similarity to Ideal Solution" (TOPSIS) which is a method to order preferences by similarity to ideal solution. ey show that their approach can support decision makers in finding the right level of subjectivity as well as objectivity for their decisions by improving visibility of the problem and flexibility in its analysis. However, they neglect interdependency of factors.
Jiang et al. [52] develop an MILP model considering optimal decisions on partner selection, technology selection, transportation mode selection, material and product supply, as well as return flows with the objective to reduce carbon footprint. ey provide a sensitivity analysis focusing on the effects of variations in the price of the carbon footprint as well as the product recovery rate and show that small increases in the carbon price have little effects on the CLSC network design, but the increase in the recovery rate greatly affects the procurement and recycling processes. Similarly, Liu and Hu [53] investigate the interaction between carbon tax policies and supply chain coordination further including customer behavior. ey use game theory to model the problem. Game theory is also used by Ma [54] to study the effects of cooperation in the green CLSC. Sadeghi et al. [55] includes location-routing decisions in their proposed MILP model for designing a CLSC. Baghizadeh et al. [56] presented a multiperiod multiobjective mixed-integer nonlinear programming (MINLP) model to maximize the profit, minimize CO 2 emission, improve social aspects, and minimize the number of lost demands. Strategic and tactical decisions have been implemented throughout the supply chain for all facilities, suppliers, and machinery. ey applied epsilon constraint and Lagrangian relaxation to solve the mathematical model.
Johnson et al. [57] present an analytical framework that considers risk mitigation strategies for robustness and resilience of a SC using optimization and simulation to quantify the impact of disruptions and the benefits from the application of risk mitigation strategies. ey provide a good overview of different strategies to gain resilience in SCs and conclude that the introduction of flexibility is best to handle disruptions and provide SC resilience.
Zahiri et al. [58] propose a stochastic multicriteria MILP for the pharmaceutical SC that encompasses aspects of sustainability as well as resilience including strategic and tactical decisions. Strategic decisions include where to open locations in the network, and tactical decisions include how much material to order, how much products to produce, put on stock, and sell based on certain technologies, how much CO 2 to emit as well as how to handle capacities, allocation of products, and their transportation. e pillars of sustainability are regarded by the networks' total costs (economical), emissions of CO 2 caused by business activities (environmental), and the social impact of business activities (social) accounted for by job creation opportunities/unemployment rate and the regional development rate where more credit is given to underdeveloped areas. Resilience is considered by introducing (1) node criticality where a node is defined as critical if total inflows and outflows exceed a certain threshold, (2) sets of technology where the utilization of new technology may be more expensive, but also more reliable, which is especially important when a certain technology is not available causing disruptions. In that case, backup technologies can ease the situation, (3) flow complexity which measures the interaction between nodes based on assignments of activities to the different neighboring center, (4) node complexity referring to the total number of active nodes, and (5) customer deservice level that expresses unmet demand as a weight function. Our approach is inspired by the modeling of node and flow complexity which we adopt for our mathematical model formulation. As Zahiri et al. [58] do not consider return flows, we do not discuss their approach in further details, but the formulations for node and flow complexity to be used to introduce resilience to our model are interesting; further details are given in Section 3.2.
Shabbir et al. [26] propose a two-stage stochastic multicriteria MILP for a CLSC network design problem including sustainability and resilience aspects. ey introduce for the first time a conditional value at risk to find robust solutions accounting for the underlying uncertainty. With respect to sustainability, their model includes all three pillars by minimizing overall network cost, energy consumption, and CO 2 emissions as well as maximizing job opportunities. Uncertainty in demand is considered, and measurements of resilience concerning capacity are taken against disruption risks due to demand variation.
Similarly, Mehrjerdi and Shafiee [29] analyze the interaction between resilience and sustainability aspects in a CLSC further examining in influence of information sharing and multiple sourcing strategies.
eir results show the benefits of an integrated view of these two concepts for decisions on the CLSC design, see also [30] for a similar approach.
Yavari and Zaker [59] develop a two-layer biobjective network MILP model for a resilient CLSC with sustainability aspects for a perishable product with a limited lifetime. Specifically, effects of electric power network disruptions are regarded. e objective is to minimize expected total network costs as well as expected total amounts of CO 2 emissions (Table 1).
Abbreviations used in the columns in Table 2: As can be seen by the literature review, resilience in CLSCs has not received enough attention of the research community. According to the importance of this subject, this study integrates flexibility measures in a CLSC. For that purpose, we develop and analyze an MINLP that includes  sustainability aspects as well as flexibility measures. We further consider the dwell time and a Jackson network simultaneously for reducing customer waiting times. is is a new aspect that has not been integrated so far into CLSCs with flexibility and sustainability aspects. We further provide numerical examples and an implementation of a case study to evaluate the proposed model. e remainder of the paper is given as follows: in Section 3, we give basics of the SC setting and our model formulation. Section 4 presents computational results. Section 5 gives the case study, while Section 6 provides a discussion. Section 7 concludes the paper. e preprint version of this research is uploaded earlier [56]. Figure 1 shows a schematic presentation of a CLSC. It shows the forward flows and the reverse flows. Forward flows start at the manufacturing centers routed through the distribution centers (DCs) to customers, whereas reverse flows start at the customers routed and operated in collection and recycling centers. In the latter, necessary disassembly, refurbishing, and repair are performed to retransform them in a "as good as new" state.

CLSC Setting and Mathematical Model
Our proposed model integrates both strategic decisions, i.e., the optimal selection of facility locations composed of DCs, recycling centers, and collection centers, and tactical decisions, i.e., production, inventory planning, capacity management, and system of transportation. Customer and factory locations are assumed as fixed.

Model Denotation.
e following notation is used to present the mathematical model: Indices: i: Index for fixed location of a factory j: Index for potential location for DCs k: Index for fix location for customer zone l: Index for potential location for collection centers m" Index for potential location for recovery centers n: Index for type of transportation mode Parameters: D kj : Demand of customer in zone k that goes to distribution center j B k : Fraction rate of return products from customer zone k fcd j : Fixed cost for opening DC j fcc l : Fixed cost for opening collection center l fcr m : Fixed cost for opening recovery center m pcc l : Processing cost at collection center l pcd j : Processing cost at distribution center j pcr m : Processing cost at recovery center j co ij : Transportation cost from plant i to distribution center j cu jk : Transportation cost from distribution j to customer zone k cq kl : Transportation cost from customer zone k to collection center l cp lm : Transportation cost from collection center l to recover center m ch mj : Transportation cost from recover center m to distribution center j px j : Maximum capacity of distribution center j pz m : Maximum capacity of recovery center m py l : Maximum capacity of collection center l pp i : Maximum capacity of plant j i : Number of created job opportunities out of opening distribution center j j l : Number of created job opportunities out of opening collection center l j m : Number of created job opportunities out of opening recovery center m ev j : Regional economic value of distribution center j ev l : Regional economic value of collection center l ev m : Regional economic value of recovery center m rd j : Regional development factor in region j rd l : Regional development factor in region l rd m : Regional development factor in region m ei ijn : Environmental impact of shipping one unit product from node i to j using transportation mode n ei jkn : Environmental impact of shipping one unit product from node j to k using transportation mode n ei kln : Environmental impact of shipping one unit product from node k to l using transportation mode n ei lmn : Environmental impact of shipping one unit product from node l to m using transportation mode n ei mjn : Environmental impact of shipping one unit product from node m to j using transportation mode n ur j : Unemployment rate in region j ur l : Unemployment rate in region l ur m : Unemployment rate in region m α 1 : Penalty coefficient for flow complexity between nodes i and j α 2 : Penalty coefficient for flow complexity between nodes j and k α 3 : Penalty coefficient for flow complexity between nodes k and m α 4 : Penalty coefficient for flow complexity between nodes m and j β 1 : Penalty coefficient for node complexity of distribution centers 6 Mathematical Problems in Engineering o ijn : Number of products transferred from plant i to distribution center j with transportation mode n u jkn : Number of products transferred from distribution center j to customer zone k with transportation mode n q kl : Number of products transferred from customer zone k to collection center l p lmn : Number of products transferred from collection center l to recovery center m with transportation mode n h mjn : Number of products transferred from recovery center m to distribution center j with transportation mode n x j : 1, if a distribution center is opened at location j; 0, otherwise y 1 : 1, if a collection center is opened at location l; 0, otherwise z m : 1, if a recovery center is opened at location m; 0, otherwise w i ′ : 1, if a plant center's node is critical; 0, otherwise w j ′ : 1, if a distribution center's node is critical; 0, otherwise w l ′ : 1, if a collection center's node is critical; 0, otherwise w m ′ : 1, if a recovery center's node is critical; 0, otherwise v ij : 1, if node i is assigned to node j; 0, otherwise v jk : 1, if node j is assigned to node k; 0, otherwise v lm : 1, if node l is assigned to node m; 0, otherwise v mj : 1, if node m is assigned to node j; 0, otherwise

Measures of Sustainability.
e environmental pillar of sustainability is considered as the reduction of harmful emissions emitted to the environment, while social aspects and impacts are regarded by the amount of job opportunities created by the opening of new facilities for distribution, collection, and recycling centers which leads to a reduction of unemployment rates. is also concerns the regional development where the focus is on less developed areas.

Measures of Resilience/Flexibility.
Network flexibility of the CLSC is integrated in the model by considering the number of notes and paths between nodes in the CLSC design selection. For instance, if the number of nodes and paths between nodes increases, so does the nodes/paths' complexity, and thus the SC's flexibility is reduced. Less flexibility leads the network to lose the ability to cope with change. us, flow and node complexity are penalized.
Given the statements in the previous section, the resilience of a network is a concept that has been less embraced in SCs, particularly in CLSCs. Hence, the proposed mathematical model bases on equations proposed by Zahiri et al. [58] to introduce resilience.  Mathematical Problems in Engineering 7 centers to the production centers. e second part deals with the relationships between the demand zone and the distribution centers. e third part allocates collection center to demand zones. e fourth part allocates collection centers to recovery centers, and finally, the last term allocates recovery centers to distribution centers [58].
According to equation (1), if overall number of connections is high, the total network flow is deemed to be utterly complex. (2) estimates the number of active nodes by calculating the number of established facilities such as factories, distribution centers, collection centers, and recovery centers. As the number of nodes increases, the communication paths between different levels of the supply chain increase sharply, which increases the complexity of the network. As mentioned in the previous section, if the number of active nodes is large, the complexity is considered as "high."

Mathematical Model
Formulation. e proposed mathematical model has the following underlying assumptions: (1) No disposal of products is considered; thus, all products that are collected are recycled (2) All customer demand must be met, and all returned products must be collected (3) ere is a pressure mechanism in the forward flow (4) ere is a tensile mechanism in the reverse flow (5) e location of factories and customer locations are predetermined and fixed e research hypothesis regarding our study is as follows: How can the proposed model and solution method support in making economically sound decisions regarding the opening of centers for distribution, collection, and recycling while considering aspects of sustainability and flexibility for a CLSC? e mathematical model is given as follows, starting with the objective functions composed of four functions where the first one (w 1 ) minimizes all costs, the second one (w 2 ) improves social aspects, the third one (w 3 ) minimizes air pollution, and the fourth one (w 4 ) improves flexibility and network resilience. e first objective function (w 1 ) tries to reduce all costs of the closed-loop supply chain.
Transportation cost, operational cost, and establishment cost of facilities are considered in the first objective. Maximizing employment rate and number of created job opportunities will be achieved by minimizing w 2 . e third objective tries to reduce the gas emission of facilities and transportation. e fourth objective function maximizes the resilience of the whole network by decreasing the complexity of supply chain.
j n u jkn ≥ l q kl · ∀k, k q kl ≤ y l py l · ∀l, l n p lmn ≤ z m pz m · ∀m, Constraint (4) ensures that customer demand is fully met. Constraint (5) states that the number of products delivered to the collection centers by customers is not greater than the number distributed to customers. It also indicates that all the return products have been collected from the customer. Constraints (6) and (7) balance the flow of products from the factory to customer area and from the costumer area to collection center. Constraint (8) indicates that all returned products are delivered to the recycling center. Constraint (9) requires that there is no waste in the recycling centers. Constraints (10)-(13) present the capacity constraint existing in the production, distribution, collection, and recycling centers. Constraints (14)-(17) are allocation constraints. ey ensure that communication between each pair of nodes only establish when the there is a flow in the given link. Finally, constraints (18)-(21) are noncritical conditions for production, distribution, collection, and recycling centers.
We use fuzzy mathematical programming (FMP) to properly manage two types of uncertainty including the inherent uncertainty of data and flexibility in the objective function and the elastic uncertainty of the constraints [60]. In the first category, the uncertain or indefinite coefficients in the objective function and related constraints are generally managed based on the available quantitative inputs and the qualitative knowledge of the decision-makers. e second category, which is flexible, is about the decision-making process under the flexible value related to the function performance and constraint elasticity [61]. Possibilistic chance-constrained programming (PCCP), one of the known methods of fuzzy programming, is used to deal with the potential data, such as the constraints that include the potential data on the left/right side (RHS) [62]. e given method provides the minimum level of confidence (Ψ) for DM to meet possible constraints. DM can determine the minimum level of confidence (Ψ) as the secure margin for the level of satisfaction of any possible chance constraint. is is considered in PCCP using the probability criterion (Pos) which represents the maximum (most optimistic) level of probability of data occurrence or the necessity (Nec) that considers the least (most optimistic) level of probability of data occurrence [63].
In real situations, DM attitudes are generally different and may differ between optimistic and pessimistic extremism. In this regard, Liu and Liu [64] proposed the criterion of mutual trust (Cr), the perspective of compromise (i.e., neutral) towards both extremes. For further explanation, Φ is considered a trapezoidal fuzzy variable. It is possible to find the apparent convergence of these measurements in the leading equations for all confidence levels:

Mathematical Problems in Engineering
e mathematical expectation value of the fuzzy variable Φ is explained in the condition of the following equation: Although reliability measurement (Cr) can be considered a relatively appropriate method to avoid the most extreme and most optimistic choices of potential data occurrence, it forces DMs to provide a flexible approach by providing a specific point between two extreme views. To cover this gap and to address inaccurate parameters of the problem, Xu and Zhou [65] proposed an advanced Me Measure method that can show optimistic and pessimistic attitudes for various decision-makers in the fuzzy environment.
e DM can choose extreme attitudes using a convex combination of any point in the range: In the above equation, c that has a value between 0 and 1 is the optimistic-pessimistic parameter and will be determined by DM's attitude and opinion. As stated before, by using the same method for Me Measure, we have e value of the mathematical expectation of the fuzzy variable for the values c 1 ≥ 0 are DMs prefer to look at the problem from the pessimistic point of view in most areas with high sensitivity to nonlinear parameters, including designing supply chain networks. Although it was previously stated that ψwould naturally take a value in the [0,1] interval, in such sensitive cases, this parameter would be less than 0.5. erefore, in this study, we obtain According to the given relationships and equations so far, we can write the possibilistic chance-constrained program (PCCP) model as follows: In the above basic model, the matrices B, C, A, G, and O are constraint coefficient matrices of the problem. Also, the vector Q is the coefficient of fixed values of the problem. Furthermore, Y and L are binary variables and N is a large value. X is also a positivity continuous variable of problem. Q, B, and d are capacity of centers, fixed costs of establishing centers, and demand in customer areas, respectively. e expected value operator is utilized for converting the probabilistic objective function to its crisp equivalent to form the PCCP model.
In this problem, a variable that follows a trapezoidal probability distribution is used to model the uncertainty parameters, which are defined by four prominent points in Figure 2. For example, Τ � χ 1 , χ 2 , χ 3 , χ 4 .
According to the given description, we have e objective function of the model in addition to the first and second constraints contains uncertainties that follow the probable trapezoidal distribution determined by experts or data history. rough above transformations, the crisp equivalent of the BPCCP model based on Me Measure and the expected value operator is given as follows: According to the inherent uncertainty of the parameters in the supply chain, dealing with them is of great importance. Demands for manufactured products were always uncertain, and frequent changes occur in them. is uncertainty poses many challenges for modeling and finding desirable solutions. To address and manage this, the model of this paper considers the demand, product return rate, operating costs at production centers, collection, and recycling, and the number of pollutant emissions emitted from the production and recycling center in an uncertain way. e final model was rewritten as follows by considering these uncertainties. e first and third objective functions are modified as follows: e resulting model is a multiobjective programming model with priorities such as cost, social dimensions, network flexibility, and environmental impacts. One of the popular ways to solve multiobjective problems is the ε-constraint method. is method has been introduced by Haimes [66]. It involves prioritizing the objective functions, Mathematical Problems in Engineering so that the objective function with the highest importance in the problem is considered as the main objective function and the other objective functions are added as constraints to the main model. In fact, by prioritizing the cost functions as the main function, the decision-maker can measure the impact of other functions on the problem. e proper application of the ε-constraint method requires determining all objectivefunction ranges, leastwise for the p-1 objective functions, which will be applied in the form of constraints. It would be challenging to calculate objective-function ranges on the efficient set. e best objective value can be effortlessly attained as the optimum value of the individual optimization; however, attaining the worst objective value over the efficient set (i.e., the nadir value) is difficult. To this aim, these ranges are frequently calculated from the payoff table, whose results are obtained from the individual optimization of the p objective functions. e approximation of the nadir value is generally carried out with the minimum value of the respective column.
Nevertheless, it should be ensured even in this context that the solutions obtained from the individual optimization of the objective functions are Pareto optimal solutions. While there exist alternative optimal solutions achieved using commercial software, the optimal solution obtained is not of a guaranteed Pareto optimal type. is ambiguity can be resolved by conducting lexicographic optimization for all the objective functions to build a payoff table with Pareto optimal solutions only.
e reservation values (which function as a lower bound or an upper bound for minimization objective functions) must be defined for all the objective functions to easily circumvent the challenges raised in estimating their nadir values. ose values that are worse than the aforesaid reservation values should be excluded.
One of the innovative additions to the algorithm is breaking out of the nested loop early when the given problem is rendered infeasible. is can dramatically speed up the algorithm when there are more than three objective functions. A more relaxed version of the constrained objective functions is initially used in this algorithm, and then bounds are gradually restricted. In maximization objective functions, the minimum is first dealt with and the right-hand side (RHS) value of the given constraint is then slowly increased; however, the reverse is true for minimization objective functions. Here, the respective objective function should not be restricted anymore when the given problem is rendered infeasible because it may lead to infeasible solutions. In this case, the algorithm will exit the innermost loop and engage in the next waiting grid point of the previous objective function corresponding to the outer loop.
In this paper, we use the improved ε-constraint method to transform the four-objective model into a single-objective model. is is done by converting all the objective functions except one to constraints. e difference between the model used in this study and the original ε-constraint method is that in the original method, all the primary objective functions that are considered as new constraints in the model are divided by equal intervals. For this reason, the model is solved using all right-valued combinations.
Previously, a model of the ε-constraint solution method is presented in which S, x, and M (x) are considered as the space of the answers of the problems, the decision vector of the problem, and the objective functions, respectively. Furthermore, equation (36) indicates the constraints of the corresponding objective functions which are the values on the right side of the constraints of the objective function with epsilon being a small value, thus a deficiency variable. e abovementioned model ensures that the obtained optimal solution is an efficient one only if all the constraints of the objective function are binding. e model presented in the above section is a large-scale MINLP model that can be solved with commercial software such as GAMS. However, since increasing the problem size of the solution would increase the diminutions of the problem sharply, it is not possible to solve large polynomial problems using standard methods and it is recommended to use more efficient and optimal approaches [67]. Hence, we use the Lagrangian relaxation method to solve the integrated optimization model. e Lagrangian relaxation method is known as one of the best approaches to solve supply chain problems aiming at obtaining robust solutions. is method is able to provide an upper and lower bound for the optimal value of the objective function, which in turn helps the quality of their solution method and to find out how far the possible solutions are from the optimal solution of the problem [68].
In this study, the Lagrangian relaxation consists of three main steps, which are as follows: in Step 1, we attempt to obtain a lower bound for the optimal solution. In Step 2, an upper bound is obtained for the optimal solution. Finally, in Step 3, the upper and lower bounds obtained from the previous steps are not close enough and the values of these upper and lower bounds should be updated. ese steps are repeated until the lower and upper bounds reach the specified limit.

Discovering the Lower Bound.
In order to achieve a lower bound of the problem, certain constraints of the problem have been relaxed in order to promote the solution of the problem even though the answer cannot be sought [68]. In this problem, the problem solving was facilitated through relaxation of the constraints (10) and (13).
According to the conducted relaxations, the first objective function changes as follows: In the abovementioned objective function, the values U 2 lmn and U 1 ijn are Lagrangian coefficients with nonnegative values. Lagrangian coefficients with fixed values were used to minimize the abovementioned relationship. Hence, the optimal value obtained from the Lagrange dual model is the lower bond of the problem.

Discovering the Upper Bound.
According to the relaxation of constraints (10) and (13), in most cases, the solution obtained for the dual Lagrangian is not feasible. Hence, the feasible solution is obtained by solving the model based on equations (38) and (39) and applying the constraint equations (7)-(20), (22)- (27), (34), and (40) and (41). is is the upper bound of the problem.
Updating upper and lower bounds, In each repetition of the Lagrangian method, the Lagrange coefficients are updated with new values, and new values are generated continuously for the upper and lower bounds. In each repetition of problem, the Lagrange coefficient is computed as follows, which was introduced by Fisher [68].
In equations (38) and (39), the value of x is the number of repetitions stp 1,x and stp 2,x which are calculated as follows: In equations (40) and (41), the values UP and LB x are the best discovered upper bound and obtained lower bound in the stp repetition, respectively. ] � 2 was given at the beginning of the solution method, if the value of LB is not improved in the five repetitions, then the value of new ] equals to the half of the previous ]. e problem-solving process ends with two conditions: (1) Obtaining a feasible solution with the appropriate tolerances (2) e value of stp has reached the minimum possible value

Computational Study for Solution Method Evaluation
In this section, 20 numerical examples with different sizes are presented to validate the model and its solution approach. In a second step, the model is solved using the Lagrangian relaxation method. e results including "objective function," "solution time," and "difference between objective functions" are presented in Table 2 and Figures 1-5.
We assume uniform distributions for all input parameters but with differing ranges. As can be easily seen in Table 2, by increasing the size of the problem and the complexity of numerical computations, a significant difference would be created in the solving time of the model with and without Lagrangian relaxation. is shows how much faster the problem is solved through Lagrangian relaxation which also indicates the effectiveness of this method. Figure 3 shows the values of the objective function obtained from solving methods with and without Lagrangian relaxation. As the dimensions of the problem increase, the objective function significantly increases in the successive repetitions. On the other hand, as shown in the chart, the difference in the objective function values of these two methods is negligible. Figure 4 also shows the difference percentage of the objective function values when solving the problem with and without Lagrangian relaxation. As can be easily seen in the smaller dimensions, this value is small and insignificant, and in the larger dimensions, the problem is acceptable.
In the same way, Figure 5 shows the solving time of the model with and without Lagrangian relaxation algorithm. Clearly, there is no significant difference between the solution finding times of these two methods when the size of the problem is small. However, in the medium and large dimensions, the Lagrangian relaxation method requires much less time for solving the problem. However, the solution finding time of the method without Lagrangian relaxation is significantly increased. erefore, it can be Mathematical Problems in Engineering 13 concluded that this solution method is more efficient and cost-effective than epsilon constraint.

Case Study
We aim at demonstrating the performance and benefits of the suggested model and solution method on a real case which is a tire factory in Iran. e focus points of this industry are high quality, reliability of availability, and costeffectiveness of resources that need to be assured by the CLSC. e design of the CLSC aims at reducing maximum network costs while meeting customer demand, thus optimal locations for DCs, collection centers, and recycling centers need to be found. Besides aiming at increasing the productivity of factories while reducing (network) costs, the company considers social issues such as job creation and unemployment reduction in less favored areas with higher unemployment rates. Because of the high pollution of the industry, the company aims at reducing the environmental damage caused by production and transportation via emissions of environmental harmful gasses. us, to reduce the volume of carbon dioxide produced by transportation, two types of transport systems have been proposed for selection, where one of them has lower transporting capacity and produces less carbon dioxide and the other one is employed for the high capacity and produce more carbon dioxide. e case company's headquarters is in Tehran. e company produces one type of tire in five different cities of Iran to meet the needs and requirements of 11 customer points. e company intends to select 12 centers out of 15 potential distribution points to establish their DCs. Because of the high costs caused by the production process and the raw material procurement where the previous is time consuming, the company intends to open seven collection centers in 10 potential points. To recycle the collected products, a budget equal to the budget needed for establishing five recycling centers in eight given cities is allocated. Accordingly, three additional constraints are added to the original proposed model: (42) Figure 6 shows the map of Iran in which the fixed location of production facilities, the demand points, and the potential locations for establishing distribution, collection, and recycling facilities are depicted. It should be noted that due to the distances between the facilities and the different transportation systems, the cost of transport between the facilities varies based on two varying parameters.
e model of the case study was solved using the solution method that was described in detail in previous sections. Figure 7 shows the selected optimal locations for the establishment of distribution, collection, and recycling facilities. Furthermore, the optimal flows are shown from the factory to distribution centers, from distribution centers to customers, the returns of second-hand products from customers to recycling centers, and the transferring of recycled products to distribution centers. e optimal solution is obtained by solving the model with the first objective function. All costs of the supply chain network are composed of transportation cost, operational cost, and establishment cost of facilities and result in 176521007.63695 toman (local currency). e optimum value of second objective function which contains job opportunity and unemployment rate is -127426.526138657 units. e best value of CO 2 emission for the third objective function is 940431.149676294 litter. e fourth objective function responses to enhancing resilience of the whole network by decreasing node and flow complexity that gains a value of 1267.45788718017 unit.

Sensitivity Analysis and Discussion
To further examine different dimensions of the models and to understand how the problem parameters are affected, a sensitivity analysis has been conducted on the performance values of the objective functions and various parameters as follows.
6.1. e Capacity of the Distribution Facility. We first analyze the effects of facility changes selected in the optimal solution based on the objective functions. e selection of facilities influences the capacity available at the DCs. e capacity of each facility is increased to 5,000 units per each test to determine the effects in a clear manner.
As shown in Figures 8 and 9, the number of selected facilities gradually decreases by increasing the capacity of DCs. It is easy to see that, by increasing the capacity up to   77 percent, the number of facilities reduces from 12 centers to 5 centers. On the other hand, after this reduction, the raising slope of the cost becomes flatter which is due to the reduction of DCs establishment costs. However, as can be seen in Figure 10, the employment rate declines because of reducing the number of established DCs which, in turn, cause increased social impact in terms of unemployment rates that rise by 11 percent compared to the initial state. Moreover, we can see in Figure 11 that these changes also increase CO 2 emissions. ese managerial insights give great support to decisionmakers in showing the effects of decisions to the overall CLSC.

e Rate of Product Return.
In this section, the relationship between the rate of the product return with the number of associated facilities and model functions is investigated. As can be seen in the diagram in Figure 12, by increasing the confidence level of the returned product rate, more facilities are needed to refurbish returned products. Hence, it is logical to increase both the collection facilities and the recycling facilities. Because of the stability of the facility, small changes, and lower costs of recycling, overall costs decrease significantly. However, if the number of these facilities increases in the objective function, the cost increases as well, and the more of these facilities are increased, the higher the slope of the costs would be. As also  shown in Figure 12, if the confidence level is increased by 5% at the beginning, then the numbers of collecting and recycling centers are 7 and 5, respectively, which shows 35% increase in the initial values. If these values are 10 and 9, then the increasing rate of 40% and 80% can be seen due to the number of these facilities. As shown in Figures 13  and 14, in the cost function, the overall cost experiences a 10 percent reduction when the confidence level is increased up to 20%. If the confidence level is further increased, then this value drops to 6%. e amount of environmental impact has also increased slightly but experienced a 2% increase when the confidence level is increased to 35%. On the other hand, by Figure 15, if these facilities are increased due to the need for more labor work, then the second function is improved significantly, indicating the increased job opportunities and reduced environmental unemployment. us, if the confidence level in increased up to 35%, then the social effects of the model would experience an 8% reduction.

Resilience Analysis.
To explore the details and impact of various parameters, we show the changes in the critical parameters of the problem and their impacts simultaneously  in Figure 16 to obtain a better understanding of the impacts of variety of changes on the flexibility of a CLSC.
As can be seen from the figure, by increasing the capacity level in the production, distribution, and the collecting centers, the flexibility is improved significantly and the complexity of the CLSC network is reduced. Accordingly, a 25% change in capacity of the collection center and a 7% improvement in flexibility can be observed. Furthermore, when the capacity of the DC is increased up to 15%, the best possible solution for the flexibility is obtained when it has a 6% reduction.
It can also be observed that an increase in the confidence level of the product returning rate greatly overshadows the flexibility of the network. e reason for this is the need to increase the nodes of distribution, collection, and recycling centers and to make decisions about the way products are transported. More precisely, a 25% increase in confidence level results in a 19% worsening of flexibility. Furthermore,         increasing the confidence level can have a similar effect on flexibility, and a negative effect of 18.5% on flexibility had a similar effect on the product returning rate.

Conclusion
is work presents an MINLP multiobjective mathematical model designed to support selecting the best possible choice for distribution, recycling, and collection centers, so that all environmental costs and impacts are minimized while social impact and network flexibility are maximized. A CLSC is developed considering resilience to determine the optimal volume of product transfer, the number of nodes, and the state and complexity of the network. Also, uncertainty factors on demand, product return rate, and costs of distribution, collecting, and recycling centers and the amount of CO 2 emissions in processes related to manufacturing and recycling centers are considered. To deal with the multiobjective function defined in the model, the improved ε-constraint method is used which performs better than the conventional ε-constraint method. Afterward, a Lagrangian relaxation method is applied. A real case study was used to evaluate the performance of the proposed model showing that the model has high credit in both the reliability and performance and can be quite effective for decision-makers and researchers. Sensitivity analysis is performed on the main parameters of the model to clarify their effect on the model. ere are some unexplored aspects for future research that can help researchers to expand the current model. For instance, the inclusion of other new methods for coping with problem complexity and the uncertain nature of input parameters are considered as fruitful directions. For instance, a scenario tree can be designed for that purpose by stochastic programming to compare the efficiency of results. Regarding complexity, developing other exact or metaheuristic algorithms to improve the solution methods' calculation time could also be pursuit to render the model more applicable for larger real-life cases than the presented case study. Notably, although this paper has considered resiliency throughout the network, considering disruptions is a trend for future research. Considering both the disruption of DCs and the routes between parties can give valuable decision support for decision-makers.

Data Availability
Data and optimal value of variables are available from the corresponding author upon reasonable request.

Disclosure
is article is available as a preprint at Arxiv.org: arXiv: 2104.06560v1.

Conflicts of Interest
e authors declare no conflicts of interest.