Designing a Scenario-Based Fuzzy Model for Sustainable Closed-Loop Supply Chain Network considering Statistical Reliability: A New Hybrid Metaheuristic Algorithm

,


Introduction
Globalization, the increased regulations of governmental and nongovernmental organizations, and the pressure and requests from customers to comply with environmental issues have led organizations to regard the necessary steps to apply sustainable closed-loop supply chain (SCLSC) management aimed at improving their environmental and economic performance [1].Te SCLSC management integrates the SCM with environmental requirements in all stages of product design, selection and supply of raw materials, production and manufacturing, distribution and transfer processes, delivery to the customer, and ultimately after consumption, recycling and reuse management aimed at maximizing the efciency of energy, and resources consumption associated with improving the performance of the entire supply chain [2,3].Te problem of vehicle routing in the supply chain distribution network is recognized as one of the subproblems of SCM, which involves selecting and allocating possible routes to available vehicles for distribution and delivery of goods to distribution centers or customers designated to minimize the relevant costs.Te optimal solving of this problem will lead to timely delivery of goods, reduce the need for warehousing and maintenance of goods, and increase customer satisfaction meanwhile reducing the distribution costs [4].One may also claim that vehicle routing is one of the most challenging issues in the context of transportation and support of the supply chain [5].A variety of products have to be made available to the customers at their request in today's global competitions.Te customers' demand for high quality and fast serviceproviding has led to enhanced pressures that have not faced before [6].In today's economic and industrial environment and given the growing trend of industries followed by the rise in environmental pollution in contrast, and most importantly, the use of limited resources, the need to recycle resources from manufactured products as well as informing the consumers about the need for a change in attitude seem to be a priority not only in the production of goods but in all the stages of production [7].Ensuring sustainable development in any country is nowadays subject to the conservation and optimal use of limited and irreplaceable resources in that country.Tus, many measures have been taken in this direction, including recycling waste in the production cycle, the reuse of consumer goods, returning the quality control returned goods to the production line, recycling, etc. Te set of these activities accompanied by applying environmental and social considerations form the concept of the SCLSC [8,9].On the other hand, aimed at accurate management and preventing the waste of resources, the management units currently have no choice other than to adopt and employ new scientifc approaches, models, processes, and techniques tailored to the current conditions to provide proper performance regarding the resources used [10].It is usually assumed in supply chain network (SCN) design studies that active utilities (facilities) are able to provide service continuously for a long period of time without any breakdown and will continue to operate without interruption.However, supply chains face a high degree of uncertainty due to their complex nature, which can adversely afect the quality of their performance [11].In the real world, utilities may sufer from disruptions and failures with possible causes of human error, natural disasters, etc. [12,13].Tus, the failure of one component of the SCN may disrupt the functioning of the entire supply chain, or in the best case, reduce the efciency of the chain.Hence, it seems essential to consider the factor of reliability in the design of the closed-loop supply chain, especially in its direct components (forward logistic stage).However, this issue has been less addressed in recent studies.In this study, a multiobjective (MO), multiperiod, multicommodity, and scenario-based fuzzy mathematical model, a new hybrid metaheuristic algorithm was proposed for locating, routing, and distributing goods in a sustainable closed loop supply chain.Te uncertainty considered in the mathematical model was fuzzy.In the proposed hybrid approach, the search mechanism and update of the solution in the basic gray wolf algorithm and the crowding distance index mechanism of the MO genetic algorithm were used to select and remove unsuccessful solutions in the external archive.

