Hybrid Harmony Search-Simulated Annealing Algorithm for Location-Inventory-Routing Problem in Supply Chain Network Design with Defect and Non-Defect Items

: This paper considers a location-inventory-routing problem (LIRP) that integrates the strategic, tactical, and operational decision planning in supply chain network design. Both defect and non-defect items of returned products are considered in the cost of reverse logistics based on the economic production quantity model. The objective of the LIRP is to minimize the total cost of location-allocation of established depots, the cost of inventory, including production setup and holding cost, as well as the cost of travelled distance by the vehicles between the open depots and assigned customers. A Hybrid Harmony Search-Simulated Annealing (HS-SA) algorithm is proposed in this paper. This algorithm incorporates the dynamic values of harmony memory considering rate and pitch adjustment rate with the local optimization techniques to hybridize with the idea of probabilistic acceptance rule from simulated annealing, to avoid the local extreme points. Computational experiments on popular benchmark data sets show that the proposed hybrid HS-SA algorithm outperforms a standard HS and an improved HS for all cases.


Introduction
The need for designing networks in a supply chain has been attracting a lot of attention in recent years; a result of the booming development of industrial commerce. A practice in traditional commerce is that customers can return the unsatisfactory products they bought within a given certain period of time. This process is called the reverse logistics by which the products are returned for proper planning of their disposal [1]. Reverse logistics can also be defined as the process of moving goods for capturing value. Efficient reverse logistics can reduce cost especially transportation cost, shipping cost, and inventory cost. Besides that, it also improves the service by offering a fast response to the demands and returns from customers. By implementing reverse logistics in supply chain network design, the expectation of a gain in both profit and increased customer retention can be expected, benefiting both from the economic and social factor of the economy respectively [2].

