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

Authors: Komeyl Baghizadeh; Julia Pahl; Guiping Hu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

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 find optimal locations for distribution centers, second-hand product collection centers, and recycling centers under the uncertainty situation alongside the factory’s fixed points. The purpose of the presented model is to minimize overall network costs including processing, establishing, and transportation of products and return flows as well as environmental impacts while maximizing social scales and network flexibility according to the presence of uncertainty parameters in the problem. To solve the proposed model with fuzzy uncertainty, first, the improved epsilon (ε)-constraints approach is used to transform a multiobjective to a single-objective problem. Afterward, the Lagrangian relaxation approach is applied to effectively 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 effects of important parameters on the optimal solution.


---

## Body

## 1. 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]. There 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]. They 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. The 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. This is done by reviewing and designing sustainable CLSCs [16]. This 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. These 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. This 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. The systemic uncertainty is associated with some processes, including production and distribution of goods. The environmental uncertainty is associated with the performance of different SC players, such as suppliers and manufacturers [19]. The 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. The 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]. This 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]. The 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.


## 2. 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]. This 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]. They 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). This 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.Özceylan 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]. Their 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. They 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., CO2 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. They 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. They 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. They 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.


### 2.2. 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. They 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. The 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. The 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. The 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. They further consider uncertainty in demand, rate of returned product, operating cost, shipping cost, and delivery time. The 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. The 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. The results show that the solution approach can find a reasonable trade-off 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. They 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. The 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. The 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. They 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. They 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. They 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 CO2 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. They 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. They 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 CO2 to emit as well as how to handle capacities, allocation of products, and their transportation. The pillars of sustainability are regarded by the networks’ total costs (economical), emissions of CO2 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. They 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 CO2 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. Their 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. The objective is to minimize expected total network costs as well as expected total amounts of CO2 emissions (Table 1).Table 1
Overview of related work including features of CLSC models with sustainability and resilience measures.
AuthorsObjectiveModel typeUncertaintySustainability featuresResilience/flexibility inclusionSolution methodApplication areaMehrjerdi and Shafiee [29]MMILPFTOPSISThree pillarsInformation sharing and multiple sourcingAugmented ECM (AUGMECON2)Tire industryMehrjerdi and Lofti [30]MMILPCVaRThree pillarsCapacities dependent on scenario, probability of disruption at center(i) LP-metric optimization(ii) SACar manufacturerShabbir et al. [26]MMILPCVaRThree pillarsManage demand uncertainty(i) LR(ii) Fix-and-optimize heuristic(iii) SAAutomobile assembler industryTalaei et al. [49]BFuzzy MILPFTCO2 emissions—ECMElectronicsYavari and Zaker [59]BMILPSTCO2 emissionsElectronic power system disruptionsLP-metric methodPerishable itemsZhu and Hu [50]MMILPFT(i) Three pillars(ii) Life cycle assessmentSustainability reporting guidelinesECMLRAUGMECONNumerical studyThis paperMMINLPFTJob creation depending on location, e.g., rural areas(i) Flow complexity(ii) Node complexity(i) ECM(ii) LRTire industryAbbreviations used in the columns in Table2:(i)
Objective:S = single criteria; B = bicriteria; M = multicriteria(ii)
Model type: GT = game theory; FP = fuzzy programming; SO = stochastic optimization(iii)
Uncertainty: RN = random numbers; FT = fuzzy theory; CVaR = conditional value at risk; ST = scenario technique; FTOPSIS = fuzzy TOPSIS(iv)
Sustainability features: three pillars include minimum network costs, minimum energy consumption/CO2 emissions, and maximun job creation(v)
Solution method: LR = Lagrangian relaxation; ECM = epsilon-constraint method; SA = sensitivity analysis; SAN = simulated annealing; MCS = Monte Carlo simulation; ACO = ant colony optimization; RO = robust optimizationTable 2
The results of the developed epsilon-constraint method and Lagrangian relaxation.
Problem numberIndicesGAMSLROptimality gapijklmnRun timeObjective valueRun timeObjective value11122212.38.13E + 074.18.13E + 070.000296%21222215.58.18E + 077.28.18E + 070.000169%31232218.982067192.819.38.21E + 070.000393%423322113.682461487.5914.38.25E + 070.000342%523422119.182970169.4218.78.30E + 070.000333%62442212680770169.42228.08E + 070.000288%724522134.382770169.4227.28.28E + 070.000400%834533149.686426698.3724.78.64E + 070.000221%935533157.688269811.2921.58.83E + 070.000370%103554326789217186.9423.98.92E + 070.000035%1136544273.493630025.62269.36E + 070.000295%1237644298.992473417.9729.59.25E + 070.000397%13489542146.895678997.5733.49.57E + 070.000386%144810542194.296060458.8436.59.61E + 070.000339%1541010542220102995413.840.71.03E + 080.000233%1641011542298.4118700925.644.81.19E + 080.002866%1741111542409.5122439701.159.41.22E + 080.029194%1841111642511.412021169774.61.20E + 080.018774%1941111652699.4128811065.986.61.29E + 080.013647%2041211652940151835127.6107.61.52E + 080.026376%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. This 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.The remainder of the paper is given as follows: in Section3, 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. The preprint version of this research is uploaded earlier [56].


## 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.Özceylan 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]. Their 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. They 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., CO2 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. They 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. They 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. They 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.