Literature Review
Te relevant literature that contributes to identify the general framework of this article is reviewed in this section.Zhang et al. [14] provided a mathematical model for a SCLSC based on economic turnaround.Te goals of this model encompass maximum proftability according to the income and costs of the entire chain, minimizing the environmental impacts due to the carbon index and maximizing the social efects according to the job opportunities created.Tey utilized the weighted sum technique to solve the model.Te sensitivity analysis results revealed that the enhanced demand rate has a remarkable impact on the goals.Hassangaviar et al. [15] provided a biobjective mathematical model for the closed-loop chain under uncertainty in prices.Te model objectives included maximizing transaction level satisfaction and the proft from the supply chain.Tey used the NSGA-II algorithm to solve the model, and the results indicated a good performance in terms of maximum expansion and distancing.Mogale et al. [16] provided a CLSCs network, in which, the demand was considered sensitive to the price, consumer motivation, and the quality level.Te core goal of the proposed model is to reduce the total cost and carbon emissions produced by the activities resulting from production, distribution, transportation, and disposal.Tey employed an NSGA-II algorithm and a cokriging approach to solve the model.Te study results revealed the positive efects of motivational pricing on the returned goods.Kazancoglu et al. [17] presented an MILP model for designing the green dual-channel and CLSCs.Te core goal of the proposed model is to optimally select echelons and transportation alternatives between these echelons in a CLSC network based on economic and environmental considerations.Te proposed model is supported by a case study in the home appliance industry in this study.Khalili Nasr et al. [18] provided a new integrated approach based on the best-worst method and MOMILP for designing an SCLSC network with the LIR problem; they applied a fuzzy method for solving their model using GP.Sadeghi Ahangar et al. [19] provided a SCLSC network to manage municipal solid waste using a MILP model based on an FPA.Tis model minimized the total cost, labor, and emission level.Goli et al. [20] examined a multiobjective, multiperiod, and multiproduct closed-loop supply chain (CLSC) model with uncertain parameters aimed at combining the fnancial cash fow as cash fow and debt constraints and employment under uncertainty.Te objectives of the proposed mathematical model in their study included increasing the cash fow, maximizing the jobs created, and maximizing the reliability of raw materials consumed.Tey developed and used MO simulated annealing, MO invasive weed optimization, and MO gray wolf metaheuristic algorithms to solve 2 Complexity the model at a large scale.Ali et al. [21] provided a novel mathematical model for reverse supply chain management of air conditioning products.Teir considered supply chain was sustainable with fuzzy demand uncertainty.Locating hub and recycling centers was among the most important goals of their research.Te case study covered the industries of Saudi Arabia and India.In their study, identifed places were prioritized with a hierarchical process analysis approach after solving the mathematical model.Yun et al. [22] mathematically modeled a MO supply chain by considering economic, environmental, and social criteria.Te objectives regarded by them included minimizing total costs, minimizing the amount of carbon dioxide released, and maximizing social impacts.Teir innovation involved considering three types of distribution channels, including normal delivery, direct delivery, and direct displacement.Te proposed model was solved using the genetic algorithm, whose results indicated the proper performance of the proposed model.Reyhani et al. [23] provided an MO and multiperiod mathematical model for inventory management for the SCLSC.Te most important decisions made in their model included determining the amount of fow at each level and locating the hubs.Te main objectives of the study contained minimizing transportation costs and the costs of carbon dioxide emissions and maximizing social responsibility.Te case study was focused on agricultural products.Rabbani et al. [1] presented a sustainable MO and multilevel mathematical model to locate distribution hubs and allocate warehouses to distribution centers.Te proposed innovation included a variety of technologies for cars, which causes the release of diferent amounts of carbon dioxide from diferent vehicles.Tey used the Epsilon constraint approach to solve the model in small and medium dimensions.Wang et al. [24] provided a mathematical model for the green inverse supply chain.One of their innovations was to consider the pricing of goods using the game theory.Tey considered three diferent prices for diferent types of goods in this article to price goods.Minimizing costs in the supply chain, including manufacturers, distributors, customers, and collection centers were some of the objectives of this study.According to their results, the chain costs were minimized to the desired level.Mohtashami et al. [25] proposed a mathematical model for reverse and green supply chains to minimize energy consumption and environmental impacts.Considering the queuing system with limited resources in hubs in this research was recognized as an innovation.Tey used the genetic algorithm to solve the proposed model in large dimensions.Sadeghi Rad et al. [26] provided a multilayer, multiperiod, multiproduct mixed-integer programming model, which included four layers in the forward (direction) fow (suppliers, production and regeneration centers, distribution center, and customers) and three layers in (backward) reverse (direction) fow (customers, inspection and collection centers, and disposal centers).Te production and inspection centers were integrated into the proposed model to reduce transportation costs.Besides economic goals, environmental aspects such as green production, technology, and transportation modes were also considered in this model.Also, the amount of raw materials purchased and the volume of greenhouse gases produced in the production process were considered dependent on the level of technology.Hajiaghaei et al. [27] provided a nonlinear mixed-integer mathematical program model to formulate a SCLSC with consideration of discounts on shipping costs.Tey suggested three hybrid RDSA, KAGA, and ICTS algorithms to solve the model, which were compared by four evaluation criteria by Pareto analysis.Te comparison result indicated that the proposed new hybrid KAGA algorithm brings better solutions compared to other algorithms but needs more time for solving.Tey fnally introduced a real example in the glass industry to confrm the proposed model and the algorithms provided.Ghomi-Avili et al. [28] presented an MO model in the green closed-loop supply chain by considering the failure of downtime of centers.Pricing of products with a collaborative approach to game theory was one of the innovations of this research.Te problem was fuzzily modeled due to the fuzzy nature of the data, whose most important goal was set to minimize the amount of pollutant gases released in the proposed chain.Rahimi and Ghezavati [29] presented an MO mathematical model for managing inventory and the fow of goods in a SCLSC.
Considering discounts for customers associated with uncertainty in demand were among the innovations in their research.Teir main goals were to minimize transportation costs and environmental costs.Te results revealed an 11% reduction in costs after implementing the model.Wang and Gunasekaran [30] provided a closed-loop supply chain where three competitive scenarios implemented by the manufacturer were studied.Tey used Stackelberg's game theory for studying this model.According to their conclusion, a producer is more inclined to perform the recycling and reproduction processes by himself and refrains from outsourcing.

Research Gap.
Te research gap can be summarized as follows according to the subject literature and Table 1: (1) Insufcient attention to the inherent uncertainty of supply chain issues together with their elements (2) Lack of attention to statistical reliability and various errors in the proposed mathematical models (3) Insufcient attention to the strengths of other metaheuristic algorithms aimed at utilizing them ogether in a new hybrid metaheuristic algorithm

