A hybrid differential evolution algorithm for a location-inventory problem in a closed-loop supply chain with product recovery

Product recovery is an important business because of its great economic, social, and environmental benefits in practice. In this paper, a location-inventory problem (LIP) in a closed-loop supply chain (CLSC) is investigated to optimize facility location and inventory control decisions by considering product recovery. The objective is to optimize facility location and inventory control decisions to minimize the total cost of business operations in a closed-loop supply chain system. We formulate this problem as a mixed-integer nonlinear programming model and design a modified hybrid differential evolution algorithm (MHDE) to solve it efficiently. Finally, numerical results are presented to validate the performance of the new algorithm. The results show that MHDE is more efficient and effective than Lingo and other algorithms for the research problem under study. Managerial insights are also derived for business managers to improve their supply chain performance.


Introduction
The economic benefits, environmental concerns, and social impacts have become an important motivation to implement a closed-loop supply chain integrated with efficient B Gang Liu lgxueshu@163.com B Ying Zhang yz29@buffalo.edu 1 School of Management, Wuhan Textile University, Wuhan 430200, Hubei, China necessary for business organizations to incorporate returned products for recovery into their traditional distribution channels.
The performance of a supply chain can be improved by better strategic and tactical decisions. Currently, there is an emerging trend to integrate strategic with tactical decisions through the development of joint location-inventory models [10]. Facility location decisions are focused on maximizing the efficiency of a supply chain, while inventory control models have focused on improving the responsiveness of the supply chain [11]. In many cases, integrating strategic decisions with inventory control policies in a closed-loop supply chain can help overcome short-sighted sub-optimal network designs [12]. Hence, it will be greatly beneficial to optimize location and inventory decisions in a closed-loop supply chain simultaneously to improve customer return and product recovery processes. This work is motivated by realworld practice, and it will address the following business questions: 1. What are the optimal number and locations of the facilities when product recovery is considered in a closed-loop supply chain? 2. In such a closed-loop supply chain, how to assign customers to different types of facilities to improve supply chain efficiency? 3. What is the impact of cost factors, facility locations, and inventory policies on the performance of the closed-loop supply chain under study?
To answer these questions, this paper incorporates product recovery into an integrated location-inventory problem in a closed-loop supply chain. With the presence of the product recovery business in a closed-loop supply chain, the solution to this problem will give (1) facility locations, (2) the assignment of customer zones to the facilities in forward and backward logistics flows, and (3) order time and size for each facility. To formulate this problem, a mixed-integer nonlinear programming model is developed to minimize the total cost in a two-stage closed-loop supply chain. Moreover, a modified hybrid differential evolution algorithm (MHDE) is designed to solve this model given the complexity of the problem. The numerical results show that MHDE can obtain good feasible solutions efficiently.
The rest of this paper is structured as follows: "Literature review" reviews research works related to this study. "The model" presents a mix-integer nonlinear programming model for the research problem under study. "Solution methodology" proposes an efficient meta-heuristic algorithm to solve this model efficiently. "Numerical experimentation" presents numerical experiments and discusses computational results. "Conclusions and future research" concludes this paper and suggests directions for future research.