## 2.2. 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. They 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. The 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. The 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. The 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. They further consider uncertainty in demand, rate of returned product, operating cost, shipping cost, and delivery time. The 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. The 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. The results show that the solution approach can find a reasonable trade-off 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. They 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. The 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. The 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. They 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. They 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. They 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 CO2 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. They 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. They 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 CO2 to emit as well as how to handle capacities, allocation of products, and their transportation. The pillars of sustainability are regarded by the networks’ total costs (economical), emissions of CO2 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. They 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 CO2 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. Their 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. The objective is to minimize expected total network costs as well as expected total amounts of CO2 emissions (Table 1).Table 1
Overview of related work including features of CLSC models with sustainability and resilience measures.
AuthorsObjectiveModel typeUncertaintySustainability featuresResilience/flexibility inclusionSolution methodApplication areaMehrjerdi and Shafiee [29]MMILPFTOPSISThree pillarsInformation sharing and multiple sourcingAugmented ECM (AUGMECON2)Tire industryMehrjerdi and Lofti [30]MMILPCVaRThree pillarsCapacities dependent on scenario, probability of disruption at center(i) LP-metric optimization(ii) SACar manufacturerShabbir et al. [26]MMILPCVaRThree pillarsManage demand uncertainty(i) LR(ii) Fix-and-optimize heuristic(iii) SAAutomobile assembler industryTalaei et al. [49]BFuzzy MILPFTCO2 emissions—ECMElectronicsYavari and Zaker [59]BMILPSTCO2 emissionsElectronic power system disruptionsLP-metric methodPerishable itemsZhu and Hu [50]MMILPFT(i) Three pillars(ii) Life cycle assessmentSustainability reporting guidelinesECMLRAUGMECONNumerical studyThis paperMMINLPFTJob creation depending on location, e.g., rural areas(i) Flow complexity(ii) Node complexity(i) ECM(ii) LRTire industryAbbreviations used in the columns in Table2:(i)
Objective:S = single criteria; B = bicriteria; M = multicriteria(ii)
Model type: GT = game theory; FP = fuzzy programming; SO = stochastic optimization(iii)
Uncertainty: RN = random numbers; FT = fuzzy theory; CVaR = conditional value at risk; ST = scenario technique; FTOPSIS = fuzzy TOPSIS(iv)
Sustainability features: three pillars include minimum network costs, minimum energy consumption/CO2 emissions, and maximun job creation(v)
Solution method: LR = Lagrangian relaxation; ECM = epsilon-constraint method; SA = sensitivity analysis; SAN = simulated annealing; MCS = Monte Carlo simulation; ACO = ant colony optimization; RO = robust optimizationTable 2
The results of the developed epsilon-constraint method and Lagrangian relaxation.
Problem numberIndicesGAMSLROptimality gapijklmnRun timeObjective valueRun timeObjective value11122212.38.13E + 074.18.13E + 070.000296%21222215.58.18E + 077.28.18E + 070.000169%31232218.982067192.819.38.21E + 070.000393%423322113.682461487.5914.38.25E + 070.000342%523422119.182970169.4218.78.30E + 070.000333%62442212680770169.42228.08E + 070.000288%724522134.382770169.4227.28.28E + 070.000400%834533149.686426698.3724.78.64E + 070.000221%935533157.688269811.2921.58.83E + 070.000370%103554326789217186.9423.98.92E + 070.000035%1136544273.493630025.62269.36E + 070.000295%1237644298.992473417.9729.59.25E + 070.000397%13489542146.895678997.5733.49.57E + 070.000386%144810542194.296060458.8436.59.61E + 070.000339%1541010542220102995413.840.71.03E + 080.000233%1641011542298.4118700925.644.81.19E + 080.002866%1741111542409.5122439701.159.41.22E + 080.029194%1841111642511.412021169774.61.20E + 080.018774%1941111652699.4128811065.986.61.29E + 080.013647%2041211652940151835127.6107.61.52E + 080.026376%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. This 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.The remainder of the paper is given as follows: in Section3, 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. The preprint version of this research is uploaded earlier [56].


## 3. CLSC Setting and Mathematical Model

Figure1 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.Figure 1
Presentation of a CLSC, own source.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.


### 3.1. Model Denotation

The following notation is used to present the mathematical model:Indices:i: Index for fixed location of a factoryj: Index for potential location for DCsk: Index for fix location for customer zonel: Index for potential location for collection centersm” Index for potential location for recovery centersn: Index for type of transportation modeParameters:Dkj: Demand of customer in zone k that goes to distribution center jBk: Fraction rate of return products from customer zone kfcdj: Fixed cost for opening DC jfccl: Fixed cost for opening collection center lfcrm: Fixed cost for opening recovery center mpccl: Processing cost at collection center lpcdj: Processing cost at distribution center jpcrm: Processing cost at recovery center jcoij: Transportation cost from plant i to distribution center jcujk: Transportation cost from distribution j to customer zone kcqkl: Transportation cost from customer zone k to collection center lcplm: Transportation cost from collection center l to recover center mchmj: Transportation cost from recover center m to distribution center jpxj: Maximum capacity of distribution center jpzm: Maximum capacity of recovery center mpyl: Maximum capacity of collection center lppi: Maximum capacity of plantji: Number of created job opportunities out of opening distribution center jjl: Number of created job opportunities out of opening collection center ljm: Number of created job opportunities out of opening recovery center mevj: Regional economic value of distribution center jevl: Regional economic value of collection center levm: Regional economic value of recovery center mrdj: Regional development factor in region jrdl: Regional development factor in region lrdm: Regional development factor in region meiijn: Environmental impact of shipping one unit product from node i to j using transportation mode neijkn: Environmental impact of shipping one unit product from node j to k using transportation mode neikln: Environmental impact of shipping one unit product from node k to l using transportation mode neilmn: Environmental impact of shipping one unit product from node l to m using transportation mode neimjn: Environmental impact of shipping one unit product from node m to j using transportation mode nurj: Unemployment rate in region jurl: Unemployment rate in region lurm: 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β2: Penalty coefficient for node complexity of collection centersβ3: Penalty coefficient for node complexity of recovery centersedi: Unit gas emission per a production in plant center iedm: Unit gas emission per recovery product in recovery center mζ1: Penalty coefficient for critical plant centersζ2: Penalty coefficient for critical distribution centersζ3: Penalty coefficient for critical collection centersDecision variables:oijn: Number of products transferred from plant i to distribution center j with transportation mode nujkn: Number of products transferred from distribution center j to customer zone k with transportation mode nqkl: Number of products transferred from customer zone k to collection center lplmn: Number of products transferred from collection center l to recovery center m with transportation mode nhmjn: Number of products transferred from recovery center m to distribution center j with transportation mode nxj: 1, if a distribution center is opened at location j; 0, otherwisey1: 1, if a collection center is opened at location l; 0, otherwisezm: 1, if a recovery center is opened at location m; 0, otherwisewi′: 1, if a plant center’s node is critical; 0, otherwisewj′: 1, if a distribution center’s node is critical; 0, otherwisewl′: 1, if a collection center’s node is critical; 0, otherwisewm′: 1, if a recovery center’s node is critical; 0, otherwisevij: 1, if node i is assigned to node j; 0, otherwisevjk: 1, if node j is assigned to node k; 0, otherwisevlm: 1, if node l is assigned to node m; 0, otherwisevmj: 1, if node m is assigned to node j; 0, otherwise


