A novel hybrid genetic algorithm for the location routing problem with tight capacity constraints

Applied Soft Computing Journal 85 (2019) 105760

Contents lists available at ScienceDirect

Applied Soft Computing Journal journal homepage: www.elsevier.com/locate/asoc

A novel hybrid genetic algorithm for the location routing problem with tight capacity constraints ∗

Xue Yu a,b , Yuren Zhou a , , Xiao-Fang Liu a a b

School of Data and Computer Science, Sun Yat-sen University, 510006 Guangzhou, China Collaborative Innovation Center of High Performance Computing, Sun Yat-sen University, China

article

info

Article history: Received 25 May 2019 Received in revised form 28 July 2019 Accepted 3 September 2019 Available online 13 September 2019 Keywords: Location routing problem Genetic algorithm Capacity constraint

a b s t r a c t Location routing problem (LRP) is a popular and challenging topic in the field of logistic systems. LRP needs to address the depot location problem and vehicle routing problem at the same time. Till now, different LRP variants have been formulated to better meet realistic requirements. In this study, we focus on the capacitated LRP (CLRP) with tight capacity constraint on both depots and vehicles. To cope with the tight constraints, a hybrid genetic algorithm (HGA) is developed to search not only feasible solution space but also infeasible solution space. The proposed HGA combines the wide exploration capacity of GA, and the fast exploitation capacity of neighbourhood local search. To evolve GA for CLRP, solutions are represented by sets and sequences, and accordingly, a multi-sequence-based crossover is designed for the offspring generation. Moreover, a population management scheme is designed to facilitate GA’s evolution. Experiments are conducted on two benchmark sets and the results show that HGA is quite competitive with existing well-known CLRP algorithms on the classical instances. Furthermore, HGA is able to obtain quite a number of new best solutions on the real-life-like instances with tighter constraints. © 2019 Elsevier B.V. All rights reserved.

1. Introduction A reasonable logistic system design should be important to reduce the total supply chain cost and to improve the supply efficiency. Generally, two-level decisions must be tackled during logistic system design: (1) a long-term decision for the facility location problem (FLP) or depot location problem (DLP); (2) a short-term decision for the vehicle routing problem (VRP) to satisfy each customer’s demand. Originally, these two problems were addressed separately, but it is very likely to obtain only suboptimal solutions in this way [1]. Therefore, the research of location routing problem (LRP) has become more and more popular to simultaneously address DLP and VRP. For details, there are research surveys on LRP till 2007 [2] and 2014 [3], 2015 [4]. In the literature, research efforts have been devoted to different variants of LRP. This work mainly focuses on the capacitated LRP (CLRP) with capacity constraint on both depots and vehicles. The CLRP is difficult since both DLP and VRP are NP-hard [5], and a tight capacity constraint further increases the problem’s challenge. The CLRP can be represented by a weighted and directed graph G = (V , A), where the node set V is composed by a set I of customer nodes and a set J of candidate depots, i.e., V = I ∪ J. ∗ Corresponding author. E-mail address: [email protected] (Y. Zhou). https://doi.org/10.1016/j.asoc.2019.105760 1568-4946/© 2019 Elsevier B.V. All rights reserved.

Each customer i ∈ I = {1, . . . , n} has a demand di . Each depot j ∈ J = {1, . . . , m} has a fixed opening cost Cj and a capacity Qj . For each arc ⟨k, l⟩ ∈ A, the arc cost ck,l indicates the travelling cost from node k to node l. A fleet of homogeneous vehicles are available to serve all customers. Each vehicle has the same capacity Q and each vehicle route consumes a fixed cost VC. Each vehicle can be assigned to only one depot. That is, the vehicle can only start and end a route at the identical depot. Each customer should be served by one and only one vehicle route. The total load of a vehicle route, i.e., the cumulative demands of customers served by the route, should not excess the vehicle capacity. The total load of one depot, i.e., the cumulative load of all vehicle routes belonging to the depot, should not excess the depot’s capacity. The objective of the CLRP is to decide which depots to be opened, and to plan vehicle routes to serve all customers, so as to minimize the total cost, including depot cost, vehicle cost, and travelling cost. In Fig. 1, we exemplify the CLRP by a simple problem instance, whose vehicle cost and capacity are 15 and 10 respectively. Also, we provide a feasible solution for the CLRP instance, in which depot 1 and depot 3 are selected and the arrow lines make up vehicle routes (values on the arrow lines indicate the travelling costs or arc costs). In reality, it is worth noting that different depots often have very different capacities. In such case, the problem becomes more challenging since a different depot location may result in a completely different plan of vehicle routing, especially when

2

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760

Fig. 1. An exemplification of CLRP and a feasible solution.

the vehicle capacity is also tightly constrained. In this work, such kind of tightly-constrained situations are also considered when to solve the CLRP. Intuitively, in CLRP, the decisions of depot location and vehicle routing are inter-dependent. That is, the depot location results will affect the following vehicle routing, and in turn, the vehicle routings results reflect the quality of depot location results. In the literature, it is popular to solve the location problem and routing problem iteratively by nested methods [3]. That is, in each iteration, the location procedure selects a set of depots and then the routing procedure tries to find the optimal routes based on the selected depots. In such case, the vehicle routing of different sets of depots is often independent and it is hard to reuse the valuable historical vehicle routing information. To overcome this defect, it should be potential to develop the CLRP algorithms that optimize the location problem and routing problem simultaneously rather than iteratively. However, the simultaneous optimization may be even more difficult for the tight CLRP since vehicle routing plans of different depots tend to become harder to reuse as discussed above. In this study, we propose a hybrid genetic algorithm (HGA), in which solutions are represented by sets and sequences, and a multi-sequence-based crossover is designed to mate depot sets and the routes of different depots at the same time. Particularly, the GA’s evolution is giant-tour-based rather than route-based to deal with the routes’ high susceptibility to depot selection. Moreover, considering the fact that the optimal solutions of constrained problems are often close to the boundaries between feasible and infeasible solution space, the HGA is designed to evolve both feasible solutions and infeasible solutions, accompanied by a suitable population management scheme. Indeed, on the other hand, the permit of infeasible solutions helps to slacken the tight capacity constraints as well. In Fig. 2, we show how the proposed HGA approaches to the optimal solution step by step by exploring both feasible and infeasible solution space. 2. Literature review 2.1. Problem definition and applications As stated in [2], the phase ‘‘location routing problem’’ (LRP) may be misleading. Nowadays, a widely-accepted definition for LRP is ‘‘location planning with tour planning aspects taken into account’’ [2]. In this respect, the depot location is the master problem in LRP, and the vehicle routing is the subordinate problem in the service of depot location. Furthermore, the inter-relation or inter-dependency of depot location and vehicle routing must be considered when to solve the problem. Researches on LRP should be quite essential since locationrouting requirements can be observed in many practical applications. For instance, some works proposed location routing methods for commercial applications such as good distribution [6–8]

