Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1598-7248 (Print)
ISSN : 2234-6473 (Online)
Industrial Engineering & Management Systems Vol.16 No.3 pp.342-362
DOI : https://doi.org/10.7232/iems.2017.16.3.342

A Bi-Objective Stochastic Closed-loop Supply Chain Network Design Problem Considering Downside Risk

Amir Mohammad Fathollahi Fard, Fatemeh Gholian-Jouybari, Mohammad Mahdi Paydar, Mostafa Hajiaghaei-Keshteli*
Department of Industrial Engineering, University of Science and Technology of Mazandaran, Behshahr, Iran
Department of Industrial Engineering, Shomal University, Amol, Iran
Department of Industrial Engineering, Babol Noshirvani University of Technology, Babol, Iran
Department of Industrial Engineering, University of Science and Technology of Mazandaran, Behshahr, Iran
Corresponding Author, mostafahaji@mazust.ac.ir
February 2, 2017 March 26, 2017 May 1, 2017

ABSTRACT

This paper deals with a closed-loop supply chain network and proposes new approaches in metaheuristics and exact methods as solution methods. Moreover, the downside risk is incorporated into the objective functions as a risk measure. Hence, the developed two-stage stochastic model aims to minimize the expected total cost and the downside risk, simultaneously. Besides, this study presents the closed-loop 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). In order to address the problem, ICA, PSO, GA, and also ε-constraint method, are utilized. 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. At the end, the results show the capability of ICA through the most of the tests problem. According to the risk management, a real data set is used to do some sensitivity analyses on the proposed model.