Research Contributions
(1) In this model, increasing the supply chain reliability, the quality of the products returned by the customers, and the routing of the goods are also considered besides the three aspects of sustainability, namely, the social efect such as job creation, customer satisfaction, and distributors, the environmental efect such as reducing air pollution, and the economic efect such as reducing the supply chain costs.(2) In the real world, all the components of a closed-loop supply chain like production centers, etc., may not fully operate and are likely to stop working due to events such as human errors, weather conditions, terrorist attacks, etc., in period t.Tus, this limitation has been considered in the form of using statistical reliability to approach the real situation.(3) Providing a hybrid metaheuristic algorithm to solve the model.

Problem Statement
Te SCN provided in this study is a closed-loop, MO, and multiperiod scenario-based network, which encompasses suppliers, manufacturers, distributors, customers, hubs, repair, recycling, landfll, and demolition centers.It should be noted that the hubs, distribution centers, recycling centers, repair, and burial centers are equipped with the capability to be reopened.In the reopening process, some places are selected as candidate points and the model chooses the optimal place from them.In the supply chain functioning process, suppliers provide raw materials to manufacturers.Te manufacturers and the warehouses under their supervision send the products to distributors.Te breakdown and failure of raw materials supply centers, warehouses, and production centers are some of the issues and problems addressed in this study.Tis focus makes the issue closer and more similar to the real world.Distributors send products to customers.Tere are hubs, recycling, landfll, and repair centers in the backward direction.Te important point is how to determine the fow of returned products is their quality.Te hubs send the returned products to repair, landfll, or production centers depending on their quality.After repairing, the repair centers send the products to distributors.Te recycling centers also send products to manufacturers.It should be noted that some of the returned goods from customers can be sent directly to the reproduction centers located in the production center and do not need recycling in the recycling centers.Figure 1 shows the proposed supply chain structure.Maximizing the social responsibility dimension is one of the goals of the proposed model, in which, the employment rate, customer satisfaction, and distribution centers are maximized by sending maximum products from raw materials supply and production centers to them.Minimizing the economic and environmental costs of the supply chain is set to be another goal of this model.Tis parameter is considered to be fuzzy due to the inherent uncertainty of demand.Locating, examining the fow rate between components, and the routing of goods are other decisions that are supposed to be made in this research.

Model Assumptions
(1) Te capacity of the centers is limited.
(2) Te demand parameter is assumed to be fuzzy.
(3) Te distances between centers is assumed to be fxed and defnite.(4) Every producer has a warehouse to store the produced goods.(5) Tere is a reproduction center in each production center.(6) Te produced goods are sent both from production centers and from their warehouses to the distribution centers.(7) Te locations of supplier centers, producers, and their warehouses are fxed and predetermined and already known.(8) Any production center and its warehouse and supply centers cannot be rehabilitated and reconstructed in the case of being destroyed by an accident.

Objective Functions.
We know that the time it takes for the warehouse of the production center j to break down over a period of T j follows an exponential distribution with a mean value of λ jt ′ .Tus, the reliability of the warehouse of the production center j in sending products to the distribution center k in the period t is equal to: Terefore, the mean value of the products sent from the warehouse of the production centers to distribution centers is equal to Terefore, the mean value of the products sent from production centers and their warehouses to distribution centers is equal to Similarly, the average of raw materials sent from supplier centers to production centers is equal to Te frst objective function represents the social responsibility dimension of the supply chain network, where J 1 : Te number of fxed jobs created.J 2 : Te number of variable jobs created.R: Te average amount of product fow sent from suppliers, production centers, and their warehouses (as the R value increases, the satisfaction level of distribution centers and customers will increase).

Complexity
In: Te work injury rate due to jobs created in the production centers and their warehouses. Min ( Te second objective function encompasses the supply chain costs as follows:  Min Te third objective function includes environmental efects and the amount of CO 2 emitted due to the establishment of potential centers (F 1 ), and transportation is in the supply chain (F 2 ).
12)-( 16) suggest that a potential center should frst be established to subsequently create a route (path) to this potential center.Te constraint (17) indicates the satisfaction of customers' demands.Te constraint (18) shows the balance between distribution, customer, and hub centers.Te constraint (29) indicates the balance between customer, hub, and production centers.Te constraint (20) shows the balance between customer, hub, and recycling centers.Te constraint (21) indicates the balance between customer, hub, and landfll and disposal centers.Te constraint (22) shows the balance between customer, hub, and repair centers.Te constraint (23) indicates the balance between manufacturer, distributor, customer, and repair centers.Te constraint (24) shows the balance between hub, recycling, and manufacturing centers.Te constraint (25) shows the balance between hub, repair, and distributor centers.Te constraints ( 26) and ( 27) determine the product fows in terms of its quality.Te constraints ( 28)- (37) suggest that there is at least one path between the supply chain elements.Te constraints (38) and (39) indicate the amount of inventory and the fnal capacity of the manufacturer.Te constraint (40) determines the amount of the producer's inventory in its

