Designing a Closedloop Supply Chain Network: A Comparison of Algorithms
Designing a closedloop supply chain network: A comparison of algorithms
Abstract
This paper attempts to design a closedloop supply chain network (CLSCN). To solve this mixed integer linear programming (MILP) model, Lagrangian relaxation (LR), Benders decomposition (BD) and Simulated Annealing (SA) algorithms are applied as well as Lshaped accelerating. The model determines the location of various types of facilities in the CLSCN while coordinating forward and reverse flows. Using algorithms results, a complete validation process is undertaken by GAMS 24.1.2 and MATLAB 2015b software. Finally, the performance of the solution methods is analyzed and compared in different sizes.
Keywords: Closedloop supply chain network; Benders decomposition; Lagrangian relaxation; Simulated Annealing.
 Introduction
Firms are steadily trying to improve the design of their supply chains among soared improvement and intense competition of contemporary global markets. Conventional supply chain practices are being constantly reevaluated as firms find ways of recovering the efficiency and accuracy of their operations. One of the traditional practices that are changing in a supply chain is the flow of products, which is no longer only in the forward direction. In order to make use of returned products, firms are adopting reverse logistics in addition to forward logistics, forming what is known as a CLSCN. One of the reasons for this is the growing interest in product recovery and material recycling, leads to extending the scope of traditional supply chains toward designing, collection, remanufacturing, and recycling centers. This increasing interest is also reflected in the number of papers published in recent years. One reverse logistics network includes inspection centers and remanufacturing facilities (Alshamsi and Diabat, 2015). The design of a reverse logistics network must take into account the coordination of the forward and reverse supply chain in terms of the type of payments, contractual incentives, and inventory risk sharing. (Govindan, Popiuc, and Diabat, 2013). Another approach is, decisions at different levels of the supply chain are now synchronous rather than consecutive. Regarding decisionmaking levels, supply chain management can be classified into strategic, tactical and operational level decisions, depending on the frequency of making them or the timeframe of impact. For instance, decisions on the strategic level have a longterm impression on the company. Tactical level decisions happen at the frequency of several times per year and comprise the selection of transportation ways as well as inventory policies. Eventually, operational level decisions have the shortestterm impact on the company; they may occur on a daily basis and involve to scheduling or routing plans. While evidently independent, these decisions are in fact closely related when it comes to selecting the optimal operations for a company. Consequently, many researchers have started to address this issue by proposing integrated models. Keyvanshokooh, M. Ryan and Kabir (2015) address a novel profit maximization model for CLSCN design as an MILP in which an accelerated stochastic Benders decomposition algorithm is proposed for solving this model. Zhang, Jiang, and Pan (2012) investigated the capacitated lot sizing problem in CLSC considering setup costs, product returns, and remanufacturing. Next, they proposed an LR approach which concentrates on minimizing the total cost of the cyclic logistic network model by determining the best locations of plants, distribution centers, and dismantlers. After that, due to the NPhard nature of the problem, they proposed a revised spanning tree and determinant encoding representation which is a hybrid of two algorithms; simulated annealing and modified genetic algorithm (Yadegari and Zandieh, 2015). In light of these facts, this paper contributes to the CLSCN design literature by presenting several solution methods and comparing them with each other.
Specifically, MILP is suggested for a singleperiod, multi–product, and capacitated CLSCN. The strategic decisions including locations and capacities of facilities as well as the tactical decisions including production amounts and shipments among the network entities are determined to minimize the overall cost of CLSCN that is consisting of investment for reverse logistics facility allocation and transportation costs. The major contributions can be summarized as follows:
First, we propose a novel mathematical modeling in CLSCN including manufacturing, collection, burning, recycling centers. Second, our model determines the location of all new reverse logistics facilities. Finally, we present three solution methods in order to compare the performance and the convergence speed of these methods in different problem sizes. The remainder of this paper is organized as follows. In the next section, we briefly review the literature on the CLSCN design problem and the relevant solution methods. The problem and its formulation are defined in Section 3. At this point, the BD algorithm with an acceleration technique for improving its convergence, the LR algorithm and finally the SA algorithm is presented in Section 4, respectively. Section 5 reports experimental design, parameter tuning, and data generation. Then, Section 6 describes numerical results and sensitivity analyses that allow us to derive managerial insights about this CLSCN. Eventually, section 7deduces this paper and offers some suggestions for future research.
 Literature review
In this section, we inquire into the literature and categorize the studies into two groups. The first category addresses the studies considering CLSCN design, and the second is the studies related to solution algorithms.
 CLSCN design problem
