CAPACITATED LOCATION ROUTING PROBLEM WITH SIMULTANEOUS PICKUP AND DELIVERY UNDER THE RISK OF DISRUPTION

This paper develops a new mathematical model to study a location-routing problem with simultaneous pickup and delivery under the risk of disruption. A remarkable number of previous studies have assumed that network components (e.g., routes, production factories, depots, etc.) are always available and can permanently serve the customers. This assumption is no longer valid when the network faces disruptions such as flood, earthquake, tsunami, terrorist attacks and workers strike. In case of any disruption in the network, tremendous cost is imposed on the stockholders. Incorporating disruption in the design phase of the network will alleviate the impact of these disasters and let the network resist disruption. In this study, a mixed integer programming (MIP) model is proposed that formulates a reliable capacitated location-routing problem with simultaneous pickup and delivery (RCLRP-SPD) services in supply chain distribution network. The objective function attempts to minimize the sum of location cost of depots, routing cost of vehicles and cost of unfulfilled demand of customers. Since the model is NP-Hard, three meta-heuristics are tailored for large-sized instances and the results show the outperformance of hybrid algorithms comparing to classic genetic algorithm. Finally, the obtained results are discussed and the paper is concluded. Mathematics Subject Classification. 90C59. Received April 25, 2019. Accepted March 30, 2021.

The rest of the paper is organized as follows. Section 2 reviews the relevant papers in the literature. The problem is stated and formulated in Section 3. Next, Section 4 proposes a classic genetic and two hybrid solution approaches to solve the proposed MILP model. Afterward, the results are provided in Section 5. Finally, the paper is concluded in Section 6.

Literature review
This section reviews the latest and most relevant RCLRP and RCLRP-SPD papers. The literature review is presented in three sections: Section 2.1 reviews the location-routing problems, Section 2.2 reviews the reliable location-routing problems and finally, Section 2.3 reviews the solution approaches, particularly the heuristic and meta-heuristic algorithms employed in LRPs.

Location-routing problem (LRP)
Location-routing problems (LRPs) have been extensively studied by several authors in the literature and the interested readers are referred to the studies dealing with LRPs and their proposed classifications by Prodhon and Prins [49] and Drexl and Schneider [12].
Karaoglany et al. [23] studied LRP with simultaneous pickup and delivery, wherein the authors developed two mixed integer formulations (i.e., a node-based and a flow-based formulation) incorporating a family of valid inequalities to strengthen the model. In order to solve the mathematical model, the authors developed a heuristic algorithm based on simulated annealing (SA) algorithm. Through the computational results, they showed that the flow-based formulation outperforms the node-based formulation in terms of the solution quality and running time for small-size instances (10-30 customers), while the node-based formulation performs better than the flow-based formulation for medium-size instances (50-100 customers). Zarandi et al. [61] examined the LRP with time windows under operational uncertainty, where customers demand and travel times are assumed to be fuzzy variables. The authors designed a fuzzy chance-constrained programming model using credibility theory, and a simulation-embedded SA algorithm was then presented to solve the problem. For initializing the solutions of the proposed SA algorithm, they used a heuristic method based on fuzzy c-means clustering with Mahalanobis distance and sweeping method.
Huang [19] presented a three-stage solution approach to deal with the multi-compartment capacitated LRP with pickup-delivery and stochastic demands. This three-stage solution approach solves the model to determine the depot locations, assign customers to the located depots and route the vehicles through customers subsequently. The objective is to minimize depot opening cost, vehicle cost, and travel costs and violation of the vehicle and depot capacity constraints. Pichka et al. [47] addressed the two-echelon open LRP (2E-OLRP). This problem is a new variant of the classical LRP that considers location and routing decisions in two-echelon supply chains in which third-party logistics providers are used. In this research, three new mathematical models and a hybrid SA algorithm are developed to solve the 2E-OLRP.
Zhang et al. [65] studied the multi-depot emergency facilities LRP with uncertain demand. To incorporate the uncertain characteristics of the emergency response system, they have developed an uncertain multi-objective programming model. Due to the computational complexity of the model, a hybrid intelligent algorithm is designed to solve the proposed model by combining uncertain simulation and a genetic algorithm. Moshref-Javadi and Lee [41] introduced the latency LRP (LLRP) aiming at the customers' waiting time minimization by optimally locating depots and routing of vehicles. To solve the model, the authors presented a memetic algorithm and a recursive granular algorithm.
Khalili-Damghani et al. [27] studied a bi-objective LRP to distribute perishable products in a supply chain. Considering the perishability of products, minimizing transportation costs is not the main objective in perishable supply chain planning. Therefore, the objectives are to minimize the total cost of the system and to balance the transportation cost of perishable products at distribution centers while the due-date of perishable products should be met. Timely delivery is an essential factor in the distribution process of perishable products. So, the authors considered a set of constraints to deliver the products no later than a predetermined time. Since the proposed model is NP-hard as well as bi-objective, they developed two solution procedures, i.e., an exact method called epsilon-constraint and an evolutionary computation, called non-dominated sorting genetic algorithm (NSGA-II) to find out the near Pareto optimal solutions. Ghaffari-Nasab et al. [17] addressed a bi-objective CLRP with probabilistic travel times, where the aim is to minimize the total cost and the maximum delivery time to the customers. The authors presented mathematical programming formulations to model the problem using two stochastic programming approaches. The deterministic equivalents of the two stochastic models are extracted and solved by a variable neighborhood descent (VND).