Converting the Indefnite Demand Constraint to Its Equivalent Defnite Constraint.
A fuzzy number is a generalization of ordinary and real numbers that refers not to a value but a connected set of values so that any possible value would have a weight between 0 and 1 [31].Tis mentioned weight is called the membership function.Hence, a fuzzy number is a certain type of the normalized convex real line fuzzy set [32].Tere are several patterns such as triangular, trapezoidal, bell-shaped patterns, etc. to describe fuzzy numbers.Triangular fuzzy numbers are the most popular type of fuzzy numbers and are widely used in representing uncertainty in applied sciences because of their ability to express the perception of experts [33].We used the triangular fuzzy number in this article to represent the fuzzy demand parameter.
We utilized the weighted-average method in this paper to convert the fuzzy demand parameter to the equivalent defnite parameter.Terefore, the fuzzy constraint [5] is determined based on equation.
In equation ( 57), w 1 + w 2 + w 3 � 1 and β are the minimum acceptable likelihood in converting the fuzzy parameter to an equivalent real number.Te symbols w 3 , w 2 , w 1 show the most pessimistic, most possible, and the most optimistic weight of fuzzy demand values, respectively.Appropriate values for these weights are usually determined based on the experience and knowledge of the decisionmaker.In this study, these weights were considered based on the proposed values of Lai and Hwang and other conducted studies [34,35] as β � 0.5, w 1 � w 3 � 1/6, w 2 � 4/6.

Solving Approaches
5.1.Epsilon Constraint.Tis method is based on turning an MO optimization problem into a single-objective optimization problem.Tus, one of the objectives of the problem is optimized as the main objective regarding other objectives as a constraint in this method [36].It is assumed that the decision to minimize the objective functions (58) is associated with the constraints (59).
One of the objective functions is selected as the main objective function according to equation (60) based on this approach.Other objective functions are considered as the constraint (61), and each time, the problem is solved according to one of the objective functions, and optimal and corresponding values of each objective function are calculated.Te range between two optimal and corresponding values of subobjective functions is subdivided into a predetermined number followed by specifying a table of values for ε j .Finally, the Pareto solutions will be obtained [37]. st:

MOPSO Algorithm.
Te MOPSO algorithm was introduced by Coello in 2004 [38], which is in fact a generalization of the PSO algorithm that is used to solve MO problems.An elitist policy is used in this algorithm to keep the superior and dominant results in the iterations of the algorithm.Te dominant solutions are stored in an external archive.Selecting Pbest and Gbest is done by a special mechanism.Pbest is only updated when a new particle dominates its previous value; then, the new particle is replaced.
If the new particle is defeated by the best previous particle, nothing occurs and the same previous particle is introduced as the solution.Otherwise, if neither of them defeats each other, one of them is randomly considered as the best particle [39].
Gbest is also selected from the nondominated solutions in the archive in each iteration.Te proper selection of optimal solutions in the MO particle swarming algorithm is known as an important step since it is subject to the convergence of the algorithm and its ability to achieve a diverse set of nondominated solutions in the decision space.If the termination conditions for the algorithm are not met, the above steps are repeated.Otherwise, the set of solutions in the archive are presented as efcient front points and fnal solutions (Figure 2).

Tuning the Coefcients.
Te equations describing the behavior of the particles in the MOPSO algorithm are as follows.Equations ( 62) and (63) determine the velocity and location of the particle i at the moment t + 1.
where x i [t]: the location of the particle i at moment t.V i [t]: the velocity of the particle i at moment t.
x Gbest [t]: the best location of the particle i at moment t. w: the coefcient of inertia r 1 , r 2 : Rand(0, 1).c 1 , c 2 : individual and collective learning coefcients.

NSGA-II Algorithm. Tis algorithm was introduced by
Deb in 2002.In this algorithm, frst, the initial population of parents P 0 is randomly generated with the size N. Te generated initial population is then sorted out by a nondominated approach and the solutions are categorized into diferent levels of degree of nondomination.Subsequently, each solution is assigned with a value or rank depending on its position on that front.Hence, the solutions of the frst front, which are at the lowest level, are ranked frst, the solutions of the second front are ranked second, and so are the rest.Ten, the population of children Q 0 is generated by the size N using the binary tournament method based on the swarm comparison operator and the intersection and jump operators.Combining two populations of parents and children generates the R 0 � P 0 ∪ Q 0 population with the size 2N.Te next generation is selected from the obtained R 0 population.Since this algorithm follows an elitist approach, the members of R 0 are categorized again by a nondominated sorting method.After creating diferent fronts of the Complexity nondomination degree, the next generation P 1 population with the size N will be flled in order from the frst front onwards by in the following approach.By creating P 1 , the same steps mentioned for the P 0 are performed and this loop continues until the algorithm is terminated.Ultimately, the frst front achieved from the R t sorting of the latest generation will be introduced as the set of Pareto solutions.We assume that the R t population resulting from combining the parents P t and their children Q t has been sorted by a nondominated approach the fronts F i have been created for i � 1, 2, • • •.Now, the solutions found on the frst front F 1 are the frst candidates to join the next generation as the best solutions available in the current generation.If the number of the F 1 members is lower than the N, all of them will be transferred to P t+1 .Te rest of the P t+1 members are selected from the F 2 and then from the F 3 , and so on.Tis procedure is continued until the F 1 would be considered as the last front from which the rest of the P t+1 members are supposed to be selected.At this time, since the total number of the F 1 members is higher than the required number of the remaining members, these members are sorted in order of decreasing the crowding distance, and then the remaining required number will be selected from the beginning of this list.
Te most signifcant distinguishing feature of the MO genetic algorithm is the concept of crowding distance, as shown in Figure 3. Tis concept is used for determining the quality of a solution in a particular Pareto frontier.In fact, NSGA-II ensures the quality of the fnal unsuccessful solutions by using the crowding distance mechanism.
In which n represents the number of objective functions.In addition, f max i , f min i represent the highest and lowest values of the i-th and k objective functions of the current solution, respectively.