To avoid suboptimality from designing forward and reverse networks separately, many researchers have integrated them into more complex CLSCN (Melo, Nickel and Saldanha 2009). Many CLSCN models are inspired by facility location theory. In this regard, Melo, Nickel, and SaldanhadaGama (2009) presented comprehensive reviews on the facility location models in supply chain planning and on supply chain network design under uncertainty sequentially. Moreover, Pokharel and Mutha (2009) summarized the current developments of reverse supply chains, while Brandenburg, Govindan, Sarkis, and Seuring (2014) and Dekker, Bloemhof, and Mallidis (2012) reviewed quantitative models that address environmental and social aspects in the supply chain.
Ramezani, Kimiagari, and Karimi (2014) suggested a financial approach to model a CLSC design in which financial features are clearly considered as exogenous variables. Hassanzadeh Amin and Zhang (2013) designed a multiple product CLSC network including multiple plants, collection centers, and demand markets using weighted sums, and εconstraint methods.
At the start, Fleischmann, Beullens, BloemhofRuwaardand Wassenhove (2001) considered the integration of forward and reverse flows as a CLSCN using some case studies. They found that this integrated approach could provide a potential for significant cost savings compared to a separate approach. The research that followed was especially carried out with simple facility location models (e.g. Aras, Aksen and Gönül Tanugur, 2008). Next, a more complex model was proposed mainly by considering the reallife case study from the agrifood industry. (e.g. CruzRivera and Ertel, 2009). The field has experienced a strong development over the last decade (e.g. Cardoso, BarbosaPóvoa and Relvas, 2013; Devika, Jafarian and Nourbakhsh, 2014; Gao and Ryan, 2014; Keyvanshokooh, Fattahi, SeyedHosseini and TavakkoliMoghaddam, 2013; Soleimaniand Govindan, 2014). Following this, Kadambalaet al. (2016) proposed a multiobjective MILP which quantitatively measures the effective responsiveness of the CLSC model in terms of time and energy efficiency. Zohal and Soleimani (2016) also designed a forward/reverse logistics network seek to address how a multiobjective logistics model in the gold industry can be created and solved through an efficient metaheuristic algorithm. Designing and optimizing a tire remanufacturing CLSC network based on tire recovery options which are included multiple products, suppliers, plants, retailers, demand markets, and dropoff depots is presented. Furthermore, the application of the model is discussed based on a realistic network in Toronto, Canada using map (Hassanzadeh, Zhang, and Akhtar, 2017).
 Solution algorithms
Designing and planning a CLSC is a complex problem(NPhard) (Krarup and Pruzan, 1983; Schrijver, 2003); therefore, utilizing heuristics and metaheuristics in solving realworld instances of the mentioned problem is unavoidable (Pétrowski and Taillard, 2006). Hence, many solution algorithms including metaheuristic, heuristic, and exact methods have been developed. Most solution methods employ standard commercial packages such as CPLEX to solve MIP formulations. However, when the number of discrete variables is large, the resulting models can be solved only by using metaheuristic or heuristic methods to obtain a near optimal solution. But, because CLSCN design involves large investment and greatly influences the operational and tactical costs as well as the efficiency of service, developing efficient exact algorithms for solving larger and more realistic cases is worthwhile (deSá, de Camargo and deMiranda 2013). Amid these exact solution approaches, branch–and–bound has been a popular methodology combined with other heuristics or LR. There are few papers that develop an exact solution scheme. As a discrete facility location problem, CLSCN design is an attractive candidate for decomposition. It involves both binary variables related to the strategic configuration and continuous variables associated with tactical and operational decisions. Among the decomposition techniques, the BD method (Benders, 1962) is a classical exact algorithm suitable for solving largescale MILP problems having a special structure in the constraint set; i.e., upon fixing the values of the complicating integer variables, the MILP problem reduces to an easy linear program.
However, classical BD and its stochastic version, called the Lshaped method, might not be efficient (Saharidis and Ierapetritou, 2010). The major issues resulting in BD slow convergence are (1) solving the relaxed master problem (RMP) which is, in fact,, an integer program or MILP, and (2) existing the quality of Benders cuts. To overcome the concerns, different acceleration techniques have been proposed to speed up BD. The Lshaped method introduced by Van Slyke and Wets (1969), they defined a cut and achieved significant improvement in convergence by applying such cuts to a problem with degenerate subproblems. Keyvanshokooh, M. Ryan and Kabir (2015) introduced an accelerated stochastic Benders decomposition algorithm for solving a MILP is developed for a multiperiod, singleproduct and capacitated CLSCN. Jeihoonian and Kazemi Zanjani (2015) also presented a BD based solution algorithm together with several algorithmic enhancements for a CLSC in the context of durable products with generic modular structures.
Other solution methods are LR and SA that they are not considered enough in CLSCN design, over the last decade. Concerning, Ferrer and Swaminathan (2006, 2010) used LR schemes in a firm that makes new products in the first period and uses returned cores to make remanufactured products (along with new products) in future periods. Zhang, Jiang, and Pan (2012) investigated A LR based approach for the capacitated lot sizing
problem in CLSC.
The proposed adaptive search approach overcomes such a drawback in solving largescale problems. Lu and Bostel (2007) also provided an MIP formulation to address a CLSC where the customers are directly served from hybrid manufacturing/remanufacturing facilities. An algorithm based on Lagrangian heuristics is developed and the model is tested on data adapted from classical test problems.
SA has the capability of escaping from local optima by accepting nonimproving solutions by certain probability at each level to reach global optimum. Pishvaee, Kianfar, and Karimi (2010) proposed an MILP model, which minimized the transportation and fixed opening costs in a multistage reverse logistics network. They applied an SA algorithm with a special neighborhood search mechanism to find nearoptimal solutions by applying spanning tree and Prufer numbers. As far as the majority of logistics network design problems are categorized as NPhard class, many powerful heuristics and metaheuristic algorithms have been developed for solving these models; nevertheless, developing efficient and effective solution methods is still a crucial need in this area. To this end, how to develop an efficient and effective solution procedure is the aim of this study to overcome the shortages of the existing models and to cope with the features of CLSC. To distinguish other research studies from this research, Table 1 reviews briefly the main studies that make decisions about CLSCN design and different solution methods.
Torabi, Namdar, Hatefi and Jolai (2016) developed a reliable CLSCN design model, which accounts for both partial and complete facility disruptions as well as the uncertainty in the critical input data. Finally, Esmaeili, Allameh, and Tajvidi (2016) proposed the short and longterm behavior of agents in implementing the appropriate collecting strategy in a twoechelon CLSC using Stackelberg game and evolutionary game theory.
Our research paper focuses on CLSCN design and allocation. It also suggests three solution methods. In figure 1 we presented the position of our paper regarding the types of solution approaches in the literature. To differentiate our efforts from those already published in this area, the main contributions of this paper can be summarized as follows:
(1) Obtaining the optimal locations for each of the centers.
(2) Developing MILP models for designing a CLSCN.
(3) Proposing a solution method based on the SA algorithm when the optimal solution cannot be obtained using GAMS software and comparing solution time and quality with LR algorithm.
(4) Representing a BD algorithm to gain the exact solution and proposing the Lshape approach in order to plunge solution time in largesize instances.
(5) Comparing the solution methods with each other in different instances.
Figure 1: Graphical position of our paper among existing methods including Benders decomposition, Lagrangian relaxation, MetaHeuristics, GAMS, CPLEX, LINGO and other methods.
 Problem definition and formulation
 Problem definition