Literature review
Nowadays, many research efforts are devoted to integrated supply chain designs because of their great economic and social benefits, and integrated location-inventory problems (LIPs) are an important topic in this research field. In this paper, we study a location-inventory problem in a closedloop supply chain by incorporating product recovery, and it is related to three research streams: closed-loop supply chain, location-inventory problems, and DE algorithms.
There are numerous works on closed-loop supply chains (CLSCs) in response to the rising awareness of environmental issues, which are reviewed thoroughly by Shekarian [13] and Peng et al. [14]. In the literature, CLSCs are extensively studied by incorporating many real-world businesses, such as take-back legislation [15], reshoring drivers [16], pricing decisions and discounts [17], social factors [18], and pricesensitive demand [19]. For CLSCs, many research efforts have been focused on the disposal of end-of-life products [20], but the consolidation of location-inventory problems and product recovery is rarely studied in the literature. To better represent the real-world practices, this paper considers both businesses under a CLSC setting, and hence fills this research gap in the literature.
In the literature, LIPs are usually studied by incorporating real-world business scenarios. This topic was first studied by Daskin, Coullard, and Shen [21], and Farahani et al. [22] present a comprehensive review of the related research efforts. Recently, location-inventory problems are generalized in several different directions. Araya-Sassi, Paredes-Belmar, and Gutiérrez-Jarpa [23] studied the distribution centers' location model that incorporates the different review inventory control policies at distribution centers. Puga and Tancrez [24] analyzed a location-inventory problem for the design of supply chain networks with uncertain demands. Also, LIPs are studied by incorporating other important topics, such as disruption risk [25], lateral transshipment policies [26], post-disaster humanitarian logistics [27], multiple customer priority classes [28], and so on. However, none of the above papers have considered reverse logistics network design and the state in which product recovery affects the facility location and inventory control decisions.
Recently, LIPs are extended with the development of closed-loop supply chains. Al-Salem et al. [7] studied an LIP with three different types of warehouses. Diabat, Abdallah, and Henschel [12] studied an LIP in an uncapacitated closedloop supply chain where returns are remanufactured as spare parts and then sold to retailers. Asl-Najafi et al. [29] formulated a bi-objective dynamic LIP by considering disruption risk in a CSLC. Zhang and Unnikrishnan [30] studied a coordinated LIP with uncertain demand. Li, Guo, and Zhang [31] studied an integrated LIP in a CLSC with third-party logistics. Guo et al. [32] proposed a new multi-commodity LIP by incorporating commercial product returns in CLSCs. As an extension, Guo et al. [33] developed a more comprehensive model by considering the secondary market of used and refurbished products, and present a more powerful algorithm to solve the extended model.
Facility location problem (FLP) is a typical NP-hard problem [34]. As an extension of FLP, LIP which integrates location and inventory decisions is more complicated. Given its complexity, many heuristic algorithms have been developed to solve LIPs, such as simulated annealing [10], hybrid genetic algorithm [25] [28], particle swarm optimization [29], differential evolution [31] [32], hybrid harmony search [35], and tabu search [36]. In addition to those algorithms, differential evolution (DE) [37] has attracted much attention because of its good capability of solving NP-hard optimization problems. For example, Li et al. [31] propose an improved hybrid differential evolution (IHDE) algorithm based on DE and GA to solve the LIP in a CLSC with third-party logistics. Guo et al. [33] design a normal distributed adaptive differential evolution (NDADE) algorithm based on a novel adaptive mutation strategy to solve LIP in presence of the secondary market. Ali et al. [38] implemented a novel discrete differential evolution (NDDE) algorithm to solve complex discrete traveling salesman problems. A hybrid differential evolution algorithm and genetic operator (HDEGO) with the fuzzy logic controller is proposed to solve the multi-trip vehicle routing problem with backhauls and heterogeneous fleet [39]. A discrete differential evolution (DDE) algorithm with new mutation and crossover operators is proposed to find near-optimal solutions for order rescheduling problem [40]. A hybrid approach (iDE-EA) is proposed for a closedloop facility layout problem [41]. Guo et al. [42] develop an improved self-adaptive differential evolution (ISaDE) algorithm by introducing an effective adaptive mechanism to solve LIP with secondary market consideration. This paper proposes a new DE-based algorithm, MHDE, to solve the research problem under study. It is well known that the control parameters of meta-heuristic algorithms play an important role in their performance. MHDE adopts selfadapting control parameters in [33] and [42] and develops a hybrid truncation-selection strategy [31] to solve locationinventory problems.
Hence, this study extends the current research from three perspectives. First, a mixed-integer nonlinear programming model is formulated to incorporate LIPs and product recovery in a CLSC. Second, a new meta-heuristics algorithm is designed to solve this model efficiently. Finally, managerial insights are provided for business managers to improve the performance of their supply chains.

