1.INTRODUCTION
Nowadays, designing an efficient supply chain is seemed essential to develop a successful business model. Researchers have defined the Supply Chain (SC) configuration as follows “an integrated system of facilities and activities consisting business functions of material procurement, material transformation to final products and also the distribution of these products to the end users” (Nourmohamadi Shalke et al., 2017). Besides, Supply Chain Management (SCM) forms some special approaches to product and to distribute the activity planning in different decision levels in order to achieve the chain goals efficiently. Such methods should integrate all the production resources of the different plants started from the suppliers by procuring of raw materials and ended with the shipping of finished products to the customer zones (Hajiaghaei Keshteli and Sajadifar, 2010; HajiaghaeiKeshteli et al., 2011). In the same way, Supply Chain Network Design (SCND) problem is defined as one of the most significant strategic decisions in subject of SCM which newly has been received the growing attention by both scholars and academia. Generally, choosing good facilities among the possible sites, determining numbers and capacities of the selected facilities as well as the material flow through the network are the main concerns and decisions in SCND (Devika et al., 2014). Besides, SCND as a key strategic issue, directly influences on both tactical and operational activities through the entire network like inventory, transportation, and production decisions or management (Aminnayeri and Hajiaghaei Keshteli, 2012; Nasiri et al., 2017).
The literature of SCND is classified into the three areas according to the type of networks. The forward, reverse and closedloop networks are the most cited classifications in the related studies. The forward one is considered the flow of products from the suppliers to the customers. To clarify the mentioned network, after purchasing from suppliers, raw materials are transformed into finished products in manufacturing centers, and then these products are delivered to customers via distribution centers to satisfy their demands (Hajiaghaei Keshteli et al., 2014, 2010). The reverse configuration is focused on the backward flow of products from customers to recovering centers. Consequently, the recovering centers divided them into three types of facilities. They are transformed to the remanufacture, recycling and disposal centers. In the Reverse Logistics (RL), the flow of returned products is processed from the customers back to the recovering centers for repair. And the last classification network is a closedloop network which considers both mentioned types on an integrated manner (Spar and La Mure, 2003). The reader can review some other comprehensive literature reviews, like Devika et al. (2014) and Govindan et al. (2015, 2016) to study the literature review.
The previous studies showed that many efforts have been performed to expand the SCND models which most of them are based on the deterministic approaches with a single objective (i.e. Jayaraman and Pirkul, 2001), while the real problems are structured by diverse source of uncertainty and numerous measures (i.e. Pishvaee et al., 2011). As pointed out by Sabri and Beamon (2000), one of the most challenging but significant problems in supply chain management is stochastic models. However, the literature shows the uncertainty which integrates all the parts in the network and allocates them is still scared (Afrouzy et al., 2016). In order to make the management considerably simpler, several firms insist on strongly the single sourcing associated with customers zones (Daskin et al., 2003). Nevertheless, researchers do not consider this feature as it makes the problem significantly more difficult.
This paper addresses the bilevel stochastic model by considering the economic and risk measure to be minimized. Due to obvious the different parts of proposed closedloop network and flow of the products, a graphical structure is provided as an example in the Figure 1.
In this work, a closedloop supply chain network is considered and solved by new approaches in metaheuristics and exact methods. The developed twostage stochastic model aims to minimize the expected total cost and the downside risk, simultaneously. Besides, this study presents the closedloop network which considers the forward and reverse networks in an integrated manner. In the forward network, this case just considers the forward network, while, the reverse logistic fully focused on the backward network by considering recovering centers (i.e. from recovering centers to remanufacture, recycling and disposal centers). Three powerful metaheuristics which used in the related studies, repeatedly, Imperialist Competitive Algorithm (ICA), Particle Swarm Optimization (PSO) and Genetic Algorithm (GA) are considered to solve the problem. Due to check the results of metaheuristics in small sizes, an exact solution approach, namely, εconstraint method is utilized in the small sizes. In addition, the parameters of algorithms are tuned by Response Surface Method (RSM) with an MODM approach. To explain the efficiency and effectiveness of methods, four assessment metrics are introduced. The other contributions of this study are investigated clearly in the following sections.
In brief, the rest of the paper is organized as follows. Section 2 explains the literature review which is described in detail. In Section 3, the proposed model in a twostage stochastic programming is presented. In Section 4, the optimization techniques and proposed algorithms are explained, respectively. Experimental evaluations are conducted to evaluate the performance of the introduced algorithms against each other in terms of different criteria and analysis in Section 5. Finally, Section 6 presents the main findings and future lines of the research.
2.LITERATURE REVIEW
As mentioned earlier, the recent studies display that the importance of reverse network has been increased. Due to the significance of collecting and treating endoflife (EOL) products, the recent researches are reviewed the various aspects of RL. For instance, Sasikumar and Kannan (2008a, 2008b, 2009), Pokharel and Mutha (2009), and also Yi et al. (2016) overviewed the different factors of RL. Accordingly, designing RL networks both in industries and academia has been one of the interesting terms in the literature in recent years (Kannan et al., 2012; Srivastava, 2008).
This paper presents a new classification in SCND models by considering the stochastic programming. As illustrated in Table 1, the expanding of stochastic models in this area is still scared. Except Yang et al. (2015), none of the previous studies have not used the metaheuristics procedures to address the problem. Hence, this work developed metaheuristics to increase the metaheuristics papers in the stochastic models in this zone.
As discussed before, the majority of network design problems are NPhard. Because of the metaheuristics are so flexible, numerous solution methods including heuristics and metaheuristics have been developed yet (Rajabi et al., 2013). For instance, Dehghanian and Mansour (2009) used the genetic algorithm to handle a sustainable supply chain network which have considered the economic, environmental and social concerns in a reverse configuration. Consequently, Devika et al. (2014) developed six hybridized metaheuristic approaches and a modified imperialist competitive algorithm to solve a closedloop supply chain with three conflict goals such as Dehghanian and Mansour (2009). Besides, Jabbarzadeh et al. (2016) and Ardalan et al. (2016) introduced models for SCND with multimode demand satisfaction having maximization objectives for the total profit of the distribution network.
In order to be simpler, most of the papers in SCND have utilized the Mixed Integer Programming (MIP) to model their problems. These models change from simple single objective forward facility location models (i.e. Pishvaee et al., 2011) up to complex multiobjective closedloop models (e.g. Ardalan et al., 2016). Elhedhli and Merrick (2012) considered emission costs as well as location and production costs in a forward SCND problem. Besides, Talaei et al. (2016) introduced a Fuzzy Mixed Integer Programming (FMIP) model for a green closedloop supply chain network. They minimized the total cost and environmental concerns in the deterministic model.
Since the closedloop SCND (CLSC) is the more realistic network, it becomes an interesting term in the literature to integrate the forward and reverse network design. The CLSC models show that this term eschews the suboptimality arising in the separated modeling of forward and reverse configurations instead of integrated ones (Soleimani et al., 2014). One of the first CLSC networks is proposed by Fleischmann et al. (2001). As the result, they stated that some of the RL operations such as remanufacturing and recycling centers can be integrated efficiently into the existing forward cycle network such as the distributing and manufacturing centers. Later, Schultmann et al. (2006) considered EOL vehicle treatment in Germany and focused on the flow of used products and reintegrated reverse flow of used products. They modeled RL with vehicle routing planning and solved their model by Tabu Search (TS).
Although the most of issues in SCND developed the deterministic models, only a few papers utilized the stochastic programming to expand the models in this area, as seen in Table 1 (Ramezani et al., 2013). In the stochastic programming, for the demands of costumers and price of products as the examples, in order to seem realistic, some different scenarios are defined. At the end, the decision maker by considering the cost and profit of each scenario decides a performance one, and the probability of the scenarios can be changed.
In order to highlight the stochastic programming, just Ramezani et al. (2013) presented the network quality and risk measure as one of the objective functions to be optimized. In addition, one of the key issues in this study is developed the metaheuristic approaches to solve the proposed Stochastic Mixed Integer Programming (SMIP) model in a closedloop SCND. As mentioned in Table 1, only a few papers use the twostage stochastic programming to address new contributions of this area like the multiobjective, multiechelon, multitype of production technology in the logistics network design problem.
In order to list the difference of this paper from the prior ones, some points should be investigated. Firstly, this paper optimizes some economic aspects and the risk measure in an easy SMIP model. Secondly, the authors consider all sections of production flow on the proposed CLSC as depicted in Figure 1. Although Yang et al. (2015) and Pishvaee et al. (2011) modeled forward and reverse networks, respectively, this essay explored the closedloop network with two different objectives. And the last but not least, a twostage stochastic programming is introduced in this work as well.
Due to the multiobjective optimization, the solutions must be satisfied and be robust in all different goals. Another contribution of this study is using the multiobjective programming and selecting the best nondominated solutions and evaluating the quality of the best front of each method by recent assessment metrics. While most of the papers in SCND, in the recent and old contents, utilize the single objective function to be minimized or maximized (e.g. Jayaraman and Pirkul, 2001; Ardalan et al., 2016).
At the end, as noticed before, the SCND models are an NPhard problem (Eckert and Gottlieb, 2002; Jo et al., 2007; Syarif et al., 2002). Subsequently, the powerful and wellknown approaches in metaheuristics are considered to find the Paretooptimal solutions. The Imperialist Competitive Algorithm (ICA) proposed by AtashpazGargari and Lucas (2007), Particle Swarm Optimization (PSO) introduced by Eberhart and Kennedy (1995) and Genetic Algorithm (GA) (Holland, 1975; Deb et al., 2000) are utilized to solve the proposed stochastic closedloop SCND problem. In addition to check the algorithm’s results in small sizes, the εconstraint method is utilized to tackle such NPhard problem, too.
3.PROBLEM DESCRIPTION
3.1.Assumptions of the Proposed Problem
This work investigates the multiobjective, multiechelon, multitype of production technology in the stochastic closedloop supply chain network which required the more efforts to analyze both the forward and the backward flow of production, simultaneously. The proposed model is mainly focused on the reverse network consisting of the treatment facilities: Remanufacturing, Recycling and Disposal centers, as shown in Figure 1. The assumptions of the developed model are on the following of the related works in this area (Hsu and Wang, 2009; Syarif et al., 2002; Wang and Hsu, 2010; Yao and Hsu, 2009; Devika et al., 2014):

