**Robust Multi-Objective Sustainable Reverse Supply Chain Planning: An Application in the Steel Industry**

**Jurgita Antucheviciene 1,\* , Ahmad Jafarnejad <sup>2</sup> , Hannan Amoozad Mahdiraji 2,3 , Seyed Hossein Razavi Hajiagha <sup>4</sup> and Amir Kargar <sup>5</sup>**


Received: 4 March 2020; Accepted: 3 April 2020; Published: 8 April 2020

**Abstract:** In the design of the supply chain, the use of the returned products and their recycling in the production and consumption network is called reverse logistics. The proposed model aims to optimize the flow of materials in the supply chain network (SCN), and determine the amount and location of facilities and the planning of transportation in conditions of demand uncertainty. Thus, maximizing the total profit of operation, minimizing adverse environmental effects, and maximizing customer and supplier service levels have been considered as the main objectives. Accordingly, finding symmetry (balance) among the profit of operation, the environmental effects and customer and supplier service levels is considered in this research. To deal with the uncertainty of the model, scenario-based robust planning is employed alongside a meta-heuristic algorithm (NSGA-II) to solve the model with actual data from a case study of the steel industry in Iran. The results obtained from the model, solving and validating, compared with actual data indicated that the model could optimize the objectives seamlessly and determine the amount and location of the necessary facilities for the steel industry more appropriately.

**Keywords:** multi-objective planning; reverse supply chain; robust optimization; uncertainty; meta-heuristic algorithm; steel making industry

### **1. Introduction**

With the expansion of the competitive environment, optimal supply chain (SC) design has become one of the fundamental issues business communities are facing [1]. This has affected all of the organization's activities to produce products, improve quality, reduce costs and provide the required services. On the other hand, with increasing greenhouse gas emissions and pollutants, managers of organizations and researchers are planning to set up networks that, in addition to considering economic optimization, have a special focus on environmental factors and the reduction of pollutants in all sectors [2].

The reverse logistics network, as part of the SC, means the accurate, correct and timely transmission of materials and the kinds of goods that are usable and unusable from the endpoint (last consumer or end-user) through the SC to the appropriate plant. In other words, reverse logistics is the process of moving and transferring goods and products that can be returned through the SC [3]. In this regard, the most important factor that is recognized in technical and economic studies of supply chains is

the demand parameter, which should be considered in the design of forward or reverse supply chain networks (SCNs) [4–7].

Moreover, many countries have an increasing interest in protecting the environment and applying environmental laws. Hence, industry owners and manufacturers have turned their attention to the design and development of the SC, taking into account environmental factors [8–12]. Green SC Design, integrating SC management with environmental requirements at all stages of product design, the selection and delivery of raw materials, production and manufacturing, distribution and transfer processes, delivery to the customer, and the management of recycling and reuse after consumption to maximize energy efficiency and the efficient use of resources are associated with improving the performance of the entire SC [8,13–20].

Several reasons justify the notion of reverse logistics and using recycled material in a reverse supply chain. The steel industry, with more than 2.5 trillion dollars worth of products, is important [21]. Usually, different economic, cost reduction, governmental regulatory, and social responsibility motivations encourage organizations to follow reverse logistic notions. Generally, the steel industry supply chain includes several stages of mining, processing, distributing and recycling. The concern of sustainability is very important in this industry. For instance, directly producing reduced iron instead of scrap requires 1120 cubic meters of water, 300,000 cubic meters of natural gas and 130,000 kilowatt-hours of electricity. This potential amount of saving has led to the introduction of a reverse logistic supply chain in the steel industry.

In this research, scenario planning was used to deal with uncertainty in the demand parameter due to unpredictable changes that have occurred during the research period in the studied case. To realize the economic, environmental and social effects of the reverse SCN and to optimize the model, three objectives were laid out, including maximizing operating profit, minimizing adverse environmental impacts and maximizing the level of service to suppliers and customers. To consider uncertainty, the multi-echelon supply chain, reverse logistics, and green supply design, the logistics network presented in this study consisted of four levels. The first level was the waste providers, considered as returning product suppliers, which could be the customers of previous periods who have returned their remaining products or can be new suppliers of scrap supplies. The second level was the gathering centers of the returned products, being responsible for supplying the scrap from the first-level suppliers of the chain, and in particular, being responsible for supplying the returned product, inspection, sorting, storage, and transferring the product to the recycling plants (product factories). The third level was the recycling plants for the production of new products, based on the received scrap from the gathering centers, which were responsible for producing new products. Finally, the fourth level was the customers. According to the given explanation, the proposed network is shown as a framework in Figure 1.

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 3 of 25

**Figure 1.** The considered supply chain (SC) scheme. **Figure 1.** The considered supply chain (SC) scheme.

The overall aim of this study was to design a sustainable reverse logistics integrated model in conditions of demand uncertainty, to optimize the flow of materials throughout the SC, and to determine the number and location of facilities and the planning of SC transportation. In this regard, the following objectives were considered in this research: The overall aim of this study was to design a sustainable reverse logistics integrated model in conditions of demand uncertainty, to optimize the flow of materials throughout the SC, and to determine the number and location of facilities and the planning of SC transportation. In this regard, the following objectives were considered in this research:


This study presents a multi-objective mathematical model for reverse SC design. The proposed model allows the NSGA-II algorithm to plan the recycling of products in the Iranian steel industry based on the modified approach of Feito Cespon et al. (2017) with their model. The proposed model has the following features [22]: 1. Using a robust optimization approach and NSGA-II algorithm for the multi-objective modeling This study presents a multi-objective mathematical model for reverse SC design. The proposed model allows the NSGA-II algorithm to plan the recycling of products in the Iranian steel industry based on the modified approach of Feito Cespon et al. (2017) with their model. The proposed model has the following features [22]:


of materials and the types of goods that are usable and unusable from the endpoint (last consumer or end-user) through the SC to the appropriate plant. In this regard, many types of research have been previously illustrated. Table 1 compares the illustrated pieces of research with the proposed method, from different perspectives. A reverse logistics network as part of the SC means the accurate, correct and timely transmission of materials and the types of goods that are usable and unusable from the endpoint (last consumer or end-user) through the SC to the appropriate plant. In this regard, many types of research have been previously illustrated. Table 1 compares the illustrated pieces of research with the proposed method, from different perspectives.



**Table 1.** Previous related research.

According to previous research, although many studies and articles have focused on the issue of sustainable SCN, there are some knowledge gaps in this area that are briefly summarized as follows:


Based on the above-mentioned gaps, the present study expresses a multi-objective mathematical model for redesigning the reverse SC network. The proposed model allows the use of a robustness approach to recycling multiple products. The proposed model has the following features:


### **2. Methods and Materials**

Based on the above mentioned theoretical background and defined problem, to resolve the uncertainty of the model parameters, a scenario-based robust optimization by the augmented epsilon constraint method was developed using the General Algebraic Modeling System (GAMS) software. An augmented epsilon constraint method was used for cases of multi-objective optimization in which one of the objective functions was more important than the other functions, and based on this, the optimization of other functions was performed. In this study, due to the greater importance of the first objective function compared to the others, this method was used. Since the main model with the network and the actual data of the case study by the GAMS cannot be solved, the objectives were defined based on an initial model on a smaller scale; they were solved by a scenario-based robust optimization, and the comparison and validation of the model were investigated by the NSGA-II algorithm in MATLAB.

### *2.1. Assumptions*


### *2.2. Model Notations*

Based on the SC structure shown in Figure 1, the problem is formulated as a multi-objective optimization model. The list of indices, parameters and decision variables of the model are presented in this section. This proposed model is based upon the work of Feito Cespon et al. 2017 [22] by using different objective functions, and case study and solving approaches, to compare the results.

The notations used in the paper are as follows.

### **Indexes**



### **Parameters**


### **Variables**


### *2.3. Model Objective Functions*

Equation (1) maximizes the operating profit of the SC.