The basis of the proposed model is based on the model proposed by Devika, Jafarian, and Nourbakhsh (2014). They proposed an MIP model which is developed for a triobjectives including minimize the total costs of the network, maximize environmental impacts of the network and the social objective such as the created fixed job
opportunities.
A schematic view of the network depicted in Figure 2, we consider a multiproduct, singleperiod, and capacitated CLSCN consisting of manufacturing, collection, recycling and burning centers as well as customers. In the forward path, manufacturing centers produce new products types. Then, these products are delivered to the customers. In the backward direction, the returned products are transferred to collection centers by customers. At this center, the end of life products is controlled and inspected and if the product can be repaired or restored, they are sent back to the manufacturing facilities, otherwise, they are sent to the recycling centers. In the recycling centers, the products are disassembled to be reused by manufacturers and the remaining parts are sent to the burning centers to be burned.
The model assumptions can be explained as follows:
 The model is considered in multiproduct and singleperiod settings.
 Customer locations are fixed.
 The potential locations of all centers in the network are given.
 All cost parameters are predetermined.
 Customers’ demands are known and must be fully fulfilled.
 The capacity of each center is limited for each product type.
 The number of return products is known.
 Problem formulation
The objective of the MILP model is to minimize the network overall costs which is consisting of transportation, centers establishing and process costs. The mathematical formulation is presented and decision variables definition is listed in Table 2. The notations used in the presented model are as follows.
Sets  Meaning 
I  Set of potential locations of manufacturing centers, i∈I 
J  Set of customer zones, j ∈J 
K  Set of potential locations of collection centers, k ∈K 
L  Set of potential locations of recycling centers, l ∈L 
N  Set of potential locations of recycling centers, n ∈N 
P  Set of product types, p ∈P 
Parameters  
fpi  Fixed cost for opening manufacturing center i 
fck  Fixed cost for opening collection center k 
frl  Fixed cost for opening recycling center l 
fbn  Fixed cost for opening burning center n 
pip  Unit production cost of product p in manufacturing center i 
bckp  Unit separation and collection cost of product p in collection center k 
brlp  Unit recycling cost of product p in recycling center l 
binp  Unit burning cost of product p in burning center n 
cpijp  Transportation cost per unit of product p transported from manufacturing center ito customer j 
ccpjk  Transportation cost per unit of product p transported from customer j to collection center k 
ciknp  Transportation cost per unit of scrapped product p transported from collection center k to burning center n 
crklp  Transportation cost per unit of recoverable product p transported from collection center k to recycling center l 
cwlip  Transportation cost per unit of product p transported from recycling center l to manufacturing center i 
djp  The amount of demand of customer j for product p 
ωjp  The amount of disposal product p of customer j (whether scrapped or recoverable products) 
πip  Maximum available capacity of manufacturing center ifor product p 
τkp  Maximum available capacity of collection center k for product p 
δlp  Maximum available capacity of recycling center l for product p 
μnp  Maximum available capacity of burning center n for product p 
Table 2. Decision variables definition.
The CLSCN design problem can be formulated as follows:
(1)  min(z)=∑iypifpi+∑kyckfck+∑nybnfbn+∑lyrlfrl+∑i∑j∑p(pip+cpijp)uijp+∑j∑k∑p(bckp+ccpjk)sjkp+∑k∑l∑p(brlp+crklp)vklp+∑k∑n∑p(binp+ciknp)zpknp+∑l∑i∑pcwlipwlip  
S.t  
(2)  ∑iuijp≥djp ∀j,p  
(3)  ∑ksjkp≥ωjp ∀j,p  
(4)  ∑jsjkp=∑lvklp+∑nzpknp ∀k,p  
(5)  ∑iuijp≥∑ksjkp ∀j,p  
(6) 