Problem description
In this paper, we study a location-inventory problem in a closed-loop logistics network including a hybrid production refurbishment center (HPRC), distribution centers (DCs), collection centers (CCs), and hybrid distribution-collection centers (HDCCs). In this network, an HPRC is a production and refurbishment facility. HDCCs are hybrid facilities that operate as distribution centers for new or refurbished products in the forward logistics flow, and collection centers for returned products in the reverse logistics flow. Customer zones with uncertain demands can be assigned to DCs or HDCCs in the forward flow, and CCs or HDCCs in the reverse flow. There are many advantages to building an HPRC and HDCCs, such as cost savings and pollution reductions as the result of sharing materials, equipment, and infrastructures [31].
As illustrated in Fig. 1, new and refurbished products are shipped from the HPRC to customer zones via DCs or HDCCs in the forward flow, and returned products are collected from customer zones to the HPRC via CCs or HDCCs for recovery or safe disposal [7]. Returned products can be recoverable or unrecoverable depending upon their conditions which are inspected in CCs or HDCCs. In this study, customer zones are predetermined and fixed. It is assumed that new and refurbished products are the same to fulfill the demands of customer zones, and the candidate locations for different types of facilities are exclusive.
The LIP under study is formulated as a mixed-integer nonlinear programming model to minimize the total cost of a CLSC which includes an HPRC and several DCs, CCs, and HDCCs. The solution to this model will provide an optimal set of facility locations, inventory policies, and the assignment between the facilities and customer zones. The notation is as follows: Sets I set of candidate locations for DCs; J set of candidate locations for CCs; K set of candidate locations for HDCCs; M set of customer zones;

Parameters
A I fixed (annual) cost of opening distribution center i, for each i ∈ I ; a j fixed (annual) cost of opening collection center j, for each j ∈ J ; for each k ∈ K ; b r fixed cost per order of new/refurbished products at DC/HDCC r, for each r ∈ R; b s fixed cost per order of returned products at CC/HDCC s, for each s ∈ S; c r fixed cost per shipment from the HPRC to DC/HDCC r, for each r ∈ R; c s fixed cost per shipment from CC/HDCC s to the HPRC, for each s ∈ S; e r shipment cost per unit of new/refurbished products between the HPRC and DC/HDCC r, for each r ∈ R; e s shipment cost per unit of returned products between the HPRC and CC/HDCC s, for each s ∈ S; d rm shipment cost per unit of new/refurbished products between DC/HDCC r and customer zone m, for each r ∈ R and m ∈ M; d sm shipment cost per unit of returned products between CC/HDCC s and customer zone m, for each s ∈ S and m ∈ M; h r (annual) holding cost per unit of new/refurbished products at DC/HDCC r, for each r ∈ R; h s (annual) holding cost per unit of returned products at CC/HDCC s, for each s ∈ S; q m quantity of returned products at customer zone m, for each m ∈ M; μ m mean (daily) demand at customer zone m, for each m ∈ M; σ 2 m variance of (daily) demand at customer zone m, for each m ∈ M; α service level; z α standard normal deviate such that P(z ≤ z α ) α; L order lead time (in days) at distribution centers; λ working days per year.
Decision variables X I 1 if a DC is opened at location i, and 0 otherwise, for each i ∈ I ; X j 1 if a CC is opened at location j, and 0 otherwise, for each j ∈ J ; X k 1 if an HDCC is opened at location k, and 0 otherwise, for each k ∈ K ; Y rm 1 if DC/HDCC r fulfills the demand at customer zone m, and 0 otherwise, for each r ∈ R and m ∈ M; Z sm 1 if CC/HDCC s collects returned products at customer zone m, and 0 otherwise, for each s ∈ S and m ∈ M.