Fig. 2. An example of the evolutionary process in the proposed HGA.

and parcel delivery [9]; some works mainly devoted efforts to public services such as waste collection [10] and telecom network design [11]. 2.2. LRP variants Originally, most researchers studied un-capacitated LRP without capacity constraint on depots or vehicles. But since the first decade of the 21st century, more and more studies are observed on capacitated LRP (CLRP) [3]. It is as expected since capacity constraints should be quite commonly-seen in lots of practical applications. In addition to these classical LRPs, in recent years, quite a number of studies have been devoted to more complicated LRP variants that often accommodate more practical requirements. Recent surveys of LRP variants can be found in [3] and [4]. Herein, we just list several typical categories of complicated LRPs as follows. Multi-echelon LRPs. Multi-echelon distribution systems are quite common in reality, e.g., a system consisting of depots, intermediate warehouses and retailers, and researches on multiechelon LRPs are hence necessary. In 2012, Nguyen et al. [12] solved the two-echelon LRP by a GRASP reinforced by a learning process and path relinking. Afterwards, the same authors further implemented a multi-start ILS+PR, strengthened by various techniques such as greedy randomized heuristics and a tabu list for short term diversification [13]. The latter work can achieve better results at the cost of higher time consumption. LRPs with special objective(s). In some special realistic applications, e.g., hazardous material disposal and emergent disaster rescue, some other objectives must also be considered besides the financial cost. For example, Wichapa et al. [14] modelled the infectious waste disposal as a multi-objective LRP and then solved

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760

it using hybrid goal programming and hybrid genetic algorithm. K. Chang et al. [15] concentrated on the multi-objective LRP with uncertain data after disasters. Rayat et al. [16] focused on the location-inventory-routing problem with partial backordering, and further considered the failure costs related to disrupted distribution centres in addition to the ordinary financial costs. Stochastic and fuzzy LRPs. Uncertain data, e.g., uncertain customer demands, is not uncommon in practical vehicle routing problems. Usually, the uncertainties can be modelled by stochastic, probabilistic, or fuzzy functions. Nowadays, some studies have been dedicated for stochastic or probabilistic LRPs. For instance, N. Ghaffari-Nasab et al. [17] modelled and solved a bi-objective LRP with probabilistic travel times; also, an improved particle swarm optimization algorithm was proposed in [18] for LRP with stochastic demands. In addition, there are researches on fuzzy LRPs. For example, a simulated annealing variant was proposed in [19] for LRP with fuzzy demands; Zarandi et al. [20] focused on the multi-depot CLRP with fuzzy travel times and then proposed a simulation-embedded simulated annealing procedure to solve the problem. Periodic LRPs. In some realistic applications, customers need to be served one or more times over multi-period horizon, e.g., working days of a week. In such case, the decision of customers’ visiting patterns is also necessary. In 2008, Prodhon presented the concept of periodic LRP (PLRP) for the first time by combining periodic VRP and LRP [21]. Afterwards, in [22], the same author provided a MIP formulation for PLRP and proposed a hybrid evolutionary algorithm composed by the evolutionary local search (ELS) and the Randomized Extended Clarke and Wright Algorithm (RECWA). Despite the diversity of existing LRP variants, the classical LRPs always play an important role in problem definition and problem solving. For example, in [23], the two-echelon LRP is viewed as two CLRPs. Thus, further research on CLRP should still be meaningful. Especially, in the past, some efforts have been devoted to the CLRP with tighter constraints so as to better simulate real-life characteristics [24].

3

Heuristics & metaheuristics. Due to the incapability of exact methods for large-scale CLRPs, in the past years, researches of heuristics and metaheuristics have recently drawn more attentions. Especially, considering the problem difficulty of CLRP, it has become popular to combine several heuristic approaches rather than use a single one. Prins et al. [29] developed a GRASP (greedy randomized adaptive search procedure) complemented a learning process for depot location and a path relinking for post-optimization. The GRASP is iteration-based and in each iteration, a trial solution is first constructed by a greedy randomized heuristic and then it is improved by local search. Subsequently, still based on GRASP, Duhamel and Prins et al. [30] presented a GRASP × ELS, in which the local search in GRASP is replaced by an evolutionary local search (ELS). ELS mainly evolves the giant tours of customers and a greedy split procedure is adopted to transform a giant tour into a legal solution composed by routes. In each iteration of ELS, several offspring solutions will be generated by perturbation and local search. Also, in 2013, Ting and Chen [5] proposed a multi-ant colony optimization algorithm (MACO), which decomposes the CLRP into three subproblems, i.e., depot selection, customer assignment, and vehicle routing, and three ant colonies are accordingly utilized. At that moment, MACO was found to perform the best or second best on some benchmark instances. As pointed out in [3], the success of MACO may come from a scheme that looks quite promising: decomposing CLRP into parts but solving these parts in a global mechanism. In fact, our proposed HGA also complies with this scheme. Recently, a simple and effective hybrid GA has also been specifically developed for CLRP [31], in which solutions are represented by routes, and a special route copy crossover and two simple mutation operators (‘‘add’’ and ‘‘swap’’) are accordingly designed for GA’s evolution. Additionally, two local search procedures are developed for depot location and route improvement respectively.

2.3. Algorithms for CLRP In the literature, existing algorithms for CLRP can be roughly categorized into three types, i.e., exact methods, heuristics, and metaheuristics. Exact methods. Recently, several exact methods have been developed for CLRP. In 2011, J. Belenguer et al. [25] proposed a branch-and-cut method for CLRP, which is based on a zero– one linear model with customized families of valid inequalities. The method is found to be able to find optimal solutions for instances with up to 40–50 customers and 5 candidate depots. Soon afterwards, R. Baldacci et al. [26] developed another exact method based on a set-partition-like formulation of CLRP, which can solve more instances to optimality. They utilized the lower bounds obtained by different bounding procedures to decompose the problem into a limited set of multi-depot VRPs (MDVRPs). Thereafter, C. Contardo et al. [27] extended the work of [26] and presented an exact method based on cut-and-column generation. They also formulated the CLRP as a set-partitioning problem and then introduced five new families of inequalities. All possible subsets of depots are first enumerated by the two-index formulation, and for each subset, the corresponding MDVRP is then solved by the column generation. More recently, researchers have put more emphasis on exact methods for CLRP’s extensions, e.g., CLRP with time windows [28]. However, as the problem complexity increases, exact methods’ capability tends to become more limited. For example, several hours may be necessary for solving 40-customer instances in [28].