(7)  ∑juijp≤πipypi ∀i,p  
(8)  ∑jsjkp≤τkpyck ∀k,p  
(9)  ∑kvklp≤δlpyrl ∀l,p  
(10)  ∑kzpknp≤μnpybn ∀n,p  
(11)  ∑lwlip≤πipypip ∀i,p  
(12)  ypi,yck,ybn,yrl∈0,1 ∀i,k,n,l  
(13)  uijp,sjkp,vklp,zpknp,wlip≥0 ∀i,j,k,n,p 
Objective (1) aims to minimize the total cost of the CLSCN, where the first four terms refer to the fixed costs of establishing centers. The rest of the terms in (1) represent the transportation costs in addition to process costs, at the network and the centers, respectively. Constraint (2) states that all the demands of customers are satisfied by manufacturers. Constraint (3) specifies that the total quantity shipped from a customer to collection must be more than the amount of disposal product of a customer. Constraint (4) ensures that the product quantity shipped out of collection centers must be equal to the amount of product which is carried to recycling centers in addition to the amount of product which is sent to burning centers. Constraint (5) assures that the quantity trucked of a customer to collection center cannot exceed the quantity of product received for the same customer. Similarly, (6) states that the quantity transported of a recycling center cannot exceed the quantity of product received. Constraints (7) to (11) express capacity restrictions for manufacturing, collection, recycling, burning centers, respectively. Constraint (12) enforces the binary variables. Finally, Constraint (13) ensures nonnegativity, integrality, and bounds on the model variables.
 Solution methodology
The small, mid, and largesized optimization problems require some optimization techniques to search the problem space profoundly and to optimize the problem effectively and efficiently. In this regard, in order to solve the mathematical model, three solution methods including SA, LR, BD algorithms are proposed.
 BD algorithm and Lshaped accelerating
BD algorithm is a classical partitioning method for MIP problems that relies on a problem manipulation using projection, followed by solution strategies of dualization, outer linearization, and relaxation. The MILP is decomposed into a master problem (MP) that involves the firststage decision variables, and Benders subproblems (BSP) to optimize the secondstage decision variables. The BSPs are scenariospecific and connected by the firststage variables. The BSP of this CLSCN can be formulated by fixing the firststage variables to the given values {
ypi=yp̅i.it, yck=yc̅k.it, yrl=yr̅l.it, ybn=yb̅n.it} at iteration it. Because our formulation possesses complete recourse, the BSP is feasible for the given values of first–stage variables, and an optimality cut (OC) may be deduced from an optimal solution to the dual of the subproblem (DSP). Let vector
w1…w10be the set of dual decision variables associated with constraints (2) to (4) and (7) to (10), (5), (6), and (11) in turn. After that the DSP which obtains an upper bound for the objective function of the original CLSN design problem at each iteration it is formulated as follows:
(14)  maxzDSP=∑j∑pw1jpdjp+∑j∑pw2jpωjp∑i∑pw4ipπipyp̅i.it∑k∑pw5kpτkpyc̅k.it∑l∑pw6lpδlpyr̅l.it
∑n∑pw7npμnpyb̅n.it∑i∑pw10ipπipyp̅i.it 
S.t  
(15)  w1jpw4ip+w8jp≤pip+cpijp ∀i,j,p 
(16)  w2jp+w3kpw5kp≤bckp+ccpjk ∀j,k,p 
(17)  w3kpw6lp+w9lp≤brlp+crklp ∀k,l,p 
(18)  w3kpw7np≤binp+ciknp ∀k,n,p 
(19)  w9lpw10ip≤cwlip ∀l,i,p 
(20)  w1jp, w2jp, w3kp,w4ip,w5kp,w6lp,w7np,w8jp,w9lp,w10ip≥0 
Following this, based on the DSP’s solution, the general MP which produces a lower bound for the objective function of original CLSN design model at each iteration can be written as:
(21)  minzMP 
S.t  
(22)  zMP≥∑iypifpi+∑kyckfck+∑nybnfbn+∑lyrlfrl 
(23)  zMP≥∑iypifpi+∑kyckfck+∑nybnfbn
+∑lyrlfrl+∑j∑p∑itww1jpitdjp +∑j∑p∑itww2jpitωjp ∑i∑p∑itww3ipitπipypi ∑k∑p∑itww5ipitτkpyck ∑l∑p∑itww6lpitδlpyrl ∑n∑p∑itww7npitμnpybn ∑i∑p∑itww10ip(it)πipypi 
(24)  ypi,yck,ybn,yrl∈0.1 ∀i,k,n,l 
In this MP, the constraint (23) represents the optimality cut where
ww1…ww10(associated with dual decision variables
w1…w10respectively) indicates the extreme point of the dual polyhedron at each iteration it that results from solving the DSP. This RMP provides a lower bound to the optimal objective value of the MP. At a given iteration of BD, the RMP is first solved to obtain the values of firststage decisions. Then, these values are used to solve DSP to obtain an extreme point and a new optimality cut (23) is included in the RMP. But this algorithm may require a large number of iterations to converge. If DSP faces unbounded solution, then extreme ray and feasibility cut must be added to MP. The constraints (25) (31) and (32) revealed extreme ray and feasibility cut, respectively. They can be written as:
(25)  ∑j∑pw1jpdjp+∑j∑pw2jpωjp
∑i∑pw4ipπipyp̅i.it∑k∑pw5kpτkpyc̅k.it ∑l∑pw6lpδlpyr̅l.it∑n∑pw7npμnpyb̅n.it ∑i∑pw10ipπipyp̅i.it=1 
(26)  w1jpw4ip+w8jp≤0 ∀i.j.p 
(27)  w2jp+w3kpw5kp≤0 ∀j.k.p 
(28)  w3kpw6lp+w9lp≤0 ∀k.l.p 
(29)  w3kpw7np≤0 ∀k.n.p 
(30)  w9lpw10ip≤0 ∀l.i.p 
(31)  w1jp. w2jp. w3kp.w4ip.w5kp.w6lp.w7np.w8jp.w9lp.w10ip≥0 
And,
(32)  ∑j∑p∑itww1jpitdjp
+∑j∑p∑itww2jpitωjp ∑i∑p∑itww4jpitπipypi ∑k∑p∑itww5jpitτkpyck ∑l∑p∑itww6jpitδlpyrl ∑n∑p∑itww7jpitμnpybn ∑i∑p∑itww10jpitπipypi≤0