### 3.2. Measures of Sustainability

The 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. This also concerns the regional development where the focus is on less developed areas.


### 3.3. 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. Thus, 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.


#### 3.3.1. Flow Complexity

Equation (1) calculates the interactions and relationships between all nodes of the network. The first part is concerned with the allocation of distribution centers to the production centers. The second part deals with the relationships between the demand zone and the distribution centers. The third part allocates collection center to demand zones. The fourth part allocates collection centers to recovery centers, and finally, the last term allocates recovery centers to distribution centers [58].(1)∑l∑jvij+∑j∑kvjk+∑k∑lvkl+∑l∑mvlm+∑m∑jvmj+vij,vlm,vmj,vkl,vjk,∈1,0,∀i,j,l,k,m.According to equation (1), if overall number of connections is high, the total network flow is deemed to be utterly complex.


#### 3.3.2. Node Complexity

Equation (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.”(2)∑jxj+∑lyl+∑mzm,xi,yl,zm∈0,1,∀j,l,m.


### 3.4. Mathematical Model Formulation

The 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)
There is a pressure mechanism in the forward flow(4)
There is a tensile mechanism in the reverse flow(5)
The location of factories and customer locations are predetermined and fixedThe 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?The mathematical model is given as follows, starting with the objective functions composed of four functions where the first one (w1) minimizes all costs, the second one (w2) improves social aspects, the third one (w3) minimizes air pollution, and the fourth one (w4) improves flexibility and network resilience. The first objective function (w1) 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 minimizingw2. The third objective tries to reduce the gas emission of facilities and transportation. The fourth objective function maximizes the resilience of the whole network by decreasing the complexity of supply chain.(3)minw1=∑ifcdixj+∑lfcclyl+∑mfcrmzm+∑j∑k∑ncujkn+pcdjujkn+∑i∑j∑ncoijn+coijn+∑k∑lcqklqkl+∑l∑m∑ncplmn+pcclplmn+∑m∑j∑nchmjn+pcrmhmjn,minw2=∑jjjurjxj+∑ljlurlyl+∑mjmurmzmmax−∑jjjurjxj+∑ljlurlyl+∑mjmurmzm∑jjjurjxj+∑ljlurlyl+∑mjmurmzmmax−∑jjjurjxj+∑ljlurlyl+∑mjmurmzmmin+∑jevjxj1−rdj+∑levlyli−rdl+∑mevmzm1−rdmmax−∑jevjxj−rdj+∑jevlyl1−rdl+∑mevmzm1−rdm∑jevjxj1−rdj+∑levlyli−rdl+∑mevmzm1−rdmmax−∑jevjxj1−rdj+∑jevlyli−rdl+∑mevmzm1−rdmmin,minw3=∑i∑j∑neiijn+edioijn+∑j∑k∑neijknujkn+∑l∑m∑neilmnplmn+∑m∑j∑neimjn+edmhmjn,minw4=∑l∑jα1¯vij+∑j∑kα2¯vjk+∑l∑mα3¯vlm+∑m∑jα4¯vmj+∑jβ1xj+∑lβ2yl+∑mβ3zm+∑iζ1wi+∑jζ2wj′+∑kζ3wl′+∑mζ4wm′,subject to the following constraints:(4)∑n∑jujkn≥Dk⋅∀k,(5)∑lqkl≤Bk×∑nujkn⋅∀k,(6)∑i∑noijn+∑m∑nhmjn=∑k∑nujkn⋅∀k,(7)∑j∑nujkn≥∑lqkl⋅∀k,(8)∑kqkl=∑m∑nplmn⋅∀l,(9)∑l∑nolmn=∑j∑nhmjn⋅∀m,(10)∑j∑noijn≤ppi⋅∀i,(11)∑i∑noijn+∑m∑nhmjn≤xjpxj⋅∀j,(12)∑kqkl≤ylpyl⋅∀l,(13)∑l∑nplmn≤zmpzm⋅∀m,(14)∑noijn≤M.vij⋅∀i,j,(15)∑nujkn≤M.vjk⋅∀k,j,(16)∑nplmn≤M.vlm⋅∀l,m,(17)∑nhmjn≤M.vmj⋅∀m,j,(18)Ui−∑j∑noijn1−wi′≥ε1−wi′⋅∀i,(19)Uj−∑i∑noijn+∑k∑nujkn+∑m∑nhmjn1−wj′≥ε⋅∀j,(20)Ul−∑m∑nplmn+∑kqkl1−wj′≥ε1−wl′⋅∀1,(21)Um−∑l∑nplmn+∑j∑nhmjn1−wm′≥ε1−wm′⋅∀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. They 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. The 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]. The 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. This 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:(22)μm=m−γ1γ2−γ1,ifγ1≤m≤γ2,1,ifγ2≤m≤γ3,γ4−mγ4−γ3,ifγ3≤m≤γ4,0,if otherwise,(23)PosΦ≤m≥Ψ⇔m−γ1γ2−γ1≥Ψ⇔m≥1−Ψγ1+Ψγ2,PosΦ≥m≥Ψ⇔γ4−mγ4−γ3≥Ψ⇔m≤1−Ψγ4+Ψγ3,NecΦ≤m≥Ψ⇔m−γ3γ4−γ3≥Ψ⇔m≥1−Ψγ3+Ψγ4,NecΦ≥m≥Ψ⇔γ2−mγ2−γ1≥Ψ⇔m≥1−Ψγ2+Ψγ1.The mathematical expectation value of the fuzzy variableΦ is explained in the condition of the following equation:(24)EPosΦ=∫0∞PosΦ≥mdx−∫−∞0PosΦ≤mdx=γ3+γ42,ENecΦ=∫0∞NecΦ≥mdx−∫−∞0NecΦ≤mdx=γ2+γ12.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. The DM can choose extreme attitudes using a convex combination of any point in the range:(25)MeF=Υ.PosF+1−Υ.NecF=NecF+Υ.PosF−NecF.In the above equation,γ 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(26)MeΦ≤m=0,ifm≤γ2,Υ×m−γ1γ2−γ1,ifγ1≤m≤γ2,1,ifγ2≤m≤γ3,Υ+1−Υm−γ3γ4−γ3,ifγ3≤m≤γ4,1,ifγ4≤m,(27)MeΦ≤m=1,ifm≤γ2,Υ+1−Υγ2−mγ2−γ1,ifγ1≤m≤γ2,1,ifγ2≤m≤γ3,Υ+1−Υ×γ4−mγ4−γ3,ifγ3≤m≤γ4,0,ifγ4≤m.The value of the mathematical expectation of the fuzzy variable for the valuesγ1≥0 are(28)EMeΦ=∫0∞MeΦ≥mdx−∫−∞0MeΦ≤mdx=1−Υ2γ3+γ4+Υ2γ1+γ2.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. Therefore, in this study, we obtain(29)MeΦ≥m≥ψ⇔1−Υγ2−γ1+Υ≥Ψ,MeΦ≤m≥ψ⇔1−Υm−γ3γ4−γ3+Υ≥Ψ.According to the given relationships and equations so far, we can write the possibilistic chance-constrained program (PCCP) model as follows:(30)minP=Q.Ys.t:M≤B.YC.X≥dA.X=0X≤N.LG.L≤1Y,L∈0,1,M≥0.In the above basic model, the matricesB, 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. The 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 Figure2. For example, Τ˜=χ1,χ2,χ3,χ4.Figure 2
The trapezoidal possibility distribution of fuzzy parameter.According to the given description, we have(31)minEP=EQ˜.Ys.t:MeM<+B.Y≥φMeC.M≥d≥σA.M=0OY≤1M<+N>LG>L<+1Y,L∈0,1,M≥0.The 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. Through above transformations, the crisp equivalent of the BPCCP model based on Me Measure and the expected value operator is given as follows:(32)minEP=1−Υ2Q1−Q2+Υ2Q3+Q4.Ys.t:M≤φ−ΥB1+1−φB21−Υ.YM>+σ−Υd3+1−σd41−ΥA.M=0O.Y≤1M≤N.LG.L≤1Y,L∈0,1,M>+0.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. This 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. The final model was rewritten as follows by considering these uncertainties.The first and third objective functions are modified as follows:(33)minΕw1=∑ifcdj.xj+∑lfccl.yl+∑mfcrm.zm+∑j∑k∑ncujkn.ujkn+∑i∑j∑ncoijn+c.oijn+∑k∑lcqkl.qkl+∑l∑m∑ncplmn.plmn+∑m∑j∑nchmjn.hmjn+∑j∑k∑n1−fzλ2Ξj1+Ξj2+fzλ2Ξj3+Ξj4.ujkn+∑l∑m∑n1−fzλ2φl1+φl2+fzλ2φl3+φl4,minEw3=∑i∑j∑neiijn.oijn+∑j∑k∑neijkn.ujkn+∑l∑m∑neilmn.plmn+∑m∑j∑neimjn.hmjn+∑i∑j∑n1−fzλ2Δ1−Δ2+fzλ2Δ3+Δ4∑n∑jujkn≥ηk−fzλdk4+1−ηkdk31−fzλ∀k,(34)∑lqkl≤ϖk−fzλdk4+1−ϖkdk21−fzλ×∑j∑nujkn∀k.The 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. This method has been introduced by Haimes [66]. It involves prioritizing the objective functions, 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. The proper application of the ε-constraint method requires determining all objective-function 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. The 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. The 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. This ambiguity can be resolved by conducting lexicographic optimization for all the objective functions to build a payoff table with Pareto optimal solutions only. The 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. Those 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. This 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. This is done by converting all the objective functions except one to constraints. The 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.(35)MaxM1x+epsilon×n2t2+10−1×n3t3+10a−2×nata,(36)s.tMbx−wb=Hb,b=2,3,…,p,x∈S;wb∈R+.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. The abovementioned model ensures that the obtained optimal solution is an efficient one only if all the constraints of the objective function are binding.The 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. The Lagrangian relaxation method is known as one of the best approaches to solve supply chain problems aiming at obtaining robust solutions. This 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. These steps are repeated until the lower and upper bounds reach the specified limit.