Model formulation
The total supply chain cost comprises two types of costs: (1) Location costs which include the fixed costs of opening facilities, and the shipping costs between facilities and customer zones; (2) inventory costs which include working inventory and safety stock costs.
(1) Location costs The fixed cost to open and run facilities is i∈I a i X i + j∈J a j X j + k∈K a k X k . Moreover, the shipping cost given facility locations is λ r ∈R m∈M μ m d rm Y rm + λ s∈S m∈M q m d sm Z sm . Thus, the total location cost is (2) Inventory costs To satisfy the demand at a customer zone, a facility in the forward flow needs to maintain two types of inventories: working inventory, which depends upon the order policy, and safety stock, which is maintained to hedge the uncertainty in the demand and delivery time [43].

(2.1) Working inventory cost
In this model, we assume that DCs and HDCCs order from the HPRC using the economic order quantity (EOQ) model, which is common to approximate the (Q, r) inventory model with type I service [44]. Thus, the EOQ model is used to optimize inventory control decisions. Let n r be the number of orders placed at DC/HDCC r to the HPRC per year. The fixed annual order cost for DC/HDCC r is b r n r , the annual shipping cost from the HPRC to DC/HDCC r is c r n r +λe r m∈M μ m Y rm , and the annual holding cost of new and refurbished products at DC/HDCC r is λh r m∈M μ m Y rm 2n r . Therefore, the total annual working inventory cost at distribution center r is given by Taking the first-order derivative of C W F (r ) with respect to n r and setting it to zero, we can obtain that n * r λh r m∈M μ m Y rm 2(b r +c r ) . Substituting this into C W F (r ), the total annual working inventory cost at DC/HDCC r can be rewritten as follows: In the reverse flow, returned products will be shipped from CCs/HDCCs to the HPRC for recovery or disposal. Similar to C * W F (r ), the total annual working inventory cost at CC/HDCC s is Therefore, the total annual working inventory cost in the CLSC is (2.2) Safety stock cost DCs and HDCCs need to keep safety stocks to maintain an appropriate service level under demand uncertainty. Given the risk pooling effects [45], the safety stock at DC/HDCC r is z α L m∈M σ 2 m Y rm to ensure that the probability of stock-outs is less than or equal to α. Therefore, the total safety stock cost is Let C T OT C L + C W I + C SS be the total cost in the CLSC. Then, the location-inventory problem under study can be formulated by the mixed-integer nonlinear programming model below Z sm ≤ X s , ∀s ∈ S, ∀m ∈ M; ( 7 ) The objective function (1) aims to minimize the total cost in the CLSC. Constraints (2) and (3) mean that at least one facility will be established in the forward and reverse flows, respectively. Constraints (4) and (5) ensure that each customer zone will be assigned to one facility in the forward and reverse networks, respectively. Constraints (6) and (7) mean that a customer can only be assigned to an existing facility. Constraints (8)- (11) mean that decision variables are binary.

Solution methodology
The facility location problem is NP-hard [34], and the joint location-inventory problem is also NP-hard given its complexity in structure. In the literature, many LIPs are solved by heuristic and meta-heuristic algorithms because of their great efficiency to solve such problems. Differential evolution [37] is an important approach that is widely applied to solve LIPs [31,32]. Although DE is a powerful tool to solve nonlinear problems, it has some drawbacks such as prematurity and local optimism. To solve the LIP under study, the MHDE algorithm is designed to enhance the performance of DE. Specifically, MHDE is improved in four ways: (1) The opposition-based learning (OBL) [46] is introduced to improve the quality of the initial population. (2) The mutation factor F is generated from a Gaussian distribution N(0,1) to enhance the diversity of the population and the global search capability. (3) The new crossover strategy is adopted to increase the population diversity and solution search ability. (4) The selection operation adopts truncation selection of GA to overcome the limitation of the greedy selection mechanism in DE. The notations used in MHDE are introduced in Table 1.
As a DE-based algorithm, MHDE uses evolutionary operators (i.e., mutation, crossover, and selection) at each generation to update the population toward the optimum after initialization, and the details are as follows.