Figure 3 depicts the flowchart for BD with the L–shaped algorithm.
 LR approach
Our LR based approach is advantageous in that it naturally gives a lower bound on the
optimal objective function value, which allows us to assess the quality of solutions found and to compare with the value obtained from SA method in optimality gap and computational time.
After applying LR to model (1) and relaxing the constraints (7), (9), (10), (11), and (5) with Lagrange multipliers
λ1, λ3, λ4, λ5,and
λ7,respectively, the mathematical model becomes:
(33)  minLR=∑iypifpi+∑kyckfck
+∑nybnfbn+∑lyrlfrl+∑i∑j∑ppip+cpijpuijp+∑j∑k∑pbckp+ccpjksjkp +∑k∑l∑pbrlp+crklpvklp +∑k∑n∑pbinp+ciknpzpknp+∑l∑i∑pcwlipwlip +(λ1(∑i∑j∑puijp∑i∑pπipypi))+(λ3(∑k∑l∑pvklp∑l∑pδlpyrl))+(λ4(∑k∑n∑pzpknp∑n∑pμnpybn))+(λ5(∑l∑i∑pwlip∑i∑pπipypip))+(λ7(∑j∑k∑psjkp∑i∑j∑puijp)) 
S.t  
(34)  ∑iuijp≥djp ∀j,p 
(35)  ∑ksjkp≥ωjp ∀j,p 
(36)  ∑jsjkp≥∑lvklp+∑nzknp ∀k,p 
(37)  ∑jsjkp≤τkpyck ∀k,p 
(38)  ∑kvklp≥∑iwlip ∀l,p 
(39)  ypi,yck,ybn,yrl∈0,1 ∀i,k,n,l 
(40)  uijp,sjkp,vklp,zknp,wlip≥0 ∀i,j,k,n,p 
(41)  λ1, λ3, λ4, λ5, λ7≥0 
Following approach helps to obtain appropriate Lagrange multipliers for the model. It is comprised of 6 steps which can be written as:
Step1. Initialization of Lagrange multipliers vector ( =0) and the selection of stepsize (k).
Step2. Solving problem which is relaxed, and calculating the optimal values of decision variables.
Step3. If constraint i which is relaxed for the optimal values of decision variables is not satisfied then,
λi=λi+k(It means which objective function is becoming worse).
Step 4. If constraint i which is relaxed for the optimal values of decision variables is satisfied, then
λi=λik.
Step 5. If after elapsed consecutive iteration m, the best lower bound value is not improved, then
k=k2.
Step 6. Return to step 2.
 SA algorithm
SA is one of the most popular iterative algorithms, which has been applied broadly to solve many combinatorial optimization problems (Jayaraman and Ross 2003; Pishvaee et al. ,2010).
According to the above references, the SA algorithm is outlined as follows:
With the aid of test problems, the parameters of the SA algorithm are tuned by changing their values and comparing the results. The number of inner loop iterations is increased dynamically with the rate of (P/2*temperature) to gain more intensification and better solutions at lower temperatures. In addition, neighborhood search is used in temperatures more than 0.001 in order to increase the speed of approaching to desirable solutions. Finally, the temperature is cooled using alpha coefficient. The other parameters of the SA algorithm and pseudocode used in this model are summarized in Fig. 4.
Figure 4. SA pseudocode.
b and alpha values are obtained for 9 instances by Taguchi approach. They are shown in Table 3.
 Experimental design
This section comparatively evaluates the performance of the proposed algorithms for the discussed integrated forward/reverse network design model. The algorithm is the SA which is compared in 9 sizes, as shown in Table 3. The algorithm was implemented 25 times in each instance. Relative percentage deviation (RPD) was used as a performance measure to compare different running. RPD was obtained by the given formula as follows:
RPD=AlgsolMinsolMinsol×100  (42) 
Where,
Algsoland
Minsolare the objective values for each replication of a trial in a given instance and the obtained best solution, respectively. After converting the objective values into RPDs, the mean RPD was calculated for each desired trial (Yadegari and Zandieh, 2015).
 Parameter tuning
It is known that different levels of the parameters clearly affect the quality of the solutions obtained by each algorithm. SA algorithm can be gained with different combinations of parameters. We try to recognize factors that are statically significant for SA algorithm in terms of objective value and CPU time. Hence, we apply parameters tuning for initial temperature and cooling rate (alpha). Table 3 illustrates the preferred parameters for 9 instances.
Table 3. Factors and their levels in the proposed algorithm.
 Data generation