1. The demand of each customer must be met and depended on the considered scenario.

2. The number of facilities and their potential sites in each echelon is predefined.

3. There are no any flow between the facilities of the same echelon.

4. As a special supposition in closedloop logistics suggested by Van Der Laan et al. (1999), it is assumed that the number of the EOL products returned to the recovering centers is a fraction of the customers’ demands (Soleimani and Kannan, 2015). In addition, they are allocated to the treatment facilities based on their qualities.

5. The rate of purchasing the production, manufacturing costs, allocating cost, demands and return rates are uncertain and are described by the set of scenarios. It should be noted that the return rate to recover centers in all customers depends on the customer’s demand.

6. The other costs (i.e., fixed costs and transportation costs) are known.
3.2.Formulation of the Proposed TwoStage Stochastic Programming
In this subsection, the notations of the proposed model are presented. The main aim of the expected total cost is to decide about the amount of the products which have to be manufactured at each plant, the allocation to the distribution centers and customers, focusing on the backward flow to the recovering centers; and also the flow of materials should be determined in a twostage stochastic programming.
The above tables which specified by Table 2, Table 3 and Table 4 are illustrated the notations of the proposed model. In the following, the proposed formulations are introduced:
The above objective minimizes expected the total costs of the network. In this regard, according to the probability of each scenario, the first and the second terms compute the fixed costs of opening facilities. The third to seventh summations have estimated the cost of purchasing, manufacturing, handling, transportation, assignment costs. The last three terms state the profit value from RL network in disposal, remanufacture and recycle centers, respectively.
In continuous, the constraints of the model are illustrated. The following avouches ensure that the flow of products is maintained and the demands are convinced.(2)
The amount of products manufactured by each manufacturing center is computed by the following constraint:(8)
At each potential site for plants, at most one technology can be established:(9)
Each customer zone should be assigned to only one distribution center and recover center:(10)
The amount of products procured from each supplier is restricted by its capacity:(11)
A manufacturing center can manufacture only when it is opened and it has idle capacity:(12)
The flow of products through a facility is allowed only when the respective facility is operating and has enough capacity as well:(13)(14)(15)
The number of facilities in each echelon is restricted by a predefined number.(16)(17)
Finally, binary and positive variables are guaranteed in stage one:(18)
According to control the unfavorable scenario, the risk management is presented to help the Decision Maker (DM). The most popular way to manage the risk is defined a risk metric in the stochastic programming model, thus leading to a multiobjective optimization model where the total cost and the risk measure are two objective functions to be optimized. The downside risk can be formulated as follows:(20)(21)
V_{cϕ} is a positive variable that measures deviation between the scenario cost value (Cost^{c}) and a target ϕ. The downside risk (DRisk_{ϕ}) is defined as the expected value of the positive variable V_{cϕ}. Since the SCND model is designed under uncertain conditions, the risk associated with the supply chain network should be measured by financial risk. This feature helps the DM in their decision about the network design. In order to integrate this metric within the mathematical model, the set of constraints are specified by Eq. (22), (23) and (24) should be added as the second stage of the developed model which were formulated by Eq. (1)(19).
4.SOLUTION APPROACH
In order to exemplify the twostage stochastic programming model, the first one is approximated the value of total costs. Then, according to the cost of each scenario, the downside risk is measured to be optimized the two objectives simultaneously. As noted in the literature, this type of problem is known an NPhard one. After the metaheuristics algorithm presented, researchers found that for solving the NPhard problems, these methods are better than the exact methods to find the solutions. Because, they need the lower procedure time and search the all feasible space carefully by two search phases: the exploration and exploitation phases. So, ICA designed intelligently and distinguished the population of solutions which can give the opportunity to the user to make a tradeoff between the exploration and exploitation phases, is employed in this paper. Besides, the PSO, GA and εconstraint method are also considered to compare the results. In the following subsections, the proposed methods are detailed to address the problem.
4.1.The Structure of the Pareto Front Optimal Solutions
Due to proposed twostage stochastic programming, the front of the solutions set is consisted on two different goals. Thus, it needs a tradeoff between the objectives. This case considers the answer as a set solutions called Pareto optimal solutions set. In nondominated sorting as point out by Deb et al. (2000), a solution dominates the other solutions, when it had better value among all objective functions. Furthermore, selecting the best solution in first front of solutions is considered by calculating the crowding distance between solutions. In fact, the optimization problem is handled for each point. Finally, the nondominated solutions are eliminated.
4.2.Encoding Scheme of Metaheuristic Solution Planning
In the engineering problems, when a metaheuristic procedure is used. The authors need some approaches to code and decode the problem in the proposed metaheuristic. The procedure is so useful, if it had no procedure time to repair after an operator by metaheuristic is used to search the feasible space. According to encoding plan, this paper illustrates a twostage technique called RandomKey (RK) to solve developed discrete problem. Researchers have used this technique repeatedly during two recent decades (Hajiaghaei Keshteli and Aminnayeri 2014; Lotfi and Tavakkoli Moghaddam, 2013). This technique helps the authors to solve the proposed problem with various methods and operators. RK is a two phases technique, in which, at first, a solution is made by random numbers and then this solution is converted to a feasible discrete solution by a procedure. For the problem, the solutions are defined in three different subsolutions. The illustration of encoding subsolutions is de picted in Figure 2. As shown in this figure, at first a matrix with p elements obtained by uniform distribution U (0, 1) is made (Figure 2(a)). This subsolution is used for binary variables such as selection the distribution centers, recover centers, remanufacturing center and recycling centers and also disposal. On the other hand, this subsolution is used in assigning the customer areas to retailers and recover centers. Furthermore, conforming to Figure 2(b), the type of technology is determined for each plant. Eventually, Figure 2(c) is determined the flow of products. In the other words, $\left{P}_{2}\right\underset{\_}{\times}\leftD\right,\text{\hspace{0.17em}}\left{P}_{2}\right\le {P}_{Max}$ random matrix is produced. Afterwards, the columns of the second matrix are normalized to specify how distribution centers and other facilities supply their demands between different selected plants.
4.3.Genetic algorithm (GA)
Holland (1975) developed GA for the first time which is inspired by genetic evolutionary in 1975. GAs are the special type of EAs and include so many methods in this classification (i.e. Genetic programming (Alemdag et al., 2016), Scatter search (Glover, 1977) and Different evolution (Storn and Price, 1997)). In 2000, Deb et al. (2000) proposed the Nondominated Sorting Genetic Algorithm II (NSGAII) to solve the multiobjective optimization problems. In GA, chromosomes are the structure of cells in animals, plants and humans. In this regard, an array of variables which called chromosome is defined. Chromosomes are altered by two operators: mutation and crossover. So, some new solutions are created by these two mentioned operators (HajiaghaeiKeshteli, 2011).
In the Genetic algorithm, like other methods, in the first step, some random solutions in feasible space are initialized. In mutation, one solution changes to a new solution by generating a neighbor of this solution. In the crossover, the two chromosomes called parents are selected. They compose together and make two new solutions named offspring. These two operators are so simple. Hence, users can utilize a creative way to do these operators. However, GA is so easy and simple to code, but it has not any special plans to explore the potential areas. As discussed earlier, the tradeoff between the two phases is so significant. GA just does these two important phases by crossover and mutation operators which are blind in the search space as mentioned before (Mahmoodjanloo et al., 2016). The employed procedures for mutation and crossover operators are illustrated as seen in Figure 3 and Figure 4. Moreover, as shown in Figure 5, steps of the algorithm are explained.
4.4.Particle Swarm Optimization (PSO)
PSO introduced by Eberhart and Kennedy (1995) as a populationbased technique in swarm intelligence is inspired by social behavior of bird flocking and fish schooling (Kennedy, 2011). Due to consider the memory of requirements and speed, it has shown a good convergence performance behavior in many fields. Particles in PSO are changed by moving according to the particle’s position and velocity. PSO keeps the best value of all particles as gbest and also saves the best particle in the neighboring of each particle as pbest. In order to balance between the phases, this algorithm utilizes the weights of gbest and lbest to move each particle (TavakkoliMoghaddam et al., 2016). Figure 6 shows the pseudocode of proposed PSO.
4.5.Imperialist Competitive Algorithm (ICA)
ICA presented by AtashpazGargari and Lucas (2007) explains a story in social developments. Actually, ICA simulates the human social evolution. In this algorithm, country is the counterpart of a solution. And the countries are divided into two types: empires and their colonies. In the next step, in original of ICA, the assimilation policy is occurred. In the real world, the imperialist countries improve their colonies in special characteristics like language, economy and culture. This fact has been simulated by moving all colonies toward the imperialist. In the other words, in this step, some new regions are discovered in the neighboring of the imperialist. After evaluating the cost of colonies, if in an empire, there is a better solution than imperialist, the positions of imperialist and a mentioned colony are exchanged. In the next step of the algorithm, the total power of each empire is calculated. The total power is relied on the cost of imperialist and percentage of related colonies. After the procedure of computing the total power of empire, the weakest empire is identified. In this regard, the one of the weakest colony is selected in the weakest empire and give it (them) to empire that has the most verisimilitude to posses it (them). This fact is called as an imperialistic competition. Consequently, the empire which has no colony is deleted. At the end, stopping condition is may be the number of iterations or just one empire exists in the world.
As mentioned before, in recent developed algorithms, the author finds the tradeoff between two phases and focuses on intelligent search (Golmohamadi et al., 2017). ICA carries the diversification by forming different empires and does the exploitation in each empire (Bagher et al., 2011). Also, in assimilating, the similarity is the sense of exploration or local search on the imperialist’s neighbor. Figure 7 shows the pseudocode of ICA.
4.6.εConstraint Method
In the twostage stochastic programming, some approaches are proposed to investigate the stochastic model by Danzig (1995) and Beale (1995). In this paper, in order to check the metaheuristic results in the small sizes, the ε constraint method is presented which was introduced for the first time by Haimes et al. (1971). This algorithm is based on the one objective function to be optimized and to consider the other objectives as the constraints with allowable bounds. Consequently, the other Pareto optimal solutions are generated by modifying the bounds consecutively. According to the mentioned method, it can be formulated as follows:
Furthermore, a set of Pareto optimal solution are provided by changing the ε. Then the methodology of the exact solution for the problem can be summarized as follows:

(1) Set a value for the objective functions and ε.

(2) Solve the model which is illustrated as Eq. (25).

(3) Need the corresponding Pareto optimal solution.

(4) Repeat the (1) to (3) in order the other Pareto optimal solutions set.

(5) Investigate the quality of Pareto optimal solutions set by proposed evaluation metrics.

(6) Select the best proposed network among all solutions set.
5.COMPUTATIONAL RESULTS
5.1.Generating Data
In order to analyze the methods, some test problems are generated. The number of potential places in each echelon directly indicates the complexity of the problem. The test problems are divided into three classifications (i.e. small, medium and large). For each classification, five random problems are initialized. In this regard, the 15 test problems are examined. In each test problem, the 10 scenarios are defined. And also the probability of each scenario is specified by DM (user). The details of the test problems are shown in Table 5, Table 6 and Table 7. In order to be faired, the searching time for metaheuristics in each test problem is defined by considering the size of problems as depicted in Table 5. And also, it should be noted that the ε constraint method is examined only on the small size.
In addition, the maximum desired number of units in each level is computed by half of its potential units. The fixed opening cost and capacity of facilities are estimated by the approach of Devika et al. (2014). The fixed opening costs of facilities are approximated as follows:(26)(27)(28)
where ∝~U(35, 45) and δ ~ U(0, 0.1). In addition, the capacities of facilities are estimated as follows:(29)(30)(31)
where γ is the summation of demands which is limited by lower hand echelon, and also β and δ are coefficients distributed as U(0, 0.1). Furthermore, sc_{d}, sc_{r} and sc_{r} and sc_{s} are approximated as the average total of cost of purchasing, manufacturing, handling a product which is taken back to the respective stage in the forward chain. At the end, the value of target ϕ is depended on the size of problem and estimated with the expected cost value by the DM.
5.2.Evaluation Metrics
As mentioned before, due to analyze the solutions of methods, four assessment metrics are defined to evaluate the quality of Pareto optimal solutions (GarciaNajera and Bullinaria, 2011) Consequently, the following is illustrated the four efficiency evaluation metrics which are utilized in this paper:

DM: the Diversification Metric (DM) is proposed to measure the spread of nondominated solution set (Govindan et al., 2015). As it is obvious by the following equation, the more of DM is preferred by the user and shown the better capability of the algorithm.(32)
$$DM=\sqrt{{\displaystyle \sum _{i=1}^{n}\text{max}\left(\Vert {x}_{i}{}^{\prime}{y}_{i}{}^{\prime}\Vert \right)}}$$(32)where n is number of Pareto fronts members and $\Vert {x}_{i}{}^{\prime}{y}_{i}{}^{\prime}\Vert $ is the Euclidean distance between the best front of nondominated solutions, x_{i} and y_{i}.

SNS: the Spread of Nondominance Solution (SNS) assesses the standard of the distance of an ideal point from a nondominated set (Behnamian and Ghomi, 2011; Govindan et al., 2015). The higher value of this metric is shown the better solution quality by the following formula:(33)
$$SNS=\sqrt{\frac{{\displaystyle {\sum}_{i=1}^{n}{\left(\overline{c}{c}_{i}\right)}^{2}}}{n}}$$(33)where ${c}_{i}=\Vert \overrightarrow{{f}_{i}}\overrightarrow{{f}_{Ideal}}\Vert ,\text{\hspace{0.17em}}\overline{c}=\frac{{c}_{i}}{n},\text{\hspace{0.17em}}\overrightarrow{{f}_{Ideal}}=\left\{\text{min}\left({f}_{1}\right),\text{\hspace{0.17em}min}\left({f}_{2}\right),\text{\hspace{0.17em}}\cdots ,\text{\hspace{0.17em}min}\left({f}_{k}\right)\right\}$