Individuals
IN the MHDE algorithm, individuals in a population are corresponding to candidate solutions to a discrete combinatorial optimization problem, and a good representation of individuals is important to the operations of a DE algorithm. For the LIP under study, an individual is composed of two parts. The first part is the assignment information in the forward logistics flow, and the second part is the assignment information in the reverse logistics flow. Hence, an individual can be represented by an array with 2M entries, which is given by ( 1 2 ) where Z g p is the pth individual in iteration g.
is the index of a DC or HDCC that is assigned to a customer zone, and y g p, m ∈ {I + 1, . . . , I + J + K }(1 ≤ p ≤ P; M + 1 ≤ m ≤ 2M) is the index of an HDCC or CC that is assigned to a customer zone. If x g p, m r (1 ≤ m ≤ M), then the demand in customer zone m will be fulfilled by DC or HDCC r in the forward flow. If y g p, m s(1 ≤ p ≤ P; M + 1 ≤ m ≤ 2M), then the returned products from customer zone m-M will be collected to HDCC or CC s in the reverse flow. Figure 2 illustrates an example individual corresponding to ten customer zones, three distribution centers, three HDCCs, and two collection centers.

Entries in an individual array
The elements in an array which represent an individual are candidate facility locations, and they are given by Eqs. (13) and (14) x 0 ( 1 4 ) where round is the rounding function, rand is a uniformly distributed random number in range [0, 1], x L and x U represent the lower and upper bound of x g p, m , respectively; y L and y U be the lower and upper bound of y g p, m . We have x L 1,

Initialization with OBL
An evolutionary algorithm will usually generate an initial population and then improve it toward an optimal solution. In this study, we can increase the opportunity of generating a better initial population by checking the opposite population of an initial population. Hence, OBL is adopted to obtain a better initial candidate population, and the steps are as follows: Step 1: Generate the initial population Z 0 p randomly.
Step 2: Obtain the opposite population Z 0 p . Step 3: Select the P best solutions from Z 0 p ∪ Z 0 p as the initial population.

Mutation operation
In  2M . In this study, the mutation operator "DE/rand/1" is adopted, which gives a mutant vector as Eqs. (15) and (16) where indices r 1 , r 2 , and r 3 are different integers that are randomly selected from [1, P] and not equal to i. round is the rounding function, and F is a mutation scaling factor that affects the differential variation between two individuals. If vx g+1 p, m or vy g+1 p, m exceeds the upper or lower bounds, then it will be randomly generated again until it is within the range.
As an improvement to DE, an adaptive mutation factor is adopted in MHDE to enhance the population diversity and global search capability, which is given by F p, g F · ξ p, g , p 1, 2, . . . , P, g 1, 2, . . . , G, where ξ p,g is generated independently from the normal distribution, and F amplifies the variation.

Crossover operation
Crossover operation is introduced in DE to increase the diversity of the perturbed parameter vectors. In each iteration, the crossover operator will be applied to each mutant vector and its corresponding target vector to yield a trial vector U Z where rand(l) ∈ [0, 1] is a random number generated from a uniform distribution, m 1 ∈ [1, 2, . . . , M] and m 2 ∈ [M + 1, M + 2, . . . , 2M] are two random numbers that ensure that trial vector UZ p g+1 gets at least one parameter from mutant vector VZ p g+1 , k is a random number whose range is [1, 2, . . . , P], and CR is a pre-defined crossover rate.
MHDE also uses an adaptive mechanism in the crossover operator to further increase the diversity of the population, for which the crossover factor CR changes adaptively by Eq. (20) where rand1 and rand2 are two random variables that are uniformly distributed on [0, 1], and τ is a constant which is the probability of changing CR. Note that τ 0.9 in this study.