To evaluate our proposed algorithms performance for solving the given model, we have generated random datasets. Table 4 includes the number of fixed locations considered in CLSCN such as the number of manufacturers, customers, collection, recycling and burning centers. Table 5 demonstrates capacities, fixed costs, demands (in unit) and disposal product of customer in the proposed model. Table 6 includes the variable costs of transportation between each consecutive level of the network depicted. It also consists of production, process, recycling and burning costs for a unit in each center. The parameters for each size of the test problems are generated randomly using uniform distributions in the range of the test problems existing in the recent literature (Keyvanshokooh, M. Ryan, and Kabir, 2015).
Table 4. Test problems characteristics.
Table 5.Capacity, fixed costs, and demand and disposal product of customer.
Table 6.Unit shipping cost between CLSCN levels and process costs in each center.
 Numerical results
In order to demonstrate the proposed methodologies application and to compare the abovementioned solution methods performances in terms of the objectivefunction values and required CPU time, various test problems of different sizes are provided in this section.
Using the above numerical data, the MILP problem modeled in Section 3.2 is solved for various instances of different sizes using all the methods on a computer with core i7CPU 2.40 GHz, RAM 8.00 GB, using the MATLAB 2015b and GAMS 24.1.2 software. For example, three graphs in Figure 5 demonstrate the values of objective, optimality gap and CPU time for 9 instances obtained employing all the methods.
Figure 5. The objective values, optimality gap and CPU time for 9 instances using all methods.
As can be seen, the objective values of all methods had a marginal increase; however, there were slight differences amid objective values. Regarding the gap optimality and Table 7, the SA gap optimality was more than LR in all instances. Concerning CPU time, the SA CPU times and LR experienced a negligible growth, while it had a considerable increase for BD with Lshaped. Furthermore, the objective, optimality gap and CPU time values for all methods are summarized in Table 7.
Table 7.The objective values, optimality gap and computational times (in seconds) for different algorithms and GAMS.
With regard to Table 7, the low convergence speed of BD algorithm in mid and large sizes leads to use of the Lshaped accelerator. For example, in Figures 6 and 7 the convergence speeds of both are compared for instance 5. The convergence speed of BD with Lshaped to an optimal solution is much more than BD algorithm; at 34.655 seconds and 200.828 seconds respectively. The BD with Lshaped after elapsed 22 iterations has converged to the optimal solution, while the BD algorithm has reached to optimality after 647 iterations. Table 7 also displays the computational times for solving each instance by solution methods and also by directly solving the extensive form by GAMS.
Figure 6. The BD algorithm convergence speed.
Figure 7. The BD with Lshaped convergence speed.
 Conclusion