#### 3.4.1. 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:(37)minWulmn2,uijn1=minw1+∑i∑j∑nuijn1×oijn−ppi+∑l∑m∑nulmn2×plmn−Zm×PZm.In the abovementioned objective function, the valuesUlmn2 and Uijn1 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.


#### 3.4.2. 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). This is the upper bound of the problem.Updating upper and lower bounds,(38)uijn1,x+1=max0,u1,x−stp1,x∑i∑j∑noijn−ppi,(39)ulmn2,x+1=max0,u2,x−stp2,x∑l∑m∑nplmn−Zm×PZm.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 stp1,x and stp2,x which are calculated as follows:(40)stp1,x=vx×UP−LBx∑i∑j∑nOijn−PPi2,(41)stp2,x=vx×UP−LBn∑l∑m∑nPijn−zm×PZm2.In equations (40) and (41), the values UP and LBx 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 ν. The problem-solving process ends with two conditions:(1)
Obtaining a feasible solution with the appropriate tolerances(2)
The value ofstp has reached the minimum possible value


## 3.1. Model Denotation

The following notation is used to present the mathematical model:Indices:i: Index for fixed location of a factoryj: Index for potential location for DCsk: Index for fix location for customer zonel: Index for potential location for collection centersm” Index for potential location for recovery centersn: Index for type of transportation modeParameters:Dkj: Demand of customer in zone k that goes to distribution center jBk: Fraction rate of return products from customer zone kfcdj: Fixed cost for opening DC jfccl: Fixed cost for opening collection center lfcrm: Fixed cost for opening recovery center mpccl: Processing cost at collection center lpcdj: Processing cost at distribution center jpcrm: Processing cost at recovery center jcoij: Transportation cost from plant i to distribution center jcujk: Transportation cost from distribution j to customer zone kcqkl: Transportation cost from customer zone k to collection center lcplm: Transportation cost from collection center l to recover center mchmj: Transportation cost from recover center m to distribution center jpxj: Maximum capacity of distribution center jpzm: Maximum capacity of recovery center mpyl: Maximum capacity of collection center lppi: Maximum capacity of plantji: Number of created job opportunities out of opening distribution center jjl: Number of created job opportunities out of opening collection center ljm: Number of created job opportunities out of opening recovery center mevj: Regional economic value of distribution center jevl: Regional economic value of collection center levm: Regional economic value of recovery center mrdj: Regional development factor in region jrdl: Regional development factor in region lrdm: Regional development factor in region meiijn: Environmental impact of shipping one unit product from node i to j using transportation mode neijkn: Environmental impact of shipping one unit product from node j to k using transportation mode neikln: Environmental impact of shipping one unit product from node k to l using transportation mode neilmn: Environmental impact of shipping one unit product from node l to m using transportation mode neimjn: Environmental impact of shipping one unit product from node m to j using transportation mode nurj: Unemployment rate in region jurl: Unemployment rate in region lurm: 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β2: Penalty coefficient for node complexity of collection centersβ3: Penalty coefficient for node complexity of recovery centersedi: Unit gas emission per a production in plant center iedm: Unit gas emission per recovery product in recovery center mζ1: Penalty coefficient for critical plant centersζ2: Penalty coefficient for critical distribution centersζ3: Penalty coefficient for critical collection centersDecision variables:oijn: Number of products transferred from plant i to distribution center j with transportation mode nujkn: Number of products transferred from distribution center j to customer zone k with transportation mode nqkl: Number of products transferred from customer zone k to collection center lplmn: Number of products transferred from collection center l to recovery center m with transportation mode nhmjn: Number of products transferred from recovery center m to distribution center j with transportation mode nxj: 1, if a distribution center is opened at location j; 0, otherwisey1: 1, if a collection center is opened at location l; 0, otherwisezm: 1, if a recovery center is opened at location m; 0, otherwisewi′: 1, if a plant center’s node is critical; 0, otherwisewj′: 1, if a distribution center’s node is critical; 0, otherwisewl′: 1, if a collection center’s node is critical; 0, otherwisewm′: 1, if a recovery center’s node is critical; 0, otherwisevij: 1, if node i is assigned to node j; 0, otherwisevjk: 1, if node j is assigned to node k; 0, otherwisevlm: 1, if node l is assigned to node m; 0, otherwisevmj: 1, if node m is assigned to node j; 0, otherwise