3. Proposed Hybrid Genetic Algorithm (HGA)

To solve the capacitated location routing problem (CLRP) with tight constraints, we propose a hybrid genetic algorithm (HGA) that allows both feasible and infeasible solutions, which are kept in two subpopulations pop_fs and pop_infs such that they two together make up the population pop. In this way, HGA would devote more efforts to the exploration of solutions near the boundaries between feasible and infeasible solutions, and we believe that these boundary solutions are more likely to spend less cost to satisfy the problem requirements while just not violating the problem constraints, i.e., better feasible solutions. To deal with both feasible and infeasible solutions, a population management mechanism is customized for HGA. Besides, for ease of population reproduction and evolution, each solution is represented by a set of selected depots and the specific customer sequence assigned to each selected depot. Accordingly, a multisequence-based crossover is designed to simultaneously consider the depot selection, customer assignment, and customer sequencing. Therein, each customer sequence is indeed a giant vehicle tour, and a splitting procedure is adopted to transform customer sequences into vehicle routes when to evaluate an individual.

4

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760

Fig. 3. An exemplification of the solution representation and transformation in HGA.

The overall framework of the proposed HGA is presented in Algorithm 1. The algorithm stops when the best-so-far solution best_sol has not been updated after a number of iterations, i.e., It NI iterations. In each iteration, np offspring individuals will be added into population, where np is designed to decrease as iteration number increases since the population tends to be less diverse as the algorithm proceeds and fewer offspring individuals are hence enough for population exploration. NP is the upper-bound value of np, which is set as the maximal subpopulation size psU. In the loop of each iteration (Lines 6–14), an offspring ind is first generated by the K -tournament offspring generation (will be detailed in Section 3.3). Afterwards, the selected offspring ind is educated or improved by the neighbourhood local search. The resultant ind will always be added into the population. Furthermore, if ind is infeasible, it will be repaired with a probability of 1/2, and the resultant feasible individual if any will also be added into the population (Lines 10–13, where rand(0,1) uniformly gets a random value in [0, 1]). The repair procedure mainly intensifies the penalty parameters first and then re-educate the target individual as Algorithm 2. 3.1. Solution representation and split-based evaluation

Solution representation. A complete solution for CLRP should indicate the selected depots, the depot assignment of customers, and the detailed vehicle routes to serve all customers. Besides, solutions should be appropriately represented or coded so as to facilitate GA’s evolution in HGA. For this purpose, each individual ind in HGA is represented as two-level chromosomes: the depot chromosome, Λ(ind), consisting of all selected depots; and the customer sequence chromosome πd (ind) for each depot d selected by ind, indicating the visit sequence of the customers assigned to the depot. Each customer sequence of a depot is indeed a giant tour without trip delimiters. Giant-tour-based representations of a similar kind have been quite popular in population metaheuristics since the development of a splitting procedure [32], which is able to transform or split a giant tour into (near) optimal vehicle routes. Splitting procedure. Given a customer sequence, a splitting procedure is aimed to obtain the corresponding (near) optimal vehicle routes by splitting the sequence into exclusive segment(s) such that each segment denotes a vehicle route. Generally, a splitting procedure works by dynamic programming, and to some degree, it is similar to the shortest path algorithm. Specifically, starting from each ith customer s in a sequence πd , the procedure will traverse each subsequent jth customer t after s, and try to update the best-so-far cumulative cost of customer t on the assumption that the sub-sequence from s to t makes up a vehicle route. During the process, the vehicle constraint, e.g., constraint on vehicle capacity, must be considered. In Algorithm 3, we provide the pseudo-code of the splitting procedure used in the proposed HGA. After the procedure, C-value of the last customer pn in πd is indeed the optimal travelling cost of πd . To obtain the vehicle routes, we can start from the end of the last route, i.e., pn, and recursively obtain the end of the route prior to each incumbent route. For the sake of better understanding, we exemplify the solution representation and the splitting-based solution transformation in Fig. 3. The depot chromosome is {1, 3}, indicating that depot 1 and depot 3 are selected. The customer sequence chromosome of depot 1 is (1, 2, 3), which means that customers 1, 2, 3 are assigned to depot 1 and they will be served following the vehicle route(s) obtained by splitting the sequence (1, 2, 3). Herein, (1, 2, 3) is split into two vehicle routes, (1, 2) and (3). That is, a vehicle will start from depot 1 to successively serve

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760

customers 1, 2, and then return; another vehicle will be sent out to serve only customer 3. In the algorithm implementation, twolevel arrays are used to represent the depot chromosome and customer chromosome, respectively. That is, an array stores all the selected depot(s), and for each selected depot, an array is used to store the customer sequence assigned to that depot. 3.2. Population management As a population-based metaheuristic, it is important for HGA to design a favourable population management mechanism to help achieve better performance. In HGA, population management mainly consists of four components: population maintenance, penalty parameter adjustment, population initialization, and population diversification. Population maintenance. As mentioned above, both feasible and infeasible solutions are permitted in HGA, and two subpopulations pop_fs and pop_infs hence need to be maintained. In this work, the size of each subpopulation is controlled in a certain range of [psL, psU]. During HGA’s evolution, once a new solution has been generated, it will always be added into one subpopulation according to its feasibility. If the subpopulation size is found to excess psU after once addition, it will be reduced to psL by removing all the other individuals except for the best psL ones. Penalty parameter adjustment. Capacity constraint is the main constraint of CLRP. In the proposed HGA, both depot capacity violation and vehicle capacity violation are permitted for infeasible solutions, and two corresponding penalty parameters, i.e., γ D and γ V , are used to help evaluate infeasible solutions. The penalized cost of a solution sol composed by D depots and R routes is defined by the weighted sum as follows

Φ (sol) = cost (sol) + γ D ×

D ∑

max (0, loadd − Qd )

as follows, where the TC (gsol) gets the travelling cost of gsol, and the denominator obtains the total demands of customers.

γ D = γ V = TC (gsol) /

+ γV ×

r =1