Literature Review
The warehouse LRP (WLRP) is first proposed by [8]. The problem is decomposed into three sub-problems: Multi-Depot Vehicle Routing Problem (MDVRP), Warehouse Location-Allocation Problem (WLAP), and Multi-Depot Routing Allocation Problem (MDRAP) and is solved using the heuristic approach. Some other studies divided the problem into two phases [9]. The first phase is a LAP and the VRP follows in the second phase. The problem is solved sequentially and iteratively by SA. Reference [10] discussed the capacitated LRP under disruption and implemented several approaches in problem-solving such as metaheuristic based on maximum-likelihood sampling method, route-allocation improvement, two-stage neighborhood search and SA.
A two-phase hybrid heuristic search approach for LRP is introduced in [11]. The problem is decomposed into LAP and VRP. For the construction of the location distribution, they proposed a Tabu Search (TS) to determine a good configuration of the facilities to be used. Meanwhile, an Ant Colony Algorithm (ACO) was used at the routing phase to obtain good routes for the given configuration. The combination of Particle Swarm Optimization (PSO) with several heuristic algorithms is called Hybrid PSO (HybPSO-LRP) and is introduced in [12]. They combined the Multiple Phase Neighborhood Search-Greedy Randomized Adaptive Search Procedure (MPNS-GRASP) algorithm, the Expending Neighborhood Search (ENS) strategy, and a Path Relinking (PR) strategy into the PSO. With the same approach used in [12], reference [13] combined the same techniques with the Honey Bee Mating Optimization (HBMO). They solved the problem using a hybrid HBMO for LRP (HBMO-LRP).
We found several studies on the integrated IRP from literature. For example, reference [14] solved a multi-objective coordination in a supply chain with three objective functions: maximizing the firm's profit, minimizing vendor's routing, and maximizing the average service level of buyers. They determined the vehicle routing, the router assigned and the amount delivered for each buyer. A HS algorithm is proposed due to the complexity of the problem. Considering the environmental impact, authors in [15] proposed a green IRP to minimize the total cost of inventory and routing, vehicle fixed cost, and emissions that are determined by load, distance, speed, and vehicle characteristics. They employed an exact method which is a Branch-and-Cut Algorithm (BCA) to add user cuts to the model and used a commercial solver to solve the linear programming relaxation problem. Similar to [16], they analyzed the benefits of collaborative energy use (CO 2 emission) and perishability in IRP with uncertain demand. The authors minimized the cost of routing, based on fuel and wage cost, waste, driving time and inventory. The problem is solved by the ILOG-OPL development studio and CPLEX 12.6 optimization.
A new bi-objective mathematical model which takes into account the economic, social, and environmental issues in the distribution of perishable products with a specified expiration date in IRPs is considered in [17]. They focused on inventory and distribution costs for the first objective and then looked at social issues to shape the second objective by considering the rate of accidents for the vehicle and the number of expired products. The Torabi-Hassani method was used and included vehicle noise emission as a constraint. The integrated Production, Inventory, and Distribution Routing Problem (PIDRP) was developed and solved heuristically with two solutions approach [18]. In the first phase, they solved the Mixed-Integer Program (MIP) with CPLEX MIP solver and in the second phase, they solved an associated consolidation problem to handle the potential inefficiency of direct shipments involving routing and inventories by Extended Optimal Partitioning (EOP) procedure. This approach can simultaneously coordinate the production, inventory, and transportation operation in PIDRP. The two-step procedure is getting more attention in solving IRP, researchers in [19] performed a comparative analysis of a series of heuristics for manufacturing supply chain with a MIP formulation and two-step procedure development. The first step is to estimate the daily delivery quantities and followed by solving a VRP in the second step. A previously developed branch-and-price algorithm is used and the production decision and inventory flow balance were taken into account in the model.
Besides IRP, the integrated LIP is also gaining popularity. An optimization model for perishable products with fuzzy capacity and carbon emission constraints for integrated LIP is proposed by [20].
They formulated a Mixed-Integer Nonlinear Programming (MINP) and is solved by using two methods which are Hybrid Genetic Algorithm (HGA) and Hybrid HS (HHS). They found that HHS gave a better solution but needed higher computing time when compared to the HGA. The LIP under vendor-managed inventory was considered by [21]. They solved the multi-objective LIP model based on the Non-dominated Sorting Genetic Algorithm (NSGA-II). This method gave promising solutions when compared with the well-known multi-objective evolutionary algorithm. By using the same approach as in [21] but hybridized with the Multi-objective Particle Swarm Optimization (MOPSO), a bi-dynamic closed-loop LIP under facility disruption risk is proposed in [22]. They considered the effectiveness of returned products in a period and retailer demand in the next period. The bi-level programming problem where the upper level is for determining the appropriate location of 3rd checking sites and the lower level is to present the coordinated inventory replenishment is studied by [23]. They considered the location and inventory policies in the supply chain when product returns are allowed.
Other related study of LIP can be found in [24]. The authors focused on three decisions which are the determination of depot locations, the allocation of flows, and the shipment sizes. Nonlinear continuous programming is formulated and an iterative heuristic is developed to estimate the depot flows, to solve the linear program, and to improve the flow depot. Meanwhile, researchers in [25] explored the continuous non-linear formulation for LIP with uncertain demand that integrates location, allocation, and inventory decisions. They solved the linear program using a heuristic algorithm and used the solution to improve the variable estimators for the next iteration. A closed-loop LIP is considered by [26] that determines which depot and remanufacturing centre need to be opened. An exact two-phase Lagrangian relaxation algorithm is developed and a mixed-integer nonlinear location-allocation model is formulated. Extension to the study in [26], an equivalent formulation with fewer nonlinear terms in the objective function is proposed in [27]. They used piecewise linearization to transform and solve it using CPLEX. Then the solution in CPLEX is compared with two previously published Lagrangian relaxation-based heuristic.
With respect to the LIRP, some of the researchers focus only on the non-defect items in reverse logistics such as in [28]. They solved the LIRP of the e-supply chain environment by using a Hybrid Genetic Simulated Annealing Algorithm (HGSAA) and compared it to a standard GA. Reference [29] used a Pseudo-Parallel Genetic Algorithm Integrating Simulated Annealing (PPGSA) to solve stochastic LIRP with continuous review inventory policy. Reverse logistics considering both defect and non-defect items of returned products has also been studied by [30]. They considered both quality of returned products in the e-commerce supply chain system. They found that the Hybrid Ant Colony Optimization algorithm (HACO) was considerably efficient and effective compared to a standard ACO. New TS (NTS) has been proposed by [31] to simultaneously integrate the LIRP and considered the minimization of manufacturing goods for new products and remanufacturing goods for damaged products in the objective function.
Based on the literature, we believe that each of the decision planning stages in the supply chain has equal significance. Recently, many industries aim to minimize the cost of suppliers by producing their own products. As the industries nowadays are trying to be greener and at the same time handle the reverse logistics effectively, they prefer to deal directly with the customers from manufacturing stage until the process of recovering the returns products. The EPQ model is gaining attention to be practiced among industries, especially among manufacturing and remanufacturing companies. With this in mind, the LIRP with EPQ model in reverse logistics is considered in this study to give new insights and solutions to the manufacturing industries in optimizing the production and inventory.