Reliable location-routing problem (RLRP)
In deterministic LRP we make decisions about depots locations, customer allocation and vehicle routing to balance the trade-off between fixed set-up costs and variable transportation costs. However, considering disruption caused by natural disasters, terrorist attacks, or labor strikes, some of the located depots may become unavailable to serve customers. When disruption occurs in a depot, its customers may have to be re-allocated from their original depot to other depots with higher transportation costs [9]. Ahmadi-Javid and Seddighi [1] studied the LRP in a supply chain network composed of a set of producer-distributors that produced and distributed a single commodity to a number of customers. In this work, the production capacity of each producer (distributor) varies randomly because of numerous sources of disruptions. In addition, the vehicles utilized the distribution system are also disrupted randomly.
Xie et al. [55] studied the RLRP wherein all located depots are independently disrupted with a same probability. One of the assumptions made in this study is that customers receive delivery service in fixed groups. In other words, a vehicle is assigned to visit the customers in each group. Under this assumption, if one depot fails, then all of its customers will be reassigned to another backup depot within the same group, which is rather restrictive. Moreover, the vehicle delivery distance is approximated by the sum of the local travel distances within the group and the line-haul distance between this group and the affected depot. Zhang et al. [64] examined an LRP that includes a set of warehouses which are randomly disrupted. A scenario dependent MILP model is proposed to optimize the location of depots, delivery routing towards faraway areas and alternative planning. The authors proposed a meta-heuristic algorithm exploiting maximum-likelihood sampling, route-reallocation improvement, two-stage neighborhood search and SA algorithm.
Mohammadi et al. [38] developed a MINLP location-allocation model in a supply chain management problem to design reliable logistics networks that perform as well as possible under normal condition, and also perform relatively well when disruptions occur. They proposed a single objective scenario based robust optimization problem to minimize the nominal cost. In addition, to solve the proposed model, a new self-adaptive metaheuristic algorithm based on genetic and imperialist competitive algorithms is developed. Li et al. [28] proposed two related models for the design of distribution networks exposed to the risk of disruption (i.e., reliable p-median problem and reliable uncapacitated fixed-charge location problem). Both models accounted for heterogeneous facility failure probabilities and one layer of supplier backup. Both models are formulated as nonlinear integer programming models which are proven to be NP-hard. To solve the proposed model, the authors developed an algorithm based on Lagrangian relaxation. Wang et al. [54] examined a facility location problem under random facility disruption and presented a MINLP model for designing a reliable supply chain network that is exposed to the risk of disruption. They solved the proposed model using a Lagrangian relaxation-based algorithm.

Heuristic/meta-heuristic approaches
Among different evolutionary algorithms, genetic algorithms (GAs) are probably the most widely used class of evolutionary algorithms and have been applied as well to the VRPs and LRPs. A comprehensive survey of different types of GAs for these problems can be found in [20].
To solve CLRP, Derbel et al. [11] developed a GA-based hybrid algorithm and an iterated local search (ILS). Since GA could fail to converge to the global optimum and ILS could fall to the local optimum too quickly, they embedded ILS into GA to refine the search through successive iterations and maximize the chance of convergence to the optimal solution. Karaoglan and Altiparmak [21] proposed a memetic algorithm for solving a CLRP with Mixed Backhauls. They conducted an experimental study and compared the results with the branch-and-cut algorithm's lower bounds. Results showed that the memetic algorithm was able to find optimal or very good quality solutions within a reasonable computation time.
Escobar et al. [14] presented a two-phase hybrid heuristic (2-Phase HGTS), wherein the first phase constructs a solution (i.e., Construction phase), and the second phase (i.e., Improvement phase) employs a modified Granular Tabu Search (GTS) approach with different diversification strategies to modify the constructed solution. Furthermore, a random perturbation procedure is considered to prevent the algorithm from remaining in a local optimum for a given number of iterations. Ardjmand et al. [3] developed a new model for the location-routing problem. They proposed a new GA that solves the model within a reasonable computational time. Yu and Lin [56] proposed an SA heuristic with a special solution encoding scheme for solving the LRP with simultaneous pickup and delivery (LRPSPD). The solution encoding scheme can broaden the search space and facilitate the search for a better solution. The results showed the efficiency of the proposed SA to solve LRPSPD and its superiority over the existing exact approaches in terms of the quality of the solution. Yu et al. [57] presented an SA algorithm in which each solution is represented as a giant tour. Move, swap and 2-opt were used as neighborhoods and they operated on the giant tour. Therefore, a neighborhood move may change the position of the customer, open or close depots or reassign customers to the locations. Ferreira and de Queiroz [16] hybridized two heuristic algorithms based on the SA algorithm combined with a diversification procedure to solve a CLRP. Table 1 summarizes the reviewed articles with their characteristics in two main aspects as problem characteristics and solution approach. In Table 1, SPD, LSbMH and PbMH stand for Simultaneous Pickup and Delivery, Local Search based Meta-Heuristics (i.e., single solution based meta-heuristics) and Population based Meta-Heuristics, respectively. This Table 1 helps to well position our work comparing to the literature.

