Location-Routing Problem of Emergency Facilities under Uncertain Demand by Branch-Price and Cut

*is paper studies the location-routing problem of emergency facilities with time window under demand uncertainty.We propose a robust mathematical model in which uncertain requirements are represented by two forms: the support set defined by cardinal constraint set. When the demand value of rescue point changes in a given definition set, the model can ensure the feasibility of each line. We propose a branch and price cutting algorithm, whose pricing problem is a robust resource-constrained shortest path problem. In addition, we take the Wenchuan Earthquake as an example to verify the practicability of the method. *e robust model is simulated under different uncertainty levels and distributions and compared with the scheme obtained by the deterministic problem. *e results show that the robust model can run successfully and maintain its robustness, and the robust model provides better protection against demand uncertainty. In addition, we find that cost is more sensitive to uncertainty level than protection level, and our proposed model also allows controlling the robustness level of the solution by adjusting the protection level. In all experiments, the cost of robustness is that the routing cost increases by an average of 13.87%.


Introduction
With the deterioration of the environment, natural disasters such as fire, earthquake, and tsunami are also endangering people's survival and development. Because of their sudden and powerful destructive power, they often lead to house collapse, road damage, and casualties, resulting in a large number of life and medical material losses in the short term. e timely delivery of postdisaster relief materials is a major guarantee to increase the survival rate. Quickly and effectively transporting emergency relief materials from emergency rescue facilities to the disaster area and meeting the needs of the affected people is the most important task. Especially in earthquake-prone China, Japan, and coastal countries with common tsunamis, the establishment of emergency rescue facilities and the transportation of relief materials are very important. For example, the Wenchuan Earthquake, a major earthquake in China in 2008, caused an big number of casualties due to the lack of timely rescue and the failure of timely delivery of relief materials. e total affected population reached 46.256 million, resulting in a direct economic loss of 845.14 billion yuan. Nowadays, COVID-19 is also continuing to intensify, and emergency drug delivery is also crucial. erefore, in the face of natural disasters, it is an important issue for us to establish emergency rescue facilities and reasonably plan emergency rescue routes.
However, the unpredictability of natural disasters and the urgency of rescue work often make it difficult for emergency departments to accurately obtain the material needs of the disaster scene. erefore, under the constraints of limited time and resources, in order to reduce the time and cost of logistics distribution and improve the robustness of facility location, it is of great significance to study the location of uncertain emergency rescue facilities and the arrangement of rescue material distribution routes.
In fact, many decision-makers usually focus on the transportation of relief materials. On the contrary, if the robustness of emergency facilities is not improved, the transportation of relief materials cannot be realized efficiently. e location of supply station is related to the operational efficiency of the whole system and also has a great impact on the subsequent material distribution path planning. e planning of material distribution path directly determines the timeliness of rescue. erefore, it is necessary to integrate emergency resource location and path planning for disaster preparedness. In recent years, scholars at home and abroad have done a lot of research on emergency resource scheduling and have achieved some results.
Research on emergency facilities can be traced back to the problem of P center of gravity and P center proposed by Hakimi [1]. Overstreet et al. [2] expounded the complexity of emergency logistics planning and the next research direction after combing a large number of emergency logistics literatures; Galindo and Batta [3] sorted out a large number of literatures on emergency logistics modeling from 2005 to 2012, found that the most widely used method of emergency logistics planning is still mathematical planning, and focused on the practicability of the assumption part of the model; Oezdamar et al. [4] discussed the emergency logistics and emergency material allocation after natural disasters, and Horner Downs [5] studied the distribution of emergency relief materials combined with GIS; Kar and Hodgson [6] constructed a linear combination model for the location of emergency shelters; Widener and Horner [7] discussed the hierarchical location model of emergency relief materials. Furthermore, many scholars combined facility location with robust optimization methods. Robust optimization makes the optimization model robust to the uncertainty of input data by defining the uncertainty set of random parameters.
Some scholars have considered the location path problem of emergency facilities under certain information with the goal of minimizing the total cost [8,9]. Some scholars have likewise studied the location path problem of emergency facilities under uncertain information. Considering the research of uncertain demand, Caunhye et al. [10] established a two-stage LRP stochastic programming model with the goal of minimizing facility construction cost and minimizing rescue time on the worst path. Bodaghi et al. [11] studied the problem of providing multiresource scheduling and routing during disasters. Considering the study of time uncertainty, Ahmadi et al. [12] constructed a model with the goal of minimizing the sum of total rescue time, penalty cost of unmet demand, and fixed cost of opening facilities LRP optimization model of scenario based mixed nonlinear programming. Zarandi et al. [13] expressed the uncertainty of time with fuzzy numbers and established a fuzzy chanceconstrained programming model with the goal of minimizing the total cost. Vahdani et al. [14] considered the uncertainty of distribution center capacity and warehouse capacity and established a multicycle and multimaterial three-level relief chain location path model with the goal of minimizing emergency rescue time and maximizing path reliability. Moreno et al. [15] assumed that the demand, supply, and road capacity were uncertain and established a scenario based two-stage stochastic programming model with the goal of minimizing the cost of emergency logistics and the cost of unmet demand loss of victims at the disaster site. From these literatures, we can conclude that there are various uncertainties in practical problems, and it is necessary to consider robust optimization methods to solve these problems. In the literature, their main methods are described by random programming, opportunity constraints, and other methods, or when using robust optimization, only one parameter is simply considered for an uncertain factor, and the use of two parameters to describe the demand uncertainty from different angles is not involved. erefore, this paper considers the variation of two parameters, that is, a cardinality constraint support to describe its demand uncertainty.
In this paper, we consider the interval uncertainty of demand; that is, the demand of each rescue point changes in a given interval. Rescue points with uncertain needs are also limited by uncertain budgets. Under these two constraints, the demand of rescue point defines uncertain support set, which is called cardinality constraint support. Uncertainty shall be applied to each rescue route to prevent the worst case. In addition, like the deterministic emergency facility location path, it is impossible to solve the robust problem with linear optimization. erefore, we use a column generation method to decompose the robust constraints into subproblems. We develop a branch-price and cut algorithm, in which the subproblem is a robust basic shortest path problem with resource constraints. In this paper, we propose a new decomposition strategy for the robust case. e method is applied to a practical case (Wenchuan Earthquake) to verify the solution in the robust case.