Mathematical Formulation of LIRP with EPQ Model
Given a set of the potential depots with fixed locations and customers with known demands and returns, the LIRP with EPQ model in reverse logistics is designed to determine the optimal location and number of open depots, the production quantity produced at each open depot as well as the routing of the assigned customers with known demands and returns to the open depots. The depots will function as a distributor as well as a collector. At the depot, the condition of the returned products will be inspected and classified into two categories: defect and non-defect items. The production setup cost will be affected where the non-defect items re-enters into the market, while the defect items will need to be counted in the process of production planning. The objective of the LIRP is to minimize the cost of establishing facilities, cost of inventory including production setup and holding, and cost of distance travelled by vehicles. Figure 1 shows an example of LIRP in the supply chain with three depots and nine customers. It can be seen from Figure 1 that customers 1, 2, and 3 are served by depot 1, while the remaining 6 customers are served by depot 2 using two vehicles. Depot 3 remains closed.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 22 products will be inspected and classified into two categories: defect and non-defect items. The production setup cost will be affected where the non-defect items re-enters into the market, while the defect items will need to be counted in the process of production planning. The objective of the LIRP is to minimize the cost of establishing facilities, cost of inventory including production setup and holding, and cost of distance travelled by vehicles. Figure 1 shows an example of LIRP in the supply chain with three depots and nine customers. It can be seen from Figure 1 that customers 1, 2, and 3 are served by depot 1, while the remaining 6 customers are served by depot 2 using two vehicles. Depot 3 remains closed. Input Parameters: D j = demand of customer j. r j = non-defect items returned by customer j. s j = defect items returned by customer j. R j = returned products by customer j (R j = r j + s j ). d ij = distance from i to j. VC k = capacity of vehicle k. N = number of customers. F i = fixed operating cost of depot i. DC = distance cost per mile. VD i = maximum through at the depot i. KC = setup cost per production. h = holding cost per unit inventory. P = production rate. TL ijk = total load of customer j at depot i by vehicle k.
Decision Variables: Q i = optimal number of production quantity for each depot i, U lk = auxiliary variable for sub-tour elimination constraint in vehicle k of customer l.
subject to: k∈K i∈U x ijk ∈ {0, 1} , ∀i ∈ I, ∀j ∈ J, ∀k ∈ K (12) The objective function of LIRP (Equation (1)) is to minimize the total fixed operating cost of depots, the total distance travelled cost by the vehicles, and the total cost of production setup and inventory holding. Constraints in Equations (2) and (3) indicate that each of the customers must be assigned in a single route and it can be served by only one vehicle. The total demands and total returns at each route cannot exceed the vehicle capacity limit and the total load on the vehicle at any arc must not exceed the vehicle capacity. These constraints are shown in Equations (4)-(6), respectively. Equation (7) states that the vehicle must start and end at the same depot. Besides the vehicle capacity limit, the capacity constraint for the depot is given in Equation (8). Equation (9) represents the new sub tour elimination constraint and Equation (10) specified that the customer will be assigned to the depot if there is a route from the depot. The binary values on the decision variable and the positive values for the auxiliary variable are defined in Equations (11)-(15), respectively.
Based on the objective function given in Equation (1), the total cost of inventory (TCI) consists of production setup cost and holding inventory cost. The optimal value of Q i can be obtained by deriving the first derivative of the cost function in Equation (16) with respect to Q i as follows: TCI = Production setup cost + Inventory holding cost Let ∂TCI The cost function of TCI is a convex function where the second derivative of the function is always non-negative (Equation (20)).
From the objective function in Equation (1), we determine three decision planning as follows: (1) the number of open depots and assigned customers at each open depot.
(2) the routing of vehicles delivering demands from the depot to the customers.
(3) the optimal number of production quantity of each open depot.
We solve the problem iteratively where the allocation of customers to each depot has been assigned at the first stage with the possible number of open depots. At this stage, the vehicle capacity limit constraint has been relaxed. At the second stage, the process of minimizing location and routing with considering vehicle capacity constraint is done. To obtain the cost of inventory as stated in the objective function, the Q i is calculated by the Equation (19) where the formulation is depending on the decision variable of y ij and the total demand and returned from each customer at each depot. In other words, the inventory cost will be calculated when the decision variable for location and routing are found throughout the process of improvisation at each iteration. The process of minimizing the overall cost depends on the total operating depots, delivery routing, and inventory planning.

