Abstract
The design for integrating the forward-reverse logistics network has attracted growing attention with the relentless pressures from environmental and social requirements to increase operational efficiency and support sustainable operations. Therefore, a dynamic facility in the closed chain is presented for the first time. In this problem, a set of facilities with determined locations are used based on whether they serve as forward-reverse, or joint facilities dealing with products from both directions. The objective is to minimize the total costs, including fixed location, fixed ordering, inventory holding, transportation, and reprocessing costs. In this paper, a mixed integer nonlinear program (MINLP) of a single type of forwarding product and a single type of returned product is proposed and a Benders decomposition method (BDM) is developed to solve the problem efficiently. The computational results on several randomly generated issues demonstrate that the proposed algorithm can efficiently solve the MINLP model compared to other algorithms.
1. Introduction
The capacitated facility location problem is a well-known discrete optimization problem, in which a manufacturer decides which facilities to open from a given set and which open facility to assign to each customer to maximize profit or capture. Logistics network design problems that consider the facility locations and product flows have been studied for decades. Due to the increasing pressures from environmental and social requirements, many manufacturers have adopted the practice of using second-hand products and incorporated used product recovery activities into the production process. Consequently, the design of reverse logistics networks has gained concern, which includes not only economic parts but also how logistics networks will affect other aspects of human life, such as the environment and sustainability of natural resources. Implementation of the reverse logistics operations requires setting up additional appropriate logistics infrastructure for the arising flows of used and recovered products. Physical location, facilities, and transportation links need to be chosen to transfer a new product from manufacturers to customers and to convey used products from users to manufacturers for recovery or safe disposal. Furthermore, the interaction between the distribution of new products and the collection of used products also needs to be considered due to the influence of the reverse logistics on forwarding logistics, such as the occupancy of storage spaces and transportation capacity.
The first model in the dynamic facility location problem was introduced in [1] in which locations of a single facility were determined over multiple periods. A comprehensive modeling framework which encompasses multicommodity and multiechelon facility location problem was presented in [2], which explicitly considered many practical aspects. Two-echelon dynamic facility location problem and its Lagrangian relaxation-based algorithm are designed in [2]. A dynamic location problem with opening, closure, and reopening of facilities was formulated in [3], and a primal-dual heuristic was developed to solve the problem. Their study has an important characteristic that distinguishes it from previous works by considering the possibility of reconfiguring one location more than once during the planning horizon. There are various applications of dynamic facility location in the literature, and many papers have been introduced in this area, for example, Fleischmann et al. [4] explained the applications in both the public and private sectors, Peeters and Antunes [5] worked on the location of production facilities, Correia and Captivo and Kim and Kim [6,7] introduced the models for healthcare facilities, and so on [8, 9].
Some of these papers have various approaches, like the temporary closing of facilities verified in [3, 10]. Also, Troncoso and Garrido [11] considered the capacity expansion and reduction and a reverse logistics network, which takes into account the optimal selection of facilities, including the capacities of inspection centers and remanufacturing facilities addressed by Alshamsi and Diabat [12]. Redesigning of the supply chain network by the focus on reviewing the coordination by contracts of the forward and reverse supply chain is addressed in [13]. Their model includes the criteria of transfer payment contractual incentives and inventory risk sharing. Recently, Salem et al. [14] considered a problem describing both the forward and reverse flow of a single product type. In this problem, a set of facilities is verified in the model whose location is determined based on whether they serve as forward-reverse or joint facilities dealing with products from both directions, inventory holding, location, reprocessing, and some other aspects. In related solution methods, because metaheuristics have a rich literature [15–28], some algorithms such as tabu search, simulated annealing, and genetic algorithms have been frequently applied to several families of location problems [29]. Also, Araya-Sassi et al. [30] used the Lagrangian relaxation heuristics-based method for single-period facility location problem and Shulman [31] proposed this algorithm for multiperiod facility location and multiechelon issues in the context of the supply chain. Salem et al. [14] successfully applied heuristic-based methods to dynamic facility location problems, and Matteo et al. [32] showed how to use the Benders method to obtain a reliable solution.
Based on [33], the Benders decomposition is preferable to metaheuristics when the considered problem can be solved within acceptable CPU time. That is why this paper uses the Benders algorithm. In this paper, we consider integrating different decision levels, inventory holding, and facility location in the context of the dynamic closed-loop supply chain. A mixed integer nonlinear program (MINLP) is proposed, and then linearization techniques are applied for getting rid of nonlinearity terms. Besides, Benders decomposition method- (BDM-) based heuristics are used which find good quality scale to solve the problem in reasonable computing times.
The remainder of the paper is organized as follows. Section 2 extends the MINLP formulation for the considered problem and linearization techniques. Section 3 explains how to use the Benders decomposition and outlines the resulting heuristics. After that, the numerical experiment results are reported in Section 4. Finally, the last section concludes the paper and provides a discussion of problem extensions and future research directions.
2. Model Description
A system of producing and remanufacturing which consisted of customers, facilities (intermediate centers), and plants (producers) was studied, as depicted in Figure 1.
In this system, a new product is delivered to several geographically dispersed customers. Also, a used product is taken back from the customers to reused facilities or safely disposed. Therefore, there are three kinds of flow. The backward flow from customers to plants is formed by the used product. The forward flow from plants to customers consists of the new product and the third flow among the facilities [34]. In this paper, a new type of intermediate depots, namely, hybrid processing facility, is also taken into account. New or/and used flows can be transferred via hybrid processing facilities. Advantages of building such hybrid processing facilities are cost-saving and pollution reduction due to sharing material handling equipment and infrastructure. All the used products are first shipped back to recovering facilities or hybrid facilities, in which some of them will be recovered and turned back to the market as a new product. Other checked return products will be then sent around to plants, in which some of them might be disposed. The customers’ demand could be met by both forward and recovered products, which meant that the renewed product from facilities is considered the same as the new product in terms of satisfying customer demand. Facilities also hold inventory, which implies that the need faced by the facilities is the overall demand among the customers in all periods. As working inventory at the facilities, we define the amount of product ordered by the facilities. This will be done for used products regarding plants, i.e., the used product may be held in the facilities and may be transported to plants in the future. These assumptions are logical because of variable pricing, including transportation costs, fixed holding costs for facilities, and other related fees in different periods. It is assumed that there are two kinds of request points: demand point for new products (by demand) and demand point for used products (by demand). A customer with both requirements can be considered as two demand points in the same places [35]. This customer is free to visit one of the following locations:(1)Visit facility type “f” for new products and type “r” for used products.(2)Visit facility type “f” for new products and type “h” for used products.(3)Visit facility type “h” for new products and type “r” for used products.(4)Visit facility type “h” for new products and used products.
Based on the notation in the appendix, the formulation for the closed-loop capacitated dynamic joint location inventory is proposed as follows:
In the proposed model, the objective function (1) is to minimize the total investment and operational costs in the dynamic logistics network, including the shipping cost of the forward products from plants to facilities and facilities to customers, the shipping cost of the returned products from customers to facilities and facilities to plants, the holding costs for forwarding products, the fixed cost associated with building forward processing facilities, collection facilities, and hybrid processing facilities, the cost associated with the relocation (closing or opening) of forward, collection and hybrid processing facilities, the costs for buying returned products from customers, the costs for shipping products among facilities, and the costs for refreshing usable the backward products, and the acquires interested from customers to sell forward product, and plants to sell backward product. Constraint (2) guarantees that each customer can select exactly one facility, and constraints (3) and (4) express that each customer can choose a facility type “f” or “h” for a new product and a facility type “r” or “h” for the used product from the leader and follower firms, respectively. Constraint (5) guarantees that each potential facility can be used for only one type of facility, and constraints (6) and (7) are about the inventory holdings for new and used products, respectively. Constraints (8) and (9) bound the number of new products shipped via facilities to customers and used products shipped via customers to facilities. Constraint (10) is about the amount of conveying new products from plants and facilities to facilities, and constraint (11) expresses the limit on the shipped new products from plant l to the firm at each period. Constraint (12) explains the number of used products sent back to plants, and constraint (13) finally describes that each variable in the network has a value of 1 if it is used and 0 otherwise.
The proposed dynamic model is distinguished from a traditional facility location model considering the three kinds of facilities, inventory holding, and the closed chain simultaneously [36]. For the solution of the above MINLP model, the linearization method is employed, which will be explained in detail in this section.
Let, , and. The following additional constraints guarantee the linearization of the model:
As such, the corresponding objective function of the MILP model is formulated in the following equation:
3. Solution Method
Benders method is one of the most robust heuristic methods for solving optimization problems that works with integer variables. In this method, the MILP model is converted to the sub-problem with continuous variables and the master problem. This iterative method aims to find an optimal issue to the considered problem [37]. Solving the master problem leads to finding a lower bound for the desired problem. Now, using the results of the master problem and not considering the integer variables, the subproblem is solved, which is regarded as an upper bound. In the next step, a cut is added to the problem as a new constraint; due to the reduction of the possible space, the new solution has the same or better quality (lower effective bound) than the previous one. The stop condition of this algorithm is satisfied if the difference between the upper and lower bounds will be less than a certain value. The flowchart for the overall procedure of the Benders decomposition is shown in Figure 2.
It is important to note that before starting the process in this algorithm, the problem should be as simple as possible as below.in which is the Benders subproblem which is developed in the following section.
3.1. Benders Subproblem
The subproblem is a minimization problem that finds the optimum values for the continuous variables for the fixed values of. This problem can be developed as follows:
In this model, constraints (5) and (6), which were in the form of equality, have been replaced with two inequality constraints to facilitate developing the dual subproblem. To generate cuts to add to the master problem, the duality of the subproblem is used. To obtain the duality of subproblem , dual variables , for constraints are, respectively, introduced. By using these dual variables, the duality of the subproblem, namely, , is developed as follows:
3.2. Benders Master Problem
The Benders master problem is obtained as follows:
In this model, (57) is the objective function of the Benders’ master problem, and constraint (58) is the optimality cut introduced to the master problem if the subproblem is solved optimally. Parameters are the values of dual variables obtained by solving Benders dual subproblem. Finally, constraint (59) is the feasibility cut added to the master problem, which is the subproblem and is infeasible.
3.3. The Proposed BDM
The pseudocode of the proposed Benders decomposition is presented in Figure 3. As it is shown in this figure, the process starts with an initial feasible solution for the master problem, which can be solved without any additional cuts. Then, the obtained solution for the master problem is given to the subproblem. If the subproblem is infeasible, i.e., the dual subproblem is unbounded, an unbounded ray is used to generate an infeasibility cut to add to the master problem in the next iteration. If the subproblem is feasible and solved to optimality, using the optimal solution obtained, an optimality cut is generated and added to the master problem for the next iteration. If the obtained solution provides a better upper bound, it is updated based on UB in the figure. Then, the master problem is solved again. This process is repeated until the gap between the lower and upper bounds is lower than a specified value.
4. Computational Experiments
In this section, the proposed model and the proposed BDM efficiency are evaluated. Firstly, the setting procedure of parameters in the algorithm is explained. Then, the computational results on some existing test problems are provided by the proposed algorithm. The computational experiments were conducted on a PC with a 2 GHz Intel CPU and 1 GB of RAM operating under the Windows XP operating system. The subproblems were implemented in AIMMS11.0 modeling software using CPLEX12 as the solver.
4.1. Test Instances and Determining the Value of Parameters
It is worth mentioning that there is no reliable input data in the literature for testing the proposed model. So, the used parameters is inspired based on some data generation techniques from the literature of facility location problems [37]. The instances are classified based on the number of potential facilities and the number of customers. Also, fifteen instances were considered and each case was run ten times. The problem tests are developed with the following procedure: some factories are in the supply chain to service without limitation. Every demand point includes both new and used products. For example, means that there are four demand points in the chain (the first two are for new and the rest are for used products). The fixed costs are generated using the below formula:
The fixed costs are doubled when where .
In this formula, capacity is generated from , capacity is generated from , capacity is generated from , and capacity is generated from Besides the points representing the facility, customers are uniformly distributed in unit square and points representing the plant are uniformly distributed in Demands are distributed from , the used product is distributed from , and the cost unit of is generated from Also, the percentage of is 35% and the value of allocation costs from facility to customer per unit for new product is calculated by multiplying the Euclidean distance by ten. Furthermore, the value of allocation costs from facility to customer per unit for used product is calculated by multiplying the Euclidean distance by eight, the value of allocation costs from plant to facility per unit is calculated by multiplying the Euclidean distance by 12, and the value of allocation costs from facility to facility per unit is calculated by multiplying the Euclidean distance by six. The inventory unit cost for holding the new and used product is generated from . The value of is generated from and the value of is generated from . The cost unit of is generated from . The value of is generated from .
4.2. Computational Results
In the case of a closed supply chain, which is a MIP problem, a comparison among the proposed BDM algorithm, genetic algorithm (GA), firefly algorithm (FA), and other four types of tuned metaheuristic algorithms (which are hybrid forms of GA and FA) including PHGF, EHGF, ISHGF, and WHGF [38], presented in Table 1, was used to provide an insightful demonstration. In this table, 15 instances with different data have been created, including small, medium, and large size examples. Therefore, the performance of the mentioned algorithm can be tested in a suitable range of samples. Also, in this table, the first column of the table demonstrates the size of each dataset, and the second column shows the CPU time of the CPLEX solver. The CPLEX solver can obtain the desired optimal solution for small and medium instances (up to size ), while for large size tests, CPLEX could not complete the solving process on time; therefore, the best value obtained from the objective function is reported as the best solution found. In fact, by increasing the problem’s size, the computing time by CPLEX increases too rapidly. Each algorithm has two subcolumns in which the column “NO” reports the number of times the best solution of CPLEX was reached or improved from 10 independent times, and the column “A.T” shows the CPU time in seconds.
The computational results show that the BDM heuristic can solve the model very efficiently, and this method solves the model faster than the CPLEX solver in all the cases. Most of the time, CPLEX is not invoked to obtain integer solutions, which shows that the model is integer-friendly. However, only 11 out of the 150 solutions obtained by the heuristic algorithm are not optimal, but the relative errors are less than 1.5%. Also, there is no significant difference between 7 algorithms in solving small size tests. Although the solution time of metaheuristic algorithms is less than the solution time of BDM, the Benders method is able to get the optimal or near-optimal solution in a reasonable time. Furthermore, the quality of the obtained solution and the number of optimal solutions from the proposed method are better than other ones based on Figure 4. In this figure, GA has the worst results in terms of finding the best solutions and only managed to get the best solutions in 93 out of 150 tests, while the two algorithms FA and WHGF have almost the same answers in this regard. More precisely, the FA algorithm at 93 and the WHGF algorithm at 98 retests could obtain the best solutions. Also, the two algorithms PHGF and EHGF, with 105 times finding the best solutions, have the same answers. In the end, the best algorithm presented, apart from the proposed algorithm, is ISHGF, which was able to achieve the best solutions in 122 of 150 test times. Hence, the result comparison of the Benders decomposition technique with the exact and other methods shows the validity and superiority of it. For better results, the comparison of the convergence of the methods in three types of test samples has been made. The computational results indicate that ISHGF almost obtains the best lower bounds in all sizes compared to other algorithms.
In Figures 5–7, the results of BDM, FA, GA, and ISHGF are compared. In all three figures, the horizontal axis represents the number of iterations of each algorithm and the vertical axis represents the value of the objective function obtained up to that time. One of each category of small, medium, and large problems has been selected to properly analyze the compared algorithms. In addition, for the small size example in Figure 5, all algorithms are executed 21 times, for the medium example in Figure 7, all algorithms are executed 80 times, and for the large example in Figure 8, all algorithms are executed 100 times. The result in these figures shows that the proposed algorithm has an excellent performance to achieve the best answers in the shortest time, so that the algorithm in the first example in the 11th iteration, in the second example in the 58th iteration, and in the third example in the 63rd iteration can achieve the best solutions.
As the problem size increases, the computation time required by CPLEX becomes very high, and the algorithm is not able to find any feasible solution. Since the GA, FA, ISHGF, and proposed BDM can obtain at least a feasible solution for large-scale instances, the performances of these method were reported in the figures. The computational results indicate that the time solutions of the metaheuristics are less than those of the BDM and slightly increase with the size of the problem; still the Benders method can obtain a better solution in fewer iterations, i.e., for small, medium, and large-scale tests, the stopping criteria are running the algorithms in 11, 57, and 68 iterations, respectively, while in other methods, to reach a solution as good as Benders’ solution, it takes at least twice the number of repetitions. It is remarkable that the existing time difference (albeit small) is related to solve the master problem in the Benders method, while in the other metaheuristics, only continuous linear problems are solved. To reduce the solution time, a heuristic or metaheuristic method can be used to solve the master problem. Generally, the time required increases as expected but is still reasonable for this kind of problem. Now, two instances of the test problems with different sizes and are considered to test the performance of the Benders method in solving the master problem and subproblem. For the first test (Figure 8), after 21 iterations (64 seconds) of the algorithms, the optimal value is obtained because the difference between the upper and lower bounds is zero.
After 38 iterations (408 seconds) of the algorithms (Figure 9), the difference between the upper and lower bounds is zero, and so the optimal value is obtained. In most cases, the upper bound quickly reaches the optimal solution than the lower bound, and it is because of generating optimal cuts in the related upper bound subproblems.
5. Conclusion
This paper studied a dynamic facility location problem in an environment where both forward and backward (reverse) logistic flows exist. To the best of our knowledge, this work is the first to present a dynamic closed-loop supply chain with the inventory holding for both forward and backward products. The paper formulates this problem by considering having inventory and its associated costs, flow capture between connected facilities, facilities and plants, facilities, and customers, and the fixed costs for opening, closing, and open out the facilities. This model was mathematically formulated as a mixed integer nonlinear programming model. After the linearization of the model by introducing additional variables and constraints, it was solved by a Benders decomposition method. The numerical experiment demonstrated that the solutions obtained by the proposed approach were close to optimal solutions. Furthermore, the solution time was better than other methods such as CPLEX solver. As a future research direction, it is reasonable to analyze an extension of the current model under uncertainty. Also, proposing bilevel and competitive models for the closed chain case can be right choices for the next works. Moreover, a combination of available algorithms such as branch and bound and Lagrangian relaxation could be helpful.
Appendix
The following notations are used for model formulation. All symbols related to the type of facilities are superscripted by “f,” “r,” and “h”. Also, the symbol related to the leader or follower firm is superscripted by “L” and “F,” respectively. Let be the indexes of potential facilities, and is the index of demand points, where is the set of demand point for new products, is the set of demand point for used products, and is the index of plants. The other parameters and variables are defined in Table 2. Decisions regarding the dynamic configuration of the logistics network include decisions about the type of facility to be built in each. Also, potential facilities of the forward and backward products should be shipped at each period link. The model decision variables are shown in Table 3.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.