This paper's contributions
Based on Table 1, the main contributions of this paper that distinguish it from the existing studies in the literature are: -Studying a reliable capacitated location-routing problem with simultaneous pickup and delivery (RCLRP-SPD) wherein depots are stochastically disrupted and become unavailable to serve the customer demands. -Proposing a new mixed-integer linear programming (MILP) model for the proposed RCLRP-SPD.
-Developing an efficient scenario-based approach to cope with the stochastic disruption of depots.
-Proposing three tailored meta-heuristic algorithms (i.e., two hybrid algorithms and a classical genetic algorithm) for solving the proposed MILP.

Problem statement and formulation
This section proposes a MILP model to formulate a RCLRP-SPD. Section 3.1 states the problem with its underlying assumptions. Afterward, Section 3.2 proposes the MILP model for the problem.

Problem statement
The goal of the proposed RCLRP-SPD is to choose and locate a set of depots (facilities) and to build vehicles' routes to meet customers delivery and pickup demands, such that the total expected cost of location, routing and disruption is minimized. In this problem, each potential depot could be stochastically disrupted, and each customer has certain demand. In Figures 1 and 2, two solutions of the CLRP and RCLRP are illustrated.
In Figure 1, three depots (i.e., depots number 1, 4 and 5) are located (opened) and the customers are allocated to the located depots and served in a particular order. Furthermore, if each of the depots fails, its allocated customers will not be served anymore and this would impose high costs on the system. On the other hand, Figure 2 illustrates the reliable version of the CLRP, wherein the located depot number 4 is disrupted and some  [16] of its customers are served by depot number 1 but with a higher transportation cost and finally some of them are served via the emergency depot.
The main assumptions of the proposed RCLRP-SPD are as follows: -Only one depot and one vehicle can be allocated to each customer.
-Route of each vehicle starts and ends at the same depot.
-A single product is delivered and picked up in this problem.
-All vehicles have a same capacity (homogeneous fleety), but depots have different capacities.  -There is a fixed cost for vehicles and depots.
-Several vehicles can be assigned to each located depot.
-Disruptions are assumed to be independent and depots can simultaneously be disrupted.
-Simultaneous pickup and delivery are assumed for each customer.
-Disruption probability of each depot is known.
-The penalty cost of the unmet demand is known.
In real-world situations, it is not possible to predict accurately the exact value of future parameters. Therefore, considering uncertainty makes the decision makers more capable to have better planning for the future.
The first step in dealing with uncertainty in modeling and optimizing supply chain related problems is to determine how to display uncertainty. Three well-known and widely-used methods for this purpose are [5]: (I) Distribution-based approach, wherein the normal distribution is typically used with a specific mean and standard deviation to model uncertain parameters [2,31,34]. (II) Fuzzy-based approach, wherein the parameters are treated as fuzzy numbers with membership functions [36,42,43,60,62]. (III) Scenario-based approach, wherein some discrete scenarios with relevant levels of probability are used to describe the expected occurrence of specific results [38,64].
In this paper, disruption of depots is modeled by a set of discrete scenarios. Each scenario specifies a subset of depots that become non-operational or disrupted with a specific probability of occurrence. In case of a disruption, the affected customers whose allocated depots are disrupted, may be served by other located depots or emergency depots. Obviously, allowing more surviving (backup) depots to cover customers that were once serviced by a disrupted depot will increase the flexibility of the distribution system. Accordingly, Section 3.2 formulates a scenario based RCLRP-SPD.

Problem formulation
This section presents a scenario-based mathematical formulation for the studied problem. This formulation is based on the model proposed by Zhang et al. [64] for the RLRP but extends it to account for the multiple vehicles. Unlike [64], which has only considered delivery decisions, the presented formulation tackles with both pickup and delivery decisions. Our model aims to minimize the total expected costs of all scenarios including fixed costs, expected travel costs and expected penalty costs. Let G = (V, A) be an undirected network, where V is the set of nodes, composed of a subset of geographically dispersed customers denoted by I and a subset of potential depots denoted by J. A = {(i, j)|i, j ∈ V, i = j} is the set of arcs that connect each pair of nodes in V . Let S be the set of scenarios, each specifying a set of disrupted depots. Before proposing the MILP model, necessary notations are explained as follows.