Contributions of is Paper
(1) Two uncertain support sets are used to define demand uncertainty (2) A branch-price and cut algorithm is proposed to solve the proposed model (3) A real case (Wenchuan Earthquake) is used to verify the effectiveness and applicability of the algorithm e rest of the paper is organized as follows: Section 2 defines the robust emergency facility location-routing problem with uncertain demand, in which the uncertainty support is based on the cardinal constraint set. In Section 3, the column generation framework and the solution to the subproblem are described in detail, and the branch-price and cut algorithm is given. In Section 4, we test our model and solution and give the comparison of robustness with the solution of deterministic model.

Robust Model
In this paper, based on the work of omas et al. [16], a new robust routing formula for disaster relief facilities (R-LRP-DR) is proposed, considering the uncertainty of the demand for rescue materials and the limitation of delivery time window.

Robust Formulation for R-LRP-DR.
We can define R-LRP-DR as a directed graph G � (V, A), and V � DA ∪ CF { } is the set of nodes and A is the set of arcs consisting of all feasible links, respectively. S � 1, 2, . . . , n { } represents sets of n suppliers, and R � n + 1, n + 2, . . . , 2n { } represents sets of n retailers. e set DA � 1, 2, . . . , M { } represents the alternative location to build some CDs, and the appropriate CDs need to be built throughout the pickup and delivery process based on the cooperating vendor. At the same time, the same types of vehicles are used in pickup and delivering to ensure its consistency, and the set is defined as K. e capacity of each vehicle is T. In the process of delivery, each customer i ∈ C has its corresponding time window [a i , b i ] and service time s i . Every customer i ∈ C has a demand q i . Here we consider that q i is uncertain. is uncertain demand parameter q i is randomly selected from the interval [q i , q i +h i ], where q i ≥ 0 denotes the nominal demand of customer i, and h i ≥ 0 denotes the maximum deviation between actual demand and nominal demand. Without losing generality, let us assume that And, Γ describes the magnitude of the deviation change. erefore, at this time, the demand uncertainty can be characterized by h i and Γ; that is, the uncertainty can be realized by adjusting the interval range of h i or value of Γ. Whether in the process of pickup or delivery, the vehicle has a cost c ij on arc ij. Among them, the cost in the pickup and delivery process is recorded as c ij . In the process of pickup and delivery, vehicles must start from it and return to it at every stage.
Robust LRP-DR (R-LRP-DR) aims to select a set of cross-docking systems and the vehicle routes associated with these cross-docking systems to minimize the total cost, including the fixed cost of establishing cross-docking systems and the cost of serving retailers. In PDPDC system, when the retailers' demand order is received, a series of homogeneous vehicles K are taken from supplier S and delivered to retailer R through cross-docking system. It should be noted that cargoes must pass through CDs and must not be allowed to be transferred directly between vehicles.
In order to clearly describe the demand uncertainty of retailers, we propose an R-LRP-PDPDC formula as follows: where objective function (1) is to minimize the cost of the whole process, including the opening cost of emergency facilities and the transportation cost of products; constraint condition (2) is the constraint on the number of open emergency facilities; in constraint (3), there is only one path through each facility; constraint (4) means that vehicles must enter and leave the supplier once in the process of picking up and delivery, which ensures the balance of traffic flow; constraint (5) means that the customer's demand should not exceed the vehicle capacity in the whole process; constraint (6) means that the service is started within the time window, whether it is pickup or delivery. By Desrochers and Laporte [17], we can define M � max i,j∈V b i − a i + t i . Constraint (7) ensures that delivery must be within the time window.
Constraints (8) and (9) are decision variables. If vehicle K goes from node i to node j in pickup process, y i,j,k ′ equals 1, and, otherwise, y i,j,k ′ equals 0. Alike, if vehicle K goes from node i to node j in delivery process, y ′ ′ i,j,k equals 1, and, otherwise, y ′ ′ i,j,k equals 0. We note that since demand q is uncertain, which belongs to the support set Q, we can further write constraint (10) as where y k ∈ 0, 1 { } |A 2 | refers to the vector set of decision variables y ′ ′ i,j,k and β(y k , where V C is a subset of customers which does not exceed Γ, so that customers in V C can be accessed by k vehicles, and customers demand deviations are accumulated and maximized. Such β(y k , Γ) ensures that the load of vehicle k can cover any demand realization and does not exceed the maximum demand of the customer.

Set Partitioning Formulation
We note that, during transportation, each vehicle meets constraint (10), even if the customer's requirements are uncertain. Here, we use a modified version of the LRSP set partition formula proposed by Akc [18]. It is considered that Ω is a feasible routing set. ese paths are related to emergency facilities. ey are transported from the supplier to the emergency facilities and then from the emergency facilities to the customers; and p j (p j ∈ Ω) is the set of routes Journal of Mathematics assigned by vehicles to emergency facilities j. en, the following parameters and decision variables are given:

Parameters:
c p � all cost of route p ∈ Ω, including transport cost and operate cost, Decision variables: e robust LRP-PDP formulation is as follows: Objective function (13) seeks to minimize the total cost, including the opening cost and transportation cost of emergency facilities. Constraint (14) ensures that each retailer location is served by only one route. Limiting condition (15) ensures the number of emergency facilities. Constraints (16) and (17) are standard binary constraints on variables.
e linear relaxation of the [P] problem is the Dantzig-Wolfe (DW) master problem, and we define a set Ω ∈ Ω so that we get the constrained DW master problem.
[P] and DW main problems have a large number of variables, and the number of possible paths is usually so large that it is impractical to explicitly enumerate all the corresponding columns. erefore, in this process, we can use the column generation (CG) method [19], which considers all possible routes through implicit expressions M(Ω). Otherwise, we solve a subproblem to find one or more routes with negative reduced costs. Let π i and μ i be the dual variables associated with constraints (15) and (16). e subproblem is an R-LRP-PDP as follows: e solution of R-ESPPRC is elementary path with minimum reduced cost. When we release the base constraint, α ip expresses the number of times path p visits customer i. e restricted DW main problem will produce a weak lower bound; R-ESPPRC differs from the deterministic ESPPRC in that constraint (18) includes internal maximization.
e shortest path problem with resource constraints is usually solved by a label-setting algorithm. Given the robustness, Lu and Gzara [20] proposed an improved label definition algorithm based on new dominance rules. Alves Pessoa [21] proposed a different label definition algorithm.
e algorithm provides a virtual solution for each initial 4 Journal of Mathematics solution; each virtual solution changes within a certain range and uses dominance rules to constantly adjust the simulated solution. Another option is to examine a series of deterministic auxiliary problems. Bertsimas and Sim [22] first proposed this solution strategy for robust combinatorial problems with uncertain objective coefficients. Alvarez Miranda et al. [23] and Goetzmann et al. [24] extended this strategy to combinatorial optimization problems with uncertain constraints. In this article, we use the strategy promoted by Alvarez Miranda et al. [23] and Goetzmann et al. [24] to convert R-ESPPRC to a series of deterministic ESPPRC and use the classical label definition algorithm to solve each. For simplicity, we first omit the vehicle index K and set Y as the feasible set satisfying constraints (17) and (19)- (22). Let C be the vector of c ij � c ij − i∈V α ij (π i − μ i ) and let Y be the vector of Y ij . Define y i � i:(i,j)∈A,j≠0 y ij for j ∈ V and h n+1 � 0. y and q are vectors of x i and d i , respectively. With w i ∈ [0, 1], ∀i ∈ V, R-ESPPRC is rewritten as Let U(y) � max 0≤w i ,∀i∈V,w i ≤ Γ i∈V h i y i w i ≤ T, let σ 0 be the dual variable of constraint i∈V w i ≤ Γ, and let σ i be the dual variable of constraint 0 ≤ w i ≤ 1, ∀i ∈ V. erefore, we can write U as follows: en, formula (26) [P'] can be rewritten as follows: Since (27) is used instead of maximizing the internal objective function U, we get constraint (34). At this time, constraints (35) and (36) are the above constraints (28) and (29). Now we can directly use the MIP solver to solve this problem. However, in order to speed up the column generation process, we consider turning the uncertainty into certainty and then choose to use dynamic programming to solve it, while generating multiple cost reduced columns.
First of all, we extend the theorem in the work of Bertsimas and Sim [22] to the case of uncertain constraint coefficients, and then, according to Alvarez Miranda et al. [23] and Goetzmann et al. [24], provide an alternative sequential approach to solve this series of N + 1 problems. By constraint (28), we can get the following: is is because of the fact that σ i ≥ h i y i − σ 0 , y i is binary and σ i ≥ 0. en, we can rewrite formula (27) as Now we need to discuss it in different situations. Since it has been assumed that h 1 , for i � 2, . . . , n + 1, max h j − σ 0 , 0 y j � h j − σ 0 y j , j � 1, . . . , i − 1. erefore, we can get the following results by sorting out: As a result, U(y) can be decomposed into n + 1 subproblems:

Journal of Mathematics
We have that n+1 i�1 max h i − σ 0 , 0 y i � 0, and, in this case, when σ 0 � h 1 , we can get the optimal solution of U(y). When i � 2, 3, . . . n + 1, Equation (42) can be written as follows: It can be seen that, for σ 0 , (44) is a first-order linear function, and σ 0 is in the interval [h i , h i−1 ]; then optimal solution for U i (y) is obtained when σ 0 � h i or σ 0 � h i−1 . So, To sum up, we can obtain a piecewise function U as follows: en, the objective function Z can be written as x ∈ X, (50) Since each subproblem is independent of the others, we set up a secondary network for each warehouse x. Each subproblem can be considered as a resource-constrained shortest path problem, and here we present the label-setting algorithm required to solve each subproblem. See Section 5 for details.

Preprocessing.
In the previous literature, there are some preprocessing schemes for subproblems. is study considers the rules corresponding to heuristic extension proposed by Desrochers, Desrosiers, and Solomon [25].
First of all, we consider the time windows of partial paths 0 ⟶ i ⟶ N + i and i ⟶ N + i ⟶ 2N + 1. Since the constraint t i ∈ [a i , b i ] of the time window must be satisfied in the whole process, we will redefine the time window as follows (see Dumas et al., 1991): (53) Due to the time window and constraints, some arcs that do not satisfy the feasible solution of the problem can be deleted. According to Dumas et al. [26], the following arcs are eliminated by using the constraints of the problem: [Priority]: none of these arcs directly cross-docking are feasible solutions, so they should be deleted: (0, N + i); (N + i, i); (2N + 1, 0); (2N + 1, i); (i, 2N + i); and (2N + 1, N + i) for i � 1, . . . , N. [Time window]: the customer's node must serve in its own time window; that is, if a i + t ij > b j , i, j ∈ 1, 2N { }, then delete arc (i, j).
[Capacity]: the capacity of the vehicle can never be exceeded; therefore, if q i + q j > T, i, j ∈ 1, N { }, i ≠ j, then remove the following arcs: [Time window with priority]: in the pickup and delivery, the travel time needs to meet the triangle inequality. If the arc does not meet this criterion (including both pickup and delivery of some customers), the arc can be removed:

Branching Rule.
At the end of the separation algorithm or in the best solution, the variable may contain a noninteger value, at which point the branch must be performed. Since there are two different binary decision variables in the formula, we adopt the following order in the branching phase: First, the algorithm attempts to branch over the most nearly 0.5 fractional Y k variable to generate two subproblems, one of which is fixed Y k0 , and the other is fixed Y k1 . If all Y k variables are integers, then the algorithm branches off from the minimum number of X ij variables.

Label-Setting Algorithm.
e objective of the pricing problem is to find a tourism tree with the greatest reduction in negative costs. Considering that the tree picks up the goods from a warehouse through the pickup point, then delivers to the pickup point and ends up in the warehouse, it is reasonable to look for it by considering each pickup process. To solve this problem, we propose an improved labeling algorithm, which was proposed by Frias and Singler [27] for resource-constrained shortest path problem. A twoway search is performed, where tags extend forward and backward. In the front, the tag extends from the warehouse to a pickup point as the first access node on the travel tree. Similarly, in the reverse direction, the label extends from the warehouse to the point of delivery and becomes the node for the last visit. Forward and reverse tags merge to form a complete travel tree. Both breadth-first and depth-first search methods are used when extending labels.

Instances Description.
is section reports a set of calculation results of the disaster related to the "5.12″ Wenchuan Earthquake. Under the existing demand, the simulation generates the relevant time window. We compare the difference between the robust solution and the deterministic solution in terms of cost, the number of vehicles used, the average percentage of infeasible scenarios, and the average capacity utilization. All algorithms are coded in Java, and the MIP and LP models are solved by CPLEX 12.5. e test is carried out on an ASUS computer with a processor of 1.8 GHz and 8 GB ram. If the solution process does not terminate within 3 hours, the program terminates.
Based on the disaster information of "5.12″ Wenchuan Earthquake, 6 candidate emergency facilities, 12 disaster sites (as shown in Figure 1), and 7 transport vehicles of the same type are selected to provide food and tent distribution services (assuming the same volume of the two types of materials). e fixed use cost of each vehicle is 1500 yuan, the transportation cost per kilometer is 5 yuan, and the vehicle capacity is 40 000. In Table 1, we give the location and capacity of the six potential emergency facilities and the different transportation cost after departure from each emergency facility. We give the nominal demand of each disaster site in Table 2, and, on the basis of the nominal values, the demand fluctuation h i at the affected point uniformly generates high uncertainty in [0, 0.5d i ] and low uncertainty in [0, 0.3d i ]. At this time, we describe the impact of the maximum deviation h i on the demand and only change its value to describe the magnitude of the demand deviation, resulting in demand uncertainty.

Numerical Experiment Results.
e following sections introduce the results of a series of experiments we conducted. e purpose of the experiment is to evaluate how a series of factors affect our proposed R-LRP-DR solution. By analyzing the influence of these factors, whether a certain factor in the uncertainty or the data change of a certain factor can produce a higher-risk solution is judged. e factors studied include the uncertainty of demand, the size of the customer's size, and some key parameters of the relevant uncertainty set.

Determining Model
Results. Firstly, we solve the deterministic model and give its optimal solution, as shown in Figure 2.
It can be seen from Figure 2 that when the material demand is determined (nominally), the facility opening points are A, D, E, and F, where facility A provides materials to demand points 4, 5, 10, 11, and 12; facility D provides materials to demand points 6 and 8; facility E provides materials to demand point 3; and facility F provides materials to demand points 1, 2, 7, and 9. At this time, the total cost � 72 737.5 yuan.

Robustness of Robust and Deterministic Solutions.
In order to understand the influence of uncertainty on the optimal solution, we compare the optimal target values of robust solution and deterministic solution. Firstly, we get the deterministic optimal solution under the nominal demand in this case (as shown in Figure 2). en, by adjusting the uncertain demand and uncertainty set, we get Figures 3 and 4. In Table 3, we give the difference of the uncertainty set and the average growth of the target value of the robust solution and the poor average optimality under the uncertainty level.
Next, we observe the robustness of the system by considering the demand for low uncertainty and the uncertainty concentration Γ � 3.
As can be seen from Figure 3, when considering the demand with low uncertainty and the uncertainty set parameter Γ � 3, there are still four facility open points, but they are changed to A, C, D, and F, in which facility A provides materials for demand points 4, 5, 11, and 12; facility C provides materials for demand points 3 and 7; facility D provides materials for demand points 6 and 8; and facility F provides materials for demand points 1, 2, 9, and 10. At this time, the total cost � 73 931.5 yuan. Figure 4 shows the opening of facilities when the demand is with high uncertainty and the uncertainty set parameter Γ � 3, and the cost � 75 046.3 yuan. Table 3 shows the difference of uncertainty sets and the average growth of the target value and average optimality of the robust solution under the uncertainty level.   Finally, we give the vehicle routing at high uncertainty level and Γ � 5, as shown in Figure 5. Table 3 reports average statistics for robust and deterministic solutions. e data represent the difference of uncertainty set, the average growth of target value, and the growth of average time. From the above chart, we can see that, under the same uncertain level set parameters, when the uncertain demand increases, the facility location will change and the cost will increase. Under the same uncertain demand, with the increase of uncertain level set parameter t, the cost increases gradually, and the facility location may change. As can be seen from Table 3, when the uncertainty level is low, even if Γ � 5, the cost increase is also less than the high uncertainty level when Γ � 3. is indicates that cost is more sensitive to the level of uncertainty than the level of protection. In all experiments, the cost of robustness is that the routing cost increases by an average of 13.87%. It is worth noting that our proposed model allows the robustness     level of the solution to be controlled by adjusting the protection level.

Conclusion
e occurrence of natural disasters will lead to serious economic losses, so the study of emergency facility location path problem is a hot spot at present. In this paper, we propose a robust model for facility location-vehiclerouting problem with uncertain demand. e motivation of this problem is the emergency material dispatching route under natural disasters. When encountering natural disasters, it is often affected by environmental factors, and its demand is usually random, which is difficult to estimate by probability distribution. erefore, robust optimization provides a more appropriate method and framework. In this paper, through formulaic processing and decomposition, we can use the column generation method to solve it.
We use cardinal constraint set to define uncertainty support and give a robust location path formula model. en, we propose a branch-price and cut algorithm, which uses the label algorithm in the cutting process, in order to generate multiple columns in each iteration of column generation, rather than a single column. Finally, we use the actual data (Wenchuan Earthquake) to verify the practicability of the model. e experimental results show that, on the one hand, the robust model solves the difficulty that the demand uncertainty problem in the linear optimization model can not be solved; on the other hand, the robust model makes up for the disadvantage that stochastic programming depends on accurate probability distribution; and the robust model still maintains certain robustness. In addition, the branching and price cutting algorithms successfully solve the optimal problem of the case in a reasonable time. We are surprised to find that when the uncertainty level is low, even Γ � 5. e cost increase is also less than the high uncertainty level with Γ � 3.
is indicates that cost is more sensitive to the level of uncertainty than the level of protection. In all experiments, the cost of robustness is that the routing cost increases by an average of 13.87%. Note that our proposed model allows the robustness level of the solution to be controlled by adjusting the protection level. e location-routing problem of emergency facilities is a hot spot in the current research. At the same time, this research application also helps to lay a theoretical foundation for emergency management and sustainable development. Although this study considers the main factor of demand uncertainty, there are still some deficiencies, for example, the application of time and cost uncertainty and facility interruption risk in emergency facility location routing. Future research can start from these aspects.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e author declares that there are no conflicts of interest.