Selection operation
The selection operation determines whether the trial vector UZ p g+1 or the target vector Z p g can get into the next generation. A greedy criterion is usually used for the selection operation in DE. In this study, a truncation selection from the genetic algorithm is adopted to overcome the limitation of the traditional greedy criterion in DE [31], and the steps are as follows: Step Step 2: Compute and sort the fitness value of the new population; Step 3: Select the first P individuals as the next generation of the population Z

Stopping criterion
The mutation, crossover, and selection operations will be repeated until a termination criterion is satisfied. In this study, MHDE will be terminated if one of the following criteria is satisfied: 1. Stop if no better solution can be generated in consecutive K generations (K is a pre-defined counter value).

Algorithm flow
In summary, MHDE consists of the following steps: Step 1: Initialization. Set k 0 and g 0, and generate an initial population with P target individuals by applying OBL.
Step 3: Generate P trial individuals. For p 1 to P, repeat Steps 3.1-3.2.
Step 3.1: Mutation operation. Update mutation factor F using Eq. (20) and generate a mutant vector VZ p g+1 by using Eqs. (18) and (19). If an infeasible solution is generated, repeat this step until a new feasible solution is generated.
Step using the truncationselection scheme.
Step 5: Set g g + 1. If f Z g best f Z g+1 best , set k k + 1. Otherwise, set k 1.
Step 6: If g G or k K, stop and output Z best and f (Z best ). Otherwise, go to Step 3.

Numerical experimentation
In this section, we present a series of experiments to evaluate the performance of MHDE and compare it with the other three algorithms in the literature, i.e., ISaDE [42], NDADE [33], and DE [37]. First, parameter analysis is conducted to find the optimal settings for the four algorithms under comparison. Furthermore, LIPs in different sizes are solved by the four algorithms to evaluate their performances. Finally, a sensitivity analysis is conducted to analyze the impact of business parameters and derive managerial insights. All the experiments are implemented by Java JDK 2020-06 on a personal computer (COU: Intel(R) Core(TM) i5-9500 T CPU @ 2.20 GHz 2.21 GHz, RAM: 4 GB, Operating System: Windows 10).

Parameter settings
In this study, the candidate facility locations and customer zones are randomly generated in a 100 × 100 square. The shipping costs between facilities and customer zones are their Euclidean distance, and the other parameters in the LIP are shown in Table 2.
In this study, parameter analysis is conducted to find the optimal settings for the four algorithms, because the quality of solutions relies heavily on their parameters. P, K, and G are the basic parameters in DE-based algorithms. Storn and Price [37] proposed that P ∈ [5M, 10M], and Li, Guo, and Zhang [31] suggested that this range can be further extended for LIPs, which are P ∈ [2M, 10M]. Hence, we are using the following parameters in this study: (1) K M in all experiments, (2) G 10 M in all experiments, (3) P 3 M in small-size problems (10 candidate DCs-5 candidate HDCCs-5 candidate CCs-100 customer zones), P 4 M in medium-size problems (15 candidate DCs-8 candidate HDCCs-7 candidate CCs-150 customer zones), and P 5 M in large-size problems (20 candidate DCs-10 candidate HDCCs-10 candidate CCs-200 customer zones).
Moreover, F and CR are chosen by tuning tests, for which an example network that consists of 10 candidate DCs, 5 candidate HDCCs, 5 candidate CCs, and 100 customer zones is used. Each combination of the two parameters was executed 30 times, and the parameters are evaluated by the average objective value and runtime. Note that the runtimes are given in seconds, and the optimal objective value from Lingo is 18372610. The parameter test results are shown in Tables 3,  4, 5, and 6.