Harmony Search
Harmony Search (HS) is a population-based metaheuristic algorithm that mimics the music improvisation of a group orchestra. The musicians will improvise their harmonies with these three options: (1) select any pitch in the recorded harmony memory; (2) select and fix any previous pitch in memory; or (3) search a new music harmonisation within the range of the pitch. The process of obtaining the optimal solution is similar to this efficient search for a perfect state of harmony [32,33]. There are five main steps in the proposed HS: Appl. Sci. 2020, 10, 6625 8 of 22 Step 1: Initialisation In HS, several solutions are generated to produce an initial population. The solutions can be generated either randomly or by the heuristics. A random initial population called harmony memory (HM) is created and sorted according to their fitness values. The number of solutions in the HM is the harmony memory size (HMS). An HM can be described in the following matrix: where, . . , n and j = 1, 2, . . . , HMS, f x j = fitness function of x j for j = 1, 2, . . . , HMS.
Step 2: Parameter Setting Harmony memory considering rate (HMCR) and pitch adjustment rate (PAR) are two main parameters used in the HS algorithm. An appropriate value of HMCR will lead to the choices of good solutions as an element of new solutions. The convergence will be slow if the value of HMCR is too high. Some potential good solutions will not be well explored. But if it is too low, only a few good solutions are selected for the next iteration. Therefore, to use memory effectively, the value of HMCR should be between 0.7 and 0.95 [34]. The speed of convergence in a standard HS algorithm may be slow if the value of PAR is too small. Besides, a high value of PAR is also not recommended. It can cause solutions to be stuck around a few potential optimal and easily get trapped in local optima. Due to this, [34] stated that the value of PAR is suggested to be between 0.1 and 0.5. However, this range may be suitable and valid for some problems only.
To balance the exploration and exploitation, instead of using a fixed value, in this paper, the proposed HS uses a dynamic value of HMCR it and PAR it introduced by [35]. The reason for reducing the HMCR it slowly is to increase the probability of exploring more solution space, not in the HM. Hence, the global optimal can be attained [36]. The dynamic value of HMCR it and PAR it can avoid the solutions getting trapped in the local optimum quickly as well. With a slight modification, the proposed HMCR it and PAR it are reduced gradually when there is no improvement found in the solution. In this paper, the range values of [HMCR min , HMCR max ] and [PAR min , PAR max ] are set to be [0.7, 0.95] and [0.3, 0.9], respectively. Equations (23) and (24) below are the formulation of dynamic HMCR it and PAR it : where, it = current iteration, MaxIt = maximum number of iterations, HMCR max = maximum value of HMCR, HMCR min = minimum value of HMCR, PAR max = maximum value of PAR, PAR min = minimum value of PAR.
Step 3: Improvisation In the proposed HS, to increase the speed of convergence, multi solutions are generated at each iteration and called HM new . Each newly generated solution is taken from the HM randomly with the probability of HMCR, otherwise, it will be generated randomly within the range of the harmony vector. To further enhance the intensification of the method, the proposed HS implements the local optimization in three different ways: within a depot, between two depots, and within the vehicle route. Hence, five multi-neighbourhood local search techniques are proposed: (1) Swap: exchange the position of two random customers within the same route or between two different routes. Swapping can be done within a depot or between the depots. (2) Insertion: insert a customer in between two other random customers within the same depot.  Step 4: Update HM To update the solutions in HM, the HM new created at each iteration is combined with current HM. The solutions are sorted according to the value of the objective function. The best solutions with the size of HMS will be kept for the next iteration.

Repeat
Step 3 and 4. The proposed algorithm will be terminated when either the MaxIt is reached or t consecutive non-improving iterations are performed.