## 3.2. Measures of Sustainability

The 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. This also concerns the regional development where the focus is on less developed areas.


## 3.3. 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. Thus, 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.


### 3.3.1. Flow Complexity

Equation (1) calculates the interactions and relationships between all nodes of the network. The first part is concerned with the allocation of distribution centers to the production centers. The second part deals with the relationships between the demand zone and the distribution centers. The third part allocates collection center to demand zones. The fourth part allocates collection centers to recovery centers, and finally, the last term allocates recovery centers to distribution centers [58].(1)∑l∑jvij+∑j∑kvjk+∑k∑lvkl+∑l∑mvlm+∑m∑jvmj+vij,vlm,vmj,vkl,vjk,∈1,0,∀i,j,l,k,m.According to equation (1), if overall number of connections is high, the total network flow is deemed to be utterly complex.


### 3.3.2. Node Complexity

Equation (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.”(2)∑jxj+∑lyl+∑mzm,xi,yl,zm∈0,1,∀j,l,m.


## 3.3.1. Flow Complexity

Equation (1) calculates the interactions and relationships between all nodes of the network. The first part is concerned with the allocation of distribution centers to the production centers. The second part deals with the relationships between the demand zone and the distribution centers. The third part allocates collection center to demand zones. The fourth part allocates collection centers to recovery centers, and finally, the last term allocates recovery centers to distribution centers [58].(1)∑l∑jvij+∑j∑kvjk+∑k∑lvkl+∑l∑mvlm+∑m∑jvmj+vij,vlm,vmj,vkl,vjk,∈1,0,∀i,j,l,k,m.According to equation (1), if overall number of connections is high, the total network flow is deemed to be utterly complex.


## 3.3.2. Node Complexity

Equation (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.”(2)∑jxj+∑lyl+∑mzm,xi,yl,zm∈0,1,∀j,l,m.


## 3.4. Mathematical Model Formulation

The 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)
There is a pressure mechanism in the forward flow(4)
There is a tensile mechanism in the reverse flow(5)
The location of factories and customer locations are predetermined and fixedThe 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?The mathematical model is given as follows, starting with the objective functions composed of four functions where the first one (w1) minimizes all costs, the second one (w2) improves social aspects, the third one (w3) minimizes air pollution, and the fourth one (w4) improves flexibility and network resilience. The first objective function (w1) 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 minimizingw2. The third objective tries to reduce the gas emission of facilities and transportation. The fourth objective function maximizes the resilience of the whole network by decreasing the complexity of supply chain.(3)minw1=∑ifcdixj+∑lfcclyl+∑mfcrmzm+∑j∑k∑ncujkn+pcdjujkn+∑i∑j∑ncoijn+coijn+∑k∑lcqklqkl+∑l∑m∑ncplmn+pcclplmn+∑m∑j∑nchmjn+pcrmhmjn,minw2=∑jjjurjxj+∑ljlurlyl+∑mjmurmzmmax−∑jjjurjxj+∑ljlurlyl+∑mjmurmzm∑jjjurjxj+∑ljlurlyl+∑mjmurmzmmax−∑jjjurjxj+∑ljlurlyl+∑mjmurmzmmin+∑jevjxj1−rdj+∑levlyli−rdl+∑mevmzm1−rdmmax−∑jevjxj−rdj+∑jevlyl1−rdl+∑mevmzm1−rdm∑jevjxj1−rdj+∑levlyli−rdl+∑mevmzm1−rdmmax−∑jevjxj1−rdj+∑jevlyli−rdl+∑mevmzm1−rdmmin,minw3=∑i∑j∑neiijn+edioijn+∑j∑k∑neijknujkn+∑l∑m∑neilmnplmn+∑m∑j∑neimjn+edmhmjn,minw4=∑l∑jα1¯vij+∑j∑kα2¯vjk+∑l∑mα3¯vlm+∑m∑jα4¯vmj+∑jβ1xj+∑lβ2yl+∑mβ3zm+∑iζ1wi+∑jζ2wj′+∑kζ3wl′+∑mζ4wm′,subject to the following constraints:(4)∑n∑jujkn≥Dk⋅∀k,(5)∑lqkl≤Bk×∑nujkn⋅∀k,(6)∑i∑noijn+∑m∑nhmjn=∑k∑nujkn⋅∀k,(7)∑j∑nujkn≥∑lqkl⋅∀k,(8)∑kqkl=∑m∑nplmn⋅∀l,(9)∑l∑nolmn=∑j∑nhmjn⋅∀m,(10)∑j∑noijn≤ppi⋅∀i,(11)∑i∑noijn+∑m∑nhmjn≤xjpxj⋅∀j,(12)∑kqkl≤ylpyl⋅∀l,(13)∑l∑nplmn≤zmpzm⋅∀m,(14)∑noijn≤M.vij⋅∀i,j,(15)∑nujkn≤M.vjk⋅∀k,j,(16)∑nplmn≤M.vlm⋅∀l,m,(17)∑nhmjn≤M.vmj⋅∀m,j,(18)Ui−∑j∑noijn1−wi′≥ε1−wi′⋅∀i,(19)Uj−∑i∑noijn+∑k∑nujkn+∑m∑nhmjn1−wj′≥ε⋅∀j,(20)Ul−∑m∑nplmn+∑kqkl1−wj′≥ε1−wl′⋅∀1,(21)Um−∑l∑nplmn+∑j∑nhmjn1−wm′≥ε1−wm′⋅∀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. They 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. The 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]. The 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. This 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:(22)μm=m−γ1γ2−γ1,ifγ1≤m≤γ2,1,ifγ2≤m≤γ3,γ4−mγ4−γ3,ifγ3≤m≤γ4,0,if otherwise,(23)PosΦ≤m≥Ψ⇔m−γ1γ2−γ1≥Ψ⇔m≥1−Ψγ1+Ψγ2,PosΦ≥m≥Ψ⇔γ4−mγ4−γ3≥Ψ⇔m≤1−Ψγ4+Ψγ3,NecΦ≤m≥Ψ⇔m−γ3γ4−γ3≥Ψ⇔m≥1−Ψγ3+Ψγ4,NecΦ≥m≥Ψ⇔γ2−mγ2−γ1≥Ψ⇔m≥1−Ψγ2+Ψγ1.The mathematical expectation value of the fuzzy variableΦ is explained in the condition of the following equation:(24)EPosΦ=∫0∞PosΦ≥mdx−∫−∞0PosΦ≤mdx=γ3+γ42,ENecΦ=∫0∞NecΦ≥mdx−∫−∞0NecΦ≤mdx=γ2+γ12.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. The DM can choose extreme attitudes using a convex combination of any point in the range:(25)MeF=Υ.PosF+1−Υ.NecF=NecF+Υ.PosF−NecF.In the above equation,γ 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(26)MeΦ≤m=0,ifm≤γ2,Υ×m−γ1γ2−γ1,ifγ1≤m≤γ2,1,ifγ2≤m≤γ3,Υ+1−Υm−γ3γ4−γ3,ifγ3≤m≤γ4,1,ifγ4≤m,(27)MeΦ≤m=1,ifm≤γ2,Υ+1−Υγ2−mγ2−γ1,ifγ1≤m≤γ2,1,ifγ2≤m≤γ3,Υ+1−Υ×γ4−mγ4−γ3,ifγ3≤m≤γ4,0,ifγ4≤m.The value of the mathematical expectation of the fuzzy variable for the valuesγ1≥0 are(28)EMeΦ=∫0∞MeΦ≥mdx−∫−∞0MeΦ≤mdx=1−Υ2γ3+γ4+Υ2γ1+γ2.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. Therefore, in this study, we obtain(29)MeΦ≥m≥ψ⇔1−Υγ2−γ1+Υ≥Ψ,MeΦ≤m≥ψ⇔1−Υm−γ3γ4−γ3+Υ≥Ψ.According to the given relationships and equations so far, we can write the possibilistic chance-constrained program (PCCP) model as follows:(30)minP=Q.Ys.t:M≤B.YC.X≥dA.X=0X≤N.LG.L≤1Y,L∈0,1,M≥0.In the above basic model, the matricesB, 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. The 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 Figure2. For example, Τ˜=χ1,χ2,χ3,χ4.Figure 2
The trapezoidal possibility distribution of fuzzy parameter.According to the given description, we have(31)minEP=EQ˜.Ys.t:MeM<+B.Y≥φMeC.M≥d≥σA.M=0OY≤1M<+N>LG>L<+1Y,L∈0,1,M≥0.The 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. Through above transformations, the crisp equivalent of the BPCCP model based on Me Measure and the expected value operator is given as follows:(32)minEP=1−Υ2Q1−Q2+Υ2Q3+Q4.Ys.t:M≤φ−ΥB1+1−φB21−Υ.YM>+σ−Υd3+1−σd41−ΥA.M=0O.Y≤1M≤N.LG.L≤1Y,L∈0,1,M>+0.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. This 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. The final model was rewritten as follows by considering these uncertainties.The first and third objective functions are modified as follows:(33)minΕw1=∑ifcdj.xj+∑lfccl.yl+∑mfcrm.zm+∑j∑k∑ncujkn.ujkn+∑i∑j∑ncoijn+c.oijn+∑k∑lcqkl.qkl+∑l∑m∑ncplmn.plmn+∑m∑j∑nchmjn.hmjn+∑j∑k∑n1−fzλ2Ξj1+Ξj2+fzλ2Ξj3+Ξj4.ujkn+∑l∑m∑n1−fzλ2φl1+φl2+fzλ2φl3+φl4,minEw3=∑i∑j∑neiijn.oijn+∑j∑k∑neijkn.ujkn+∑l∑m∑neilmn.plmn+∑m∑j∑neimjn.hmjn+∑i∑j∑n1−fzλ2Δ1−Δ2+fzλ2Δ3+Δ4∑n∑jujkn≥ηk−fzλdk4+1−ηkdk31−fzλ∀k,(34)∑lqkl≤ϖk−fzλdk4+1−ϖkdk21−fzλ×∑j∑nujkn∀k.The 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. This method has been introduced by Haimes [66]. It involves prioritizing the objective functions, 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. The proper application of the ε-constraint method requires determining all objective-function 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. The 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. The 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. This ambiguity can be resolved by conducting lexicographic optimization for all the objective functions to build a payoff table with Pareto optimal solutions only. The 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. Those 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. This 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. This is done by converting all the objective functions except one to constraints. The 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.(35)MaxM1x+epsilon×n2t2+10−1×n3t3+10a−2×nata,(36)s.tMbx−wb=Hb,b=2,3,…,p,x∈S;wb∈R+.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. The abovementioned model ensures that the obtained optimal solution is an efficient one only if all the constraints of the objective function are binding.The 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. The Lagrangian relaxation method is known as one of the best approaches to solve supply chain problems aiming at obtaining robust solutions. This 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. These steps are repeated until the lower and upper bounds reach the specified limit.


