Specialist Dissertation Services
Dissertation Examples

Designing a closed-loop supply chain network: A comparison of algorithms

Abstract

This paper attempts to design a closed-loop 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 L-shaped 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: Closed-loop supply chain network; Benders decomposition; Lagrangian relaxation; Simulated Annealing.

1. 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 re-evaluated 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 decision-making levels, supply chain management can be classified into strategic, tactical and operational level decisions, depending on the frequency of making them or the time-frame of impact. For instance, decisions on the strategic level have a long-term 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 shortest-term 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 NP-hard 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 single-period, 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.

2. 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.

1. CLSCN design problem

To avoid sub-optimality 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 Saldanha-da-Gama (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, Bloemhof-Ruwaardand 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 real-life case study from the agri-food industry. (e.g. Cruz-Rivera 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, Seyed-Hosseini and TavakkoliMoghaddam, 2013; Soleimaniand Govindan, 2014). Following this, Kadambalaet al. (2016) proposed a multi-objective 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 multi-objective logistics model in the gold industry can be created and solved through an efficient meta-heuristic 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 drop-off 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).

1.  Solution algorithms

Designing and planning a CLSC is a complex problem(NP-hard) (Krarup and Pruzan, 1983; Schrijver, 2003); therefore, utilizing heuristics and metaheuristics in solving real-world 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 (de-Sá, de Camargo and de-Miranda 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 large-scale 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 L-shaped 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 L-shaped 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 sub-problems. Keyvanshokooh, M. Ryan and Kabir (2015) introduced an accelerated stochastic Benders decomposition algorithm for solving a MILP is developed for a multi-period, single-product 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 large-scale 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 non-improving 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 near-optimal solutions by applying spanning tree and Prufer numbers. As far as the majority of logistics network design problems are categorized as NP-hard 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.

Table 1: literature review.

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 long-term behavior of agents in implementing the appropriate collecting strategy in a two-echelon 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 L-shape approach in order to plunge solution time in large-size instances.

(5)  Comparing the solution methods with each other in different instances.

3. Problem definition and formulation
   1.  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 tri-objectives 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 multi-product, single-period, 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 multi-product and single-period 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.

  Figure 2. The proposed CLSCN.

2.  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.

Table 2. Decision variables definition.

The CLSCN design problem can be formulated as follows:

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 non-negativity, integrality, and bounds on the model variables.

4. Solution methodology

The small, mid, and large-sized 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.

1.  BD algorithm and L-shaped 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 first-stage decision variables, and Benders sub-problems (BSP) to optimize the second-stage decision variables. The BSPs are scenario-specific and connected by the first-stage variables. The BSP of this CLSCN can be formulated by fixing the first-stage 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 sub-problem (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:

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:

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 first-stage 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:

And,

Figure 3 depicts the flowchart for BD with the L–shaped algorithm.

2. 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:

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 step-size (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=λi-k.

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.

3.  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 pseudo-code used in this model are summarized in Fig. 4.

Figure 4. SA pseudo-code.

b and alpha values are obtained for 9 instances by Taguchi approach. They are shown in Table 3.

5. 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:

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).

1. 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.

2. 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.

6. Numerical results

In order to demonstrate the proposed methodologies application and to compare the above-mentioned solution methods performances in terms of the objective-function 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 i7-CPU 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 L-shaped. 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 L-shaped accelerator. For example, in Figures 6 and 7 the convergence speeds of both are compared for instance 5. The convergence speed of BD with L-shaped to an optimal solution is much more than BD algorithm; at 34.655 seconds and 200.828 seconds respectively. The BD with L-shaped 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 L-shaped convergence speed.

7. 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 L-shaped 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 Dantzig-Wolfe decomposition, Particle Swarm Optimization, etc. Moreover, to solve this large-scale problem, other versions of the BD approach such as a Benders-based branch-and-cut approach, where a single branch-and-cut tree is constructed and then the Benders cuts are added during the exploration of this tree can be applied for performance comparisons.

References

Al-Salem, M., Diabat, A., Dalalah, D., andAlrefaei, M. (2016) ‘A closed-loop supply chain management problem: Reformulation and piecewise linearization’, Journal of Manufacturing Systems 40: 1-8.

Alshamsi, A., and Diabat, A. (2015) ‘A reverse logistics network design’, Journal of Manufacturing Systems 37: 589-598.

Amin, S. H., Zhang, G., & Akhtar, P. (2017) ‘Effects of uncertainty on a tire closed-loop supply chain network’, Expert Systems with Applications, 73: 82-91.

Amin, S. H., and Zhang, G. (2013) ‘A multi-objective facility location model for closed-loop supply chain network under uncertain demand and return’, Applied Mathematical Modelling 37(6): 4165-4176.

Aras, N., Aksen, D., and GönülTanugur, A. (2008) ‘Locating collection centers for incentive-dependent returns under a pick-up 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., Barbosa-Pó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.

Cruz-Rivera, R., and Ertel, J. (2009) ‘Reverse logistics network design for the collection of End-of-Life 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): 185-202.

Devika, K., Jafarian, A., and Nourbakhsh, V. (2014) ‘Designing a sustainable closed-loop 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 closed-loop supply chain from short-and long-term perspectives’, International Journal of Production Research 54(7): 2152-2169.

Fleischmann, M., Beullens, P., Bloemhof-Ruwaard, 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 closed-loop 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: 319-334.

Jayaraman, V., and Ross, A. (2003) ‘A simulated annealing methodology to distribution network design and management’, European Journal of Operational Research 144(3): 629-645.

Jeihoonian, M., Zanjani, M. K., and Gendreau, M. (2016) ‘Accelerating Benders decomposition for closed-loop supply chain network design: Case of used durable products with different quality levels’, European Journal of Operational Research 251(3): 830-845.

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., Seyed-Hosseini, S. M., and Tavakkoli-Moghaddam, 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 closed-loop 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): 299-323.

Melo, M. T., Nickel, S., and Saldanha-da-Gama, 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(1-4): 269-281.

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) ‘Closed-loop supply chain network design: A financial approach’, Applied Mathematical Modelling38(15):4099-4119.

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 closed-loop supply chain network design’, International Journal of Production Research 54(5): 1358-1387.

Zhang, Z. H., Jiang, H., andPan, X. (2012) ‘A Lagrangian relaxation based approach for the capacitated lot sizing problem in closed-loop supply chain’, International Journal of Production Economics 140(1): 249-255.

Zohal, M., and Soleimani, H. (2016) ‘Developing an Ant Colony Approach for Green Closed-Loop Supply Chain Network Design: A Case Study in Gold Industry’, Journal of Cleaner Production.

Table 2. Decision variables definition.

Table 3. Factors and their levels in the proposed algorithm.

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.

Table 6. The objective values, optimality gap and computational times (in seconds) for different algorithms and GAMS.

Figure 2. The proposed CLSCN.

Figure 3. BD with an L-shaped flowchart.

Figure 4. SA pseudo-code.

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 L-shaped convergence speed.