$$\begin{array}{ll} \max f\_{1} &= \sum\_{\mathbf{s}} P \mathbf{S}\_{\mathbf{s}} \big( \sum\_{k} P R \mathbf{I}\_{p} \sum\_{l} \sum\_{m} Q \mathbf{P} \mathbf{C}\_{klmps} \\ &- \left( \sum\_{m} \operatorname{CUT}\_{m} \cdot \left( \sum\_{j} \sum\_{k} V \mathbf{P} \mathbf{P}\_{jkmps} d\_{jk}^{\mathbf{P} \mathbf{C}} \right. \\ &+ \sum\_{k} \sum\_{l} V \mathbf{P} \mathbf{C}\_{klmps} d\_{kl}^{\mathbf{P} \mathbf{C}} \\ &+ \sum\_{i} \sum\_{j} V \mathbf{S} \mathbf{R}\_{ijms} d\_{ij}^{\mathbf{S} \mathbf{R}})) \\ &- \sum\_{k} \operatorname{CFP}\_{k} \mathbf{P}\_{k} - \sum\_{k} \sum\_{p} \operatorname{CLP}\_{kp} \sum\_{l} \operatorname{QPC}\_{klmps} \\ &- \sum\_{j} \operatorname{CFN}\_{j} \mathbf{R}\_{j} - \sum\_{j} \sum\_{p} \operatorname{CILR}\_{jp} \sum\_{k} \operatorname{QRP}\_{jkms} \end{array} \tag{1}$$

Equation (2) minimizes the adverse environmental impacts of the SC.

$$\begin{array}{ll} \min f\_{2} &= \sum\_{s} P S\_{s} (\\ & \sum\_{m} \varPi\_{m} (\sum\_{j} \sum\_{p} \bigQ S R\_{ijmps} d\_{ij}^{SR} \\ &+ \sum\_{j} \sum\_{k} \sum\_{p} \bigQ P P\_{jkmps} d\_{jk}^{RP} \\ &+ \sum\_{k} \sum\_{l} \sum\_{p} Q P C\_{klmps} d\_{kl}^{PC}) \\ &+ \operatorname{IE} \left( \sum\_{k} \operatorname{Cfe}\_{k} P\_{k} + \sum\_{j} \operatorname{Cfe}\_{j} R\_{j} + \sum\_{p} \operatorname{Cve}\_{p} \sum\_{k} \sum\_{m} Q P C\_{klmps} \right) \\ &+ \operatorname{IP} \left( \sum\_{k} \alpha\_{k} P\_{k} + \sum\_{j} \beta\_{j} R\_{j} \right) + \operatorname{IA} \sum\_{p} \operatorname{Cve}\_{k} \sum\_{l} \sum\_{m} Q P C\_{klmps} \right) \end{array} \tag{2}$$

*Symmetry* **2020**, *12*, 594

Equation (3) maximizes the supplier's and customers' service levels, being different from that in the work of Feito Cespon et al. 2017 [22]:

$$\max f\_3 = \sum\_s P\_s(\frac{\left[\sum\_k \sum\_l \sum\_m QPC\_{klmps} + \sum\_i \sum\_j \sum\_m QSS\_{ijmps}\right]}{\left[\sum\_l \sum\_p D\_{lps} + \sum\_i \sum\_p G\_{ip}\right]}\tag{3}$$

### *2.4. Model Constraints*

The model constraints are shown in Equations (4) to (19). Each constraint has been discussed below [22]:

• Equations (4) to (6) guarantee the flow of materials through the SCN. The output from each center is, at most, equal to the inputs from different centers at the previous level of the SC;

$$\sum\_{j} \sum\_{m} QSR\_{ij\text{mpps}} \le G\_{ip} \,\forall i, p, s \tag{4}$$

$$\sum\_{k} \sum\_{m} QRP\_{jkmps} \le \sum\_{i} \sum\_{m} QSR\_{ijmps} \,\forall j, p\_\prime s \tag{5}$$

$$\sum\_{l} \sum\_{m} QPC\_{klmps} \le \sum\_{j} \sum\_{m} QRP\_{jkmps} \,\forall k, p\_\prime s \tag{6}$$

• Equations (7) and (8) respectively guarantee that the flow of materials rate does not exceed the maximum capacity of the recycling plants and the product demand;

$$\sum\_{l} \sum\_{m} QPC\_{klmps} \le \mathbb{C}\_{kp} \,\forall k, p, s \tag{7}$$

$$\sum\_{k} \sum\_{m} \text{QPC}\_{\text{klmps}} \le D\_{\text{lps}} \,\forall l, p, s \tag{8}$$

• Equations (9) to (11) maintain the balance between two facilities concerning the number of transportation. Since the number of transport must be an integer value, a series of inactive variables have been suggested to maintain the model's probability;

$$\sum\_{p} \frac{\text{QSR}\_{\text{ijmsps}}}{\text{CT}\_{\text{mp}}} + \text{HSR}\_{\text{ijms}} = \text{VSR}\_{\text{ijms}} \,\forall i, j, m, \text{s} \tag{9}$$

$$\sum\_{p} \frac{\mathcal{Q}RP\_{jkmps}}{\mathcal{C}T\_{mp}} + HRP\_{jkms} = VRP\_{jkms} \,\forall \, j, k, m, s \tag{10}$$

$$\sum\_{p} \frac{Q \text{PC}\_{\text{klmps}}}{\text{CT}\_{\text{mp}}} + H \text{PC}\_{\text{klms}} = V \text{PC}\_{\text{klms}} \,\forall \text{k}, l, m, \text{s} \tag{11}$$

• Equations (12) to (14) guarantee that ineffective variables focus only on differences in the number of transport;

$$VSR\_{ijms} + HSR\_{ijms} \ge 0 \quad \forall i, j, m, s,\tag{12}$$

*Symmetry* **2020**, *12*, 594

$$VRP\_{jkms} + HRP\_{jkms} \ge 0 \quad \begin{array}{l} \forall j, k, m, s, \\ \left(-1 < HRP\_{jkms} < 1\right) \end{array} \tag{13}$$

$$\text{VPC}\_{klms} + \text{HPC}\_{klms} \ge 0 \quad \begin{array}{l} \forall k, l, m, s, \\ (-1 < \text{HPC}\_{klm} < 1) \end{array} \tag{14}$$

• Equation (15) limits the number of trips per transportation mode [21];

$$\sum\_{i} \sum\_{j} VSR\_{ij\text{ms}} + \sum\_{j} \sum\_{k} VRP\_{jk\text{ms}} + \sum\_{k} \sum\_{l} VPC\_{kl\text{ms}} \le NV\_m \,\forall m,s \tag{15}$$

• According to Equations (16) and (17), binary variables should be assumed, such that if a gathering center or recycling plant is used in the model, then the value is 1, and otherwise it is zero;

$$\sum\_{k} \sum\_{m} \sum\_{p} \mathcal{Q} \mathcal{R} P\_{jkmps} \le M \, R\_j \, \forall j, s \tag{16}$$

$$\sum\_{l} \sum\_{m} \sum\_{p} QPC\_{klmps} \le MP\_k \,\,\forall k,s\tag{17}$$

• Finally, Equations (18) and (19) show the nature of the variables.

*QRPjkmps*, *QPCklmps*, *QSRijmps*, *VPCklms*, *VRPjkms*, *VSRijms* ≥ 0 (18)

$$R\_{j\prime}P\_k \in \{0, 1\} \tag{19}$$

### **3. Results**

This proposed model is based upon the work of Feito Cespon et al. 2017 [22] by using a different objective function and case study and different solving approaches to compare the results. This is due to the main model (real-world case study) not being solvable with mathematical programming methods. At first, a smaller scale research model with fewer data called the "Initial Model" can be solved by an augmented epsilon constraint method, and subsequently, it is solved by the NSGA-II algorithm [46,47]. Furthermore, the performance of the NSGA-II algorithm is evaluated by solving several examples in the proposed research model, and the corresponding criteria are calculated. Afterwards, with the confidence of the model's validity, and since the main problem is Np-hard, the model of the real-world case study called "Main Model" is solved by the NSGA-II algorithm and will be analyzed at the end. Note that sensitivity analysis will be implemented on the initial model. Before solving the model, the method for determining the chromosome in the NSGA-II algorithm is presented. The algorithmic scheme of this section is illustrated in Figure 2.

**Figure 2.** The algorithmic scheme of empirical research. **Figure 2.** The algorithmic scheme of empirical research.

### *3.1. Definition of Chromosomes in the NSGA-II Algorithm. 3.1. Definition of Chromosomes in the NSGA-II Algorithm*

The matrix of the answer in the model has two sections, called allocation and assignment. The allocation section has two parts, the first part of which is the location of the gathering centers (J), and the second part is the location of the recycling plants (K). The cells of the matrix are filled with numbers 0 or 1, for example, if J = 5 and K = 3; an example of this matrix is given in Table 2. The matrix of the answer in the model has two sections, called allocation and assignment. The allocation section has two parts, the first part of which is the location of the gathering centers (J), and the second part is the location of the recycling plants (K). The cells of the matrix are filled with numbers 0 or 1, for example, if J = 5 and K = 3; an example of this matrix is given in Table 2.

**Table 2.** The allocation matrix. **J1 J2 J3 J4 J5** According to Table 2, gathering center No. 5 has been constructed, and gathering centers No. 1, 2, 3 and 4 have not been constructed. Moreover, the recycling plant No.1 has been constructed and the recycling centers No. 2 and 3 have not been constructed. Each set is the answer that is called a chromosome. The assignment section shows the flow rate from the waste supplier to the gathering centers, from the gathering centers to the recycling plants and from the recycling plants to the customers. In Figure 3, the assignment section is shown. As a case in point, in the first part of the following table, which is a K ∗ L dimensional matrix, the flow rates from the recycling plants (k) to the customer (l) for the transportation mode (m), that is m = 1, p = 1 and s = 1, are shown.

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 15 of 25

Tables 3 and 4 show that when the values of the first objective function deteriorate, the values of the other objectives function do not. In other words, these values remain either constant or close to their optimal values, which is the expected process that the multi-objective models suggest. Based on Table 3, it is obvious that the flow rate of product 1 from the recycling plant 1 to customer 1 by transportation mode 1 in scenario 1 is 3. In this regard, each answer set is called a chromosome, and each cell is a gene.



### *3.2. NSGA-II Operator Selection*

Achieving a high performance of genetic algorithms is highly dependent on the performance of the genetic operators. One of the main operators in genetic algorithms is crossover. The crossover operator is used to generate a new chromosome by crossing over two selected chromosomes. Different

**Figure 4.** Pareto points obtained by the NSGA-II algorithm in the initial model. **Figure 3.** Pareto points obtained by the NSGA-II algorithm in the initial model.

crossover operators are represented in previous studies. Here in this paper, the single point crossover is used. The next important operator is a mutation to assure diversity. Beyond the mutation probability that is tuned in Section 3.6, in this study, the reverse and replace operators are used randomly to mutate the selected chromosomes. In reverse mutation, two genes are selected in the considered chromosome, and the values of remaining genes between these two selected genes are reversed from right to left. In the replacement mutation, two genes are selected, and their positions are swapped with one another's.

### *3.3. Initial Model Solving Results*

Table 4 indicates SC characteristics in the initial model. Furthermore, the probability of the occurrence of each scenario is obtained using the analytical hierarchical process (hereafter AHP) method, which for scenarios 1 to 2, is 52.4% and 47.6%, respectively. Moreover, in the initial model, a big M value is 10,000.


**Table 4.** The specifications and parameters of the initial model.

The Pareto result according to the GAMS and NSGA-II algorithm of five sets of answers derived from the initial model solving is shown in Tables 5 and 6. Figures 3 and 4 show these Pareto points. It is conceivable that the results obtained for the small scale problems by GAMS will outperform the NSGA-II results, as this can be seen in similar studies [48–50]; however, as it is illustrated in the next sections, the main advantage of NSGA-II is in its ability to solve large scale and real-world problems. The time of the GAMS solving in this model, although the problem dimensions are low, is 326 seconds, which is increased sharply by increasing the dimensions of the problem.


**Table 5.** Pareto points set by the GAMS for the initial model.


*Symmetry* **2020**, *12*, x FOR PEER REVIEW 15 of 25

**Figure 3.** Pareto points obtained by solving the GAMS in the initial model. **Figure 4.** Pareto points obtained by solving the GAMS in the initial model.

Tables 3 and 4 show that when the values of the first objective function deteriorate, the values of the other objectives function do not. In other words, these values remain either constant or close to their optimal values, which is the expected process that the multi-objective models suggest.

### *3.4. Model Validation*

Constraint method.

To evaluate the performance of the model and to compare the performance of the NSGA-II algorithm with the augmented epsilon constraint method, five examples with different dimensions randomly compiled on the research model, and the criteria for comparing the efficiency of the multi-objective algorithms, are calculated, the results of which are shown in Table 5. The results are obtained by running the proposed algorithm in a single trial with a population size of 10,000 and 250 repetitions. As the results indicate, it can be seen that using the NSGA-II algorithm has the necessary validity to solve the main model.

**Figure 4.** Pareto points obtained by the NSGA-II algorithm in the initial model. Tables 3 and 4 show that when the values of the first objective function deteriorate, the values of the other objectives function do not. In other words, these values remain either constant or close to their optimal values, which is the expected process that the multi-objective models suggest. In this table, five measures are reported. Mean Idear Distance (MID) measures the convergence of an algorithm by averaging the distances of solutions from the best feasible solution [51,52]. Spacing measures the standard deviation of the distances among the Pareto front solution [52]. Diversity evaluates the spread of the Pareto front [52]. The Number of solutions (NoS) is the number of different Pareto solutions [53,54]. Time(s) is the time for which the algorithm needs to be run to reach the near-optimal solution [53,54].

*3.4. Model Validation* To evaluate the performance of the model and to compare the performance of the NSGA-II algorithm with the augmented epsilon constraint method, five examples with different dimensions randomly compiled on the research model, and the criteria for comparing the efficiency of the multi-The lower the index MID, the better the research results. As can be seen, the performance of the Epsilon Constraint (E.C). method in two sets of responses is better than that of the NSGA-II algorithm; however, with increasing dimensions of the problem, the method of E.C. loses its effectiveness. Since this difference is not very high, both indicators have shown good performance.

objective algorithms, are calculated, the results of which are shown in Table 5. The results are obtained by running the proposed algorithm in a single trial with a population size of 10,000 and 250 repetitions. As the results indicate, it can be seen that using the NSGA-II algorithm has the necessary validity to solve the main model. The lower the index spacing, the better the research results. According to Table 7, the performance of the NSGA-II algorithm is better than that of the E.C. method. The higher the index (diversity), the better the research results. Based on the results, it can be seen that the proposed E.C method is better than the NSGA-II algorithm.

**Table 7.** Comparison of indices for five examples with NSGA-II algorithms and the Epsilon

 8556.59 51.73 2562.22 14 8 7283.27 13.05 2858.35 96 55.96 7574.8 73.37 3820.57 16 13 7749.36 79.91 2247.99 97 63.73 5383.46 181.73 6709.34 23 48 6720.7 30.82 3616.23 98 58.43 6317.78 181.41 5209.74 16 93 7150.89 28.77 2423.52 99 57.68 7109.71 242.98 5219.62 18 407 7903.34 13.88 1650.52 95 58.53

**MID Spacing Diversity NoS Time(s) MID Spacing Diversity NoS Time(s)**

**Item Epsilon Constraint NSGA-II**


**Table 7.** Comparison of indices for five examples with NSGA-II algorithms and the Epsilon Constraint method.

The higher the index NoS, the better the research results. Based on the results, the NSGA-II algorithm obtained a greater number of Pareto members. It is logical to increase the solution time of the algorithms by increasing the dimensions of the problem.

Therefore, according to the results, the same trend is observed; with an increase in the dimensions of the problems, the time taken to solve by the method of E.C. increases exponentially, and this method loses its efficiency in high-dimensional issues. However, it is almost constant for the NSGA-II algorithm.

### *3.5. The Parameter Adjustment of the NSGA-II Algorithm*

Under the meta-heuristic algorithms that do not guarantee an exact optimal solution, the algorithm may be followed by a different response at any time by solving it. Therefore, a meta-heuristic algorithm is good when used, with almost identical answers each time. The most influential parameters in the NSGA-II algorithm are the number of initial population (nPop), the number of repetitions (MaxIt), the intersection rate (Pc) and the rate of mutation (Pm). With using the Taguchi design of experiments method, the parameters of this algorithm are based on comparative criteria for nine exams that have been determined under the following steps.

### *3.6. Taguchi Design of Experiment*

In the NSGA-II algorithm, the four factors/parameters MaxIt, nPop, Pc, and Pm should be set to optimal levels. For this purpose, at first, for each parameter, three levels of low (1), medium (2) and high (3) are considered, as shown in Table 8. The proposed Taguchi experiments for four factors at three levels are shown in Table 9 for nine experiments. These experiments are designed based on Taguchi methods [55].


**Table 8.** The setting up NSGA-II parameters at three levels.


**Table 9.** Taguchi designed experiments to adjust the parameters of the NSGA-II algorithm.

Table 8 reveals the levels of the NSGA-II parameters, for each parameter considering three different levels. Table 9 demonstrates the experiments designed to adjust the parameters of the NSGA-II algorithm. In Table 10, the results of the NSGA-II algorithm for nine independent experiments are presented.


**Table 10.** Results from the experiments of the NSGA-II algorithm.

According to Table 10, to create an output from each test and for five criteria, using the fuzzy unambiguous technique and the ideal planning approach, all indicators become responses after normalization. The normalization of the results and the calculation of the response variable are shown in Table 11.

**Table 11.** Normalized results and the calculation of responses for setting the parameter of the NSGA-II algorithm.


In the last step, based on the calculated response variable in the previous step, the S/N rate is calculated, and the optimal levels of the input parameters are determined. This operation is performed by the MINITAB software, and the results are illustrated in Figure 5. This figure illustrates the main effect plot of different algorithm parameters. The plots are plotted by fixing parameters at their three levels, and then comparing the means of the S/N ratios against those at different levels.

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 18 of 25

**Figure 5.** MINITAB output to adjust the parameter of the NSGA-II algorithm. **Figure 5.** MINITAB output to adjust the parameter of the NSGA-II algorithm.

The optimal levels for the parameters of the algorithm examined according to Figure 5 and the above tables are shown in Table 12. The optimal levels for the parameters of the algorithm examined according to Figure 5 and the above tables are shown in Table 12.


**Table 12.** Optimized levels for the NSGA-II algorithm. **Table 12.** Optimized levels for the NSGA-II algorithm.

### *3.7. Results of Solving the Main Model (Case Study) 3.7. Results of Solving the Main Model (Case Study)*

Number of Suppliers 5

Table 13 shows the main model of the case study. The case study is related to a steel production company with 8.1 million tons of yearly capacity in Iran. The company produced a set of intermediate and final products. The considered network includes suppliers, recycling plants, transportation modes and gathering centers. Additionally, five types of product are identified to be delivered to three types of customer. The problem parameters are gathered from the company's databases, which are not presented due to their large magnitude. The problem specification and its parameters are illustrated in Table 13. Table 13 shows the main model of the case study. The case study is related to a steel production company with 8.1 million tons of yearly capacity in Iran. The company produced a set of intermediate and final products. The considered network includes suppliers, recycling plants, transportation modes and gathering centers. Additionally, five types of product are identified to be delivered to three types of customer. The problem parameters are gathered from the company's databases, which are not presented due to their large magnitude. The problem specification and its parameters are illustrated in Table 13.

conditions are defined for these conditions, which are considered as a scenario for each mode, **Table 13.** The specifications and parameters of the model in a case study.

According to the requirements of the considered organization as the case study, five incidental

Number of Gathering


Products <sup>5</sup> Centers Number of Recycling Plants 3 Number of Customers 3 225 Number of Scenarios 5 Number of Transportation Modes As stated above, since the size of the case study is high and the GAMS software is not able to According to the requirements of the considered organization as the case study, five incidental conditions are defined for these conditions, which are considered as a scenario for each mode, providing different data. Moreover, the probabilities of occurrence of each scenario are obtained using the AHP method, which are—for scenarios 1 to 5—16%, 22%, 48%, 9% and 5%, respectively. Furthermore, in the main model, a big M value of 10,000 is considered.

5

Number of

NSGA-II algorithm, the Pareto points from which are presented in Table 14 and Figure 6.

solve it at an acceptable time, the original model in the MATLAB software is solved based on the

**Answer** 

**The Value of the First** 

As stated above, since the size of the case study is high and the GAMS software is not able to solve it at an acceptable time, the original model in the MATLAB software is solved based on the NSGA-II algorithm, the Pareto points from which are presented in Table 14 and Figure 6. **The Value of the Third Objective Function Objective Function Objective Function No.** 1 1,168,678,032,301 17,570,971,633 0.4738 2 1,535,360,428,421 18,258,944,734 0.5864

**The Value of the Second** 

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 19 of 25


**Table 14.** A set of the Pareto points generated by algorithm NSGA-II for the case study. 3 1,532,610,197,259 18,179,813,506 0.5912 4 1,226,611,751,948 17,585,580,845 0.3662

**Figure 6.** The reverse supply chain (SC) after a solution. **Figure 6.** The reverse supply chain (SC) after a solution.

Besides, the values of the parameters of the NSGA-II algorithm are described in Table 15 in the main model. **Table 15.** The values of the algorithm NSGA-II operators for the case study. **MaxIt nPop Pc Pm** As expected, due to the high cost of constructing facilities in the steel industry, one center (center number 5) of the five candidates was selected, and one factory (factory number 1) of the three nominated recycling/production plants for construction is calculated, and it provides an excellent answer to the case study. In this way, the schematic of the reverse SC model is changed after the solution, as shown in Figure 6.

NSGA-II 100 100 0.6 0.3 Besides, the values of the parameters of the NSGA-II algorithm are described in Table 15 in the main model.


**Table 15.** The values of the algorithm NSGA-II operators for the case study.

For any Pareto optimal solution, the optimal values of decision variables are obtained(see Figure 7). For instance, the magnitudes of some decision variables for the first Pareto solution in Table 14 are represented in Table 16. According to this table, under the first scenario, a magnitude of 2.4605 of the first product type should be transported from the first supplier to the fifth gathering center. Other values can be interpreted similarly. For each Pareto optimal solution, a similar set of optimal decision variables is obtained.

*3.8. Sensitivity Analysis*

are shown in Table 17 and Figure 9.

**Amount of 1st O.F.**

**Amount of 2nd O.F.**

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 20 of 25

**Figure 7.** A set of the Pareto points for the case study. **Figure 7.** A set of the Pareto points for the case study.

For any Pareto optimal solution, the optimal values of decision variables are obtained(see figure **Table 16.** A sample of optimal decision variables for the first Pareto solution of Table 14.


**Variable Value Variable Value** QSR1,5,1,1,1 2.4605 QSR1,5,2,1,1 7.9593 QSR2,5,1,1,1 2.9064 QSR2,5,2,1,1 7.7476 QSR3,5,1,1,1 6.1607 QSR3,5,2,1,1 11.3636 QSR4,5,1,1,1 8.0460 QSR4,5,2,1,1 7.8387 To assure the convergence of the proposed algorithm, with its tuned parameters, the above problem is repetitively solved 100 times. Figure 8 illustrates the obtained results of the solutions for different objectives. According to this figure, it is clear that the algorithm-proposed solution in all of the objectives represents a controllable variance, and it can be considered as the proposed algorithm reliability. No outlier solutions can be detected in repetitions. *Symmetry* **2020**, *12*, x FOR PEER REVIEW 21 of 25

**Figure 8.** Algorithm consistency over repetitions. **Figure 8.** Algorithm consistency over repetitions.

parameters. In the present study, this type of analysis—as compared to the change in the demand parameter in the initial model, which decreased by 10%—and the results of this sensitivity analysis

**Table 17.** Changes in Pareto points for the initial model induced by changing the value of demand.

**Amount of 3rd. O.F.**

 −128.17 2271.64 0.8 −170.57 2117.32 0.8 −247.05 1135.82 0.6 −92.47 801.17 0.41 −92.47 801.17 0.41 42.57 1058.66 0.2 67.78 1135.82 0.2 49.15 2117.32 0.6 87.77 2271.64 0.6 158.82 1058.66 0.08

**Figure 9.** A set of the Pareto points for the initial model from changing the value of demand.

**Before Changing the Demand Parameter After Changing the Demand Parameter No. Amount of**

**Amount of 1st O.F.**

**3rd. O.F.**

**Amount of 2nd O.F.**

To investigate how the values of the objective functions varied, sensitivity analysis should be

0 2E+11 4E+11 6E+11 8E+11 1E+12 1.2E+12 1.4E+12 1.6E+12 1.8E+12

### *3.8. Sensitivity Analysis 3.8. Sensitivity Analysis*

1

6

11

16

21

26

31

36

41

46

51

56

61

To investigate how the values of the objective functions varied, sensitivity analysis should be performed on some of the parameters. Regarding the multiplicity of the model, two types of analysis are performed. The first type is the change in Pareto's values relative to the change in one of the parameters. In the present study, this type of analysis—as compared to the change in the demand parameter in the initial model, which decreased by 10%—and the results of this sensitivity analysis are shown in Table 17 and Figure 9. To investigate how the values of the objective functions varied, sensitivity analysis should be performed on some of the parameters. Regarding the multiplicity of the model, two types of analysis are performed. The first type is the change in Pareto's values relative to the change in one of the parameters. In the present study, this type of analysis—as compared to the change in the demand parameter in the initial model, which decreased by 10%—and the results of this sensitivity analysis are shown in Table 17 and Figure 9.

**Figure 8.** Algorithm consistency over repetitions.

66

71

76

81

86

91

96

0

0.1

0.2

0.3

0.4

2nd objective 1st objective 3rd objective

0.5

0.6

0.7

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 21 of 25


**Table 17.** Changes in Pareto points for the initial model induced by changing the value of demand. **Table 17.** Changes in Pareto points for the initial model induced by changing the value of demand.

**Figure 9.** A set of the Pareto points for the initial model from changing the value of demand. **Figure 9.** A set of the Pareto points for the initial model from changing the value of demand.

As demonstrated, with the decrease in the average demand, the value of the objective function decreased. By virtue of shipping costs, other items are reduced by decreasing demand. In Figure 4, the stars represent the Pareto before the change, and the circles represent the Pareto after the change.

The second type of sensitivity analysis is employed for one Pareto point (here, for the tenth Pareto point), which is done for the demand parameter. The variability of the value of the objective function concerning demand in the initial model is shown in Table 18 and Figure 10.

**Table 18.** Changes to the first objective function for the initial model induced by changing the amount of demand. 2 17.76 423.237 3 30.765 369.45

**Table 18.** Changes to the first objective function for the initial model induced by changing the amount

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 22 of 25

function concerning demand in the initial model is shown in Table 18 and Figure 10.

1 9.57 459.08

As demonstrated, with the decrease in the average demand, the value of the objective function decreased. By virtue of shipping costs, other items are reduced by decreasing demand. In Figure 4, the stars represent the Pareto before the change, and the circles represent the Pareto after the change. The second type of sensitivity analysis is employed for one Pareto point (here, for the tenth Pareto point), which is done for the demand parameter. The variability of the value of the objective

**Figure 10.** Sensitivity analysis of the first objective function in terms of the demand average. **Figure 10.** Sensitivity analysis of the first objective function in terms of the demand average.

**4. Discussion and Conclusions** In the present study, the model was defined as multi-objective functions based on the conditions of the uncertainty of demand and five scenarios; the model was solved by an augmented epsilon constraint method and the NSGA-II algorithm, and finally analyzed. This proposed model is based In Table 18, by assuming that the second and third objective functions are fixed, and that only the changes in the first objective function have been investigated, it is seen that with increasing average demand, the first objective function is reduced. In fact, with increasing demand, the amount of cost is higher than the amount of income.

### upon the work of Feito Cespon et al. 2017 [22] by using a different objective function and case study **4. Discussion and Conclusions**

and different solving approaches to compare the results. Because the number of levels and actual data of the model would be Np-hard by solving the GAMS, and it would not be able to achieve the optimal response, model validation and sensitivity analysis were done on a smaller scale. The results of the comparative indices showed that solving the model with the NSGA-II algorithm yielded acceptable results, and the main model was solved accordingly. Additionally, the optimal levels of the NSGA-II algorithm parameters were adjusted in the original model, based on the Taguchi design of experiments method. In analyzing the results, as expected, in the locating facility, one gathering center was selected from five candidates, and one recycling plant was selected from three candidate plants. The number of objective functions in different Pareto points has been obtained with a suitable In the present study, the model was defined as multi-objective functions based on the conditions of the uncertainty of demand and five scenarios; the model was solved by an augmented epsilon constraint method and the NSGA-II algorithm, and finally analyzed. This proposed model is based upon the work of Feito Cespon et al. 2017 [22] by using a different objective function and case study and different solving approaches to compare the results. Because the number of levels and actual data of the model would be Np-hard by solving the GAMS, and it would not be able to achieve the optimal response, model validation and sensitivity analysis were done on a smaller scale. The results of the comparative indices showed that solving the model with the NSGA-II algorithm yielded acceptable results, and the main model was solved accordingly. Additionally, the optimal levels of the NSGA-II algorithm parameters were adjusted in the original model, based on the Taguchi design of experiments method. In analyzing the results, as expected, in the locating facility, one gathering center was selected from five candidates, and one recycling plant was selected from three candidate plants. The number of objective functions in different Pareto points has been obtained with a suitable and acceptable dispersion criterion. The model showed that it could be integrated into optimizing the objectives, determining the number and location of necessary facilities and planning the transportation between different levels of the steel industry. Some of the assumptions implied in the current study can be adjusted in future studies. First, there are some general assumptions including the uncertainty of demand parameter, the determinedness of costs, the number and capacity

of transportation modes and the number of supply chain levels. These assumptions can be generalized straightforwardly, by considering the uncertainty of other parameters and extending the model into more levels. Additionally, in the current paper, it is assumed that the number of suppliers, customers, gathering centers and recycling plants are determined. However, in some cases, the problem can be formulated to select different markets, suppliers, and gathering and recycling centers in a broader scope. These extensions will not change the structure of the proposed method drastically. However, altering some assumptions requires fundamental changes in the proposed model. Among these assumptions, reference to transshipment among levels (ignoring the sixth assumption) and limiting the capacity of gathering centers and recycling plants can be made. Further research can be done on the above-mentioned problems by changing the current model's assumptions.

**Author Contributions:** J.A.—conceptualization, formal analysis, writing—review and editing; A.J.—conceptualization, methodology, writing—original draft preparation, supervision; H.A.M.—conceptualization, formal analysis, resources, writing—review and editing, visualization; S.H.R.H.—conceptualization, methodology, software, validation, resources, writing—original draft preparation; A.K.—conceptualization, methodology, software, validation, investigation, data curation, writing—original draft preparation, project administration. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