Gray Wolf Algorithm.
Te gray wolf optimization algorithm is one of the new metaheuristic algorithms proposed by Mirjalili et al. [40]; being inspired by the hunting behavior of gray wolves in the wild.In order to model the mathematical relations of the hierarchical structure of the wolf community for designing the GWO, the fttest solution is considered as alpha (α) wolf.Furthermore, the second and third optimal solutions are considered beta (β) and delta (δ), respectively.Other candidate solution are called omega (ω).In the GWO algorithm, hunting is guided by α, β, and δ.Wolves ω follow these three wolves in the hope of fnding an optimal solution.In addition to leading the wolf community, the following relations state the simulation of the siege behavior of gray wolves while hunting prey.
In these relations, a → decreases linearly from 2 to 0.
Te GWO algorithm uses community leadership simulation, as well as the siege mechanism to fnd the optimal solution to optimization problems.Tis algorithm selects the frst three solutions from the desired solution and requires other search factors including omega wolves to improve their position over the top solutions.Te following relations are applied to each factor during the optimization process to perform the hunting process and thus fnd the appropriate areas in the search space.p g (t) x i (t) Te mechanism of action of the particle swarm optimization algorithm.
12 Complexity Te search in the problem space is guaranteed by the vector A → with random values between [-1,1], necessitating the factors to distance themselves from the prey.Another GWO parameter which helps to search further is the vector C → .Te vector C → selects random values between [0, 2], leading to random weight values for the prey.Tese values show the efect of prey on the defnition of distance in (64).Te GWO algorithm begins the optimization process by generating a set of initial random solutions.During the optimization process, three of the best solutions are saved (α, β, δ).Ten, the wolves' situation is updated through equations ( 67)-(69).Meanwhile, the values of A and a increase linearly with increasing repetitions.Tus, wolves will tend to distance from the prey and get closer to the prey.Ultimately, the position and value of the alpha solution are returned to the algorithm as the best solution found in terms of all constraints.A distinctive feature of the gray wolf algorithm is the use of several guides and the avoidance of falling into the optimal local trap.
Figure 4 shows the multisegment chromosome presented in this study.For example, the fgure on the right shows the chromosome related to the Inv ct j,s variable.Te value within each gene represents the amount of residual inventory of the product c in the warehouse of the production center j in the period t in scenario s.Te fgure on the left also shows the chromosome corresponding to the variable x t p .If the value within each gene is equal to 1, the recycling center at location p is established in period t.
Figure 5 illustrates the intersection operator.As seen in the fgure, the single-point intersection has been used in this research.In this type of intersection, the two sides of the chromosome will be displaced together.
As shown in Figure 6, the mutation operator considered in this study is a reverse mutation type.In this type of mutation, a row of chromosomes is selected and will be then reversed.

Proposed Hybrid Metaheuristic Algorithm GW-NS.
In order to implement the optimization process of the proposed method based on two algorithms NSGA-II and GWO, two ideas were hybridised with each other.
Te frst idea is the strategy of selecting a leader from the external archive of nondominated solutions.
In the proposed algorithm, the value of the crowding distance is calculated for each member in the external archive and using the roulette wheel mechanism, we select the three members alpha (α), beta (β), and delta (δ).Members of the external archive that have more crowding distance should have a more chance of being selected.In other words: p i : probability of selecting the i-th element d i : crowding distance of the i-th element in the external archive.
Te second idea is the control process, when a nondominated solution aims to enter the external archive or when the number of archive members is higher than its capacity (it should be noted that there is always a certain number of members for the external archive).
In the proposed algorithm for deleting a nondominated solution, when the number of archive members exceeds its capacity, some members of the archive should be selected for deletion.Te members with less crowding distance should have a more chance of being selected for deletion.Te selection process is done by roulette wheel mechanism.
In general, the steps of the proposed algorithm are as follows (Figure 7): (1) Defning a random primary population in the problem search space.(2) Defning an external archive of nondominated members of the primary population.(3) Selecting three members alpha (α), beta (β), and delta (δ) from the members in the external archive randomly, using the roulette wheel mechanism, based on the crowding distance In which:  8) If the number of members of the external archive is higher than the determined capacity, extra members will be removed.Note: Deletion is conducted by the roulette wheel mechanism and based on the crowding distance index, i.e., the smaller the crowding distance, the more likely it is to be eliminated.
In which: H: A total of 50% better total external archive members based on crowding distance index.
(9) If the termination conditions are met, go to the next stage, otherwise return to stage 3. (10) Te End.