Sets
I Set of customers. J Set of potential depots. S Set of scenarios. H Set of vehicles. e Emergency depot.
Delivery goods demand of customer i. p i Pickup goods demand of customer i. θ i Penalty of serving customer i using the emergency depot. cv Capacity of vehicles. f v Fixed cost of each vehicle. tc ij Traveling cost from node i to node j. a js If depot j in scenario s is not disrupted 1, otherwise 0. Disrupted depots are known a priori ; therefore, a js is a binary parameter. pr s Probability that scenario s occurs. M A big number.

Decision variables
If customer i is served by depot j in scenario s 1, otherwise 0. X ikhjs If node k is visited immediately after node i on a route originating from depot j by vehicle h in scenario s 1, otherwise 0. T ijs Delivery goods transported between node i and node j in scenario s. W ijs Pickup goods transported between node i and node j in scenario s. U ihj An auxiliary integer variable to eliminate the subtours.
Emergency depot e, has no fixed cost, a failure probability of 0 and infinite capacity. If a customer i is not served in a scenario, we assign it to the emergency depot e, which induces an unmet demand penalty θ i .
In certain scenarios, if a depot survives, then it is "normal" (non-disrupted); if it is both normal and open, we call it "available". Considering this point, customers can only obtain service from available depots in each scenario. Consider J ns = {j|j ∈ J, a js = 1} as the set of normal depots in scenario s. In the strategic decision, assume that the set of depots selected to be opened is denoted by J o . Then, the set of available depots in scenario s is J as = J ns ∩ J o . Since each scenario specifies a subset of depots which become non-operational or disrupted (a js , ∀j ∈ J, ∀s ∈ S), the probability that scenario s occurs is represented using equation (3.1).
A total of 2 |J| possible scenarios exist, which leads to an exponential growth in problem size and huge computational burden to solve the problem. By increasing the number of depots by one, we double the number of scenarios.
Based on the provided notations, the proposed MILP model for the RCLRP-SPD is as follow: subject to: The objective function (3.2) minimizes sum of the depots fixed costs (scenario-independent) and the expected costs of each scenario including routing costs, penalty costs and vehicles fixed costs (scenario-dependent) respectively. Constraint (3.3) ensures that a customer in any scenario could be assigned to a depot if the depot is opened. Constraint (3.4) ensures that each customer is assigned to exactly one depot in each scenario. Constraints (3.5) and (3.6) are depots capacity constraints in each scenario. Constraint (3.7) links the allocation and routing in each scenario. Constraint (3.8) ensures that, in each scenario, each customer is visited only once. Constraint (3.9) is the flow conservation constraint which forces that, in each scenario, the input flow to each customer/depot is equal to the output flow from that customer/depot. Constraint (3.10) checks the feasibility of the route in each scenario when starting the routes from the depots. Constraint (3.11) ensures that each vehicle is assigned at most to a single depot in each scenario. Constraints (3.12)-(3.14) control and impose the vehicles capacities in each scenario. Constraint (3.15) ensures that the amount of delivery demand when vehicles return to depots is zero. Constraint (3.16) ensures the amount of pickup demand when vehicles start routes from depots is zero. Constraint (3.17) is the sub-tour elimination constraint in each scenario. Constraints

