Abstract
In this study, a new nonlinear mathematical programming model of mixed integer was presented to formulate the problem of designing a sustainable closed loop supply chain, in which the three aspects of sustainability, i.e., social effect such as job creation, customer satisfaction, and distributors, environmental effects such as reducing air pollution, and economic effects such as reducing supply chain costs, increasing supply chain reliability, quality of returned products by customers, and product routing were considered. In order to solve the proposed model, a new hybrid metaheuristic algorithm based on the distinctive features of gray wolf algorithm and genetic algorithm was proposed in addition to MOPSO and NSGA-II algorithms. After tuning their parameters by the Taguchi method, their performance in problems with different dimensions was tested and evaluated by MID, DM, and SM criteria. The results of statistical analysis of indices indicated that no significant difference 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. In addition, the variability of DM index in all three algorithms was almost the same, but in MID index, GW-NS algorithm, and in SM index, MOPSO algorithm had the highest change and less sustainability.
1. 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]. The 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 efficiency of energy, and resources consumption associated with improving the performance of the entire supply chain [2, 3]. The 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. The 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. The customers’ demand for high quality and fast service-providing 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. Thus, 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. The 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 scientific 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 affect the quality of their performance [11]. In the real world, utilities may suffer from disruptions and failures with possible causes of human error, natural disasters, etc. [12, 13]. Thus, the failure of one component of the SCN may disrupt the functioning of the entire supply chain, or in the best case, reduce the efficiency 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. The 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.
2. Literature Review
The 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. The goals of this model encompass maximum profitability according to the income and costs of the entire chain, minimizing the environmental impacts due to the carbon index and maximizing the social effects according to the job opportunities created. They utilized the weighted sum technique to solve the model. The 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. The model objectives included maximizing transaction level satisfaction and the profit from the supply chain. They 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. The 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. They employed an NSGA-II algorithm and a cokriging approach to solve the model. The study results revealed the positive effects of motivational pricing on the returned goods. Kazancoglu et al. [17] presented an MILP model for designing the green dual-channel and CLSCs. The 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. The 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. This 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 financial cash flow as cash flow and debt constraints and employment under uncertainty. The objectives of the proposed mathematical model in their study included increasing the cash flow, maximizing the jobs created, and maximizing the reliability of raw materials consumed. They developed and used MO simulated annealing, MO invasive weed optimization, and MO gray wolf metaheuristic algorithms to solve the model at a large scale. Ali et al. [21] provided a novel mathematical model for reverse supply chain management of air conditioning products. Their considered supply chain was sustainable with fuzzy demand uncertainty. Locating hub and recycling centers was among the most important goals of their research. The case study covered the industries of Saudi Arabia and India. In their study, identified 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. The objectives regarded by them included minimizing total costs, minimizing the amount of carbon dioxide released, and maximizing social impacts. Their innovation involved considering three types of distribution channels, including normal delivery, direct delivery, and direct displacement. The 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. The most important decisions made in their model included determining the amount of flow at each level and locating the hubs. The main objectives of the study contained minimizing transportation costs and the costs of carbon dioxide emissions and maximizing social responsibility. The 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. The proposed innovation included a variety of technologies for cars, which causes the release of different amounts of carbon dioxide from different vehicles. They 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. They considered three different prices for different 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. They 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) flow (suppliers, production and regeneration centers, distribution center, and customers) and three layers in (backward) reverse (direction) flow (customers, inspection and collection centers, and disposal centers). The 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. They suggested three hybrid RDSA, KAGA, and ICTS algorithms to solve the model, which were compared by four evaluation criteria by Pareto analysis. The comparison result indicated that the proposed new hybrid KAGA algorithm brings better solutions compared to other algorithms but needs more time for solving. They finally introduced a real example in the glass industry to confirm 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. The 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 flow of goods in a SCLSC. Considering discounts for customers associated with uncertainty in demand were among the innovations in their research. Their main goals were to minimize transportation costs and environmental costs. The 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. They 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.
2.1. Research Gap
The research gap can be summarized as follows according to the subject literature and Table 1:(1)Insufficient 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)Insufficient attention to the strengths of other metaheuristic algorithms aimed at utilizing them ogether in a new hybrid metaheuristic algorithm
2.2. 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 effect such as job creation, customer satisfaction, and distributors, the environmental effect such as reducing air pollution, and the economic effect 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. Thus, 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.3. Problem Statement
The SCN provided in this study is a closed-loop, MO, and multiperiod scenario-based network, which encompasses suppliers, manufacturers, distributors, customers, hubs, repair, recycling, landfill, 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. The manufacturers and the warehouses under their supervision send the products to distributors. The breakdown and failure of raw materials supply centers, warehouses, and production centers are some of the issues and problems addressed in this study. This focus makes the issue closer and more similar to the real world. Distributors send products to customers. There are hubs, recycling, landfill, and repair centers in the backward direction. The important point is how to determine the flow of returned products is their quality. The hubs send the returned products to repair, landfill, or production centers depending on their quality. After repairing, the repair centers send the products to distributors. The 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. This parameter is considered to be fuzzy due to the inherent uncertainty of demand. Locating, examining the flow rate between components, and the routing of goods are other decisions that are supposed to be made in this research.
4. Formulating the Model
4.1. Model Assumptions
(1)The capacity of the centers is limited.(2)The demand parameter is assumed to be fuzzy.(3)The distances between centers is assumed to be fixed and definite.(4)Every producer has a warehouse to store the produced goods.(5)There is a reproduction center in each production center.(6)The produced goods are sent both from production centers and from their warehouses to the distribution centers.(7)The locations of supplier centers, producers, and their warehouses are fixed 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.4.2. Objective Functions
We know that the time it takes for the warehouse of the production center j to break down over a period of follows an exponential distribution with a mean value of . Thus, 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:
Therefore, the mean value of the products sent from the warehouse of the production centers to distribution centers is equal to
Therefore, 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
The first objective function represents the social responsibility dimension of the supply chain network, where : The number of fixed jobs created. : The number of variable jobs created. : The average amount of product flow sent from suppliers, production centers, and their warehouses (as the R value increases, the satisfaction level of distribution centers and customers will increase). : The work injury rate due to jobs created in the production centers and their warehouses.
The second objective function encompasses the supply chain costs as follows: : The costs of establishing hub, repair, distribution, demolition, landfill, and recycling centers. : The costs of production and maintenance of products in the production sector. : The costs of recycling, demolition, etc.
The third objective function includes environmental effects and the amount of CO2 emitted due to the establishment of potential centers (), and transportation is in the supply chain ().
4.3. Constraints
The constraints (7)–(11) suggest that at least one of the potential distribution, repair, hub, recycling, and disposal centers in the supply chain is established. The constraints (12)–(16) suggest that a potential center should first be established to subsequently create a route (path) to this potential center. The constraint (17) indicates the satisfaction of customers’ demands. The constraint (18) shows the balance between distribution, customer, and hub centers. The constraint (29) indicates the balance between customer, hub, and production centers. The constraint (20) shows the balance between customer, hub, and recycling centers. The constraint (21) indicates the balance between customer, hub, and landfill and disposal centers. The constraint (22) shows the balance between customer, hub, and repair centers. The constraint (23) indicates the balance between manufacturer, distributor, customer, and repair centers. The constraint (24) shows the balance between hub, recycling, and manufacturing centers. The constraint (25) shows the balance between hub, repair, and distributor centers. The constraints (26) and (27) determine the product flows in terms of its quality. The constraints (28)–(37) suggest that there is at least one path between the supply chain elements. The constraints (38) and (39) indicate the amount of inventory and the final capacity of the manufacturer. The constraint (40) determines the amount of the producer’s inventory in its warehouse. The constraints (41)–(48) indicate the capacity of the fixed and potential centers. The constraints (49)–(54) explain the limitations of vehicle routing. Finally, the constraint (55) represents the decision variables of the proposed model.
4.4. Converting the Indefinite Demand Constraint to Its Equivalent Definite 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]. This mentioned weight is called the membership function. Hence, a fuzzy number is a certain type of the normalized convex real line fuzzy set [32]. There 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. Thus, we defined the demand fuzzy parameter as , in which the upper indices o, m, and p represent the most pessimistic, most possible, and the most optimistic values for the parameter, respectively. Therefore, the demand membership function would be as follows:
We utilized the weighted-average method in this paper to convert the fuzzy demand parameter to the equivalent definite parameter. Therefore, the fuzzy constraint [5] is determined based on equation.
In equation (57), and are the minimum acceptable likelihood in converting the fuzzy parameter to an equivalent real number. The symbols 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 decision-maker. In this study, these weights were considered based on the proposed values of Lai and Hwang and other conducted studies [34, 35] as .
5. Solving Approaches
5.1. Epsilon Constraint
This method is based on turning an MO optimization problem into a single-objective optimization problem. Thus, 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. The range between two optimal and corresponding values of subobjective functions is subdivided into a predetermined number followed by specifying a table of values for . Finally, the Pareto solutions will be obtained [37].
5.2. MOPSO Algorithm
The 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. The dominant solutions are stored in an external archive. Selecting and is done by a special mechanism. 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]. is also selected from the nondominated solutions in the archive in each iteration. The 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 efficient front points and final solutions (Figure 2).
5.2.1. Tuning the Coefficients
The 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 : the location of the particle i at moment t. : the velocity of the particle i at moment t. : the best location of the particle i at moment t. w: the coefficient of inertia , : . , : individual and collective learning coefficients.
5.3. NSGA-II Algorithm
This algorithm was introduced by Deb in 2002. In this algorithm, first, the initial population of parents is randomly generated with the size . The generated initial population is then sorted out by a nondominated approach and the solutions are categorized into different 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 first front, which are at the lowest level, are ranked first, the solutions of the second front are ranked second, and so are the rest. Then, the population of children is generated by the size 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 population with the size . The next generation is selected from the obtained population. Since this algorithm follows an elitist approach, the members of are categorized again by a nondominated sorting method. After creating different fronts of the nondomination degree, the next generation population with the size will be filled in order from the first front onwards by in the following approach. By creating , the same steps mentioned for the are performed and this loop continues until the algorithm is terminated. Ultimately, the first front achieved from the sorting of the latest generation will be introduced as the set of Pareto solutions.
We assume that the population resulting from combining the parents and their children has been sorted by a nondominated approach the fronts have been created for . Now, the solutions found on the first front are the first candidates to join the next generation as the best solutions available in the current generation. If the number of the members is lower than the , all of them will be transferred to . The rest of the members are selected from the and then from the , and so on. This procedure is continued until the would be considered as the last front from which the rest of the members are supposed to be selected. At this time, since the total number of the 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.
The most significant distinguishing feature of the MO genetic algorithm is the concept of crowding distance, as shown in Figure 3. This concept is used for determining the quality of a solution in a particular Pareto frontier. In fact, NSGA-II ensures the quality of the final unsuccessful solutions by using the crowding distance mechanism.
In which n represents the number of objective functions. In addition, represent the highest and lowest values of the i-th and k objective functions of the current solution, respectively.
5.4. Gray Wolf Algorithm
The 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 fittest 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 finding 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.where t represents iterations, shows coefficient vectors, represents prey position vectors, and shows gray wolf position vectors. The vectors and are calculated as follows:
In these relations, decreases linearly from 2 to 0. The GWO algorithm uses community leadership simulation, as well as the siege mechanism to find the optimal solution to optimization problems. This algorithm selects the first three solutions from the desired solution and requires other search factors including omega wolves to improve their position over the top solutions. The following relations are applied to each factor during the optimization process to perform the hunting process and thus find the appropriate areas in the search space.
The search in the problem space is guaranteed by the vector 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 . The vector selects random values between [0, 2], leading to random weight values for the prey. These values show the effect of prey on the definition of distance in (64). The 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 (α, β, δ). Then, the wolves' situation is updated through equations (67)–(69). Meanwhile, the values of A and a increase linearly with increasing repetitions. Thus, 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 figure on the right shows the chromosome related to the variable. The 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. The figure on the left also shows the chromosome corresponding to the variable . 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 figure, 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.
5.5. 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.
The first 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: : probability of selecting the i-th element : crowding distance of the i-th element in the external archive.
The 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. The members with less crowding distance should have a more chance of being selected for deletion. The selection process is done by roulette wheel mechanism.
In general, the steps of the proposed algorithm are as follows (Figure 7):(1)Defining a random primary population in the problem search space.(2)Defining 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: : A total of 50% better total external archive members based on crowding distance index.(4)Updating the position of each current member based on equations (67)–(69) and forming a new population.(5)Selecting a percentage of the new population of the previous stage (stage 4) for mutation and creation of a new population.(6)Determining the nondominated members of the new population produced in the previous stage (stage 5).(7)Adding the nondominated population created in the previous stage (stage 6) to the external archive and updating it.(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: : 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)The End.
5.6. 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. The 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 effects of factors on the mean values of the solution and variations. The most appropriate test design in this research was found to be the three-level experiments at three low, moderate, and high levels. Then, the arrays and were chosen as the suitable test design to tune the parameter of the proposed algorithms due to the Taguchi standard orthogonal arrays. The 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. This ratio encompasses mean values and variations, and it would be more desirable at higher levels. The solution variable considered in this study is the ratio of the two indices of the Mean Ideal Distance (MID, equation (77)) and the Diversification 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). The 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.
6. Results Analysis and Comparisons
The model is first solved in small and medium dimensions aimed at evaluating the accuracy and precision of the 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 landfill center in sample number 1.
Figure 8 shows the model solution time for NSGA-II, MOPSO, GW-NS algorithms, and the Epsilon constraint 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.
6.1. Multiobjective Performance Metrics
The performance of the algorithms used to solve the proposed model was evaluated by the following criteria, which are described below.
6.1.1. Mean Ideal Distance (MID)
This criterion measures the degree of closeness between the solutions found on the Pareto front and the ideal points , which is calculated by equation (77). The lower the value of this criterion for an algorithm, the better the performance of that algorithm would be.
In equation (77), is the number of nondominated solutions, while and are the best and worst values of the ith objective function, respectively [42].
6.1.2. Diversification Metric (DM)
This 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].
6.1.3. Spacing Metric (SM)
This criterion measures the scattering pattern of the nondominated solutions, which is calculated by the following equation.
In equation (79), determines the Euclidean distance between successive solutions in the set of the nondominated solutions obtained by the algorithm and is the average of these distances. According to the definition 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. There is a significant difference between the DM of the solutions generated by the three algorithms GW-NS, MOPSO, and NSGA-II.
Hypothesis 2. There is a significant difference between the MID of the solutions generated by the three algorithms GW-NS, MOPSO, and NSGA-II.
Hypothesis 3. There is a significant difference 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.
The 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 efficiency 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 efficiency 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.
6.2. Statistical Analysis Comparisons
For the purpose of statistical analysis, one-way variance analysis technique and SPSS software are applied. As well, to confirm 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 value of ANOVA and Kruskal—Wallis tests are all more than 0.05, Thus, there is no significant difference between the performances of the three algorithms at the 5% error level. Figures 12–14 indicate that GW-NS, NSGA-II, and MOPSO algorithms have better performance in terms of MID index, respectively. In other words, this algorithm exhibited higher potentials for convergence to the ideal solution compared to other algorithms. The GW-NS algorithm also offered 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 15–17 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.
7. 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 flow 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 different dimensions was tested and evaluated by MID, DM, and SM criteria. The results of statistical analysis of indices indicated that no significant difference 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. The 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. The robust planning approach can be also used to deal with the supply chain uncertainty, making the model more powerful and flexible 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 financial, environmental, and social goals combined with dynamic constraints can provide more effective 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
Indices| R: | Routes |
| L: | Customers |
| J: | Manufacturers |
| : | The set of quality levels , |
| : | Quality level of products that are sent from hub stations to repair stations |
| : | Quality level of products that are sent directly from the hub stations to the production stations |
| : | Quality level of products that are sent from hub stations to recycling stations |
| : | Quality level of products that are sent from hub stations to disposal stations, (where is the quality level of the comparison operator) |
| I: | Suppliers |
| O: | The set of raw materials |
| P: | Candidate points for recycling stations |
| M: | Candidate points for hubs |
| W: | Candidate points for repair stations |
| : | The set of means of transportation |
| C: | Products set index |
| N: | Candidate points for landfill and disposal stations |
| K: | Candidate points of distribution stations |
| : | Set of scenarios |
| T: | Term index |
| : | Distance of the production station j from its own warehouse in scenario s |
| : | Distance from the customer station l to the hub m in scenario s |
| : | Distance of the hub station m from the production station j in scenario s |
| : | Distance of the distribution station k from the customer station l in scenario s |
| : | Distance of the supplier i from the production station of j in scenario s |
| : | Distance of the production station j from the distribution station k in scenario s |
| : | Distance of the hub m from the landfill canter n in scenario s |
| : | Distance of the recycling station p from the production station j in scenario s |
| : | Distance of the hub m from the recycling station p in scenario s |
| : | Distance of the hub m from the repair station in scenario s |
| : | Distance of the repair station from the distribution station k in scenario s |
| : | Distance of the producer warehouse j from the distribution station k in scenario s |
| : | CO2 emitted during the construction of the recycling station p |
| : | CO2 emitted during the construction of the distribution station k |
| : | CO2 emitted during the construction of the hub m |
| : | The CO2 emitted during the construction of the landfill and disposal station n |
| : | CO2 emitted during the construction of the repair station |
| : | CO2 emitted from the shipment of a yield unit by a type vehicle over one kilometer |
| : | CO2 emitted from transporting a unit of product from the production station j to its own warehouse |
| : | Valency of the recycling station p in the term t |
| : | Valency of the landfill and disposal station n in the term t |
| : | Valency of the distribution station k in the term t |
| : | Valency of the supply station i in the term t |
| : | Valency of the repair station in the term t |
| : | Valency of the producer’s j warehouse in the term t |
| : | Valency of the production station j in the term t |
| : | Valency of the hub m in the term t |
| : | Valency of remanufacturing products in the production station j in the term t |
| : | Fixed cost of constructing the landfill and disposal (demolition) station n in the term t |
| : | Fixed cost of constructing the hub station m in the term t |
| : | Fixed cost of constructing the distribution station k in the term t |
| : | Fixed cost of constructing the recycling station p in the term t |
| : | Fixed cost of constructing the repair station in the term t |
| : | Amount of demand for the yield c by the customer l in the term t in scenario s |
| : | Return rate of the yield c with the quality q from the customer station l in the term t in scenario s |
| : | Return rate of the yield c from the hub station m to the production station j in the term t |
| : | Return rate of the yield c from the hub station m to the landfill and disposal station n in the term t |
| : | Return rate of the yield c from the hub station m to the recycling station in the term t |
| : | Return rate of the yield c from the hub station m to the repair station in the term t |
| : | 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 |
| : | Maintenance cost of each unit of the yield c in the producer’s warehouse at location j in the term t in scenario s |
| : | Transfer fee per unit of the returned yield c with quality q from the hub station m to the landfill station n in scenario s |
| : | Transfer fee per unit of the yield c from the distribution station k to the customer station l in scenario s |
| : | Transfer fee per unit of the raw materials o from the supply station i to the production station j in scenario s |
| : | Transfer fee per unit of the yield c from the production station j to the distribution station k in scenario s |
| : | Transfer fee per unit of the comeback yield c with quality q from the hub station m to the production station j in scenario 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 |
| : | Transfer fee per unit of the comeback yield c from the recycling station p to the production station j in scenario s |
| : | Transfer fee per unit of the comeback yield c with quality q from the hub station m to the repair station in scenario s |
| : | Transfer fee per unit of the yield c from the repair station to the distribution station k in scenario s |
| : | Transfer fee per unit of the yield c from the warehouse of the producer j to the distribution station k in scenario s |
| : | Transfer fee per unit of the yield c from the production station j to its own warehouse in scenario s |
| : | Transfer fee per unit of the returned yield c from the customer l to the hub station m in scenario s |
| : | Transfer fee per unit of the yield c in the production station j in the term t in scenario s |
| : | Transfer fee per unit of the yield c in the reproduction station j in the term t in scenario s |
| : | Transfer fee per unit of the yield c at the demolition station n in the term t in scenario s |
| : | Cost of collecting and sending a unit of the yield c to the hub station m in the term t in scenario s |
| : | Cost of recycling a unit of the yield c at the recycling station in the term t in scenario s |
| : | Cost of repairing a unit of the yield c in the repair station in the term t in scenario s |
| : | Cost of distributing a unit of the yield c in the distribution station k in the term t in scenario s |
| : | Fixed number of permanent job founded if established the hub station m in scenario s |
| : | Fixed number of permanent job founded if established the recycling station p in scenario s |
| : | Fixed number of permanent job founded if established the repair station in scenario s |
| : | Fixed number of permanent job founded if established the demolition station n in scenario s |
| : | Fixed number of permanent job founded if established the distribution station k in scenario s |
| : | Variable number of permanent job founded if established the hub station m in scenario s |
| : | Variable number of permanent job founded if established the recycling station in scenario s |
| : | Variable number of permanent job founded if established the repair station in scenario s |
| : | Variable number of permanent job founded the demolition station n in scenario s |
| : | Variable number of permanent job founded if established the distribution station k in scenario s |
| : | Missed working days due to work injury in the production station j in the term t |
| : | Equivalent to 1 if the vehicle of type goes from the supplier i to the manufacturer j from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle of type goes from the manufacturer j to the distributor k from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle of type 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 |
| : | Equivalent to 1 if the vehicle of type goes from the distributor k to the customer l from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle of type goes from the customer l to the hub m from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle of type goes from the hub m to the recycling station p from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle of type goes from the recycling station p to the manufacturer j from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle of type goes from the hub m to the manufacturer j from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle of type goes from the hub m to the demolition station n from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle type goes from the hub m to the repair station from the route r in the term t in scenario s, otherwise 0 |
| : | Equivalent to 1 if the vehicle type goes from the repair station to the distributor k from the route r in the term t in scenario s, otherwise 0 |
| : | 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 |
| : | 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 |
| : | Flow amount of the raw materials o from the supply station i to the production station j in the term t in scenario 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 |
| : | 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 |
| : | 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 |
| : | Flow amount of the returned yield c with the quality of the type q from the hub m to the repair station in the term t in scenario s |
| : | Flow amount of the yield c from the repair station to the distribution station k in the term t in scenario s |
| : | Flow amount of the yield c from the production station j to the distribution station k in the term t in scenario s |
| : | Flow amount of the yield c from the distribution station k to the customer l in the term t in scenario 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 |
| : | Flow amount of the yield c from the production station j to its own warehouse in the term t in scenario s |
| : | If the distribution station is established at the location k in the term t, its value would be equal to 1, otherwise 0 |
| : | If the recycling station is established at the location p in the term t, its value would be equal to 1, otherwise 0 |
| : | If the hub is established at the location m in the term t, its value would be equal to 1, otherwise 0 |
| : | If the landfill and demolition station is established at the location n in the term t, its value would be equal to 1, otherwise 0 |
| : | If the repair station is established at the location in the term t, its value would be equal to 1, otherwise 0. |
Data Availability
The data used to support the findings of this study will be made available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.