Tuning the Parameters.
Since the output of the problems strongly depends on the parameters of the proposed algorithms, thus, we used the Taguchi method to tune their parameters.Te advantage of the Taguchi method over other tests design methods in addition to cost is to obtain the optimal levels of parameters in less time [41].Choosing an orthogonal array appears to be one of the most important steps, which estimates the efects of factors on the mean values of the solution and variations.Te most appropriate test design in this research was found to be the three-level experiments at three low, moderate, and high levels.Ten, the arrays L 9 and L 27 were chosen as the suitable test design to tune the parameter of the proposed algorithms due to the Taguchi standard orthogonal arrays.Te levels of the parameters of NSGA-II and MOPSO algorithms are given in Table 2 for their tuning.
A statistical measure of performance, known as the signal to noise ratio (S/N) is considered in the Taguchi method to tune the optimal parameters.Tis ratio encompasses mean values and variations, and it would be more desirable at higher levels.Te solution variable considered in this study is the ratio of the two indices of the Mean Ideal Distance (MID, equation ( 77)) and the Diversifcation Metric (DM, equation ( 78)) for MO algorithms.Since this solution variable is of the type "the lower, the better", then, its corresponding S/N ratio is considered as equation (76).Te proposed metaheuristic algorithms are implemented for each Taguchi experiment, and then, the S/N ratios will be calculated by the Minitab 19.2020.1 software.

Results Analysis and Comparisons
Te model is frst solved in small and medium dimensions aimed at evaluating the accuracy and precision of the 14 Complexity proposed model.Table 3 shows the problems with small dimensions (samples 1-5) and medium dimensions (samples 6-10).For example, there is a customer, a manufacturer, a supplier, a recycling center, a hub, a repair center, a distributor, and a landfll center in sample number 1.
Figure 8 shows the model solution time for NSGA-II, MOPSO, GW-NS algorithms, and the Epsilon constraint For each member of the current population, three alpha members are selected randomly using the Roulette wheel mechanism based on crowding distance and then updating its situation using Eqs.78-80 and creating a new population Complexity method.As can be seen, the solution time by the Epsilon constraint method is increasing exponentially as the dimension of the problem increases.Algorithms are programmed using GAMS and MATLAB (R2020 a) and implemented on a PC under Windows 7, Intel Core i3, 3.3 GHz, 4 GB RAM.

Multiobjective Performance Metrics.
Te performance of the algorithms used to solve the proposed model was evaluated by the following criteria, which are described below.

Mean Ideal Distance (MID). Tis criterion measures
the degree of closeness between the solutions found on the Pareto front and the ideal points (f best 1 , f best 2 , f best 3 ), which is calculated by equation (77).Te lower the value of this criterion for an algorithm, the better the performance of that algorithm would be.16 Complexity

MID
In equation (77), n is the number of nondominated solutions, while f best i and f nadir i are the best and worst values of the i th objective function, respectively [42].

Diversifcation Metric (DM)
. Tis criterion measures the scattering of the Pareto solutions, which is calculated by the following equation.
According to equation (78), a higher value of DM indicates the better performance of the algorithm [42].

Spacing Metric (SM)
. Tis criterion measures the scattering pattern of the nondominated solutions, which is calculated by the following equation.
In equation (79), d i determines the Euclidean distance between successive solutions in the set of the nondominated solutions obtained by the algorithm and d is the average of these distances.According to the defnition of SM, the lower the value of this index, the better the algorithm would be [42].
In order to compare the results of the algorithms, three metaheuristics are applied to solve 10 problems in various dimensions, and the obtained results are reported in Table 4.For a closer look at the results of the three algorithms, the following hypothesis was tested according to DM, SM, and MID indexes.

Hypothesis 1.
Tere is a signifcant diference between the DM of the solutions generated by the three algorithms GW-NS, MOPSO, and NSGA-II.

Hypothesis 2.
Tere is a signifcant diference between the MID of the solutions generated by the three algorithms GW-NS, MOPSO, and NSGA-II.Hypothesis 3. Tere is a signifcant diference between the SM of the solutions generated by the three algorithms GW-NS, MOPSO, and NSGA-II.
Table 4 shows the values of MID, DM, and SM indices for the problems listed in Table 3.
Te MID and SM indices, respectively, increase with low and high slopes as the size of the problem enhances according to Figures 9 and 10.To put it better, increasing the dimensions of the problem reduces the efciency of algorithms in terms of MID and SM indices, while the MID index shows less sensitivity to the increase in dimensions.According to Figure 11, the increased size of the problem enhances the efciency of the algorithms based on the DM index.In other words, increasing the dimensions of the problem exhibited a higher potential for exploring and extracting the region of the solution.

Statistical Analysis Comparisons.
For the purpose of statistical analysis, one-way variance analysis technique and SPSS software are applied.As well, to confrm the parametric results, a nonparametric test called Kruskal-Wallis test was used.If data are suitable for variance analysis, nonparametric test is used where there is no precondition for uniformity of the variance or normal distribution.Results of one-way variance analysis and nonparametric test for three measures are provided in Tables 5-8.
Based on Tables 5-8 for all three indices SM, DM, and MID, the P value of ANOVA and Kruskal-Wallis tests are all more than 0.05, Tus, there is no signifcant diference    .17Complexity algorithm also ofered higher-diversity solutions(DM) than the NSGA-II and MOPSO algorithms.In other words, this algorithm exhibited higher potentials for exploring and extracting the region of the solution as compared with the NSGA-II and MOPSO algorithms.In addition, in terms of SM index, NSGA-II, MOPSO, and GW-NS algorithms indicate better performance, respectively.In other words, the NSGA-II algorithm generated more uniform solutions as compared with the GW-NS and MOPSO algorithms.Box plots Figures show that the variability of the DM index is almost the same in all three algorithms.But in the MID index, the GW-NS algorithm and in the SM index, the MOPSO algorithms have the highest changes and less stability.