POD: Percentage of Domination (POD) is firstly proposed on basis of the Coverage Metric (CM) and developed by Zitzler and Theile (1998). CM measures the extent to which one Pareto set X′′ is covered by another solution set X′ by comparing the number of solutions in X′′ covered by solutions in X′ . POD is presented to measure the ability of an algorithm to dominate the solutions of other algorithms. For more than two algorithms, all pair of combinations should be analyzed. To compute the POD, all of the nondominated solutions required by algorithms are composed in one Pareto set, then the percentage of solutions belonging to each method is computed. The algorithm with higher POD is better. It is calculated as follows:(34)
$$\begin{array}{l}POD\left({X}^{\prime},\text{\hspace{0.17em}}{X}^{\u2033},\text{\hspace{0.17em}}\cdots ,\text{\hspace{0.17em}}{X}^{n}\right)\\ \text{\hspace{1em}}=\frac{\left\left({x}_{ki}\in {X}_{i}\right)\right\exists {x}_{ki}\in {X}_{j},\text{\hspace{0.17em}}i\in j:\text{\hspace{0.17em}}{x}_{ki}\le {x}_{kj}}{\left{X}_{j}\right}\end{array}$$(34) 
DEA: the Data Envelopment Analysis (DEA) determines the efficiency of solutions (Wu et al., 2017). DEA is proposed by Charnes et al. (1978) for the first time and it is usually used to measure the performance of some selections with some attributes (Govindan et al., 2015). In this paper, this method is used to decide the efficiency of nondominated solutions considered by each method. To compute DEA, each Pareto optimal solution can be considered as a DMU. Furthermore, an objective function such as total costs may be considered as input, while the other objective is assumed as output (Pakkar, 2016). In addition, all nondominated solutions used by algorithms are composed and the efficiency of these points is computed by DEA (Amin and Toloo, 2007). Like other above metrics, the higher value is better and preferable.
5.3.Parameters Setting
To evaluate the performance of any metaheuristic, it should be used a significant section to set the parameters (Vahdani and Zandie, 2010). As discussed earlier, in each presented metaheuristic, it is so necessary to balance between the phases and tune the parameters. Here we utilize RSM introduced by Box and Wilson (Mayers et al., 2009). In this method, we define for each factor (X_{i}) which is measured at two levels that can be coded as 1 to 1 the low level (x_{l}) and high level (x_{h}) of each variable, respectively (Bezerra et al., 2008). The following equation shows the independent variables which are used in our study:(35)
where K is the number of variables. The response surface when variables are assumed measurable can be declared as $y={f}_{\left({x}_{1},{x}_{2},\text{\hspace{0.17em}}\cdots ,\text{\hspace{0.17em}}{x}_{n}\right)}$. And the goal is to optimize response variables. In this way, we define y which is utilized to describe the variation in response variables, as follows:
where y, β_{0}, β_{j}, β_{ij}, are an estimated response, a constant, the linear coefficient, and the interaction coefficient, respectively. β_{jj} value, linear coefficient, interaction coefficient and quadratic coefficient. In addition, equation (36) should be assumed while there is curvature in the system. The factors and their levels as well as the number of experiments are shown in Table 8. In this table, the lower limit and upper limit are written on the left and right, respectively.
Since some performance evaluation metrics are utilized together to make quantitative comparisons in multiobjective metaheuristics, such methods will be tuned by a medium problem (i.e. P8) to set the input value of factors based on the represented evaluation metrics as illustrated in Table 8 simultaneously, which leads to desired response in multiobjective decision making (MODM). In this way, Derringer and Suich (1980) presented a utility function according to Eq. (37) to optimize the multiple responses, as following equations:
where response functions in form of minimum y_{i} are transformed to utility function d_{i} (y_{i}). l_{i} and h_{i} are lower and upper bounds, respectively. s measures the severity of d_{i} (y_{i}) by which the shape of the desirability can be arrayed for each goal by severity value. Severity is utilized to put emphasis on the lower/upper bounds. If s = 1, d_{i} (y_{i}) will array from 0 to 1 in a linear fashion. When s is between 0.1 and 1, less emphasis is placed on the goal. Conversely, when s is between 1 and 10, the more emphasis is on the goal function. Hence, parameters of s are set to 1, 1, 2 and 3 for DM, DEA, SNS and POD, due to their importance, respectively. Besides, desirability is estimated through the following equation, where m is the number of objectives:(38)
Also, it is obvious the higher value of D displays the better capability for algorithms. The tuned values for parameters, Rsquared (R^{2}) and desirability (D), are approximated as displayed in Table 9.
5.4.Implementation and Computational Results
5.4.1.Comparing the Effectiveness and Efficiency of Proposed Methods
This subsection aims to probe the effectiveness and efficiency of the presented algorithms. Due to it, each metaheuristic is performed in all test problems for 30 times runs. In this case, the behavior of algorithms in the two objective functions during the thirty run times is considered in Figure 8. In this figure, a test problem such as P7 is selected. In addition, the best solution is considered by computing the crowding distance as mentioned earlier. Finally, the average of outputs is saved. First of all, the proposed evaluation metrics are calculated as given in Table 10 and Table 11.
Later, the obtained results for each problem are converted to Relative Percentage Deviation (RPD). The RPD is considered by this formula:(39)
where Alg_{sol} is the output of algorithm and Best_{sol} is the best value ever found in size of problem. It should be noted that the lower value for RPD is preferred.
Figure 9 illustrates the nondominated solutions for one run which are reached by algorithms in P8 and P13 as instances. It should be noted that the value of the target φ is estimated 5,000 and 13,000 to calculate the downside risk in the second, respectively. As depicted in this figure, the number of solutions in best Pareto front is set to 10. As it can be seen, the solutions of ICA approximately overcome the other algorithms solutions. Furthermore, to verify the statistical validity of the results, authors have performed an analysis of variance (ANOVA) to accurately analyze the results (as seen in Table 10 and Table 11 based on the RPD). The results demonstrate that there is a clear statistically significant difference between the performances of the algorithms. The means plot and LSD intervals (at the 95% confidence level) for all methods are shown in Figure 10 and Figure 11. The first figure explores the metaheuristics and the εconstraint method. Then, Figure 11, the algorithms are explored for the medium and large of the size of test problems based on the evaluation metrics. The results show the superiority of the ICA in this study.
5.4.2.Sensitivity Analyses of the Downside Risk Management
In this section, three sample test problems such as P3, P8 and P13 are selected. The data of this sample problems are similar to a real case study involves from a Food company in Mazandaran province, Iran. The aim of this part shows that the authors only focused on the downside risk to do some sensitivity analyses. Due to risk management, we set two phases. In the first step, the value of risk and the target of φ are assumed as 0, respectively. In this way, the total cost is computed. Consequently, by considering the risk management, we calculate the objectives again. As mentioned earlier, the proposed powerful algorithm in this paper is known as ICA. Hence, this algorithm is utilized in this section to perform the analysis. In addition, in each front, the best solution is selected.
The outputs of proposed algorithm before and after the downside risk are given in Table 12. In order to emphasize the changes of cost for the downside risk model, a metric is presumed as Eq (40) in the following:
As can be seen in the above table, after the risk management, most of the expected value of the cost is increased. So, the decision maker should take a suitable decision according to the downsie risk value and the other objectives. And the probability of each scenario can be changed by the value of downside risk via decision maker's preference as seen in Figure 12. This figure shows the decision to the change the probability for a large scale test problem numbered as P13.
6.CONCLUSION AND FUTURE WORKS
This work addresses the metaheuristics procedures and εconstraint method the twostage stochastic strategically CLSC configuration. In this paper, the authors had some special plans to handle this problem. ICA, PSO and GA as three powerful metaheuristics in the related researchers are considered in this study. In the literature, a classification for the stochastic models in SCND is provided. In addition, the authors highlight the need of utilizing the metaheuristic approach to solve the stochastic SCND models which show still scare in the literature. Furthermore, Response Surface Method (RSM) with an MODM approach was used to adjust the parameters of proposed methods and better performance of them. In this study by utilizing the evaluation metrics which introduced in the related researches the effectiveness and efficiency of methods are measured. In order to be faired, the searching time for each metaheuristic is specified and is based on the size of problems as long as increased. As a result, ICA had shown better capability to PSO and GA. At the end of the paper, some sensitivity analyses are obtained to explore the characteristics of the proposed model.
For future works, more comprehensive analysis on the proposed model may be required to be explored. Furthermore, some other real study cases such as: mine or gold or glass industry can be tested on the model. Besides, to expand the mathematical model and the suggestions for new works, the model can be added by more reallife constraints such as: crossdocking operations and vehicle routing practicing which can be diminished the transportation total cost by utilizing the same facility in the forward and reverse logistics and can also be abated the capacity of vehicles which is unused, all of them by considering a stochastic model. In addition, researchers can add the environmental and social aspects to develop the proposed model. Also, some other recent powerful metaheuristics can be utilized to probe the results of this paper.