Solution approaches
Metaheuristic algorithms have been widely used in the literature to solve optimization problems [13,24,40]. In this regard, the papers in the literature use either population-based evolutionary algorithms or single solutionbased local search algorithms. The former is powerful in diversification and the latter is powerful in intensification. In this paper, we attempt to propose a hybrid algorithm which is powerful in both diversification and intensification [13,24]. Indeed, hybrid algorithms attempts to benefit from the strengths of different algorithm in a single algorithm [25,26,37,39,59,63]. Accordingly, we proposed a hybrid algorithm by the combination of genetic algorithm (GA) (as a powerful algorithm in diversification) and two local search algorithms (as powerful algorithms in intensification).
Genetic algorithm (GA), which was first developed by Holland [18], is one of the most well-known evolutionary optimization algorithms. GA evolves a population of randomly generated solutions with size P size using crossover and mutation operators. At each iteration, a number of solutions are generated through crossover operator with rate C R and a set of other solution are generated through mutation operator with rate M R . In crossover, a pair of selected parents (using tournament selection mechanism) are crossed with the hope of creating better children. Despite crossover, mutation operator is employed on an individual solution in order to diversify the solution space and escape from the local optimum. These operators are employed to generate new solutions as a next offspring so that the current population and generated offspring are combined together for creating next generation. In this paper, we do not consider the elitism operator. GA can be used in various types of problems with appropriate genetic representation, fitness function, and evolution operators.
For solving the proposed RCLRP-SPD, we propose two efficient hybrid GA by hybridizing the GA with variable neighborhood descent (VND) (HGAVND) and the GA with a local search (HGALS) [4]. In addition, the proposed model aims at locating depots which are less sensitive to disruption and lead to the minimum cost no matter which disruption scenario occurs. Indeed, the final solution would be the most reliable solution when exposed to any disruption scenario. Therefore, the located depot(s) is (are) the same for all scenarios. Since in the RCLRP-SPD the open depots are the same for all scenarios, a two-stage mechanism is developed in the algorithms. The first and the second stages are called "Location-Allocation" stage and "Vehicle Routing" stage respectively. In the Location-Allocation stage, the depots are located and the customers are allocated to the opened depots. Once the first stage is solved, its results are fed to the second stage and the routing decisions are made in the second stage. Indeed, the routing is solved for each scenario independently. Besides to the proposed HGAVND and HGALS algorithms, we also employ classical GA to better show how these hybridizations improve the performance of the classical GA. Therefore, we consider three combinations of algorithms to solve the two "Location-Allocation" and "Vehicle Routing" stages as GA→GA, GA→HGAVND, and GA→HGALS. For instance, the combination GA→GA means that the GA is used for solving both stages; and the combination GA→HGAVND states that stage 1 is solved using GA and stage 2 is solved using HGAVND. As it can be seen, for solving the "Location-Allocation" stage (i.e., stage 1), we only use classical GA but the "Vehicle Routing" stage (i.e., stage 2) is solved using all three algorithms (i.e., GA, HGAVND, and HGALS).

Solution representation
Any evolutionary algorithm's performance remarkably depends on its representation of the solution. Since GA is a population-based algorithm, the solution representation should not consume much memory. The solution representation contains two different parts and each part corresponds to a stage of the proposed algorithms. These parts are described in the following.
(i) Location-Allocation: Figure 3a shows the location-allocation decisions in a network with 8 customers.
Numbers 1-8 represent the ycustomers, and each "D" shows a depot. The number of "D"s plus 1 demonstrates the number of located depots. In Figure 3, there are totally three located depots. Starting from the left, the customers are allocated to the first located depot until they will not exceed the capacity of the corresponding depot. For instance, customer number 1 has been allocated to the first depot. Then, customers number 3, 5, 8 and 4 should be allocated to the second located depot but only customers 3, 5 and 8 can be allocated since adding customer 4 will exceed the capacity of the second depot. Therefore, the demand of customer number 4 is unmet. This procedure continues until the most possible customers are allocated to the located depots. Figure 3b illustrates the scheme of the location-allocation decisions. (ii) Vehicle Routing: Figure 4a shows the solution representation regarding the vehicle assignments and routing.
Numbers 1-12 represent the customers. In this figure, each "V" signifies that the route is terminated and another route from the same depot is started, even though the accumulated demand has not exceeded vehicle capacity, while each "D" terminates the assignments of routes to the current depot and starts the   assignment to the next depot. Note that every "D" represents a new depot and a new vehicle. A vehicle has been assigned to the first and the third depots, but two vehicles have been assigned to the second located depot. The route of the third located depot serves customers 9, 11 and 3. Since adding customer 6 exceeds the vehicle capacity, the route is terminated. Note that the capacity of depots is not taken into consideration during the decoding process. Therefore, a per unit penalty cost P is added to the objective function value whenever the total demand served by a depot exceeds its capacity.

Initial population construction
Each stage of the proposed algorithms requires an initial population. The initial population of the first stage is generated in a greedy manner. We generate an initial population including greedy heuristic solution as well as a random solution. The Location-Allocation part of the solution representation is then constructed by a greedy method with the hope that a good initial solution can be found within a reasonable time. In this greedy algorithm, all potential depots are opened and each customer is allocated to the closest located depot according to an increasing order of D i + P i until the depot capacity is not exceeded; if exceeded, the next closest depot is selected. This greedy solution is then improved through the algorithm by closing unnecessarily opened depots. The population of the second stage is generated randomly.

Fitness function
The fitness function is defined by the solution representation and measures the quality of each solution. The fitness of each stage is calculated as follow: (i) Location-Allocation: in this stage, the fitness function includes depots fixed cost, direct distance cost between the depots and customers location, and penalty cost of unmet customers. Since we want to minimize the cost, lower value of the fitness function implies higher fitness. (ii) Vehicle Routing: in the second stage, the fitness function includes vehicles fixed cost, vehicles routing cost, constant P if the depot capacity is exceeded and the penalty cost of unmet customers. Since we want to minimize the cost, lower value of the fitness function implies higher fitness.