Conclusion and Future Works
In this study, a new multiobjective, multiperiod, multiproduct scenario-based fuzzy mathematical model for CLSC was presented.In the proposed model, in addition to the three aspects of sustainability, namely social, economic and environmental impact, reliability, quality of returned products by customers, and routing of goods fow in the supply chain were considered.In this network, the demand parameter was considered as fuzzy.In order to solve the proposed model, in addition to the two metaheuristic algorithms NSGA-II and MOPSO, a new hybrid metaheuristic algorithm using the strengths of the GWO algorithms (using multiple guides and avoiding falling into the optimal local trap) and NSGA-II (guaranteeing the quality of unused solutions crowding distance mechanism) was presented.After tuning their parameters by Taguchi method, their performance in problems with diferent dimensions was tested and evaluated by MID, DM, and SM criteria.Te results of statistical analysis of indices indicated that no signifcant diference between the performance of the three algorithms at 5% error level.In general, GW-NS, NSGA-II, and MOPSO algorithms had better performance in terms of MID index, respectively.In addition, GW-NS, NSGA-II, and MOPSO algorithms performed better in terms of DM index.NSGA-II, MOPSO, and GW-NS algorithms performed better in terms of SM index, respectively.Te model proposed in this paper can be developed in future research by considering several decision-makers in the closed loop network and the use of the concept of the game.Te robust planning approach can be also used to deal with the supply chain uncertainty, making the model more powerful and fexible in the face of uncertainty.Also, instead of the exponential distribution assumed to model the reliability of the direct logistics elements, other probability distributions such as Erlang or Weibull can be considered.Multiple fnancial, environmental, and social goals combined with dynamic constraints can provide more efective and practical solutions.At the end, identify the strengths of other metaheuristic algorithms with the aim of using them to propose new hybrid algorithms.Combining two or more performance indicators of metaheuristic algorithms with each other and using them to compare algorithms can be considered for future study.

Abbreviations
Quality level of products that are sent from hub stations to repair stations Q 2 : Quality level of products that are sent directly from the hub stations to the production stations Missed working days due to work injury in the production station j in the term t Variables ro rt vij,s : Equivalent to 1 if the vehicle of type v goes from the supplier i to the manufacturer j from the route r in the term t in scenario s, otherwise 0 ro rt vjk,s : Equivalent to 1 if the vehicle of type v goes from the manufacturer j to the distributor k from the route r in the term t in scenario s, otherwise 0 Ro rt vjk,s : Equivalent to 1 if the vehicle of type v goes from the warehouse of the manufacturer j to the distributor k of the route r in the term t in scenario s, otherwise 0 ro rt vkl,s : Equivalent to 1 if the vehicle of type v goes from the distributor k to the customer l from the route r in the term t in scenario s, otherwise 0 ro rt vlm,s : Equivalent to 1 if the vehicle of type v goes from the customer l to the hub m from the route r in the term t in scenario s, otherwise 0 ro rt vmp,s : Equivalent to 1 if the vehicle of type v goes from the hub m to the recycling station p from the route r in the term t in scenario s, otherwise 0 ro rt vpj,s : Equivalent to 1 if the vehicle of type v goes from the recycling station p to the manufacturer j from the route r in the term t in scenario s, otherwise 0 ro rt vmj,s : Equivalent to 1 if the vehicle of type v goes from the hub m to the manufacturer j from the route r in the term t in scenario s, otherwise 0 ro rt vmn,s : Equivalent to 1 if the vehicle of type v goes from the hub m to the demolition station n from the route r in the term t in scenario s, otherwise 0 ro rt vmw,s : Equivalent to 1 if the vehicle type v goes from the hub m to the repair station w from the route r in the term t in scenario s, otherwise 0 ro rt vwk,s : Equivalent to 1 if the vehicle type v goes from the repair station w to the distributor k from the route r in the term t in scenario s, otherwise 0 z ct lmq,s : Flow amount of the returned yield c with the quality of the type q from the customer l to the hub m in the term t in scenario s z ct mpq,s : Flow amount of the returned yield c with the quality of the type q from the hub m to the recycling station p in the term t in scenario s z ot ij,s : Flow amount of the raw materials o from the supply station i to the production station j in the term t in scenario s z ct pj,s : Flow amount of the reused yield c from the recycling station p to the production station j in the term t in scenario s z ct mnq,s : Flow amount of the returned yield c with the quality of the type q from the hub m to the demolition station n in the term t in scenario s z ct mjq,s : Flow amount of the returned yield c with the quality of the type q from the hub m to the production station j in the term t in scenario s z ct mwq,s : Flow amount of the returned yield c with the quality of the type q from the hub m to the repair station w in the term t in scenario s z ct wk,s : Flow amount of the yield c from the repair station w to the distribution station k in the term t in scenario s z ct jk,s : Flow amount of the yield c from the production station j to the distribution station k in the term t in scenario s z ct kl,s : Flow amount of the yield c from the distribution station k to the customer l in the term t in scenario s Q ct jk,s : Flow amount of the yield c from the warehouse of the producer j to the distribution station k in the term t in scenario s Q ct jj,s : Flow amount of the yield c from the production station j to its own warehouse in term t in scenario s x t k : If the distribution station is established at the location k in the term t, its value would be equal to 1, otherwise 0 x t p : If the recycling station is established at the location p in the term t, its value would be equal to 1, otherwise 0 x t m : If the hub is established at the location m in the term t, its value would be equal to 1, otherwise 0 x t n : If the landfll and demolition station is established at the location n in the term t, its value would be equal to 1, otherwise 0 x t w : If the repair station is established at the location w in the term t, its value would be equal to 1, otherwise 0.