where the first term cost(sol) is the total cost of sol, including travelling cost, vehicle cost and depot cost, the second term denotes the penalty of depot capacity violation, the third term denotes the penalty of vehicle capacity violation. Furthermore, in order to explore more the boundary between feasible and infeasible solution space, the population is expected to always consist of a reasonable number of both feasible and infeasible individuals. For this purpose, the penalty parameters are self-adjusted during HGA’s evolution according to an expected feasible ratio ρ Exp , which is the expected ratio of feasible individuals in the population. Specifically, at the end of each iteration of HGA, a penalty parameter will be increased or decreased when the real-time feasible ratio ρ X is too small or too large: if ρ X ≤ 0.5 × ρ Exp − 0.05, then γ X = γ X × 2.0; else if ρ X ≤ ρ Exp − 0.05, then γ X = γ X × 1.2; else if ρ X ≥ ρ Exp + 0.05, then γ X = γ X × 0.8; where X is D or V corresponding to the depot constraint or vehicle constraint, ρ X is the ratio of potentially feasible individuals in terms of constraint X. It is noteworthy that a potentially feasible individual of one constraint X may still violate other constraint(s). In the beginning, the penalty parameter initialization is based on a greedy solution gsol, which is constructed in three steps: firstly, repeatedly select a depot with higher capacity until the cumulative depot capacities equal or excess the total demands; secondly, assign each customer to its nearest selected depot while ignoring the depot capacity constraint; thirdly, greedily generate routes for each depot. The two penalty parameters are initialized

di .

Population initialization. In order to quickly locate a promising solution space at the beginning of HGA, we introduce randomness into the classical Clarke-Wright (CW) heuristic [33] for population initialization. The CW is a well-known heuristic for the vehicle routing problem (VRP) and it has always been a popular ingredient when to develop more complicated algorithms for many VRP variants. For example, the GRASP in [29] is indeed an extended and randomized version of CW. The CW is a typical ‘‘route-first, cluster-second’’ method. At the beginning of CW, each customer is assigned to a different route. Then, two routes with the least merging cost or the largest saving will be merged. The merging process will be repeated until no routes can be merged. In the proposed randomized CW, the randomness mainly manifests in the assignment of single-customer routes to depots, i.e., customer assignment to depots. The assignment is completed by successively selecting a depot for each customer, following a random permutation of customers. Herein, all depots’ probabilities are first obtained and then the roulette wheel selection is adopted for depot selection. For each customer i, the probability of each depot d is calculated as follows Prob (d) =

Prob_g(d) (cd,i + ci,d ) + γ D × depot_cap_v io

(1)

where the summation in parentheses gets the travelling cost of the single-customer route of i if it is assigned to depot d, γ D is the penalty parameter of depot capacity violation, depot_cap_vio is the expected capacity violation of d if i is assigned to d, Prob_g(d) is the global probability of depot d and it is computed as follows at the beginning of HGA, Prob_g(d) =

max(0, loadr − Q )

∑ i∈I

d=1 R ∑

5

Qd Cd

+

Qd TCd

+

nd n

(2)

where Qd and Cd are the capacity and opening cost of depot d, n is the total number of customers. The computation is based on a greedy sequence construction. Specifically, for each depot d, starting from the closest customer to d, a greedy sequence πd is constructed for depot d by repeatedly selecting the closest yet-not-selected customer to append at the end of πd until the capacity of d is used up. Then, nd is the number of customers in πd , TC d is travelling cost to sequentially visit customers in πd without considering vehicle constraints. As a result, a depot will have a larger probability value if it is able to satisfy more customer demands while spending less cost. After customer assignment, a single-customer route will be constructed for each customer and the route is assigned to the customer’s depot. Thereafter, routes will be merged in the classic CW way. Notably, we only merge the routes belonging to the same depot, but after a merging operation, all the other depots will be re-tested for the resultant route in case that there is a better depot. Population diversification. Population diversification strategies are often essential for a population-based metaheuristic to jump out of local optima. In the proposed HGA, if the best-so-far solution has not been updated for It RI iterations, both subpopulations, pop_fs and pop_infs, will be shrunk to the minimal acceptable size, i.e., psL/3, by removing the worse individuals. Then, psU individuals will be randomly generated and added into the population. In this way, we introduce diversity into the population and meanwhile preserve the promising historical information. To randomly generate an individual, firstly, a set of depots is selected by roulette wheel selection according to the depots’ global

6

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760

probabilities defined by Eq. (2). These selected depots should have enough capacity to serve all customer demands. Secondly, following a random permutation of customers, we assign each customer i to one of afore-selected depots by roulette wheel selection based on depot probability, which is similarly defined as Eq. (1) except that the numerator is a constant value. Specifically, customer i will be appended to the tail of customer sequence of the selected depot. 3.3. Offspring generation and multi-sequence-based crossover

The population evolution, especially the offspring generation, is the main part of HGA. In this work, a K -tournament offspring generation procedure (Algorithm 4) is customized for HGA in order to generate high-quality offspring individuals. In the procedure, at most K trail individuals are generated and the best one will be selected as the offspring individual and subsequently added into the population. As shown in Algorithm 4, to generate each kth trail individual, two parents P1, P2 are first selected and then a multi-sequence-based crossover is applied to the parents. Particularly, P1 is set by the kth best individual in the population (herein, a feasible individual is thought to be better than an infeasible one), P2 is selected by a binary tournament selection, which randomly picks two individuals from the population and then selects the better one. Herein, by performing crossover on well-behaved individuals (P1) and ordinary individuals (P2), we would like to further exploit the so-far-found promising solution space and meanwhile, explore new solution space.