### 3.4.1. 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:(37)minWulmn2,uijn1=minw1+∑i∑j∑nuijn1×oijn−ppi+∑l∑m∑nulmn2×plmn−Zm×PZm.In the abovementioned objective function, the valuesUlmn2 and Uijn1 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.


### 3.4.2. 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). This is the upper bound of the problem.Updating upper and lower bounds,(38)uijn1,x+1=max0,u1,x−stp1,x∑i∑j∑noijn−ppi,(39)ulmn2,x+1=max0,u2,x−stp2,x∑l∑m∑nplmn−Zm×PZm.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 stp1,x and stp2,x which are calculated as follows:(40)stp1,x=vx×UP−LBx∑i∑j∑nOijn−PPi2,(41)stp2,x=vx×UP−LBn∑l∑m∑nPijn−zm×PZm2.In equations (40) and (41), the values UP and LBx 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 ν. The problem-solving process ends with two conditions:(1)
Obtaining a feasible solution with the appropriate tolerances(2)
The value ofstp has reached the minimum possible value


## 3.4.1. 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:(37)minWulmn2,uijn1=minw1+∑i∑j∑nuijn1×oijn−ppi+∑l∑m∑nulmn2×plmn−Zm×PZm.In the abovementioned objective function, the valuesUlmn2 and Uijn1 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.