Figure 1 :
Figure 1: Te structure of the proposed model.

C 2 :
Te costs of production and maintenance of products in the production sector.C 3 : Te costs of recycling, demolition, etc.
warehouse.Te constraints (41)-(48) indicate the capacity of the fxed and potential centers.Te constraints (49)-(54) explain the limitations of vehicle routing.Finally, the constraint (55) represents the decision variables of the proposed model.

Figure 8 :Figure 10 :
Figure 8: Te model solution time in terms of increasing the problem dimension.

Figure 9 :
Figure 9: Te result of comparing the algorithms in terms of MID.

Figure 12 :
Figure 12: Mean changes in the performance index of SM algorithms.

Figure 13 :Figure 14 :Figure 11 :
Figure 13: Mean changes in the MID performance index of algorithms.

Table 1 :
Literature review table.
and p represent the most pessimistic, most possible, and the most optimistic values for the parameter, respectively.Terefore, the demand membership function would be as follows: Tus, we defned the demand fuzzy parameter as ((d ct ls ) p , (d ct ls ) m , (d ct ls ) o ), in which the upper indices o, m,

Table 2 :
Te algorithm parameter table.

Table 3 :
Te dimensions of the problem.

Table 4 :
Te computational results of comparison measurement criteria of MOPSO, NSGA-II, and GW-NS algorithms.

Table 5 :
Results of one-way variance analysis for MID, DM, and SM.
Quality level of products that are sent from hub stations to recycling stations Q 4 : Quality level of products that are sent from hub stations to disposal stations, (where Q 1 ≻Q 2 ≻Q 3 ≻Q 4 ≻ is the quality level of the comparison operator) : Transfer fee per unit of the comeback yield c with quality q from the hub station m to the repair station w in scenario s h c wk,s : Transfer fee per unit of the yield c from the repair station w to the distribution k in scenario s cq c jk,s : Transfer fee per unit of the yield c from the warehouse of the producer j to the distribution station k in scenario s cq c jj,s : Transfer fee per unit of the yield c from the production station j to its own warehouse in scenario s h c lmq,s : Transfer fee per unit of the returned yield c from the customer l to the hub station m in scenario s cj ct,s : Transfer fee per unit of the yield c in the production station j in the term t in scenario s rcj ct,s : Transfer fee per unit of the yield c in the reproduction station j in the term t in scenario s cn ct,s : Transfer fee per unit of the yield c at the demolition station n in the term t in scenario s cm ct,s : Cost of collecting and sending a unit of the yield c to the hub station m in the term t in scenario s cp ct,s : Cost of recycling a unit of the yield c at the recycling station p in the term t in scenario s cw ct,s : Cost of repairing a unit of the yield c in the repair station w in the term t in scenario s ck ct,s : Cost of distributing a unit of the yield c in the distribution station k in the term t in scenario s αj m,s : Fixed number of permanent job founded if established the hub station m in scenario s αj p,s : Fixed number of permanent job founded if established the recycling station p in scenario s αj w,s : Fixed number of permanent job founded if established the repair station w in scenario s αj n,s : Fixed number of permanent job founded if established the demolition station n in scenario s αj k,s : Fixed number of permanent job founded if established the distribution station k in scenario s βj m,s : Variable number of permanent job founded if established the hub station m in scenario s βj p,s : Variable number of permanent job founded if established the recycling station p in scenario s βj w,s : Variable number of permanent job founded if established the repair station w in scenario s βj n,s : Variable number of permanent job founded the demolition station n in scenario s βj k,s : Variable number of permanent job founded if established the distribution station k in scenario s dl t j : s : Return rate of the yield c with the quality q from the customer station l in the term t in scenario s rb ct mj : Return rate of the yield c from the hub station m to the production station j in the term t rb ct mn : Return rate of the yield c from the hub station m to the landfll and disposal station n in the term t rb ct mp : Return rate of the yield c from the hub station m to the recycling station p in the term t rb ct mw : Return rate of the yield c from the hub station m to the repair station w in the term t fp ct lmq,s : Return cost of each unit of the returned yield c with quality q from the customer station l to the hub station m in the term t in scenario s s : Transfer fee per unit of the comeback yield c with quality q from the hub station m to the recycling station p in scenario s