In this paper, an MILP model for a single–period, multi–product, and capacitated CLSCN design problem is formulated to minimize the overall cost of CLSCN. As the major contribution, we presented different solution methods including BD, LR, and SA for solving the proposed mathematical model. Then, some numerical instances were created to analyze and validate the formulation. After that, an accelerated BD algorithm was proposed so as to be accelerated convergence. Subsequently, the SA algorithm optimality gap is more than LR algorithm. Finally, BD with Lshaped performance is much better than BD about convergence speed, at mid and large sizes. However, the LR algorithm is better than the SA algorithm because it makes it possible to gain the best solution.
As this paper introduces several solution methods for CLSCN design, there are some opportunities for future research such as applying other approaches including DantzigWolfe decomposition, Particle Swarm Optimization, etc. Moreover, to solve this largescale problem, other versions of the BD approach such as a Bendersbased branchandcut approach, where a single branchandcut tree is constructed and then the Benders cuts are added during the exploration of this tree can be applied for performance comparisons.
References
AlSalem, M., Diabat, A., Dalalah, D., andAlrefaei, M. (2016) ‘A closedloop supply chain management problem: Reformulation and piecewise linearization’, Journal of Manufacturing Systems 40: 18.
Alshamsi, A., and Diabat, A. (2015) ‘A reverse logistics network design’, Journal of Manufacturing Systems 37: 589598.
Amin, S. H., Zhang, G., & Akhtar, P. (2017) ‘Effects of uncertainty on a tire closedloop supply chain network’, Expert Systems with Applications, 73: 8291.
Amin, S. H., and Zhang, G. (2013) ‘A multiobjective facility location model for closedloop supply chain network under uncertain demand and return’, Applied Mathematical Modelling 37(6): 41654176.
Aras, N., Aksen, D., and GönülTanugur, A. (2008) ‘Locating collection centers for incentivedependent returns under a pickup policy with capacitated vehicles’, European Journal of Operational Research 191: 1223–1240.
Brandenburg, M., Govindan, K., Sarkis, J., andSeuring, S. (2014) ‘Quantitative models for sustainable supply chain management: Developments and directions’, EuropeanJournal of Operational Research233: 299–312.
Cardoso, S. R., BarbosaPóvoa, A. P. F. D., andRelvas, S. (2013) ‘Design and planning of supply chains with integration of reverse logistics activities under demand uncertainty’, EuropeanJournal of Operational Research226: 436–451.
CruzRivera, R., and Ertel, J. (2009) ‘Reverse logistics network design for the collection of EndofLife Vehicles in Mexico’, EuropeanJournal of Operational Research196: 930–939.
Dekker, R., Bloemhof, J., and Mallidis, I. (2012) ‘Operations Research for green logistics – An overview of aspects, issues, contributions and challenges’, European Journal ofOperational Research 219: 671–679.
deSá, E. M., de Camargo, R. S., andde Miranda, G. (2013) ‘An improved Benders decomposition algorithm for the tree of hubs location problem’, European Journal of Operational Research 226(2): 185202.
Devika, K., Jafarian, A., and Nourbakhsh, V. (2014) ‘Designing a sustainable closedloop supply chain network based on triple bottom line approach: A comparison of metaheuristics hybridization techniques’, EuropeanJournal of Operational Research235: 594–615.
Esmaeili, M., Allameh, G., and Tajvidi, T. (2016) ‘Using game theory for analysing pricing models in closedloop supply chain from shortand longterm perspectives’, International Journal of Production Research 54(7): 21522169.
Fleischmann, M., Beullens, P., BloemhofRuwaard, J. M., and Wassenhove, L. N. (2001) ‘Theimpact of product recovery on logistics network design’, Production and Operations Management 10(2): 156–173.
Gao, N., and Ryan, S. (2014) ‘Robust design of a closedloop supply chain network for uncertain carbon regulations and random product flows’, EURO Journal on Transportation and Logistics3: 5–34.
Govindan, K., Popiuc, M. N., andDiabat, A. (2013) ‘Overview of coordination contracts within forward and reverse supply chains’, Journal of Cleaner Production 47: 319334.
Jayaraman, V., and Ross, A. (2003) ‘A simulated annealing methodology to distribution network design and management’, European Journal of Operational Research 144(3): 629645.
Jeihoonian, M., Zanjani, M. K., and Gendreau, M. (2016) ‘Accelerating Benders decomposition for closedloop supply chain network design: Case of used durable products with different quality levels’, European Journal of Operational Research 251(3): 830845.
Kadambala, D. K., Subramanian, N., Tiwari, M. K., Abdulrahman, M., and Liu, C. (2016) ‘Closed loop supply chain networks: Designs for energy and time value efficiency’, International Journal of Production Economics.
Keyvanshokooh, E., Fattahi, M., SeyedHosseini, S. M., and TavakkoliMoghaddam, R. (2013) ‘A dynamic pricing approach for returned products in integrated forward/reverse logistics network design’, Applied Mathematical Modelling37: 10182–10202.
Keyvanshokooh, E. (2015) ‘Hybrid robust and stochastic optimization for closedloop supply chain network design using accelerated Benders decomposition’, Master.Thesis, Iowa State University, USA.
Lu, Z., and Bostel, N. (2007) ‘A facility location model for logistics systems including reverse flows: The case of remanufacturing activities’, Computers & Operations Research 34(2): 299323.
Melo, M. T., Nickel, S., and SaldanhadaGama, F. (2009) ‘Facility location and supply chain management – A review’, EuropeanJournal of Operational Research 196: 401–412.
Pishvaee, M. S., Kianfar, K., and Karimi, B. (2010) ‘Reverse logistics network design using simulated annealing’, The International Journal of Advanced Manufacturing Technology 47(14): 269281.
Pokharel, S., and Mutha, A. (2009) ‘Perspectives in reverse logistics: A review’, Resources, Conservation and Recycling 53: 175–182.
Ramezani, M., Kimiagari, A. M., and Karimi, B. (2014) ‘Closedloop supply chain network design: A financial approach’, Applied Mathematical Modelling38(15):40994119.
Soleimani, H., and Govindan, K. (2014) ‘Reverse logistics network design and planning utilizing conditional value at risk’, European Journal of Operational Research 237:487–497.
Torabi, S. A., Namdar, J., Hatefi, S. M., and Jolai, F. (2016) ‘An enhanced possibilistic programming approach for reliable closedloop supply chain network design’, International Journal of Production Research 54(5): 13581387.
Zhang, Z. H., Jiang, H., andPan, X. (2012) ‘A Lagrangian relaxation based approach for the capacitated lot sizing problem in closedloop supply chain’, International Journal of Production Economics 140(1): 249255.
Zohal, M., and Soleimani, H. (2016) ‘Developing an Ant Colony Approach for Green ClosedLoop Supply Chain Network Design: A Case Study in Gold Industry’, Journal of Cleaner Production.
Research study 
Demand  No. of SC echelon 
Mathematical model 
Solution approach 
Supply chain entities  
Stochastic 
Deterministic 

(Pishvaee and Kianfar, 2010)   4  MILP  SA and LINGO  Collection, recovery, disposal centers customer zone  
(Zhang, Jiang and Pan, 2012)   2  MILP  LR  Manufacturer and remanufacturing centers 

(Keyvanshokooh, Fattahi, SeyedHosseini and TavakkoliMoghaddam, 2013)   6  MILP  CPLEX  Collection, production, distribution and disposal centers and customer zones and a hybrid processing facility  
(Hassanzadeh Amin and Zhang, 2013)   3  MILP  Econstraint and weighted sums  Multiple plants, collection centers and demand markets  
(Devika, Jafarian and Nourbakhsh, 2014)   5  MILP  Metaheuristic  Collection, disposal, distribution, remanufacturing, recycling and manufacturing centers and customer zone  
(Fattahi, Mahootchi, Husseini, Keyvanshokooh and Alborzi, 2014)   4  MILP  Metaheuristic  Supplier, manufacturer, retailer and customers  
(Gao and Ryan, 2014)   3  MIP  CPLEX  Factory, warehouse and retailer  
(Ramezani, Kimiagari and Karimi, 2014)   4  MILP  Financial  Supplier, plant, distribution, collection, repair, disposal and customer  
(Giovanni and Zaccour, 2014)   3  Analytical  Manufacturer ,remanufacturing and retailer 

(Jeihoonian and KazemiZanjani, 2015)   5  MIP  BD  Disassembly, bulk recycling, collection, disposal, manufacturer, supplier and distribution centers  
(Keyvanshokooh, M. Ryan and Kabir, 2015)   4  MILP  BD  Manufacturer ,remanufacturing, distribution and disposal centers and retailers 