Performance comparison
In this section, we compare the performance of MHDE with the other three DE-based algorithms in the literature (i.e., ISaDE, NDADE, and DE) and a commercial solver (i.e., Lingo). The experiments are conducted on nine example problems in three different sizes, which are 3 small-size, 3 medium-size, and 3 large-size problems, respectively. For each problem, MHDE, ISaDE, NDADE, and DE were executed 30 times because of the random nature of DE-based algorithms, while the problem is solved only once in Lingo.
The numerical results are shown in Tables 7, 8  Since Lingo is not able to solve large-size problems after excessively long runtimes, the numerical results from Lingo are not provided in Table  9.
The following conclusions can be made from Tables 7, 8  and 9: 1. The optimal solutions from MHDE and Lingo are almost identical for small-size and medium-size problems, but MHDE is more efficient than Lingo.  and algorithm efficiency. Moreover, gaps are larger when the problem size increases. This indicates that MHDE is a more effective approach and it will have more advantages when the model is bigger. 3. MHDE has less standard deviation than ISaDE, NDADE, and DE in each instance, which means that MHDE is more stable and hence more likely to get good solutions with fewer executions of the algorithm.
To justify the performance of MHDE, non-parametric statistical tests for independent populations can be useful for pairwise comparison between solutions of the algorithms. The Mann-Whitney U test is an appropriate candidate for this situation. This test with a 95% confidence interval is used to verify whether there is a significant difference in the total cost. Tables 10, 11 and 12 show the test results for smallsize, medium-size, and large-size problems, respectively. We did all statistical computations in the IBM SPSS statistics 25 software.
From Tables 10, 11 and 12, we can see that p values are less than 0.05 in most problems. MHDE is statistically significantly better than other approaches. Therefore, it is concluded that MHDE is superior to these other algorithms for solving the proposed problem. Hence, MHDE is a more effective approach that can provide better feasible solutions.

Managerial implications
In this section, sensitivity analysis is performed to analyze the impact of business factors on supply chain performance. Figures 3, 4, 5, 6, 7 and 8 show the changes of the parameters being tested and the changes of supply chain costs accordingly. For example, Fig. 3 shows that C L , C I , and C T will increase as q m increases. Figure 4 shows that C L and C T will increase, but C I will remain the same when d rm or d ms increases. Figures 5, 6, 7 and 8 show that C I and C T will increase but C L will remain the same when any of b r , b s , c r , c s , e r , e s , h r , h s increases. According to the numerical results above, some managerial insights are obtained as follows: 1. Considering the complex environment of product recovery, numerous items are involved in the recovery operations. Managers should pay attention to product recovery to increase sales and profitability and reduce the environmental impact. 2. From our experiments, the more increase in the parameter q m , the more increase in the total cost of the supply chain system. The parameter q m has a significant impact on supply chain costs, location decisions, and inventory policies. Therefore, managers should pay attention to the related cost factors. 3. The corresponding adjusting direction of cost parameters is provided to control the total cost for enterprises according to the results of the cost sensitivity analysis. 4. The numerical results show that MHDE is a fast algorithm to find more feasible solutions and save the cost of the supply chain system. Therefore, MHDE can help recovery businesses to solve complex location-inventory problems.

Conclusions and future research
In this paper, an integrated location-inventory problem in a closed-loop supply chain with product recovery is studied. This problem is formulated as a mix-integer nonlinear programming model and solved by a novel heuristic algorithm which guarantees a good solution accuracy with great time  efficiency. This study presents a new methodology to optimize facility location, inventory replenishment, and demand planning decisions, and it also provides a powerful tool to evaluate the business impact of real-world factors and identify opportunities to improve supply chain performance.
This study can be extended in several directions. For example, this problem can be extended to a more complicated location-inventory-routing problem which considers vehicle routing between facilities and customer zones. It can also be extended by considering more business processes and scenarios related to product recovery such as product disassembly. The algorithm can also be improved for better solution quality and time efficiency.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.