After parent selection, a multi-sequence-based crossover (MSIX) as shown in Algorithm 5 is performed to obtain a new individual C from two parents P1, P2. Firstly, MSIX classifies all candidate depots into three disjoint sets (Lines 1–3): the common set Λ1 composed by the depots selected by both P1 and P2, the exclusive set Λ2 consisting of each depot d selected only one parent Pd (Pd may be P1 or P2), and the non-used set Λ3 composed by the remain depots not used by any parent. Depots in Λ1 , Λ2 , and Λ3 will have degressive priorities during crossover. Secondly, depot selection and customer sequence inheritance are completed at the same time (Lines 4–14). Specifically, each depot d in the common set Λ1 will always be selected by C and the customer sequence πd (C ) will be obtained by crossover on πd (P1) and πd (P2) (Lines 4–9). If the total capacity of the alreadyselected depots is not enough, MSIX continues to select depots from the exclusive set Λ2 until the total customer demands can be satisfied (Lines 10–14), where each depot is selected by the binary-tournament selection considering the global depot probability, and customer sequence inheritance is performed after each depot selection. During the process, the inheritance of a customer sequence (segment) is completed by copying the sequence but skipping any visited customer. Thus, the order of customer sequence inheritance is important since an earlier-inherited sequence tends to encounter fewer visited customers. Also, during customer inheritance, we skip any demand-unsatisfiable customer whose demand is larger than the depot’s yet-not-occupied capacity. Besides, the sequence inheritance for a specific depot stops once the depot’s remaining capacity is too small to afford the customer with the least demand. Fig. 4 illustrates an example of the multi-sequence-based crossover. Looking at Depot1’s sequence inheritance, customers 6 and 2 are selected as the start point and the end point of the two-point crossover, thus the part of sequence (6, 7, 1, 2) of parent1 is inherited into the child. Then, to inherit customer sequence from parent2, customers 2, 1, 6 are skipped since they have been visited and customer 8 is skipped for the lack of capacity. After the sequence inheritance from two parents, the unvisited customer 5 is inserted into the child. After the completion of depot selection, there may still exist unvisited customers. If so, each unvisited customer i will be assigned to the nearest depot with enough capacity to hold i (Lines 15–21). Finally, some customers may be transferred from their pre-assigned depots to empty depot(s) such that each depot serves at least one customer (Lines 22–29). 3.4. Education by neighbourhood local search Local search (LS) is often necessary for route improvement when solving various VRPs. As mentioned above, when to evaluate an offspring, it will be first transformed into routes by the splitting procedure. Thereafter, a neighbourhood local search (NLS) is adopted to further improve the routes. The neighbourhood N(u) of each customer u is defined as the set of the nsize closest customers to u, where nsize is the neighbourhood size. The NLS mainly utilizes three classic LS operators (as exemplified in Fig. 5), i.e., 2-opt, swap, and relocate, where at most two customers will be swapped or relocated to another position. As given in Algorithm 6, to achieve fast route improvement, NLS always performs the best move on each customer in an iterative way. The iteration is controlled by the local search size lsize. That is, a new iteration starts only when the routes have been improved more than lsize times in the last iteration. For each customer u, each kind of LS move is evaluated on u and its each neighbour v. Then, the best move is performed for u if the routes can be improved. Notably, penalized cost saving is considered during LS evaluation, i.e., evaluate_LS, since infeasible solutions are permitted in HGA. Generally, nsize and lsize should be set according to the total

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760

7

Fig. 4. An example of the multi-sequence-based crossover.

Fig. 5. Exemplifications of three local search operators (swap, relocate, and 2opt).

number of customers, i.e., n. A recommended setting is, nsize=n/4 and lsize = n/10.

constructed in 2004 [34]. The DLPP’s set is a newer benchmark set proposed in 2009 [24], which considers more realistic problem characteristics such as asymmetric cost graph, non-euclidean distances, and strongly heterogeneous depots. Recall that this work focuses on the CLRP with both capacitated depots and vehicles. Although there are also other widely-used LRP instances such as the Tuzun and Burke’s instances (1999) and the Barreto’s instances (2004), these instances have no constraint on either depots or vehicles and will hence not be included in our experiments. Parameter setting. HGA is tested with the same parameter setting: the maximum permitted unimproved iteration number It NI = 25, the number of unimproved iterations for population diversification It RI = 5, the range of subpopulation size [psL, psU] is [25, 100], and K is set as 10 for K -tournament offspring generation. For penalty parameter adjustment, the expected feasible ratio ρ Exp is set as 0.2. As for neighbourhood local search, the neighbourhood size nsize and local search size lsize are set as n/4 and n/10 respectively, where n is the customer number.

4. Computational results

4.1. Prodhon’s instances