(Torabi, Namdar, Hatefi and Jolai, 2016)   5  MILP  GAMS  productionrecovery, distributioncollection, disposal centers and customer zone  
(AlSalem, Diabat, Dalalah and Alrefaei, 2016)   3  MINLP  GAMS  Warehouse, retailer and customer  
(Zohal and Soleimani, 2016)   5  MILP  LINGO and Metaheuristic  Supplier, manufacturer, distribution, collection, recycling and disposal centers  
(Kadambala, Subramanian, Tiwari, Abdulrahman, Liu, 2016)   3  MILP  Metaheuristic  Manufacturer and distributor centers and customers  
(Hassanzadeh, Zhang and Akhtar, 2017)   5  MILP  A new decision treebased  multiple products, suppliers, plants, retailers, demand markets, and dropoff depots  
Our work   5  MILP  SA, LR, BD and GAMS 
Manufacturing, collection, recycling and burning centers and customers 
Table 8: literature review.
Table 2. Decision variables definition.
Decision variables  
ypi  Binary variable equal to 1 if a manufacturing center i is opened, 0 otherwise 
yck  Binary variable equal to 1 if a collection center k is opened, 0 otherwise 
yrl  Binary variable equal to 1 if a recycling center l is opened, 0 otherwise 
ybn  Binary variable equal to 1 if a burning center n is opened, 0 otherwise 
uijp  Quantity of products p transported from manufacturing center ito customer j 
sjkp  Quantity of returned products p transported from customer jto collection center k 
vklp  Quantity of recoverable products p transported from collection center kto recycling center l 
zpknp  Quantity of scrapped products p transported from collection center kto burning center n 
wlip  Quantity of products p transported from recycling center lto manufacturing center i 
Table 3. Factors and their levels in the proposed algorithm.
Instances  Parameter  Level 
1  b  1000 
alpha  0.8  
2  b  600 
alpha  0.8  
3  b  500 
alpha  0.7  
4  b  500 
alpha  0.6  
5  b  300 
alpha  0.95  
6  b  1000 
alpha  0.7  
7  b  4000 
alpha  0.95  
8  b  750 
alpha  0.95  
9  b  3000 
alpha  0.9 
Table 4. Test problems characteristics.
No. instance  No. products P  Manufacturers I  Customers J  Collection centers K  Recycling centers L  Burning centers N 
Small size  
1  2  3  4  2  2  2 
2  4  5  5  3  3  3 
3  4  6  7  4  4  4 
Midsize  
4  6  8  10  6  6  6 
5  7  10  12  8  8  8 
6  8  12  15  9  9  9 
Large size  
7  9  13  17  11  11  11 
8  10  15  20  12  12  12 
9  12  18  25  15  15  15 
Table 5. Capacity, fixed costs, and demand and disposal product of customer.
Manufacturers, collection, recycling and burning centers  Capacity  ~uniform (1000, 1500) 
Fixed costs  ~uniform (10, 40)  
Customer  Demand  ~uniform (250, 400) 
Customer  ω  ~uniform (150, 200) 
Table 6. Unit shipping cost between CLSCN levels and process costs in each center.
Flow  Distribution 
Manufacturer customers, customer collection center, collection center recycling and burning centers, recycling center manufacturer  ~uniform (1, 10) 
Production, process, recycling and burning costs for a unit  ~uniform (1, 10) 
Table 6. The objective values, optimality gap and computational times (in seconds) for different algorithms and GAMS.
Instance no.  BD  BD with L shaped  LR  GAMS  
Objective  Gap (%)  Time  Objective  Gap (%)  Time  Objective  Gap (%)  Time  objective  
1  27055  0  5.641  27055  0  5.541  26664.45  1.44  5.533  27055 
2  67880  0  33.252  67880  0  15.075  67256.963  0.92  5.631  67880 
3  72575  0  40.696  72575  0  17.196  72070.152  0.7  6.236  72575 
4  147449  0  129.907  147449  0  29.002  143177.443  2.9  6.747  147449 
5  270235  0  200.828  270235  0  34.655  264713.623  2.04  7.61  270235 
6  359808  0  263.042  359808  0  42.485  348468.25  3.15  8.259  359808 
7  395232  0  10983.998  395232  0  62.508  389761  1.38  9.826  395232 
8  514610  0  14023  514610  0  72.954  504943  1.88  9.919  514610 
9  807346  0  25287  807346  0  81.798  762328  5.576  11.591  807346 
Instance no.  SA  
Objective  Gap (%)  Time  
1  27805  2.77  0.12042  
2  69700  2.83  0.124  
3  74214  2.26  0.13785  
4  155911  5.74  0.163  
5  280379  3.75  1.05477  
6  377383  4.88  1.47  
7  411675  4.16  2.071907  
8  574521  6.4  2.98  
9  862215  6.796  3.528  
Figure 1: Graphical position of our paper among existing methods including Benders decomposition, Lagrangian relaxation, MetaHeuristics, GAMS, CPLEX, LINGO and other methods.
Figure 2. The proposed CLSCN.
Figure 3. BD with an Lshaped flowchart.
Figure 4. SA pseudocode.
Figure 5. The objective values, optimality gap and CPU time for 9 instances using all methods.
Figure 6. The BD algorithm convergence speed.
Figure 7. The BD with Lshaped convergence speed.