## 3.4.2. 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). This is the upper bound of the problem.Updating upper and lower bounds,(38)uijn1,x+1=max0,u1,x−stp1,x∑i∑j∑noijn−ppi,(39)ulmn2,x+1=max0,u2,x−stp2,x∑l∑m∑nplmn−Zm×PZm.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 stp1,x and stp2,x which are calculated as follows:(40)stp1,x=vx×UP−LBx∑i∑j∑nOijn−PPi2,(41)stp2,x=vx×UP−LBn∑l∑m∑nPijn−zm×PZm2.In equations (40) and (41), the values UP and LBx 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 ν. The problem-solving process ends with two conditions:(1)
Obtaining a feasible solution with the appropriate tolerances(2)
The value ofstp has reached the minimum possible value


## 4. 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. Table3 provides the input parameters. All numerical examples are calculated using GAMS Windows application on a laptop with Windows 10, Intel Core i7 CPU, and 8 GB of RAM. At a first stage, the presented numerical examples are solved using GAMS without applying the Lagrangian relaxation method. In a second step, the model is solved using the Lagrangian relaxation method. The results including “objective function,” “solution time,” and “difference between objective functions” are presented in Table 2 and Figures 1–5.Table 3
Values of input parameters.
DkjUniform16000,17000BkUniform0.03,0.1fcdjUniform18000,19000fcdjUniform18000,19000fcc1Uniform17500,18500fcrmUniform22000,25000pcc1Uniform1,2pcdjUniform1.2,1.9pcrmUniform.9,1.3coijUniform10,20cujkUniform4,7chmjUniform5,11pxjUniform25000,36000pzmUniform33000,36000py1Uniform21000,26000ppiUniform42000,46000jjUniform120,160jmUniform220,260evjUniform1.2,1.6ev1Uniform0.6,1.1evjUniform1.3,1.7cqklUniform2,6cp1mUniform1,8Figure 3
Values of the objective function obtained from solving methods with and without LR.Figure 4
Difference percentage of the objective function values with and without LR method.Figure 5
Solving time of the model with and without LR.We assume uniform distributions for all input parameters but with differing ranges. As can be easily seen in Table2, 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. This shows how much faster the problem is solved through Lagrangian relaxation which also indicates the effectiveness of this method.Figure3 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.Figure4 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, Figure5 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. Therefore, it can be concluded that this solution method is more efficient and cost-effective than epsilon constraint.