Crossover operator
To obtain a near -optimum solution, genetic evolutionary operators are used to create a better solution and replace them with those that existed in the initial population. Crossover is the GA's main operator, sharing information between two randomly selected parents in the hope of better offspring. These offspring are compared with respect to the fitness function and passed onto the next generation. In this paper, partially mapped crossover (PMX) is employed in both parts of "Location-Allocation" and "Vehicle Routing". The PMX indexes in the middle part of the solution representation are exchanged and the missing stops replace their matching index in the second parent. Stop indexes are exchanged based on the mapping set, which is constructed by comparing the exchanging parts of two parents [6]. This popular crossover method is illustrated in Figure 5. Note that each "V" and "D" can be replaced by numbers.

Mutation operator
Mutation is a strategy of diversification that aims to explore new areas of solution space, avoid early convergence, and escape from local optima. Through the mutation operator, three "Swap", "2-Opt-a", and "2-Opt-b" moves in the both "Location-Allocation" and "Vehicle Routing" parts. Using the "Swap" operator, two genes from the chromosome are selected and swapped. Using the "2-Opt-a" move, two points in the chromosome are selected and the order of the in-between genes is reversed. Using the "2-Opt-b" move, two points in the chromosome are selected and the order of the in-between alleles is randomly changed. These mutation operators are illustrated in Figure 6

Hybrid genetic algorithm
In this paper, we propose two hybrid GA by hybridizing the GA with VND (HGAVND) and the GA with a local search (HGALS). Both HGAVND and HGALS are employed to jointly solve routing in the "Vehicle Routing" part. Since the GA may fail to converge to a global optimum, we use VND and a local search to refine Figure 6. Examples of the "swap", "2-Opt-a" and "2-Opt-b" mutation operators.
the GA search through successive iterations and maximize the chance of convergence to an optimal solution [7,32,35,48].
The VND is a deterministic version of the variable neighborhood search (VNS) algorithm proposed by Mladenović and Hansen [35]. VND's basic idea is to explore a set of predefined neighborhood structures (N (l)(l = 1, 2, . . ., l max )) successively in order to obtain a better solution. It explores a set of neighborhoods systematically in order to obtain various local optima and escape from local optima. The overview of the proposed hybrid GA-VND (HGAVND) algorithm is depicted in Figure 7. In the proposed HGAVND, we incorporate VND into the GA's general scheme. This allows us to take advantage of VND features in order to improve the population generated by the GA and thus to complement the genetic search. In Figure 7, MaxIt GA is the maximum number of iterations of the GA (number of generations). In addition, p hybrid is a factor to calculate the probability of doing VND in the proposed HGAVND.
It is noteworthy that implementing the VND at each iteration of the GA increases the computational time. Therefore, the VND is involved in the GA with a probability equal to prob(i). This value increases by the number of iterations. Accordingly, we give more chance to the VND in the later iterations. Then, as we find the promising regions of the feasible solution space, we gradually increase the probability of employment of VND [48]. The pseudo code of the VND proposed for HGAVND algorithm is depicted in Figure 8, wherein Itr H is the maximum iteration for the VND.
In the VND used to complement the GA in our hybrid approach, we use seven neighborhoods (l max = 7) in both hybridization algorithm, namely N (1), N (2), N (3), N (4), N (5), N (6) and N (7), which are described in the following [14,21,32]: -N (1): randomly selects two customers assigned to two different routes and interchanges them. -N (7): randomly selects two sequential customers assigned to a route and adds them to the random position in another random route with same sequential.
HGALS is similar to HGAVND, with difference in local search. In HGAVND we obtain local optima in every neighborhood structure, but in HGALS, all the neighborhood structures are applied together and we obtain local optima in all neighborhood structures. The pseudo code of the local search proposed for HGALS algorithm   is depicted in Figure 9, wherein Itr H is the maximum number of iterations of the proposed local search in the HGALS.

Computational study
In this section, we present extensive computational results in order to assess the effectiveness of the proposed GA, HGAVND and HGALS algorithms for solving the RCLRP-SPD in two stages. The GA, HGAVND, and HGALS algorithms are coded using MATLAB R2014a and the mathematical model is solved with CPLEX 12.7.1.0 (with a time limit of 3 h) solver in GAMS 24.9.1 software on Intel Core TM i7-720QM (1.6 GHz, 4 GB RAM).