In this section, we present the computational results of our proposed HGA. The HGA is coded in Visual Studio C++ IDE 2017 R and runs on a PC with Intel⃝ CoreTM i7-8750H CPU (2.20 GHz), 8G RAM, Windows 10 operating system. Benchmarks. In this section, the proposed HGA is tested on two sets of benchmark instances (available at http://prodhonc.free. fr/): (1) the Prodhon’s instances, and (2) the DLPP’s instances. The Prodhon’s set is a popular and classical benchmark set originally

The Prodhon’s benchmark set consists of 30 instances as shown in Table 1, where the number of customers ranges from 20 to 200 (column ‘‘n’’), the number of depots ranges from 5 to 10 (column ‘‘m’’), and the vehicle capacity is 70 or 150 (column ‘‘Q ’’). Column ‘‘BKS’’ lists the best-known solution cost of each instance, where a cost value with an asterisk indicates it is a proven optimal solution cost. Table 1 also presents the computational results obtained by our HGA based on 10 independent runs,

8

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760 Table 1 Computational results of HGA on the Prodhon’s instances.

a

Instance

n

m

Q

BKS

C

bestIter

CPU

Gap

20-5-1a 20-5-1b 20-5-2a 20-5-2b 50-5-1 50-5-1b 50-5-2 50-5-2b 50-5-2bis 50-5-2bbis 50-5-3 50-5-3b 100-5-1 100-5-1b 100-5-2 100-5-2b 100-5-3 100-5-3b 100-10-1 100-10-1b 100-10-2 100-10-2b 100-10-3 100-10-3b 200-10-1 200-10-1b 200-10-2 200-10-2b 200-10-3 200-10-3b

20 20 20 20 50 50 50 50 50 50 50 50 100 100 100 100 100 100 100 100 100 100 100 100 200 200 200 200 200 200

5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 10 10 10 10 10 10 10 10 10 10 10 10

70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150 70 150

54793a 39104a 48908a 37542a 90111a 63242a 88298a 67308a 84055a 51822a 86203a 61830a 274814a 213615 193671a 157095a 200079a 152441a 287983 231763 243590a 203988a 250882 204317 477248 378351 449571 374330 469433 362817

54793 39104 48908 37542 90111 63242 88298 67308 84055 51822 86203 61830 276846 214250 194059 157178 200635 152441 291709 231049 244938 203988 252363 204664 481154 378648 452043 375232 474971 363971

2 1 1 1 1 8 9 15 5 13 11 7 12 8 12 20 9 15 9 7 11 13 10 11 18 15 12 15 12 15

1.839 1.27 0 1.226 7.413 25.398 31.191 106.519 21.031 39.149 37.258 20.307 310.827 285.589 143.976 222.388 302.402 169.706 272.956 361.904 380.699 162.364 152.807 176.899 1281.01 869.477 699.779 738.66 738.494 827.535

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.74 0.30 0.20 0.05 0.28 0.00 1.29 −0.31 0.55 0.00 0.59 0.17 0.82 0.08 0.55 0.24 1.18 0.32

Indicates a proved optimal solution. A solution better than or equal to the BKS is in boldface.

where the column ‘‘C’’ lists best solution cost found by HGA, the column ‘‘CPU’’ lists the spent time to reach the best solution (in seconds), the column ‘‘bestIter’’ lists the number of iteration when the best solution is obtained by HGA on each instance, the column ‘‘Gap’’ gives each solution cost C’s relative deviation over the best-known solution cost BKS, computed by Gap (C) = 100% × (C − BKS)/BKS. Hence, a smaller Gap value means a better performance, the Gap value 0 indicates the BKS can be found, and a negative Gap value means that the BKS can be updated. From Table 1, we observe that HGA is able to find the optimal solution on 14 out of 30 instances. On the other 16 instances, the best solutions obtained by HGA are also very approximate to the best-known solutions. In particular, except for two instances ‘‘100-10-1’’ and ‘‘200-10-3’’, the HGA can achieves a Gap less than 1% on each of the remain instances. As for time consumption, from column ‘‘CPU’’, the running time of HGA becomes longer as the problem scale increases as expected. HGA takes less than 1 min to find the optimal solutions on almost all small-scale instances with no more than 50 customers, except for instance ‘‘50-5-2b’’ which consumes about 107 s. On medium-scale instances with 100 customers, about 2 to 6 min seems to be enough to find promising solutions, and on the remain largescale instances, it takes about 10 to 20 min. Actually, the time consumption of HGA should still be acceptable considering that in reality, a CLRP cares most a long-term decision and relatively little timeliness is required by the problem. Moreover, from column ‘‘bestIter’’, we can see that HGA usually finds the best solution within 20 iterations, and on some small-scale instances, the best solutions are obtained within 5 iterations. Thus, it should be quite reasonable to diversify the population after 5 iterations without improvement and stop the algorithm after 25 iterations without improvement. Furthermore, Table 2 presents the comparison results (on Gap values) of HGA with the classical and the state-of-the-art algorithms in Table 2, including the GRASP [29], MA | PM [35],

LRGTS [36], GRASP × ELS (2010 [30]), SALRP (2010 [37]), MACO (2013 [5]), hybridGA (2016 [31]). Actually, there are also some more recent algorithms, e.g., a PSO variant (2017, [38]), which are not introduced into the comparison due to their moderate performance. Among the compared algorithms, hybridGA seems to have the best overall performance, and SALRP seems to win more times on larger scale instances. From Table 2, as shown in the row ‘‘avg.’’, our HGA can achieves the best average Gap, i.e., 0.24%, over the best-known solutions, and our HGA finds the largest number of optimal solutions, i.e., 14 out 30, as shown in row ‘‘#opt’’. The best Gap value(s) on each instance has/have been marked in bold and the times of each algorithm to obtain the best Gap have been summarized in row ‘‘#best’’. Therein, we can see that HGA obtains the best Gap on more instances, i.e., 19 out of 30 instances, than all the other compared algorithms. Indeed, as shown in row ‘‘#≤’’, HGA is able to achieve the same good or better performance than each compared algorithm on most even all instances. In conclusion, the proposed HGA is quite promising and competitive compared with existing well-known algorithms for CLRP. 4.2. DLPP’s instances As mentioned above, the DLPP’s instances is a newer set of CLRP instances aimed to simulate more realistic characteristics [24]. Taking the instance ‘‘DLPP_40’’ for an example, the instance’s graphic representation is shown in Fig. 6(a) with 14 customers (grey nodes) and 3 depots (black nodes), and Fig. 6(b) details a solution (with two routes) for the instance obtained by HGA. From Fig. 6, we can see that the instance map is quite similar to a real-life road net, where nodes are connected by directional arcs and there may be zero, one or two unidirectional arc(s) between two nodes. Each arc has a specific cost. Note that not every node is a customer and the best path from one customer

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760

9

Table 2 Comparative results (of Gap value) on the Prodhon’s instances. Instance

GRASP

MA|PM

LRGTS

GRASP × ELS

SALRP

MACO

hybridGA

HGA

20-5-1aa 20-5-1ba 20-5-2aa 20-5-2ba 50-5-1a 50-5-1ba 50-5-2a 50-5-2ba 50-5-2bisa 50-5-2bbisa 50-5-3a 50-5-3ba 100-5-1a 100-5-1b 100-5-2a 100-5-2ba 100-5-3a 100-5-3ba 100-10-1 100-10-1b 100-10-2a 100-10-2ba 100-10-3 100-10-3b 200-10-1 200-10-1b 200-10-2 200-10-2b 200-10-3 200-10-3b

0.42 0.00 0.00 0.00 0.58 2.40 0.55 1.09 0.00 0.46 1.37 0.10 1.68 1.19 3.02 1.56 1.96 1.41 12.22 17.14 4.31 1.26 7.95 5.80 2.84 10.15 14.04 1.51 5.81 7.22

0.00 0.00 0.00 0.00 0.05 0.00 0.00 0.87 0.00 0.00 0.00 0.00 2.59 1.42 0.98 0.15 0.83 0.58 9.93 16.61 0.63 0.52 1.11 0.24 1.31 0.45 0.50 0.18 1.85 0.56

0.62 0.00 0.00 0.00 0.05 0.02 0.47 0.58 0.15 0.33 0.00 0.00 1.14 0.59 1.48 0.44 0.94 1.49 1.36 1.63 1.28 0.22 3.10 0.77 0.93 0.60 0.84 0.81 1.54 0.67

0.00 0.00 0.00 0.00 0.00 0.00 0.39 0.00 0.00 0.00 0.00 0.00 0.78 1.05 0.31 0.18 0.13 0.06 4.67 16.32 0.08 0.00 1.05 0.38 1.93 1.05 0.45 0.22 1.91 0.65

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.00 0.00 0.29 1.41 0.81 1.12 0.23 0.04 0.08 0.02 1.06 1.06 0.91 0.65 0.00 0.34 0.79 1.38 0.28 0.63 0.95 0.24

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.51 0.33 0.40 0.08 0.48 0.18 1.09 1.55 0.69 0.75 1.36 0.23 0.33 0.14 0.42 0.17 1.22 0.71

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.61 0.59 0.45 0.08 0.14 0.15 1.38 0.75 0.11 0.00 1.15 0.22 0.70 0.73 0.44 0.10 1.37 0.77

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.74 0.30 0.20 0.05 0.28 0.00 1.29 −0.31 0.55 0.00 0.59 0.17 0.82 0.08 0.55 0.24 1.18 0.32

avg. #opt #best #≤

3.60 4 4 30/30

1.38 10 10 30/30

0.74 5 5 30/30

1.05 12 13 26/30

0.41 9 16 22/30

0.35 12 14 25/30

0.32 13 14 24/30

0.24 14 19

a

Indicates an optimal solution is already known for the instance. The best Gap value(s) on each instance is/are in boldface. Row ‘‘avg.’’ gives the average Gap obtained by each algorithm; row ‘‘#opt’’ counts the number of the optimal solutions found by each algorithm; row ‘‘#best’’ summarizes the times of each algorithm obtaining the best Gap among all compared algorithms; row ‘‘#≤’’ lists the times of each algorithm obtaining the same or worse Gap compared with HGA.

Fig. 6. Instance ‘‘DLPP_40’’: (a) graphic representation of the instance, and (b) a solution with two routes for the instance obtained by HGA (in each route, the red marks the selected depot, the served customers, and the detailed path).

10

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760 Table 3 Computational results of HGA on the DLPP’s instances. Instance DLPP_16 DLPP_19 DLPP_21 DLPP_25 DLPP_28 DLPP_32 DLPP_34 DLPP_36 DLPP_38 DLPP_40 DLPP_42 DLPP_47 DLPP_51 DLPP_52 DLPP_55 DLPP_68 DLPP_76 DLPP_87 DLPP_90 DLPP_100 DLPP_110 DLPP_126 DLPP_130 DLPP_137 DLPP_170 DLPP_210 DLPP_224 DLPP_260 DLPP_285 DLPP_341 avg.

#nodes 16 19 21 25 28 32 34 36 38 40 42 47 51 52 55 68 76 87 90 100 110 126 130 137 170 210 224 260 285 341

#arcs 30 52 66 69 76 94 103 92 112 107 106 135 138 159 162 228 244 268 292 319 350 403 401 433 534 659 710 863 907 1124

n 5 8 10 11 10 15 13 14 13 14 18 24 24 29 26 31 30 39 35 50 49 59 68 55 72 100 96 113 130 164

m 2 3 2 2 1 2 2 2 2 3 3 3 3 3 5 6 6 5 5 7 5 7 10 6 10 13 11 14 17 23

BKS

HGA

C

CPU

C

bestIter

CPU

Gap

8163.35 8782.99 9242.07 9255.75 15042.5 7383.21 8018.75 9877.2 6376.11 11328.4 19475.6 19029.7 20225.7 19869.2 24978.8 19183.2 19438.6 31626.7 24337.9 51475.6 47444.6 35052.4 36021.6 50233.2 60786.5 65201.9 81068.6 89787.1 98500.6 145951

1 0 1 1 0 0 2 0 0 0 1 1 3 4 39 76 45 19 44 23 106 49 71 30 154 11 393 580 133 532

8163.35 8782.99 9242.07 9255.75 15042.5 7383.21 8018.75 9663.36 6376.11 11328.4 19337.6 18303 19855.2 19869.1 25028.2 19138.1 19438.6 30942.9 23990.2 51475.6 47560.4 35051.4 35925 49119.5 60743.2 64811.1 78232.3 88762.4 96474.3 140307

0 1 0 0 0 0 3 0 0 0 0 1 0 0 3 2 13 10 2 4 5 4 7 8 4 16 6 4 8 5

0 0 0 0 0 0 1 0 0 0 2 2 0 3 6 4 14 22 5 18 20 125 53 41 43 240 109 108 229 227

0.000 0.000 0.000 0.000 0.000 0.000 0.000 −2.165 0.000 0.000 −0.709 −3.819 −1.832 −0.001 0.198 −0.235 0.000 −2.162 −1.429 0.000 0.244 −0.003 −0.268 −2.217 −0.071 −0.599 −3.499 −1.141 −2.057 −3.867

35105.29

77.3

34587.39

3.5

42.4

−0.854

The composite column ‘‘BKS’’ presents the total cost of BKS (column ‘‘C’’) and CPU time in seconds to obtain the BKS (column ‘‘CPU’’). Column ‘‘Gap’’ lists the gap of HGA’s solution costs over the best-known solutions.

to another customer may need to pass by other nodes or customers. Besides, depots in the DLPP’s instances have much more heterogeneous capacities than most classical CLRP instances [24], and the corresponding capacity constraint hence becomes tighter and the depot location becomes more challenging. For example, in instance ‘‘DLPP_40’’, the three candidate depots’ capacities are 1871, 1569, and 4786, and the vehicle capacity is 1547. To serve all customers, selecting the former two depots or selecting the third depot will lead to very different routing plans. Also, due to the vehicle capacity constraint, two customers (node 5, and node 33 in Fig. 6) cannot be served by route 1 even though the route needs to pass by the two customers as shown in Fig. 6(b1). The information of DLPP’s instances and the corresponding computational results obtained by HGA have been provided in Table 3, where an obvious superiority of HGA can be observed compared with the best-known solutions. In the table, the four columns ‘‘#nodes’’, ‘‘#arcs’’, ‘‘n’’, ‘‘m’’ list the number of nodes, arcs, customers, and depots respectively. The best-known solutions (BKS) are drawn from [24]. From Table 3, as shown in bold, HGA is able to find the same good or better solutions than BKS on 28 out of 30 instances. Particularly, HGA obtains a negative Gap over the BKS on 16 instances, and the Gap is less than −1.0% on 10 instances. That is, the best-known solutions can be updated by our HGA on these instances. On the whole, as shown in row ‘‘avg.’’, HGA can achieve an average Gap of −0.854% over the bestknown solutions, and its average running time is 42.7s, less than the 77.3 s of BKS. Although the CPU statistics are not obtained on the same machine, the time comparison at least shows that the proposed HGA is not very time-expensive compared with other algorithms and the good performance of HGA should not be achieved at the expense of a higher time complexity.

5. Conclusion In this study, focusing on the capacitated location routing problem, we propose a hybrid genetic algorithm (HGA) that is able to deal with tight capacity constraints, which explores both feasible solution space and infeasible solution space. Experiments have been conducted on two sets of CLRP benchmark instances, i.e., a classical set and a real-life-like set with tighter capacity constraints. The computational results show that the proposed HGA is competitive with other well-known algorithms for CLRP in terms of both solution quality and time efficiency. In particular, HGA is able to find the optimal solutions on quite a number of instances. Moreover, for the benchmark set with tighter constraints, the best-known solutions of many instances can be updated by HGA. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work is supported by the National Natural Science Foundation of China (61773410, 61673403). References [1] S. Salhi, G.K. Rand, The effect of ignoring routes when locating depots, European J. Oper. Res. 39 (2) (1989) 150–156. [2] G. Nagy, S. Salhi, Location-routing: Issues models methods, European J. Oper. Res. 177 (2) (2007) 649–672.

X. Yu, Y. Zhou and X.-F. Liu / Applied Soft Computing Journal 85 (2019) 105760 [3] C. Prodhon, C. Prins, A survey of recent research on location-routing problems, European J. Oper. Res. 238 (1) (2014) 1–17. [4] M. Drexl, M. Schneider, A survey of variants and extensions of the location-routing problem, European J. Oper. Res. 241 (2) (2015) 283–308. [5] C. Ting, C. Chen, A multiple ant colony optimization algorithm for the capacitated location routing problem, Int. J. Prod. Econ. 141 (1) (2013) 34–44. [6] C. Watson-Gandy, P.J. Dohrn, Depot location with van salesmen—a practical approach, Omega 1 (3) (1973) 321–329. [7] J. Perl, M.S. Daskin, A unified warehouse location-routing methodology, J. Bus. Logist. 5 (1) (1984) 92–111. [8] J. Perl, M.S. Daskin, A warehouse location-routing problem, Transp. Res. B 19 (5) (1985) 381–396. [9] A. Bruns, A. Klose, P. Stähly, Restructuring of swiss parcel delivery services, OR Spektrum 22 (2) (2000) 285–302. [10] H. Farrokhi-Asl, R. Tavakkoli-Moghaddam, B. Asgarian, E. Sangari, Metaheuristics for a bi-objective location-routing-problem in waste collection management, J. Ind. Prod. Eng. 34 (4) (2017) 239–252. [11] A. Billionneta, S. Elloumib, L.G. Djerbi, Designing radio-mobile access networks based on synchronous digital hierarchy rings, Comput. Oper. Res. 32 (2) (2005) 379–394. [12] V. Nguyen, C. Prins, C. Prodhon, Solving the two-echelon location routing problem by a GRASP reinforced by a learning process and path relinking, European J. Oper. Res. 216 (1) (2012) 113–126. [13] V. Nguyen, C. Prins, C. Prodhon, A multi-start iterated local search with tabu list and path relinking for the two-echelon location-routing problem, Eng. Appl. Artif. Intell. 25 (1) (2012) 56–71. [14] N. Wichapa, P. Khokhajaikiat, Solving a multi-objective location routing problem for infectious waste disposal using hybrid goal programming and hybrid genetic algorithm, Int. J. Ind. Eng. Comput. 9 (1) (2018) 75–98. [15] K. Chang, H. Zhou, G. Chen, H. Chen, Multiobjective location routing problem considering uncertain data after disasters, Discrete Dyn. Nat. Soc. 2017 (2017). [16] F. Rayat, M. Musavi, A. Bozorgi-Amiri, Bi-objective reliable locationinventory-routing problem with partial backordering under disruption risks: A modified AMOSA approach, Appl. Soft Comput. 59 (2017) 622–643. [17] N. Ghaffari-Nasab, M.S. Jabalameli, M.B. Aryanezhad, A. Makui, Modeling and solving the bi-objective capacitated location-routing problem with probabilistic travel times, Int. J. Adv. Manuf. Technol. 67 (9–12) (2013) 2007–2019. [18] Y. Marinakis, An improved particle swarm optimization algorithm for the capacitated location routing problem and for the location routing problem with stochastic demands, Appl. Soft Comput. 37 (2015) 680–701. [19] N. Ghaffari-Nasab, S.G. Ahari, M. Ghazanfari, A hybrid simulated annealing based heuristic for solving the location-routing problem with fuzzy demands, Sci. Iran. 20 (3) (2013) 919–930. [20] M.H.F. Zarandi, A. Hemmati, S. Davari, The multi-depot capacitated location-routing problem with fuzzy travel times, Expert Syst. Appl. 38 (8) (2011) 10075–10084.

11

[21] C. Prodhon, A metaheuristic for the periodic location-routing problem, 4OR (2008) 159–164. [22] C. Prodhon, A hybrid evolutionary algorithm for the periodic locationrouting problem, European J. Oper. Res. 210 (2) (2011) 204–212. [23] V.C. Hemmelmayr, J. Cordeau, T.G. Crainic, An adaptive large neighborhood search heuristic for two-echelon vehicle routing problems arising in city logistics, Comput. Oper. Res. 39 (12) (2012) 3215–3228. [24] C. Duhamel, P. Lacomme, C. Prodhon, C. Prins, A GRASP × ELS approach for real-life location routing problems, in: 2009 International Conference on Computers & Industrial Engineering, IEEE, 2009, pp. 1082–1087. [25] J. Belenguer, E. Benavent, C. Prins, C. Prodhon, R.W. Calvo, A branch-andcut method for the capacitated location-routing problem, Comput. Oper. Res. 38 (6) (2011) 931–941. [26] R. Baldacci, A. Mingozzi, R. Wolfler Calvo, An exact method for the capacitated location-routing problem, Oper. Res. 59 (5) (2011) 1284–1296. [27] C. Contardo, J. Cordeau, B. Gendron, An exact algorithm based on cut-andcolumn generation for the capacitated location-routing problem, INFORMS J. Comput. 26 (1) (2013) 88–102. [28] S. Ponboon, A.G. Qureshi, E. Taniguchi, Branch-and-price algorithm for the location-routing problem with time windows, Transp. Res. E 86 (2016) 1–19. [29] C. Prins, C. Prodhon, R.W. Calvo, Solving the capacitated location-routing problem by a GRASP complemented by a learning process and a path relinking, 4OR 4 (3) (2006) 221–238. [30] C. Duhamel, P. Lacomme, C. Prins, C. Prodhon, A GRASP × ELS approach for the capacitated location-routing problem, Comput. Oper. Res. 37 (11) (2010) 1912–1923. [31] R.B. Lopes, C. Ferreira, B.S. Santos, A simple and effective evolutionary algorithm for the capacitated location–routing problem, Comput. Oper. Res. 70 (2016) 155–162. [32] C. Prins, A simple and effective evolutionary algorithm for the vehicle routing problem, Comput. Oper. Res. 31 (12) (2004) 1985–2002. [33] G. Clarke, J.W. Wright, Scheduling of vehicles from a central depot to a number of delivery points, Oper. Res. 12 (4) (1964) 568–581. [34] C. Prins, C. Prodhon, R. Wolfler-Calvo, Nouveaux algorithmes pour le problème de localisation et routage sous contraintes de capacité, Mosim, 2004, pp. 1115–1122. [35] C. Prins, C. Prodhon, R.W. Calvo, A memetic algorithm with population management (MA|PM) for the capacitated location-routing problem, in: European Conference on Evolutionary Computation in Combinatorial Optimization, Springer, 2006, pp. 183–194. [36] C. Prins, C. Prodhon, A. Ruiz, P. Soriano, R. Wolfler Calvo, Solving the capacitated location-routing problem by a cooperative Lagrangean relaxation-granular tabu search heuristic, Transp. Sci. 41 (4) (2007) 470–483. [37] F.Y. Vincent, S. Lin, W. Lee, C. Ting, A simulated annealing heuristic for the capacitated location routing problem, Comput. Ind. Eng. 58 (2) (2010) 288–299. [38] Z. Peng, H. Manier, M. Manier, Particle swarm optimization for capacitated location-routing problem, IFAC-PapersOnLine 50 (1) (2017) 14668–14673.