Hybrid Harmony Search-Simulated Annealing
To prevent the possibility of the HM population been trapped into local optima due to premature convergency, the proposed HS algorithm is further enhanced by hybridizing a SA with the implementation of probabilistic acceptance rule based on the Boltzmann factor, e − ∆ T it during the improvisation procedure of the HS. In SA, the worst solution will be accepted within a probability and allows the search to proceed to its neighbourhood (Algorithm 1). By preserving the worst solution in a population, it provides practical randomness into the search to avoid the local extreme points. The probability of acceptance in SA is depending on the value of annealing temperature and the difference in the value of an objective function. At each iteration, the annealing temperature will be reduced according to the value of the cooling rate, which is a function of temperature. The higher value in the cooling rate, the slower the temperature reduction [37]. The procedures of probabilistic acceptance rule in SA described in Algorithm 1 are of a standard SA that has been adapted into the proposed hybrid HS-SA as shown in Algorithm 2. if f (x ) < f (x) 8: x ← x 9: else 10:
1: begin 2: define the parameter setting: HMCR max , HMCR min , PAR max , PAR min , HMS, T 0 , α, MaxIt, t 3: generate solution vectors and calculate fitness function for the initial population of HM. 4: while ((it < MaxIt) OR (non_improving < t)) do 5: while (no. of harmony vectors < HM new ) do 6: do location-allocation and route allocation 7: calculate HMCR it and PAR it 8: if rand < HMCR it 9: select solution vectors from the HM randomly 10: if rand < PAR it 11: implement one of the multi-neighbourhood search techniques randomly to create new neighbourhood solution, x 12: if f (x ) < f (x) 13: x ← x 14: else 15:

Computational Experiments and Discussion
In this section, computational experiments are performed to illustrate the performance of the proposed hybrid HS-SA compared with other approaches in the literature. The proposed algorithm is coded in MATLAB software R2017b on a laptop computer with 1.60 GHz Intel ® Core™ i5-4200U CPU with 8 GB of RAM. The parameters configuration of the proposed algorithm is provided in Table 1. The test problem instances are three popular benchmark datasets of Perl, Gaskell, and Christofides with additional simulated data of returned products and production rate. The additional data were generated randomly by using a uniform distribution. The production rate is set to be greater than the total demand and returns, P > j∈J D j + r j + s j while the returned products, R j are generated by the uniform distribution of R j ∼ U 0, D j . We set 70% of the returned products as the non-defects (r j ) items while the remaining are the defect items (s j ). The characteristics of the dataset are given in Table 2. To solve the LIRP problem, the location and allocation of the depots and customers are solved first. In the location-allocation problem, each of the customers is initially assigned to the nearest depots based on the Euclidean distance formulation. The customers will be moved around to the possible depots during the process of improvisation. The vehicle capacity limit constraint has been relaxed to minimize the number of open depots and it is assumed that only one vehicle is being used at each depot. However, the depot capacity limit should not be violated during the process. Then, the allocation of customers at each depot needs to be sequenced and divided into vehicles according to the vehicle capacity limit. This process is performed among the open depots only. The process of local optimization is performed within the open depots. Both depot capacity limits and vehicle capacity limits should not be violated. Finally, the inventory part is included during the process of minimizing the cost.
Since there is no study on the LIRP that utilizes the EPQ model for defect and non-defect items in reverse logistics from the literature, in order to validate the performance of the proposed hybrid HS-SA algorithm, the LRP will be solved first and compared with other heuristics and metaheuristics found in the literature that used the same benchmark datasets. Tables 3-5 present the comparative results of the three benchmark datasets on the proposed hybrid HS-SA with the heuristic method (HR) of [8], simulated annealing in [9] (SA-W) and [10] (SA-Z), hybrid tabu search and ant colony (TACO) from [11], particle swarm optimization (PSO) with different variants of particle swarm optimization which are: standard PSO, PSO with MPNS-GRASP (PMG), PSO with MPNS-GRASP and ENS (PMGE) and hybrid PSO (HLRP) in [12] and honey bee mating optimization (HBMO) in [13]. A standard HS (SHS) and the proposed HS (PHS) from [7] are also included for comparison in all experiments. For each problem instances, the proposed algorithm as well as SHS and PHS are performed for 10 independent runs. In each table, note that the numerical results of the selected algorithms are extracted from the original papers found in the literature and the best result of each problem instances is highlighted in bold.     As compared to the SHS and PHS in [7], the hybrid HS-SA performed significantly better for all instances in Perl, Gaskell, and Christofides, except Perl 1. However, for Perl 3, results in SA-Z [10] and TACO [11] are slightly better than the hybrid HS-SA algorithm. All algorithms managed to find the optimal solution in Perl 1 since the number of customers and depots are small. The SHS did not perform well for medium and large size. This indicates that the SHS should be improved and modified to get better solutions [7]. From the numerical experiments in LRP, hybrid HS-SA has successfully solved the problem with comparable results. Therefore, the proposed algorithm will be used to solve the LIRP that utilizes the EPQ model for defect and non-defect items in reverse logistics.
The computational results for all test instances in LIRP are shown in Tables 6-8   As compared to the SHS and PHS in [7], the hybrid HS-SA performed significantly better for all instances in Perl, Gaskell, and Christofides, except Perl 1. However, for Perl 3, results in SA-Z [10] and TACO [11] are slightly better than the hybrid HS-SA algorithm. All algorithms managed to find the optimal solution in Perl 1 since the number of customers and depots are small. The SHS did not perform well for medium and large size. This indicates that the SHS should be improved and modified to get better solutions [7]. From the numerical experiments in LRP, hybrid HS-SA has successfully solved the problem with comparable results. Therefore, the proposed algorithm will be used to solve the LIRP that utilizes the EPQ model for defect and non-defect items in reverse logistics. As compared to the SHS and PHS in [7], the hybrid HS-SA performed significantly better for all instances in Perl, Gaskell, and Christofides, except Perl 1. However, for Perl 3, results in SA-Z [10] and TACO [11] are slightly better than the hybrid HS-SA algorithm. All algorithms managed to find the optimal solution in Perl 1 since the number of customers and depots are small. The SHS did not perform well for medium and large size. This indicates that the SHS should be improved and modified to get better solutions [7]. From the numerical experiments in LRP, hybrid HS-SA has successfully solved the problem with comparable results. Therefore, the proposed algorithm will be used to solve the LIRP that utilizes the EPQ model for defect and non-defect items in reverse logistics.
The computational results for all test instances in LIRP are shown in Tables 6-8. Since there is no benchmark for the EPQ model in LIRP, the comparison is assessed between the hybrid HS-SA algorithm with SHS and PHS only. All algorithms are performed for 10 independent runs and the mean, standard deviation, coefficient of variation (C.V.), and the best result of each algorithm are reported in column 3-6 respectively. The Wilcoxon signed-rank test with α = 0.05 is performed to test the significance of the results. As the Wilcoxon signed-rank test does not assume normality in the data, it can be used when this assumption has been violated and the use of the dependent t-test is inappropriate. The p-value of the test is reported in the last column of the tables. The best result of each problem instances is highlighted in bold. In addition, the best result of the all problem instances are summarized as a bar chart in Figure 6.  The Wilcoxon signed-rank test elicits a statistical significant difference in the performance of finding the minimum cost of the LIRP between the hybrid HS-SA and the SHS and PHS across all the datasets (p-value < 0.005) based on the difference in the median cost. These tests are used to show that the results are significant, supporting that the proposed algorithm will produce results that are different from the results of the SHS and PHS. Apart from the Perl 1 instance where the proposed algorithm produced the same results as SHS and PHS, all other datasets showed p-values of less than 0.05. This means the results generated by the hybrid HS-SA is significant and different from the values from the SHA and PHS. So, the results from the proposed algorithm is significant from the other works, but is it better? From Tables 6-8, the hybrid HS-SA outperforms the solutions in SHS and PHS for all cases of datasets of Perl, Gaskell, and Christofides. For this, we use the C.V. values. The lower the value of the C.V., the more precise the estimate. As reported in Tables 6-8, all the C.V. values are low which supports the high predictiveness of the proposed algorithm.
Tables 9-11 present the detailed results of all problem instances produced by the proposed hybrid HS-SA for LIRP. The total cost for each problem instance is highlighted in bold. These tables served as the benchmark for the LIRP. For instance, Gaskell 1 required two depots with a fixed cost of 50 each to serve all the 21 customers with four vehicles (see , Table 10). In depot 1, with the inventory cost of 573.05, two vehicles are deployed to serve the customers (19→21→20→17) and (18→15→12→14→16) with costs of 59.45 and 86.90 respectively. While in depot 2, with the inventory cost of 577.74, two vehicles are used to serve the customers (9→7→5→2→1→6) and (8→3→4→11→13→10) with costs of 83.01 and 95.55 respectively. In this result, the best solution obtained suggest to close Depot 3-5. The overall minimum cost for Gaskell 1 is 1575.7. Figure 7 illustrate the topological layout of the Gaskell 1 produced by the proposed hybrid HS-SA for LIRP. The Wilcoxon signed-rank test elicits a statistical significant difference in the performance of finding the minimum cost of the LIRP between the hybrid HS-SA and the SHS and PHS across all the datasets (p-value < 0.005) based on the difference in the median cost. These tests are used to show that the results are significant, supporting that the proposed algorithm will produce results that are different from the results of the SHS and PHS. Apart from the Perl 1 instance where the proposed algorithm produced the same results as SHS and PHS, all other datasets showed p-values of less than 0.05. This means the results generated by the hybrid HS-SA is significant and different from the values from the SHA and PHS. So, the results from the proposed algorithm is significant from the other works, but is it better? From Tables 6-8, the hybrid HS-SA outperforms the solutions in SHS and PHS for all cases of datasets of Perl, Gaskell, and Christofides. For this, we use the C.V. values. The lower the value of the C.V., the more precise the estimate. As reported in Tables 6-8, all the C.V. values are low which supports the high predictiveness of the proposed algorithm.
Tables 9-11 present the detailed results of all problem instances produced by the proposed hybrid HS-SA for LIRP. The total cost for each problem instance is highlighted in bold. These tables served as the benchmark for the LIRP. For instance, Gaskell 1 required two depots with a fixed cost of 50 each to serve all the 21 customers with four vehicles (see , Table 10). In depot 1, with the inventory cost of 573.05, two vehicles are deployed to serve the customers (19→21→20→17) and (18→15→12→14→16) with costs of 59.45 and 86.90 respectively. While in depot 2, with the inventory cost of 577.74, two vehicles are used to serve the customers (9→7→5→2→1→6) and (8→3→4→11→13→10) with costs of 83.01 and 95.55 respectively. In this result, the best solution obtained suggest to close Depot 3-5. The overall minimum cost for Gaskell 1 is 1575.7. Figure 7 illustrate the topological layout of the Gaskell 1 produced by the proposed hybrid HS-SA for LIRP.