Parameter setting
The efficiency and effectiveness of all meta-heuristic algorithms are significantly influenced by its parameters. In this study, the parameters optimization of the proposed GA, HGAVND and HGALS are performed based on Taguchi methods. The experimental design of Taguchi has been widely used for optimization problems. Taguchi design is based on two major tools: orthogonal array (OA) and signal-to-noise (S/N) ratio. The OA is a matrix of numbers containing experimental schemes based on different levels of factors. The S/N ratio is the measure of variation and guarantees the robustness of this kind of experimental design. The term "signal" represents "mean response variable" as the desired value and "noise" represents "standard deviation" as undesirable value [7]. Taguchi method is utilized via Minitab 18 software for design of experiments (DOE) and analyzing their results. The aim of the Taguchi method is to maximize the S/N ratio which is calculated by equation (5.1) for minimization problems for each parameter i on its related level j.
where z ij is the objective function value using parameter i on level j and n is the number of repetitions for the corresponding experiment.
Since the proposed algorithms are executed in two stages, the parameters of each stage are tuned separately as follows: Table 2. Different levels of GA parameters in "Location-Allocation" parameters.  (i) Location-Allocation: since GA solves the first stage, different levels of the "Location-Allocation" parameters for tuning process is shown in Table 2. According to the S/N ratio plot shown in Figure 10, it is inferred that the number of generations, P size , C R and M R are tuned at 100, 200, 0.8 and 0.5, respectively. (ii) Vehicle Routing: since hybrid algorithms solve the second stage, different levels of the "Vehicle Routing" parameters for tuning process is shown in Table 3. According to the S/N ratio plot shown in Figure 11, it is inferred that the number of generation, population size, crossover rate, mutation rate, number of generation in hybridization algorithms and parameter that affects the probability of starting the hybridization algorithms, are tuned at 150, 200, 0.9, 0.4, 100 and 1.15, respectively.