## 5. 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. The focus points of this industry are high quality, reliability of availability, and cost-effectiveness of resources that need to be assured by the CLSC. The 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. Thus, 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.The case company’s headquarters is in Tehran. The company produces one type of tire in five different cities of Iran to meet the needs and requirements of 11 customer points. The 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)∑jxj=12,∑lyl=7,∑mzm=6.Figure6 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.Figure 6
Geographical location of candidate and exist facilities.The model of the case study was solved using the solution method that was described in detail in previous sections. Figure7 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.Figure 7
Optimal facilities’ locations and transfer lines.The 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). The optimum value of second objective function which contains job opportunity and unemployment rate is -127426.526138657 units. The best value of CO2 emission for the third objective function is 940431.149676294 litter. The 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.


## 6. 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. The Capacity of the Distribution Facility

We first analyze the effects of facility changes selected in the optimal solution based on the objective functions. The selection of facilities influences the capacity available at the DCs. The capacity of each facility is increased to 5,000 units per each test to determine the effects in a clear manner.As shown in Figures8 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 CO2 emissions. These managerial insights give great support to decision-makers in showing the effects of decisions to the overall CLSC.Figure 8
The impact of distribution centers’ capacity change on the number of distribution centers.Figure 9
The impact of distribution centers’ capacity changes on the total cost.Figure 10
The impact of distribution centers’ capacity changes on the social effects.Figure 11
The impact of distribution centers’ capacity changes on the environmental effects.


### 6.2. The 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 Figure12, 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%. The 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. Thus, if the confidence level in increased up to 35%, then the social effects of the model would experience an 8% reduction.Figure 12
The impact of the confidence level of returned product changes on the number of facilities.Figure 13
The impact of the confidence level of returned product changes on total costs.Figure 14
The impact of the confidence level of returned product changes on gas emission.Figure 15
The impact of the confidence level of returned product changes on social effects.


### 6.3. 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 Figure16 to obtain a better understanding of the impacts of variety of changes on the flexibility of a CLSC.Figure 16
The impact of confidence important parameter changes on the resilience concept.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. The 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.


## 6.1. The Capacity of the Distribution Facility

We first analyze the effects of facility changes selected in the optimal solution based on the objective functions. The selection of facilities influences the capacity available at the DCs. The capacity of each facility is increased to 5,000 units per each test to determine the effects in a clear manner.As shown in Figures8 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 CO2 emissions. These managerial insights give great support to decision-makers in showing the effects of decisions to the overall CLSC.Figure 8
The impact of distribution centers’ capacity change on the number of distribution centers.Figure 9
The impact of distribution centers’ capacity changes on the total cost.Figure 10
The impact of distribution centers’ capacity changes on the social effects.Figure 11
The impact of distribution centers’ capacity changes on the environmental effects.


## 6.2. The 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 Figure12, 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%. The 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. Thus, if the confidence level in increased up to 35%, then the social effects of the model would experience an 8% reduction.Figure 12
The impact of the confidence level of returned product changes on the number of facilities.Figure 13
The impact of the confidence level of returned product changes on total costs.Figure 14
The impact of the confidence level of returned product changes on gas emission.Figure 15
The impact of the confidence level of returned product changes on social effects.


## 6.3. 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 Figure16 to obtain a better understanding of the impacts of variety of changes on the flexibility of a CLSC.Figure 16
The impact of confidence important parameter changes on the resilience concept.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. The 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.


## 7. Conclusion

This 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 CO2 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.There 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.



---

*Source: 9951220-2021-09-27.xml*