Conclusions and Future Work
The integration of facility location, inventory planning, and vehicle routing is one of the most challenging problems in the supply chain network design. In this study, the Economic Production Quantity (EPQ) model of reverse logistics that considers the returned products from the customers is included during the process of optimization. Since the problem is NP-hard, a hybrid populationbased metaheuristic called hybrid Harmony Search-Simulated Annealing (HS-SA) algorithm is proposed. Computational experiments on three sets of popular benchmark instances show that the proposed hybrid HS-SA algorithm is successful in finding better solutions as compared to a standard HS, as well as other approaches from the literature of Location-Routing Problem (LRP) and Location-Inventory-Routing Problem (LIRP). For future work, the HS can be hybridized with other metaheuristic methods such as differential evolution, genetic algorithm, and tabu search to further enhanced the efficiency of the approach. Environmental issues such as the CO2 emission and the use of electric vehicles can also be integrated into the model to be solved.

Conclusions and Future Work
The integration of facility location, inventory planning, and vehicle routing is one of the most challenging problems in the supply chain network design. In this study, the Economic Production Quantity (EPQ) model of reverse logistics that considers the returned products from the customers is included during the process of optimization. Since the problem is NP-hard, a hybrid population-based metaheuristic called hybrid Harmony Search-Simulated Annealing (HS-SA) algorithm is proposed. Computational experiments on three sets of popular benchmark instances show that the proposed hybrid HS-SA algorithm is successful in finding better solutions as compared to a standard HS, as well as other approaches from the literature of Location-Routing Problem (LRP) and Location-Inventory-Routing Problem (LIRP). For future work, the HS can be hybridized with other metaheuristic methods such as differential evolution, genetic algorithm, and tabu search to further enhanced the efficiency of the approach. Environmental issues such as the CO 2 emission and the use of electric vehicles can also be integrated into the model to be solved.