초록


    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; Hajiaghaei-Keshteli 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 closed-loop 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 closed-loop 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 bi-level stochastic model by considering the economic and risk measure to be minimized. Due to obvious the different parts of proposed closed-loop network and flow of the products, a graphical structure is provided as an example in the Figure 1.

    In this work, a closed-loop supply chain network is considered and solved by new approaches in metaheuristics and exact methods. The developed two-stage stochastic model aims to minimize the expected total cost and the downside risk, simultaneously. Besides, this study presents the closed-loop 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 two-stage 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 end-oflife (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 NP-hard. 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 closed-loop 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 multi-mode 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 multi-objective closed-loop 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 closed-loop supply chain network. They minimized the total cost and environmental concerns in the deterministic model.

    Since the closed-loop 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 closed-loop SCND. As mentioned in Table 1, only a few papers use the two-stage stochastic programming to address new contributions of this area like the multi-objective, multi-echelon, multi-type 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 closed-loop network with two different objectives. And the last but not least, a two-stage stochastic programming is introduced in this work as well.

    Due to the multi-objective 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 NP-hard problem (Eckert and Gottlieb, 2002; Jo et al., 2007; Syarif et al., 2002). Subsequently, the powerful and well-known approaches in metaheuristics are considered to find the Pareto-optimal solutions. The Imperialist Competitive Algorithm (ICA) proposed by Atashpaz-Gargari 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 closed-loop SCND problem. In addition to check the algorithm’s results in small sizes, the ε-constraint method is utilized to tackle such NP-hard problem, too.

    3.PROBLEM DESCRIPTION

    3.1.Assumptions of the Proposed Problem

    This work investigates the multi-objective, multiechelon, multi-type of production technology in the stochastic closed-loop 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 closed-loop 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 Two-Stage Stochastic Programming

    In this sub-section, 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 two-stage 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:

    M i n E ( C o s t ) = c p b c ( p t f c p t Y p t + e f c f c f e Y f e ) + m p p c m c X m p c + p t m c p t c H p t c + e e f e f e t c f e f e X f e f e c + i n a c n i c d i c + i j z i j c v c j i d i c s c s ( s n X s n c ) s c r ( r n X r n c ) s c d ( d p X d p c )
    (1)

    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)

    m X m p c c , p
    (2)
    (3)(4)(5)(6)(7)

    p X p n c = i Z n i c d i c c , n
    (3)

    f e X j f e c = i α i d i c Z i j c c , j , e { d , r , s }
    (4)

    m X s m C = j X j s C c , s
    (5)

    p X d p c = j X j d c c , d
    (6)

    n X r n c = j X j r c c , r
    (7)

    The amount of products manufactured by each manufacturing center is computed by the following constraint:(8)

    t H p t c = m X m p c c , p
    (8)

    At each potential site for plants, at most one technology can be established:(9)

    t Y p t 1 p
    (9)

    Each customer zone should be assigned to only one distribution center and recover center:(10)

    n Z n i c = j Z i j c c , i
    (10)

    The amount of products procured from each supplier is restricted by its capacity:(11)

    n X m p c p m c , p
    (11)

    A manufacturing center can manufacture only when it is opened and it has idle capacity:(12)

    H p t c p p t Y p t p , t , c
    (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)

    p X p n c p n Y n c , k
    (13)

    i i d i c Z i j c p j X j c , j
    (14)

    j X j f e c p f e Y f e e { d , r , s }
    (15)

    The number of facilities in each echelon is restricted by a pre-defined number.(16)(17)

    p t Y p t max p
    (16)

    f e Y f e max e e { d , m , n , p , r , s }
    (17)

    Finally, binary and positive variables are guaranteed in stage one:(18)

    Y p t , Y f e , Z n i c , Z i j c { 0 , 1 }
    (18)

    X f e f e c , H p t c 0
    (19)

    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 multi-objective 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)

    D R i s k φ = E ( V c φ )
    (20)

    w h e r e V c φ = { C o s t c φ i f C o s t c > φ c 0     o t h e r w i s e }
    (21)

    V is a positive variable that measures deviation between the scenario cost value (Costc) and a target ϕ. The downside risk (DRiskϕ) is defined as the expected value of the positive variable V. 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).

    M i n D R i s k φ = c p b c V c φ
    (22)

    V c φ C o s t c φ , c
    (23)

    V c φ 0 ,    c
    (24)

    4.SOLUTION APPROACH

    In order to exemplify the two-stage 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 NP-hard one. After the metaheuristics algorithm presented, researchers found that for solving the NP-hard 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 trade-off 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 sub-sections, the proposed methods are detailed to address the problem.

    4.1.The Structure of the Pareto Front Optimal Solutions

    Due to proposed two-stage stochastic programming, the front of the solutions set is consisted on two different goals. Thus, it needs a trade-off between the objectives. This case considers the answer as a set solutions called Pareto optimal solutions set. In non-dominated 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 two-stage technique called Random-Key (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 sub-solutions 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 sub-solution 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 sub-solution 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, | P 2 | × _ | D | , | P 2 | P M a x 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 Non-dominated Sorting Genetic Algorithm II (NSGA-II) to solve the multi-objective 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 (Hajiaghaei-Keshteli, 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 trade-off 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 population-based 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 (Tavakkoli-Moghaddam et al., 2016). Figure 6 shows the pseudo-code of proposed PSO.

    4.5.Imperialist Competitive Algorithm (ICA)

    ICA presented by Atashpaz-Gargari 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 trade-off 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 pseudo-code of ICA.

    4.6.ε-Constraint Method

    In the two-stage 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:

    min  f 1 s . t .    E q . ( 2 ) ( 19 )  and  ( 23 ) ( 24 )    f 2 ε
    (25)

    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)

    f c f m e a n = × max f F × f = 1 F f = 1 F t c f f F × F     f = { d , j , n , p , r , s }
    (26)

    f c f min = [ ( 1 δ ) × f c f m e a n ]  and f c f max = [ ( 1 + δ ) × f c f m e a n ] , f = d j n p r s }
    (27)

    f c f = U ( f c f min , f c f max )
    (28)

    where ∝~U(35, 45) and δ ~ U(0, 0.1). In addition, the capacities of facilities are estimated as follows:(29)(30)(31)

    p f m e a n = γ × ( 1 + β ) max f , f = { d , j , n , p , r , s }
    (29)

    p f min = [ ( 1 δ ) × p f m e a n ]  and p f max = [ ( 1 + δ ) ] × p f m e a n , f = { d , j , n , p , r , s }
    (30)

    p f = U ( p f min , p f max ) , f = { d , j , n , p , r , s }
    (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, scd, scr and scr and scs 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 (Garcia-Najera 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 non-dominated 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)

      D M = i = 1 n max ( x i y i )
      (32)

      where n is number of Pareto fronts members and x i y i is the Euclidean distance between the best front of non-dominated solutions, xi and yi.

    • SNS: the Spread of Non-dominance Solution (SNS) assesses the standard of the distance of an ideal point from a non-dominated 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)

      S N S = i = 1 n ( c ¯ c i ) 2 n
      (33)

      where c i = f i f I d e a l , c ¯ = c i n , f I d e a l = { min ( f 1 ) ,  min ( f 2 ) , ,  min ( f k ) }

    • 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 non-dominated 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)

      P O D ( X , X , , X n ) = | ( x k i X i ) | x k i X j , i j : x k i x k j | | X j |
      (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 non-dominated 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 (Xi) which is measured at two levels that can be coded as -1 to 1 the low level (xl) and high level (xh) of each variable, respectively (Bezerra et al., 2008). The following equation shows the independent variables which are used in our study:(35)

    x i = 2 X i x h x l x h x l i = { 1 , ​  2 , ,    K }
    (35)

    where K is the number of variables. The response surface when variables are assumed measurable can be declared as y = f ( x 1 , x 2 , , x n ) . 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:

    y = β 0 + j = 1 k β j x j + j = 1 k i < j k β i j x i x j + j = 1 k β j x j j 2 + ε
    (36)

    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 multi-objective 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:

    d i ( y i ) = ( h i y i h i l i ) ,    l i < y i < h i
    (37)

    where response functions in form of minimum yi are transformed to utility function di (yi). li and hi are lower and upper bounds, respectively. s measures the severity of di (yi) 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, di (yi) 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)

    D = d 1 ( y 1 ) × d 2 ( y 2 ) × × d m ( y m ) m
    (38)

    Also, it is obvious the higher value of D displays the better capability for algorithms. The tuned values for parameters, R-squared (R2) 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 sub-section 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)

    R P D = | A lg s o l B e s t s o l | B e s t s o l
    (39)

    where Algsol is the output of algorithm and Bestsol 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 non-dominated 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:

    G a p = 100 × | V a l u e a f t e r t h e r i s k V a l u e b e f o r e t h e r i s k V a l u e b e f o r e t h e r i s k |
    (40)

    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 two-stage 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 real-life constraints such as: cross-docking 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.

    Figure

    IEMS-16-342_F1.gif

    The graphical structure of a simple closed-loop network utilizing in this paper.

    IEMS-16-342_F2.gif

    The graphical illustration of structure solution with three parts, (a) Facilities selection sub-solution, (b) Technology levels selection and customer zones allocation sub-solution, (c) Flow of products sub-solution.

    IEMS-16-342_F3.gif

    The crossover operators in this paper.

    IEMS-16-342_F4.gif

    The mutation operators in this study.

    IEMS-16-342_F5.gif

    The pseudo-code of GA.

    IEMS-16-342_F6.gif

    The pseudo-code of proposed PSO.

    IEMS-16-342_F7.gif

    The pseudo-code of ICA.

    IEMS-16-342_F8.gif

    The behavior of algorithms in terms of expected total cost (a) and downside risk (b) among thirty runs.

    IEMS-16-342_F9.gif

    The dispersion of non-dominated solutions obtained by different algorithms in samples problems (i.e. P8(a) and P13(b)).

    IEMS-16-342_F10.gif

    Means plot and LSD intervals to specify RPD for evaluation parameters in the small sizes (i.e. DM (a), SNS (b), DEA (c) and POD (d)).

    IEMS-16-342_F11.gif

    Means plot and LSD intervals to specify RPD for evaluation parameters in the medium and large sizes (i.e. DM (a), SNS (b), DEA (c) and POD (d)).

    IEMS-16-342_F12.gif

    The probability of the cost for each scenario after and before the risk management in P13 by DM.

    Table

    The literature review for stochastic models in SCND

    Model Notation

    Model Notation

    Model Notation

    Level, number and size of problems.

    Parameters and their surfaces

    Parameters and their surface

    Algorithms and factor levels with number of experiments

    nPop = number of population, Nemp = number of empires, ε = colonies mean cost coefficient, PAS = rate of assimilation, PR= probability of revolution, W = inertia weight, C1 = acceleration coefficient of local optimum, C2 = acceleration coefficient of global optimum, PC = probability of crossover, PM = probability of mutation, TC = Type of crossover (ceil (0, 1) = single-point; ceil (1, 2) = doublepoint; ceil (2, 3) = uniform), TM = Type of mutation (ceil (0, 1) = swap; ceil (1, 2) = reversion; ceil (2, 3) = insertion).

    Final values of algorithms parameters, R-squared (R2) and desirability (D)

    The evaluation metrics to algorithms performance (DM, SNS, DEA and POD) for test problems in small size

    The evaluation metrics to algorithms performance (DM, SNS, DEA and POD) for test problems in medium and large sizes

    Cost comparison before and after the downside risk

    REFERENCES

    1. Afrouzy Z.A. , Nasseri S.H. , Mahdavi I. (2016) A genetic algorithm for supply chain configuration with new product development , Comput. Ind. Eng, Vol.101 ; pp.440-454
    2. Alemdag S. , Gurocak Z. , Cevik A. , Cabalar A.F. , Gokceoglu C. (2016) Modeling deformation modulus of a stratified sedimentary rock mass using neural network, fuzzy inference and genetic programming , Eng. Geol, Vol.203 ; pp.70-82
    3. Amin G.R. , Toloo M. (2007) Finding the most efficient DMUs in DEA: An improved integrated model , Comput. Ind. Eng, Vol.52 (1) ; pp.71-77
    4. Aminnayeri M. , Hajiaghaei-Keshteli M. (2012) The cost function for a two-warehouse and three-echelon inventory system , International Journal of Industrial Engineering & Production Research, Vol.23 (4) ; pp.285-290
    5. Ardalan Z. , Karimi S. , Naderi B. , Khamseh A.A. (2016) Supply chain networks design with multimode demand satisfaction policy , Comput. Ind. Eng, Vol.96 ; pp.108-117
    6. Atashpaz-Gargari E. , Lucas C. (2007) Imperialist competitive algorithm: An algorithm for optimization inspired by imperialistic competition , IEEE Congress on Evolutionary Computation, ; pp.4661-4667
    7. Bagher M. , Zandieh M. , Farsijani H. (2011) Balancing of stochastic U-type assembly lines: An imperialist competitive algorithm , Int. J. Adv. Manuf. Technol, Vol.54 (1-4) ; pp.271-285
    8. Beale E.M. (1955) On minimizing a convex function subject to linear inequalities , J. R. Stat. Soc. B, Vol.17 ; pp.173-184
    9. Behnamian J. , Ghomi S.M. (2011) Hybrid flowshop scheduling with machine and resourcedependent processing times , Appl. Math. Model, Vol.35 (3) ; pp.1107-1123
    10. Bezerra M.A. , Santelli R.E. , Oliveira E.P. , Villar L.S. , Escaleira L.A. (2008) Response surface methodology (RSM) as a tool for optimization in analytical chemistry , Talanta, Vol.76 (5) ; pp.965-977
    11. Charnes A. , Cooper W.W. , Rhodes E. (1978) Measuring the efficiency of decision making units , Eur. J. Oper. Res, Vol.2 (6) ; pp.429-444
    12. Danzig G.B. (1955) Linear programming under uncertainty , Manage. Sci, Vol.1 (3/4) ; pp.197-206
    13. Daskin M.S. , Sriyder L.V. , Berger R.T. Langevin A. , Riopel D. (2003) Logistics Systems: Design and Optimization, Kluwer,
    14. Deb K. , Agrawal S. , Pratap A. , Meyarivan T. (2000) A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II , International Conference on Parallel Problem Solving From Nature, ; pp.849-858
    15. Dehghanian F. , Mansour S. (2009) Designing sustainable recovery network of end-of-life products using genetic algorithm , Resour. Conserv. Recycling, Vol.53 (10) ; pp.559-570
    16. Derringer G. , Suich R. (1980) Simultaneous optimization of several response variables , J. Qual. Technol, Vol.12 (4) ; pp.214-219
    17. Devika K. , Jafarian A. , Nourbakhsh V. (2014) Designing a sustainable closed-loop supply chain network based on triple bottom line approach: A comparison of metaheuristics hybridization techniques , Eur. J. Oper. Res, Vol.235 (3) ; pp.594-615
    18. Eberhart R. , Kennedy J. (1995) A new optimizer using particle swarm theory , Proceedings of the Sixth International Symposium on, ; pp.39-43
    19. Eckert C. , Gottlieb J. Guervós J. , Adamidis P. , Beyer H-G. , Schwefel H.P. , Fernández-Villacañas J-L. (2002) Parallel problem solving from nature, - PPSN VII, Springer, ; pp.77-87
    20. Elhedhli S. , Merrick R. (2012) Green supply chain network design to reduce carbon emissions , Transp. Res. Part D Transp. Environ, Vol.17 (5) ; pp.370-379
    21. Fleischmann M. , Beullens P. , Bloemhof-Ruwaard J.M. , Van Wassenhove L.N. (2001) The impact of product recovery on logistics network design , Prod. Oper. Manag, Vol.10 (2) ; pp.156-173
    22. Fonseca M.C. , García-Sánchez Á. , Ortega-Mier M. , Saldanha-da-Gama F. (2010) A stochastic biobjective location model for strategic reverse logistics , Top (Madr.), Vol.18 (1) ; pp.158-184
    23. Garcia-Najera A. , Bullinaria J.A. (2011) An improved multi-objective evolutionary algorithm for the vehicle routing problem with time windows , Comput. Oper. Res, Vol.38 (1) ; pp.287-300
    24. Glover F. (1977) Heuristics for integer programming using surrogate constraints , Decis. Sci, Vol.8 (1) ; pp.156-166
    25. Golmohamadi S. , Tavakkoli-Moghaddam R. (2017) Solving a fuzzy fixed charge solid transportation problem using batch transferring by new approaches in meta-heuristic , Electron. Notes Discrete Math, Vol.58 ; pp.143-150
    26. Govindan K. , Jafarian A. , Nourbakhsh V. (2015) Biobjective integrating sustainable order allocation and sustainable supply chain network strategic design with stochastic demand using a novel robust hybrid multi-objective metaheuristic , Comput. Oper. Res, Vol.62 ; pp.112-130
    27. Govindan K. , Paam P. , Abtahi A. (2016) A fuzzy multi-objective optimization model for sustainable reverse logistics network design , Ecol. Indic, Vol.67 ; pp.753-768
    28. Haimes Y.Y. , Lasdon L.S. , Wismer D.A. (1971) On a bicriterion formulation of the problems of integrated system identification and system optimization , IEEE Trans. Syst. Man Cybern, Vol.1 (3) ; pp.296-297
    29. Hajiaghaei-Keshteli M. (2011) The allocation of customers to potential distribution centers in supply chain network; GA and AIA approaches , Appl. Soft Comput, Vol.11 (2) ; pp.2069-2078
    30. Hajiaghaei-Keshteli M. , Aminnayeri M. (2014) Solving the integrated scheduling of production rail transportation problem by Keshtel algorithm , Appl. Soft Comput, Vol.25 ; pp.184-203
    31. Hajiaghaei-Keshteli M. , Sajadifar S.M. (2010) Deriving the cost function for a class of three-echelon inventory system with N-retailers and one-for-one ordering policy , Int. J. Adv. Manuf. Technol, Vol.50 (1-4) ; pp.343-351
    32. Hajiaghaei-Keshteli M. , Aminnayeri M. , Fatemi Ghomi S.M. (2014) Integrated scheduling of production and rail transportation , Comput. Ind. Eng, Vol.74 ; pp.240-256
    33. Hajiaghaei-Keshteli M. , Molla-Alizadeh-Zavardehi S. , Tavakkoli-Moghaddam R. (2010) Addressing a nonlinear fixed-charge transportation problem using a spanning tree-based genetic algorithm , Comput. Ind. Eng, Vol.59 (2) ; pp.259-271
    34. Hajiaghaei-Keshteli M. , Sajadifar S.M. , Haji R. (2011) Determination of the economical policy of a three-echelon inventory system with (R, Q) ordering policy and information sharing , Int. J. Adv. Manuf. Technol, Vol.55 (5-8) ; pp.831-841
    35. Holland J.H. (1975) Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence, University of Michigan Press,
    36. Hsu H-W. , Wang H-F. Wang H-F. (2009) Web-based green products life cycle management systems: Reverse supply chain utilization, IGI Global Publication, ; pp.268-282
    37. Jabbarzadeh A. , Pishvaee M.S. , Papi A. (2016) A multi-period fuzzy mathematical programming model for crude oil supply chain network design considering budget and equipment limitations , Journal of Industrial and Systems Engineering, Vol.9 ; pp.88-107
    38. Jayaraman V. , Pirkul H. (2001) Planning and coordination of production and distribution facilities for multiple commodities , Eur. J. Oper. Res, Vol.133 (2) ; pp.394-408
    39. Jo J-B. , Li Y. , Gen M. (2007) Nonlinear fixed charge transportation problem by spanning treebased genetic algorithm , Comput. Ind. Eng, Vol.53 (2) ; pp.290-298
    40. Kannan D. , Diabat A. , Alrefaei M. , Govindan K. , Yong G. (2012) A carbon footprint based reverse logistics network design model , Resour. Conserv. Recycling, Vol.67 ; pp.75-79
    41. Kennedy J. (2011) Encyclopedia of Machine Learning, Springer,
    42. Listeş O. , Dekker R. (2005) A stochastic approach to a case study for product recovery network design , Eur. J. Oper. Res, Vol.160 (1) ; pp.268-287
    43. Lotfi M.M. , Tavakkoli-Moghaddam R. (2013) A genetic algorithm using priority-based encoding with new operators for fixed charge transportation problems , Appl. Soft Comput, Vol.13 (5) ; pp.2711-2726
    44. Mahmoodjanloo M. , Parvasi S.P. , Ramezanian R. (2016) A tri-level covering fortification model for facility protection against disturbance in r-interdiction median problem , Comput. Ind. Eng, Vol.102 ; pp.219-232
    45. Mayers R.H. , Montgomery C.D. , Anderson-Cook C.M. (2009) Response surface methodology, hoboken, New Jersey: John Wiley & Sons , Inc, Vol.2 ; pp.38
    46. Nasiri E. , Afshari A.J. , Hajiaghaei-Keshteli M. (2017) Addressing the freight consolidation and containerization problem by recent and hybridized meta-heuristic algorithms, International Journal of Engineering-Transactions C , Aspects, Vol.30 (3) ; pp.403-410
    47. Nourmohamadi Shalke P. , Paydar M. M. , Hajiaghaei-Keshteli M. (2017) Sustainable supplier selection and order allocation through quantity discounts , International Journal of Management Science and Engineering Management, ; pp.1-13
    48. Pakkar M.S. (2016) Using DEA and AHP for hierarchical structures of data , Industrial Engineering and Management Systems, Vol.15 (1) ; pp.49-62
    49. Pishvaee M. , Rabbani M. , Torabi S.A. (2011) A robust optimization approach to closed-loop supply chain network design under uncertainty , Appl. Math. Model, Vol.35 (2) ; pp.637-649
    50. Pokharel S. , Mutha A. (2009) Perspectives in reverse logistics: A review , Resour. Conserv. Recycling, Vol.53 (4) ; pp.175-182
    51. Rajabi F. , Najafi S.E. , Hajiaghaei-Keshteli M. , Molla-Alizadeh-Zavardehi S. (2013) Solving fuzzy step fixed charge transportation problems via metaheuristics , International Journal of Research in Industrial Engineering, Vol.2 (3) ; pp.24-34
    52. Ramezani M. , Bashiri M. , Tavakkoli-Moghaddam R. (2013) A new multi-objective stochastic model for a forward/reverse logistic network design with responsiveness and quality level , Appl. Math. Model, Vol.37 (1) ; pp.328-344
    53. Rezaee A. , Dehghanian F. , Fahimnia B. , Beamon B. (2015) Green supply chain network design with stochastic demand and carbon price , Ann. Oper. Res, Vol.••• ; pp.1-23
    54. Sabri E.H. , Beamon B.M. (2000) A multi-objective approach to simultaneous strategic and operational planning in supply chain design , Omega, Vol.28 (1) ; pp.581-598
    55. Sasikumar P. , Kannan G. (2008) Issues in reverse supply chains, Part I: End-of-life product recovery and inventory management - An overview , Int. J. Sustain. Eng, Vol.1 (3) ; pp.154-172a
    56. Sasikumar P. , Kannan G. (2008) Issues in reverse supply chain, Part II: Reverse distribution issues - An overview , Int. J. Sustain. Eng, Vol.1 (4) ; pp.234-249b
    57. Sasikumar P. , Kannan G. (2009) Issues in reverse supply chain, Part III: Classification and simple analysis , Int. J. Sustain. Eng, Vol.2 (1) ; pp.2-27
    58. Schultmann F. , Zumkeller M. , Rentz O. (2006) Modeling reverse logistic tasks within closed-loop supply chains: An example from the automotive industry , Eur. J. Oper. Res, Vol.171 (3) ; pp.1033-1050
    59. Soleimani H. , Kannan G. (2015) A Hybrid particle swarm optimization and genetic algorithm for closedloop supply chain network design in large-scale networks , Appl. Math. Model, Vol.39 (14) ; pp.3990-4012
    60. Soleimani H. , Esfahani M.S. , Govindan K. (2014) Incorporating risk measures in closed-loop supply chain network design , Int. J. Prod. Res, Vol.52 (6) ; pp.1843-1867
    61. Spar D.L. , La Mure L.T. (2003) The power of activism: Assessing the impact of NGOs on global business , Calif. Manage. Rev, Vol.45 (3) ; pp.78-101
    62. Srivastava S.K. (2008) Network design for reverse logistics , Omega, Vol.36 (4) ; pp.535-548
    63. Storn R. , Price K.V. (1997) Differential evolution: A simple evolution strategy for fast optimization , Dr. Dobb’s J, Vol.22 (4) ; pp.18-24
    64. Syarif A. , Yun Y. , Gen M. (2002) Study on multistage logistic chain network: A spanning tree-based genetic algorithm approach , Comput. Ind. Eng, Vol.43 (1-2) ; pp.299-314
    65. Talaei M. , Moghaddam B. F. , Pishvaee M. S. (2016) A robust fuzzy optimization model for carbon-efficient closedloop supply chain network design problem: A numerical illustration in electronics industry , J. Clean. Prod, Vol.113 ; pp.662-673
    66. Tavakkoli-Moghaddam R. , Hajiaghaei-Keshteli M. , Mousavi M. , Ranjbar-Bourani M. (2016) Two meta-heuristics to solve a coordinated air transportation and production scheduling problem with time windows for the due date , Systems, Man, and Cybernetics SMC), 2016 IEEE International Conference on, IEEE,
    67. Vahdani B. , Zandieh M. (2010) Scheduling trucks in cross-docking systems: Robust meta-heuristics , Comput. Ind. Eng, Vol.58 (1) ; pp.12-24
    68. Van Der Laan E. , Salomon M. , Dekker R. , Van Wassenhove L. (1999) Inventory control in hybrid systems with remanufacturing , Manage. Sci, Vol.45 (5) ; pp.733-747
    69. Wang H.F. , Hsu H.W. (2010) A closed-loop logistic model with a spanning-tree based genetic algorithm , Comput. Oper. Res, Vol.37 (2) ; pp.376-389
    70. Wu D. , Ding W.D. , Koubaa A. , Chaala A. , Luo C. (2017) Robust DEA to assess the reliability of methyl methacrylate-hardened hybrid poplar wood , Ann. Oper. Res, Vol.248 (1-2) ; pp.515-529
    71. Yang G.Q. , Liu Y.K. , Yang K. (2015) Multiobjective biogeography-based optimization for supply chain network design under uncertainty , Comput. Ind. Eng, Vol.85 ; pp.145-146
    72. Yao M.J. , Hsu H.W. (2009) A new spanning tree-based genetic algorithm for the design of multi- stage supply chain networks with nonlinear transportation costs , Optim. Eng, Vol.10 (2) ; pp.219-237
    73. Yi P. , Huang M. , Guo L. , Shi T. (2016) A retailer oriented closed-loop supply chain network design for end of life construction machinery remanufacturing , Journal of Cleaner Production,
    74. Zitzler E. , Thiele L. (1998) An evolutionary algorithm for multi-objective optimization: The strength Pareto approach , Swiss Federal Institute of Technology, TIK-Report,