Computational results
This section provides comprehensive numerical results to demonstrate the efficacy of the proposed MILP model and the solution approaches. Numerous instances are generated based on [64] and solved in this section. which are expressed as "a" to "e", respectively. These instances are labeled as |I| − |J| − |H|-a, where |I| is the number of customers, |J| is the number of potential depots, |H| is the number of vehicles, and finally "a" denotes the level of failure probability. Each instance of the problem is executed 30 times and the results, including the best, the worst, the average and the standard deviation (σ) values are reported afterwards.
The potential location of customers and depots are randomly generated in a continuous space between [1,10]. The cost between two nodes is calculated as the Euclidean distance between them. The delivery demand and pickup demand of each customer are randomly distributed in [10,100]. The unmet-demand penalty ( i is chosen from [100, 120]. The fixed setup cost of each depot and vehicle are drawn uniformly from [100, 500] and [50,100], respectively. We first calculate the C and D as: Each capacity cd j andcv are then drawn uniformly from [1.5C, 2.5C] and [D, 1.5D], respectively. In medium and large instances, all parameters are generated like small instances, except cd j , cv and θ i that are drawn uniformly from [1.5C, 1.8C], 2/3[D, 1.5D] and [50,100], respectively.

Analysis of small-size instances
Three sets of instances are solved with different failure probability from "a" to "e", where the failure probability increases from "a" to "e". Two number of potential depots are considered for small-size instances. Therefore, the number of all scenarios is then 2 2 and instances are solved in four scenarios. The results of CPLEX, GA, HGAVND, and HGALS have been provided in Table 4. Since for all 30 executions, the optimal solution is found, Table 4 contains only the best values (best, worst and average are the same and σ = 0). The results are compared in terms of objective function value, CPU time and the number of opened depots.
It can be seen that for all the small-size instances, all three GA, HGAVND, and HGALS reach the optimal solution; however, GA has lower CPU time since no extra local search operator is employed. In addition, the proposed HGALS leads to smoothly lower CPU time comparing to HGAVND. In addition, in the solution of instances 7-2-3-c,d,e, 8-2-3-e and 9-2-3-e, no depot has been opened since the failure probability is high and the model prefer not to allocate customers to unreliable depots. Two analyses are illustrated on the result of the CPLEX algorithm in Figure 12. Figure 12a illustrates that the total cost increases with increasing the failure probability, while there is no change in 7-2-3 ("c", "d" and "e") instances because in these instances the penalty cost is paid to all customers and no location and routing decisions are performed. Figure 12b connects the running times resulting from running the CPLEX on different sized instances, there is no growth in "e" instances because in these instances the penalty cost is paid to all customers and location and routing is not performed. Note that all the meta-heuristic algorithms have reached the optimal solution. In the majority of instances (11 out of 15), metaheuristic algorithms reach the optimal solution in significantly less CPU time comparing to CPLEX.

Analysis of medium instances
Four instances are solved with "a" failure probability because we mentioned in the introduction section that the occurrence probability of disruption is small. In addition, the same observation and results were obtained for medium and large instances comparing to small instances. Therefore, the results of medium and large instances only contain the failure probability "a" which is low and more close to reality. Furthermore, three potential depots are considered that leads to the total 2 3 number of different scenarios. Among them, we choose those scenarios that have the largest probability of occurrence. Here, instances are solved in four (|J|+1) scenarios that have the largest probability of occurrence. The results of the best, the worst, the average and σ values are given in Table 5 and CPLEX time limit is 10 800 s. The column "DIFF (%)" shows the gap between the metaheuristic cost and CPLEX cost and is calculated as DIFF = OFV meta-heuristic −OFVCPLEX

OFVCPLEX
. Figure 13a shows that only in the third instance (20-3-4-a), CPLEX has obtained a better result than all meta-heuristic algorithms. According to the Table 5 and "Open depot" column, we found out that CPLEX and solution algorithms have led to the different location of depots because in the solution algorithms, the "Location-Allocation" stage is solved first and this lead to a local optimum (this is a heuristic approach); but still closed-to-optimal solutions are found. Figure 13b connects the running times resulting from running the meta-heuristic algorithms on different instances and shows that hybrid algorithms require more time. Moreover, the results show that GA, HGAVND, and HGALS algorithms have satisfactory performance in terms of solution quality and CPU time.

Analysis of large instances and benchmark
Twelve instances with |J|+1 scenarios that have the largest probability of occurrence are solved in "a" failure probability and results are given in Table 6, in which the number and probability of scenarios are written below the instances as well. Figure 14 shows that HGAVND and HGALS algorithms significantly outperform GA in terms of solution quality. In the sixth instance (40-5-6-a), HGALS has obtained a disappointing result, but according to Table 6 and "Open depot" columns, we found out that they have a different location, because in the meta-heuristic algorithms, first the "Location-Allocation" is solved and this is a heuristic approach. Figure 14 shows that HGAVND is between the GA and HGALS in term of CPU time for small and medium instances. But this algorithm performs faster when dealing with large instances. Finally, the performance of the proposed algorithms is compared against the instances derived from Barreto's test set by Angelelli and Mansini's    separation approach (i.e., BAM) benchmark in the literature [22] as reported in Table 7. This comparison shows the superiority of the proposed algorithms compared to those of in the literature.

Analysis in the different number of scenarios
According to the previous section, we choose HGALS and all results are obtained from HGALS. The reason for selecting HGALS is that the HGALS obtains better solutions in terms of the objective function value and Table 6. Computational results of the large instances.   the computational time. In addition, the HGALS has had also higher convergence rate (i.e., lower standard deviation when being executed several times). Therefore, the analyses obtained from the HGALS algorithm would be more validated. Since the number of scenarios is different in Four instances are solved in "a" failure probability. Four potential depots are considered and the number of all scenarios is 2 4 , so we choose the different number of scenarios that have the largest probability of occurrence. The results are given in Table 8. Figure 16a shows that the smallest cost of each instance is in the one scenario, which does not take into account the disruptions. With increasing the number of scenarios, costs are rising, but there is not much difference between costs in 5, 11 and 16 scenarios. As shown in Figures 16a and 16b, we conclude that the number of scenarios that previously considered (|J| + 1) was correct and suitable in terms of time.
The last column "Cost increase (%)" in Table 8 is calculated as OFV16 scenarios−OFV1 scenario OFV1 scenario that OFV 1 scenario represents the total cost if disruptions are not assumed and the OFV 16 scenarios represent the total cost if all disruptions are assumed. The results show that percentage of increased costs are between one to six percent that each manager in the supply chain is willing to accept it, because not considering the disruptions, imposed the staggering cost on the system. Figures 17 and 18 illustrate the sensitivity of the total cost vs. the capacity of depots and vehicles, respectively. It can be seen in both Figures 17 and 18 that increasing the capacity will decreases the total cost. Actually, when the capacity of depots and vehicles increases, two event happen as: (1) less number of depots are opened and consequently lower fixed cost is imposed to the solution and (2) customers are allocated to closer depots and consequently lower transportation cost is paid in the solution.

Conclusion
This paper develops a new mixed-integer mathematical model for a reliable capacitated location-routing problem with simultaneous pickup and delivery, wherein the depots are exposed to the risk of disruption. Any disruption in the location-routing network imposes a huge cost on the company and its stockholders. Considering disruption in the design phase of the network will alleviate this impact and let the network to resist disruption. The objective function of the proposed model minimizes the total expected cost of the network including the sum of location cost of depots, routing cost of vehicles and cost of unfulfilled demand of customers. The proposed model is solved for small-sized instances and for large-sized problems, three efficient meta-heuristics are tailored and their performance is investigated through a comprehensive experimental research. All the meta-heuristics algorithms reach the optimal solution in the small instances and better solution comparing to CPLEX in the most cases in the medium instances, wherein the CPLEX is stopped after an enough huge computational time. The computational results demonstrate that the hybrid algorithms significantly outperform the classical genetic algorithm in terms of solution quality.
Possible future extensions include using multiple products, using time window for products, considering the location-routing-inventory problem and accounting for the cost of carbon emissions and developing other meta-heuristics for the problem. In addition, comparing the proposed hybrid meta-heuristic algorithms with other existing meta-heuristic algorithms such as Tabu Search (TS), Particle Swarm Optimization (PSO) and Differential Evolution (DE) algorithm would be an important future research direction.