Abstract
Facility location, inventory management, and vehicle routing are three important decisions in supply chain management, and location-inventory-routing problems consider them jointly to improve the performance and efficiency of today’s supply chain networks. In this paper, we study a location-inventory-routing problem to minimize the total cost in a closed-loop supply chain that has forward and reverse logistics flows. First, we formulate this problem as a nonlinear integer programming model to optimize facility location, inventory control, and vehicle routing decisions simultaneously in such a system. Second, we develop a novel heuristic approach that incorporates simulated annealing into adaptive genetic algorithm to solve the model efficiently. Last, numerical analysis is presented to validate our solution approach, and it also provides meaningful managerial insight into how to improve the closed-loop supply chain under study.
1. Introduction
Supply chain management is critical for many business organizations to gain advantage in a competitive environment, and its impact has increased steadily in the past decades [1]. Although most practices in supply chain management were focused on forward logistics in early days, reverse logistics flows that are caused by consumer returns have gained a lot of attention recently and hence are considered by many firms to improve their business. According to a National Retail Federation report, the total merchandise returns accounted for $260.5 billion and $28.3 billion for the loss of the U.S. retailers and Canadian retailers in 2015, respectively [2]. Consumer returns also have a significant impact on e-commerce, and it is shown that at least 30% of all the products ordered online are returned as compared to 8.89% in traditional offline stores [3]. Particularly, for fashion products such as fashion apparel, the return rate can be as high as 75% [4]. Therefore, consumer returns represent a growing financial and operational concern for many firms in different industries, and they also have a significant impact on their supply chains.
Closed-loop supply chains (CLSCs) [5, 6] are an emerging topic in supply chain management because of the growing concern about consumer returns and environmental sustainability. Unlike traditional supply chains that only consider forward logistics flows directed from manufacturers to consumers, CLSCs also consist of reverse flows of new or used products that are directed from consumers to manufacturers. In practice, business managers need to make many strategic, tactical, and operational decisions such as facility locations, inventory control, and vehicle routing decisions to improve the efficiency and sustainability of their supply chains. In this paper, we study a location-inventory-routing problem (LIRP) that integrates those three decisions in a multiechelon closed-loop supply chain network for a manufacturer. This network comprises a manufacturing factory, multiple hybrid distribution-collection centers (HDCCs), and several retailers, where HDCCs will operate as warehouses and collection centers in the forward and reverse flows, respectively. From a practical perspective, the research questions that motivate this study are summarized as follows: (1)For a manufacturer, how to decide HDCC locations in a supply chain network when forward and reverse logistics flows are both considered, and how to use those HDCCs to fulfill the demands and collect returns from retailers?(2)What is the optimal stock replenishment policy for those HDCCs?(3)How to optimize vehicle routes in the forward and reverse flows when retailers are served by those HDCCs?
The rest of this paper is organized as follows: Section 2 reviews the related literature. Section 3 describes the research problem under study and formulates it as a nonlinear integer programming model. Section 4 proposes an adaptive hybrid simulated annealing genetic algorithm (AHSAGA) to solve the model efficiently. Section 4.1 presents the numerical study and computational results. Section 5 concludes this paper and provides directions for future research.
2. Literature Review
The location-inventory-routing problem (LIRP) comprises three subproblems: facility location, inventory control, and vehicle routing. Since they are highly correlated in the real-world business, many research efforts have been conducted to study those problems jointly. The location inventory problem (LIP) is the integrated form of the first two problems, and they are first proposed by Daskin et al. [7] and Shen et al. [8]. LIPs have been extended with many business scenarios such as lateral transshipment [9], perishable products [10], correlated demands [11], disruption risk [12], and inventory control strategies [13], and most of those works are reviewed by Farahani et al. [14]. Recently, LIPs are also studied by incorporating CLSCs. For example, Diabat et al. [15] study a LIP by considering spare parts in a closed-loop system, Guo et al. [16] study location-inventory decisions for closed-loop supply chain management with secondary market consideration, and Li et al. [17] present an important and meaningful work by studying LIP and CLSC with third-party logistics (3PL) because it is a fundamental logistics strategy that has been adopted by many firms in practice. The location routing problems (LRPs) integrate facility location and vehicle routing problems but ignore inventory management decisions. Min et al. [18] and Nagy and Salhi [19] review the research works related to LRPs in early days, and Schneider and Drexl [20] examine the most recent works that are published in the literature since the survey by Nagy and Salhi [19].
LIRPs incorporate all three decisions above, and hence they are a more comprehensive form. In early days, Shen and Qi [21] develop a location-allocation model which approximates routing costs according to the locations of opened depots, and then Javid and Azad [22] study such a problem without any approximation. Moreover, LIRPs are extensively studied under many practical settings such as perishable products [23, 24], deterministic or stochastic demand [25–27], and disruption risks [28]. Closed-loop supply chains (CLSCs) have attracted considerable attention from researchers and practitioners because of the significant impact of consumer returns, and it is emergent to study LIRPs under a CLSC setting. For example, Li et al. [29] study a LIRP by considering returns in an electronic supply chain environment, and Deng et al. [30] develop and solve a model when returned products can be either defective or nondefective. From the perspective of sustainability, Zhalechian et al. [31] design a closed-loop system with location routing inventory decisions under mixed uncertainty.
In this paper, a nonlinear integer program model is formulated to study a LIRP in a closed-loop logistics system by considering many real-world business scenarios such as vehicle capacity and the disposal of different types of returned products. To solve this model efficiently, we develop a novel solution approach that extends the power of the adaptive genetic algorithm by incorporating simulated annealing, and numerical study shows that it is more powerful and efficient than other similar heuristics in the literature.
3. The Model
3.1. Problem Description
In this paper, we study a closed-loop supply chain network that comprises a manufacturing factory, multiple hybrid distribution-collection centers (HDCCs), and several retailers. This network can be represented by a directed graph in which vertices are the factory, HDCCs, and retailers, and the edges can be directed from the factory to retailers via HDCCs, or vice versa. More specifically, in the forward flow, new products are first shipped from the factory to HDCCs and then from HDCCs to retailers by vehicles on certain routes. In the reverse flow, returned products are sent from retailers to HDCCs for inspection first. A returned product will be disposed immediately at a HDCC if it cannot be refurbished. Otherwise, it will be sent from HDCCs to the factory for repair. In this system, HDCCs operate as warehouses and return collection centers on working days in the forward and reverse flows, respectively. Vehicles are used to deliver new products from HDCCs to retailers as well as to collect returned products from retailers to HDCCs, and a vehicle must return to the same HDCC after it visits all retailers on a route. Figure 1 illustrates the closed-loop supply chain network under study.
For simplicity, we consider a single type of products and vehicles and assume that a retailer will be assigned to a same HDCC in the forward and reverse flows. Given the locations of the factory and retailers, HDCCs will be built at selected locations, and a HDCC will order new products from the factory and serve at least one retailer in the forward and reverse flows. To minimize the total cost in this system, the following decisions will be optimized: (1)HDCC location and retailer assignment: selecting locations to build HDCCs and assigning retailers to those HDCCs(2)Inventory replenishment: deciding the optimal order frequency and quantity for each HDCC(3)Vehicle routing: designing circular vehicle routes starting from and ended by each HDCC
3.2. Objective Function
In the closed-loop supply chain under study, the total cost is composed by the following: (1) location cost which is the fixed cost of building and operating HDCCs; (2) working inventory cost including order, holding, and shipping costs; (3) routing cost between HDCCs and retailers; (4) return cost. The individual costs per year are calculated as follows: (1)Location cost: (2)Working inventory cost The working inventory cost comprises three individual terms. The first term is the order cost that is incurred when placing orders to the factory at HDCCs, the second term is the holding cost of new products in inventory, and the third term is the shipping cost of new products from the factory to HDCCs. Similar to [22], we adopt a (Q, r) inventory model with type I service to manage inventories at HDCCs, and the holding cost is adapted from the standard form in the economic order quantity (EOQ) model. Obviously, the order frequency and quantity at a HDCC is determined by the expected demands of the retailers that are served by the HDCC. Therefore, the individual terms of the working inventory cost can be written as follows: (i)Order cost: (ii)Holding cost of new products: (iii)Shipping cost from the factory to HDCCs: Consequently, the total working inventory cost per year is given as follows: (3)Vehicle routing cost (i)Forward logistics: (ii)Reverse logistics: Therefore, the total annual routing cost is given as follows: (4)Return cost (i)Inspection cost: (ii)Disposal cost at HDCCs: (iii)Cost of refurbish returned products at the factory: (iv)Shipment cost from HDCCs to the factory: (v)Holding cost of returned products:
For simplicity, we assume that the holding cost of a returned product is independent of how long it stays in inventory. Therefore, the total return cost per year is given as follows:
According to the individual costs above, the total annual cost in the CLSC is calculated as follows:
Therefore, the location-inventory-routing problem under study can be formulated as follows:
The constraints of this model are explained as follows. Constraint (5) means that at least one HDCC will be built. Constraint (6) means that a vehicle can be assigned to at most one HDCC. Constraint (7) means that a retailer will be served by exactly one HDCC. Constraint (8) means that a retailer will be placed in exactly one route. Constraint (9) is the vehicle capacity constraint which means that the total demand in a route cannot exceed the capacity of a vehicle. Constraint (10) means that vehicles can be assigned to a HDCC only if the HDCC has been built. Constraint (11) means that at least one vehicle will be assigned to a HDCC. Constraint (12) means that retailers can be served by a HDCC only if it has been built. Constraint (13) means that at least one retailer will be served by a HDCC. Constraint (14) means that for each HDCC, the number of vehicles or routes is less than the number of retailers. Constraints (15) and (16) enforce the relationship between and , and they mean that a retailer can be placed in a route only if the route exists as well as that at least one retailer will be included in a route. Constraints (17) and (18) enforce the relationship between and , and they mean that a retailer will be placed in exactly one route which belongs to a HDCC if and only if it is served by the same HDCC. Constraints (19) and (20) enforce the relationship between and for retailers, and they mean that a retailer cannot have neighbor locations in a route if it is not in the route as well as that no closed subloop will present in a route. Constraints (21), (22), (23), and (24) enforce the relationship between and for HDCCs, and they mean that exactly one HDCC will be placed in each route. Constraints (25) and (26) mean that in a route, a HDCC can be directed to/from at most one retailer. Constraint (27) is the flow conservation constraint, and it means that a vehicle must leave a retailer after it arrives at this retailer and hence the route is circular. Constraint (28) means that a route cannot be directed from a retailer to a HDCC if the retailer is not served by this HDCC. Constraint (29) guarantees that a route will be in one direction but not in two directions. Constraints (30), (31), (32), (33), and (34) specify that , , , , and are binary variables.
It is obvious that the objective function is convex with respect to . To calculate the optimal number of orders placed at HDCC annually, let , then we have
By substituting into (4), the objective function can be rewritten as
4. Solution Approach
Facility location and vehicle routing problems are NP-hard in general [22], and LIRPs can be more complex to solve because of the integration of those problems. In this paper, we propose a two-phase heuristic method that incorporates simulated annealing (SA) into adaptive genetic algorithm (AGA). More specifically, facility location and vehicle routing decisions will be encoded as chromosomes in AGA, and then the two decisions will be optimized by an evolution process. Thereafter, the optimal inventory replenishment decision will be determined accordingly.
GA is a popular search technique to solve optimization problems based on the principles of natural selection and genetics [32]. In practice, a GA process may converge prematurely or do not converge at all, both of which will lead to bad solutions, and hence adaptive coefficients are usually used to compensate those shortcomings. SA is a probabilistic method which was first proposed to find the global minimum of a cost function that may possess several local minima [33], and it has been widely used to solve many research problems. In this paper, we propose a novel heuristic algorithm, that is, adaptive hybrid simulated annealing genetic algorithm (AHSAGA), to solve the nonlinear integer programming model presented in Section 3. AHSAGA is an improved form of traditional AGAs by adopting the great local search capacity of SA, and our numerical experiments show that it is an effective approach in terms of both solution accuracy and time efficiency.
4.1. Basis
4.1.1. Encoding and Decoding
When GA is applied to solve an optimization problem, chromosomes are usually used to represent the candidate solutions to this problem, and they will evolve to better solutions iteratively. In this study, the solutions to the location and routing problems will be first encoded as chromosomes and then solved by AHSAGA. Once the location and routing problems are solved, inventory decisions can be easily optimized by solving (35).
The length of a chromosome is R + S, where and are the number of candidate HDCC locations and retailers, respectively. Let be the population size, then an initial population can be created by randomly choosing chromosomes that satisfy (5), (6), (7), (8), (9), (10), (11), (12), (13), (14), (15), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), (29), (30), (31), (32), (33), and (34). In a chromosome, a HDCC and its subsequent retailers comprise a circular route. Therefore, each chromosome will start with a candidate HDCC location. If a chromosome starts with a retailer, then the first allele will be replaced by a candidate HDCC location because the retailers before the first candidate HDCC location in the chromosome will not be assigned to any HDCC. Moreover, if there are consecutive candidate HDCC locations in a chromosome, then a HDCC will be built at the location represented by the last allele of this string.
4.1.2. Fitness Function and Selection Method
Once a population is created, chromosomes or candidate solutions will be evaluated by their fitness to decide whether they will be kept in its offspring population. In this study, the fitness of an individual is measured as follows: where is the objective function given by (36).
In AHSAGA, roulette-wheel selection is adopted to select and copy solutions with higher fitness values into new populations. Let be the population size, then chromosome will be reproduced in the next generation if it satisfies the equation below: where is the fitness value of chromosome , is a random number that follows the uniform distribution.
4.1.3. Crossover Operator
In general, a GA process will start with an initial population that is generated randomly, and the fitness of solutions will be improved iteratively by applying selection, crossover, mutation, and replacement operators. In AHSAGA, crossover operator will be applied in an iteration by the following three steps to recombine individuals for a better offspring: (1)Choose two parents from a population randomly and decide two crossover points arbitrarily.(2)Generate two intermediate chromosomes by moving all the alleles positioned between the crossover lines in a parent to the beginning of the other.(3)In each intermediate chromosome, remove the same alleles which appear in the string moved from the other parent.
An example of this procedure is illustrated in Figure 2.
Usually, fitness values of chromosomes will be significantly different at the beginning of a GA process, and hence crossover is greatly beneficial to speed up the evolution. In AHSAGA, the probability of crossover is given by (39), which is similar to [34] in spirit. where is the initial probability of crossover, and are the two adaptive coefficients, , , and are the maximal, average, and minimal fitness values in a population, respectively.
4.1.4. Mutation Operator
Crossover operators cannot work effectively if individuals have similar fitness values in a population. For example, in some cases, new chromosomes cannot be generated by crossover if two parents have the same allele at a given gene. To solve this problem, mutation is designed to add diversity to the population and make it possible to explore the entire search space [32]. AHSAGA uses an inverse function as the mutation operator to select two points in a parent chromosome randomly and invert the order of the alleles between the two points. For example, if (7 8 5 |3 6 9 1| 4 10 2) is a parent, then the first and second split points are located after the third and seventh genes, respectively, and hence the offspring will be (7 8 5|1 9 6 3|4 10 2) after inverse. Mutation will occur randomly, and the probability of mutation is given as follows: where is the initial probability of mutation, and are the two adaptive coefficients, , , and are the maximal, average, and minimal fitness values in a population, respectively.
If a chromosome starts with a retailer, then the initial allele will be inverted with the first allele that represents a candidate HDCC location.
4.1.5. Simulated Annealing and Individual Replacement
In AGAs, individuals will be replaced by new ones for evolution. AHSAGA adopts SA as the steady-state technique [32], and the probability that a chromosome will be replaced is given as follows: where and are the fitness values of the new and old individuals, respectively, and is the temperature given as follows: where is the temperature at time , and is the change rate of temperature.
4.2. Algorithm
The pseudocode of AHSAGA is shown in Algorithm 1, and the steps in this algorithm are briefly explained as follows:
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Step 1. Initialize parameters such as population size pop_size, iteration number , crossover factors , , , mutation factors , , , initial temperature , and cooling factor .
Step 2. Create an initial population randomly.
Step 3. Generate an offspring population by applying selection, crossover, and mutation operators.
Step 4. Calculate the fitness values of the individuals in a new population, and identify those with the maximal and minimal fitness values.
Step 5. Check whether the maximal fitness value in an offspring population is greater than that in its parent population. If yes, go to Step 6. Otherwise, apply SA to decide whether the best solution in the parent population will be introduced into the offspring population, then update temperature in SA.
Step 6. Check whether the termination condition is satisfied. If yes, return the chromosome with the maximal fitness value. Otherwise, go to Step 3.
5. Numerical Study
In this study, AHSAGA is implemented by Matlab R2014a and all numerical experiments are conducted on a workstation equipped with an Intel Core i7-4790 CPU at 3.60 GHz and 8.0 GB of RAM under Windows 7.
To validate its performance, AHSAGA has been tested on five data sets that are adapted from LRP files provided by the University of Aveiro [35] for locations, fixed costs, and demands. For example, the input data adapted from the Gaskell67-21 × 5 files for HDCCs and retailers are shown in Tables 1 and 2, respectively. All other parameters are provided in Table 3.
5.1. Sensitivity Analysis
Since the performance of AHSAGA can be affected significantly by its parameters, a sensitivity analysis is conducted on the parameters shown in Table 4 to identify the optimal setting in this study. To eliminate the excessive number of combinations, other parameters will be set to their median values when a parameter is tested. The numerical results are shown in Figure 3, where red and blue lines represent the mean objective values and average computational times, respectively. Figure 3 shows that solution accuracy and computational times are both affected by those parameters, and the best result can be archived when , , , , , , , , , and . The optimal setting is presented under the “Experiment setting” column in Table 4, and it will be used in the subsequent experiments.
From Figure 3, we can see that AHSAGA can be affected by those parameters in the following way: (1) and When or increases, the optimal value and computational time will decrease and increase, respectively. This indicates that more iterations or greater population diversity will lead to a better optimal solution with an additional time cost.(2) and The optimal value will always decrease when increases, but it will not always decrease when increases. This indicates that a larger probability of mutation is always helpful to get a better optimal solution, but the probability of crossover should be moderately large. Moreover, we can see that when and increase, the computational time will increase and decrease, respectively. This indicates that a larger and will slow down and speed up the convergence of AHSAGA, respectively.(3) and When increases, the optimal solution can be always improved with an additional time cost. But should have a moderate value to achieve the best performance in terms of optimal solution and computational time.(4) and needs to be moderate to get the best optimal solution, but a larger can always improve the convergence and reduce the computational time. When increases, the optimal solution can always be improved with an additional time cost.(5) and A higher initial temperature in SA can always reduce computational time, but it is not always helpful to get a better optimal solution. However, a larger cooling factor will always improve the optimal solution and computational time.
5.2. Illustrative Example
In this section, Gaskell67-21 × 5 files [35] are used as an example to show the application and performance of AHSAGA. To get started, an initial chromosome is generated randomly, which is {5, 9, 10, 8, 26, 7, 14, 21, 15, 12, 13, 25, 1, 2, 11, 16, 3, 4, 18, 20, 6, 17, 24, 22, 19, 23}. According to the encoding-decoding scheme in Section 4.1, HDCC locations and vehicle routes can be decoded as that in Table 5, and then the corresponding optimal number of orders per year can be calculated by (36). When the algorithm is executed iteratively, objective values and adaptive probabilities change monotonically as shown in Figures 4 and 5, respectively.
(a) Crossover probability
(b) Mutation probability
To validate its performance, AHSAGA is compared with other two heuristics in the literature, which are adaptive annealing genetic algorithm (IAGA) [34] and hybrid genetic simulated annealing algorithm (HGSAA) [29]. To avoid any bias, each algorithm is replicated 50 times by using a same data set, and the mean objective values and computational times are compared. To illustrate the stochastic nature of the three algorithms, Figure 6 shows the 50 objective values from AHSAGA, HGSAA, and IAGA in a descending order by using the data adapted from Gaskell67-21 × 5 files, and Table 6 presents a more thorough comparison between the three algorithms on this data set, which shows that AHSAGA is more effective than HGSAA and IAGA from the perspectives of robustness, solution quality, and time efficiency.
5.3. Performance Comparison
The section presents a comprehensive comparison between AHSAGA, HGSAA, and IAGA on three types of problems by the number of retailers. More specifically, the number of retailers is less than 50 in small-size problems, between 50 and 100 in medium-size problems, and more than 100 in large-size problems. The numerical results on small-size, medium-size, and large-size problems are shown in Tables 7–9, respectively, from which we can make the following conclusions: (i)The mean objective values from AHSAGA are significantly lower than those from IAGA and HGSAA for most problems. This indicates that AHSAGA has a great capability to search global optimums and hence can provide better solutions.(ii)AHSAGA takes less computational times and convergence generations to find the optimal solution than IAGA and HGSAA for all problems. This indicates that AHSAGA is the most efficient approach.(iii)The variation of the optimal values from AHSAGA, which is measured by the coefficient of variation, is lower than that from IAGA and HGSAA for all problems. This indicates that AHSAGA is more robust and consistent than the other two algorithms.
6. Conclusions and Future Study
Closed-loop supply chains are an emerging and important topic due to the tremendous economic and environmental impact of consumer returns. In this paper, we study a location-inventory-routing problem in a closed-loop supply chain by formulating it as a nonlinear integer programming model. Since the problem is NP-hard, we also design a novel adaptive genetic algorithm by incorporating simulated annealing to solve this model efficiently. To make this study more practical, many real-world business scenarios such as vehicle capacity and the disposal of different types of returned products are also considered and modeled precisely.
This study can be extended in several directions in the future: first, this problem will be more practical and flexible if some assumptions are relaxed. For example, it will be flexible to allow a many-to-many relationship between vehicles and retailers, and it will also be more practical to relax the assumption that a retailer will be visited by a vehicle every working day. Second, since secondary markets have become an important channel to sell used products, it will be greatly beneficial to study LIRPs in a CLSC by considering those markets. Third, our model will be more valuable if it incorporates more business scenarios such as supply risk and multiple sourcing.
Sets
| : | set of candidate HDCC locations, where |
| : | set of vehicles, where |
| : | set of retailers, where |
| : | set of locations, which is the union of HDCCs and retailers (i.e., ). |
| : | fixed cost of building and operating a HDCC at location |
| : | shipping cost per unit of product between the factory and HDCC |
| : | vehicle capacity |
| : | daily demand of retailer |
| : | fixed cost per shipment from the factory to HDCC |
| : | fixed administrative and handling cost of placing an order to the factory from HDCC |
| : | disposal cost per unit of returned product which cannot be refurbished at HDCC |
| : | holding cost per unit of new product per year at HDCC |
| : | holding cost per unit of returned product at HDCC |
| : | fixed cost of repairing and repacking one unit of returned product at the factory |
| : | inspection cost per unit of returned product at HDCC |
| : | daily returns from retailer , where |
| : | distance from HDCC to retailer in a route (forward logistics) |
| : | distance from retailer to HDCC in a route (reverse logistics) |
| : | shipping cost per unit of product and distance |
| : | probability that a returned product cannot be refurbished |
| : | workdays per year (remark: similar to [21], we assume that a retailer will be visited by a vehicle every workday. Hence, is also the number of road trips for a vehicle per year). |
| Nr: | number of orders placed at HDCC r per year |
Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was supported by the National Natural Science Foundation of China under Grant nos. 71672074 and 71772075.