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

Authors: Xuchen Deng
Journal: Journal of Mathematics
 (2021)

---

## Abstract

This 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. The robust model is simulated under different uncertainty levels and distributions and compared with the scheme obtained by the deterministic problem. The 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%.


---

## Body

## 1. 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. The 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. The 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. Therefore, 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. Therefore, 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. The 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. The planning of material distribution path directly determines the timeliness of rescue. Therefore, 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 chance-constrained 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. Therefore, 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. Therefore, 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. The method is applied to a practical case (Wenchuan Earthquake) to verify the solution in the robust case.


### 1.1. Contributions of This 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 algorithmThe rest of the paper is organized as follows: Section2 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.


## 1.1. Contributions of This 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 algorithmThe rest of the paper is organized as follows: Section2 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.


## 2. Robust Model

In this paper, based on the work of Thomas 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.


### 2.1. Robust Formulation for R-LRP-DR

We can define R-LRP-DR as a directed graphG=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. The 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. The capacity of each vehicle is T. In the process of delivery, each customer i∈C has its corresponding time window ai,bi and service time si. Every customer i∈C has a demand qi. Here we consider that qi is uncertain. This uncertain demand parameter qi is randomly selected from the interval qi¯,qi¯+hi, where qi¯≥0 denotes the nominal demand of customer i, and hi≥0 denotes the maximum deviation between actual demand and nominal demand. Without losing generality, let us assume that h1≥h2≥…≥hn. Therefore, we define an uncertain support set Q=qi=qi¯+hiwi,q∈RC,∀i∈C as customer demands, where wi∈0,1, and ∑i=1nwi≤Γ. And, Γ describes the magnitude of the deviation change. Therefore, at this time, the demand uncertainty can be characterized by hi and Γ; that is, the uncertainty can be realized by adjusting the interval range of hi or value of Γ. Whether in the process of pickup or delivery, the vehicle has a cost cij on arc ij. Among them, the cost in the pickup and delivery process is recorded as cij. 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 vehiclesK 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:(1)min∑j∈CDfjXj+∑k∈K∑ij∈Acijyijk,(2)s.t.∑i,j∈Ayi,j,k≤Xj,(3)∑i,j∈Ayi,j,k=1,∀i∈V,k∈K,(4)∑i,j∈Ayi,j,k=∑i,j∈Ayj,i,k,∀j∈V,k∈K,(5)∑i,j∈Adiyi,j,k≤T,∀d∈D,k∈K,(6)ti,k+si+ti,j−tj,k≤M1−yi,j,k,∀i,j∈A,k∈K,(7)ai≤ti,k≤bi,∀i∈V,(8)Xi∈0,1,∀i∈VD,(9)yi,j,k∈0,1,∀i,j∈VS∪V0,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=maxi,j∈Vbi−ai+ti. 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, yi,j,k′ equals 1, and, otherwise, yi,j,k′ equals 0. Alike, if vehicle K goes from node i to node j in delivery process, yi,j,k′′ equals 1, and, otherwise, yi,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(10)∑i,j∈A2qi¯yi,j,k+βyk,Γ≤T,∀q∈Q,k∈K,where yk∈0,1A2 refers to the vector set of decision variables yi,j,k′′ and βyk,Γ=maxVC¯|VC¯∈VC,|VC¯|≤Γ∑i∈VC¯wi∑j∈A2yi,j,k′′, where VC is a subset of customers which does not exceed Γ, so that customers in VC can be accessed by k vehicles, and customers demand deviations are accumulated and maximized. Such βyk,Γ ensures that the load of vehicle k can cover any demand realization and does not exceed the maximum demand of the customer.


## 2.1. Robust Formulation for R-LRP-DR

We can define R-LRP-DR as a directed graphG=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. The 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. The capacity of each vehicle is T. In the process of delivery, each customer i∈C has its corresponding time window ai,bi and service time si. Every customer i∈C has a demand qi. Here we consider that qi is uncertain. This uncertain demand parameter qi is randomly selected from the interval qi¯,qi¯+hi, where qi¯≥0 denotes the nominal demand of customer i, and hi≥0 denotes the maximum deviation between actual demand and nominal demand. Without losing generality, let us assume that h1≥h2≥…≥hn. Therefore, we define an uncertain support set Q=qi=qi¯+hiwi,q∈RC,∀i∈C as customer demands, where wi∈0,1, and ∑i=1nwi≤Γ. And, Γ describes the magnitude of the deviation change. Therefore, at this time, the demand uncertainty can be characterized by hi and Γ; that is, the uncertainty can be realized by adjusting the interval range of hi or value of Γ. Whether in the process of pickup or delivery, the vehicle has a cost cij on arc ij. Among them, the cost in the pickup and delivery process is recorded as cij. 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 vehiclesK 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:(1)min∑j∈CDfjXj+∑k∈K∑ij∈Acijyijk,(2)s.t.∑i,j∈Ayi,j,k≤Xj,(3)∑i,j∈Ayi,j,k=1,∀i∈V,k∈K,(4)∑i,j∈Ayi,j,k=∑i,j∈Ayj,i,k,∀j∈V,k∈K,(5)∑i,j∈Adiyi,j,k≤T,∀d∈D,k∈K,(6)ti,k+si+ti,j−tj,k≤M1−yi,j,k,∀i,j∈A,k∈K,(7)ai≤ti,k≤bi,∀i∈V,(8)Xi∈0,1,∀i∈VD,(9)yi,j,k∈0,1,∀i,j∈VS∪V0,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=maxi,j∈Vbi−ai+ti. 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, yi,j,k′ equals 1, and, otherwise, yi,j,k′ equals 0. Alike, if vehicle K goes from node i to node j in delivery process, yi,j,k′′ equals 1, and, otherwise, yi,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(10)∑i,j∈A2qi¯yi,j,k+βyk,Γ≤T,∀q∈Q,k∈K,where yk∈0,1A2 refers to the vector set of decision variables yi,j,k′′ and βyk,Γ=maxVC¯|VC¯∈VC,|VC¯|≤Γ∑i∈VC¯wi∑j∈A2yi,j,k′′, where VC is a subset of customers which does not exceed Γ, so that customers in VC can be accessed by k vehicles, and customers demand deviations are accumulated and maximized. Such βyk,Γ ensures that the load of vehicle k can cover any demand realization and does not exceed the maximum demand of the customer.


## 3. 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. These paths are related to emergency facilities. They are transported from the supplier to the emergency facilities and then from the emergency facilities to the customers; and pjpj∈Ω is the set of routes assigned by vehicles to emergency facilities j. Then, the following parameters and decision variables are given:Parameters:(11)cp=all cost of routep∈Ω, including transport cost and operate cost,αip=1,if customeri∈Ris assigned to routep∈Ω, text and,0,otherwise.Decision variables:(12)θp=1,if routepis selected in the solution, and,0,otherwise,Xj=1,if emergency facilitiesjis selected and used, and,0,otherwise.The robust LRP-PDP formulation is as follows:(13)Pmin∑j∈CDfjXj+∑p∈Ωcpθp,(14)s.t.∑p∈Ωαipθp=1,∀i∈VC,(15)Xj−∑p∈Ωθp≥0,∀j∈CD,p∈Ω,(16)Xj=0,1,∀j∈CD,(17)θp=0,1,∀p∈Ω.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.The 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. Therefore, 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). The subproblem is an R-LRP-PDP as follows:(18)min∑j∈CDfjXj+∑ij∈Acijyijk−∑i∈V∑ij∈Aαijkπi−μijyijk,(19)s.t.∑i,j∈Ayi,j,k≤Xj,(20)∑i,j∈Ayi,j,k=∑i,j∈Ayj,i,k,∀j∈V,k∈K,(21)∑i,j∈A2qi¯yi,j,k+βyk,Γ≤T,∀q∈Q,(22)ti,k+si+ti,j−tj,k≤M1−yi,j,k,∀i,j∈A,k∈K,(23)ai≤ti,k≤bi,∀i∈V,(24)Xi∈0,1,∀i∈VD,(25)yi,j,k∈0,1,∀i,j∈VS∪V0.The 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. The restricted DW main problem will produce a weak lower bound; R-ESPPRC differs from the deterministic ESPPRC in that constraint (18) includes internal maximization.The 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. The algorithm provides a virtual solution for each initial 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 indexK and set Y as the feasible set satisfying constraints (17) and (19)–(22). Let C˜ be the vector of cij˜=cij−∑i∈Vαijπi−μi and let Y be the vector of Yij. Define yi=∑i:i,j∈A,j≠0yij for j∈V and hn+1=0. y and q are vectors of xi and di, respectively. With wi∈0,1, ∀i∈V, R-ESPPRC is rewritten as(26)P′minFx+cys.t.x∈Xy∈Yy≤xqTy+max0≤wi,∀i∈V,wi≤Γ∑i∈Vhixiwi≤T.LetUy=max0≤wi,∀i∈V,wi≤Γ∑i∈Vhiyiwi≤T, let σ0 be the dual variable of constraint ∑i∈Vwi≤Γ, and let σi be the dual variable of constraint 0≤wi≤1,∀i∈V. Therefore, we can write U as follows:(27)Uy=minΓσ0+∑i∈VCσi,(28)s.t.σ0+σi≥hiyi,∀i∈VC,(29)σ0,σi≥0,∀i∈VC.Then, formula (26) [P’] can be rewritten as follows:(30)p″minFx+cy,(31)s.t.x∈X,(32)y∈Y,(33)y≤x,(34)qTy+Γσ0+∑i∈VCσi≤T,(35)σ0+σi≥hiyi,∀i∈VC,(36)σ0,σi≥0,∀i∈VC.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: σi=maxhi−σ0,0yi,∀i∈Vc. This is because of the fact that σi≥hiyi−σ0,yi is binary and σi≥0. Then, we can rewrite formula (27) as(37)Uy=minΓσ0+∑i∈VCmaxhi−σ0,0yi,(38)s.t.σ0≥0.Now we need to discuss it in different situations. Since it has been assumed thath1≥h2≥⋯≥hn≥hn+1=0, if σ0∈h1,+∞,maxhi−σ0,0yi=0,i=1,…,n+1; and if σ0∈hi,hi−1, for i=2,…,n+1,maxhj−σ0,0yj=hj−σ0yj,j=1,…,i−1. Therefore, we can get the following results by sorting out:(39)∑i=1nmaxhi−σ0,0xi=∑j=1i−1hj−σ0yj,ifσ0∈hi,hi−1,i=2,…,n+1,0,ifσ0∈h1,+∞.As a result,Uy can be decomposed into n+1 subproblems:(40)U1y=minΓσ0,(41)s.t.σ0∈h1,+∞.We have that ∑i=1n+1maxhi−σ0,0yi=0, and, in this case, when σ0=h1, we can get the optimal solution of Uy. When i=2,3,…n+1,(42)Uiy=minΓσ0+∑j=1i−1hj−σ0yj,(43)s.t.σ0∈hi,hi−1.Equation (42) can be written as follows:(44)Uiy=minΓ−∑j=1i−1yjσ0+∑j=1i−1hjyj,(45)s.t.σ0∈hi,hi−1.It can be seen that, forσ0, (44) is a first-order linear function, and σ0 is in the interval hi,hi−1; then optimal solution for Uiy is obtained when σ0=hi or σ0=hi−1. So,(46)Uiy=minΓhi−1+∑j=1i−1hj−hi−1yj,Γhi+∑j=1ihj−hiyj.To sum up, we can obtain a piecewise functionU as follows:(47)Uiy=Γh1,ifi=1,Γhi+∑j=1ihj−hiyj,ifi=2,…,n,∑j=1nhjyj,ifi=n+1.Then, the objective functionZ can be written as(48)Z=minFx+c,(49)s.t.qTy+Uiy≤T,i=1,2,…,n+1,(50)x∈X,(51)y∈Y,(52)y≤x.Since each subproblem is independent of the others, we set up a secondary network for each warehousex. 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.


### 3.1. Preprocessing

In the previous literature, there are some preprocessing schemes for subproblems. This 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 paths0⟶i⟶N+i and i⟶N+i⟶2N+1. Since the constraint ti∈ai,bi 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)bN+i=minbN+i,b2N+i−ti,2N+i,bi=minbi,bN+i−ti,N+i,ai=maxai,a0+t0,i,aN+i=maxaN+i,ai+ti,N+i.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, ifai+tij>bj,i,j∈1,2N, then delete arc i,j.[Capacity]: the capacity of the vehicle can never be exceeded; therefore, ifqi+qj>T,i,j∈1,N,i≠j, then remove the following arcs: i,j;j,i;i,N+j;j,N+i;N+i,N+j;N+j,N+i.[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:If pathj⟶i⟶N+j⟶N+i is not feasible for tj=aj, then arc i,N+j is deletedIf pathi⟶N+i⟶j⟶N+j is not feasible for ti=ai, then arc N+i,j is deletedIf pathi⟶j⟶N+i⟶N+j is not feasible for ti=ai, then arc i,j is deletedIf pathi⟶j⟶N+i⟶N+j and j⟶i⟶N+i⟶N+i⟶N+j⟶N+j is not feasible for ti=ai and tj=aj, remove arc N+i,N+j


### 3.2. 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 fractionalYk variable to generate two subproblems, one of which is fixed Yk0, and the other is fixed Yk1. If all Yk variables are integers, then the algorithm branches off from the minimum number of Xij variables.


### 3.3. Label-Setting Algorithm

The 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 two-way 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.


## 3.1. Preprocessing

In the previous literature, there are some preprocessing schemes for subproblems. This 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 paths0⟶i⟶N+i and i⟶N+i⟶2N+1. Since the constraint ti∈ai,bi 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)bN+i=minbN+i,b2N+i−ti,2N+i,bi=minbi,bN+i−ti,N+i,ai=maxai,a0+t0,i,aN+i=maxaN+i,ai+ti,N+i.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, ifai+tij>bj,i,j∈1,2N, then delete arc i,j.[Capacity]: the capacity of the vehicle can never be exceeded; therefore, ifqi+qj>T,i,j∈1,N,i≠j, then remove the following arcs: i,j;j,i;i,N+j;j,N+i;N+i,N+j;N+j,N+i.[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:If pathj⟶i⟶N+j⟶N+i is not feasible for tj=aj, then arc i,N+j is deletedIf pathi⟶N+i⟶j⟶N+j is not feasible for ti=ai, then arc N+i,j is deletedIf pathi⟶j⟶N+i⟶N+j is not feasible for ti=ai, then arc i,j is deletedIf pathi⟶j⟶N+i⟶N+j and j⟶i⟶N+i⟶N+i⟶N+j⟶N+j is not feasible for ti=ai and tj=aj, remove arc N+i,N+j


## 3.2. 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 fractionalYk variable to generate two subproblems, one of which is fixed Yk0, and the other is fixed Yk1. If all Yk variables are integers, then the algorithm branches off from the minimum number of Xij variables.


## 3.3. Label-Setting Algorithm

The 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 two-way 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.


## 4. Computational Experiments

### 4.1. Instances Description

This 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. The 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 Figure1), 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). The 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 hi at the affected point uniformly generates high uncertainty in 0,0.5di and low uncertainty in 0,0.3di. At this time, we describe the impact of the maximum deviation hi on the demand and only change its value to describe the magnitude of the demand deviation, resulting in demand uncertainty.Figure 1
6 emergency facilities and 12 disaster sites.Table 1
Relevant data of candidate emergency facilities.
Candidate facilitiesLocation coordinatesCapacityTransportation costA(96.855 675, 33.507 342)122 00012.2B(97.213 021, 32.711 051)100 00010.0C(95.818 005, 32.604 522)140 00014.0D(96.328 710, 34.087 447)160 00016.0E(95.966 188, 33.411 089)185 00018.5F(96.548 062, 32.983 319)100 00010.0Table 2
Nominal demand for materials in disaster areas.
Disaster areaLocation coordinatesNominal demandLongbao(96.423 140, 33.268 340)23 950Xialaxiu(96.610 170, 32.663 040)3560Sahuteng(95.624 890, 33.043 210)4200Chengwen(97.095 760, 33.366 309)2280Qingshuihe(97.137 670, 33.803 491)4800Jiajibo(95.630 240, 33.844 551)29 050Xiangda(96.475 460, 32.203 150)3120Yuegai(95.811 340, 34.131 201)3880Jiegu(96.978 617, 33.001 271)4200Xiewu(97.355 080, 33.130 181)8700Zhaduo(96.743 690, 33.767 039)10 600Zhenqin(97.301 630, 33.407 661)13 700


### 4.2. Numerical Experiment Results

The following sections introduce the results of a series of experiments we conducted. The 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. The factors studied include the uncertainty of demand, the size of the customer’s size, and some key parameters of the relevant uncertainty set.


#### 4.2.1. Determining Model Results

Firstly, we solve the deterministic model and give its optimal solution, as shown in Figure2.Figure 2
Facility location under nominal demand.It can be seen from Figure2 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.


#### 4.2.2. 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 Figure2). Then, 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.Figure 3
Facility location under low uncertainty andΓ=3.Figure 4
Facility location under high uncertainty andΓ=3.Table 3
Comparison of results under uncertain demand and uncertain set.
UncertaintyCostOpen facility numberTime∗Gap∗Γ=3Γ=4Γ=5Γ=3Γ=4Γ=5072 737.573 061.473 528.3A, D, E, FA, D, E, FA, C, D, F34.520.28Low73 931.574 325.674 984.3A, C, D, FA, D, E, FA, D, E, F38.070.34High75 046.376 316.478 297.5A, B, D, E, FA, B, D, E, FA, B, D, E, F45.630.53∗Average results of 10 calculations.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 Figure3, 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.Figure4 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.Table3 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.Figure 5
Vehicle routing with high uncertainty level andΓ=5.Table3 reports average statistics for robust and deterministic solutions. The 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. This 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.


## 4.1. Instances Description

This 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. The 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 Figure1), 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). The 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 hi at the affected point uniformly generates high uncertainty in 0,0.5di and low uncertainty in 0,0.3di. At this time, we describe the impact of the maximum deviation hi on the demand and only change its value to describe the magnitude of the demand deviation, resulting in demand uncertainty.Figure 1
6 emergency facilities and 12 disaster sites.Table 1
Relevant data of candidate emergency facilities.
Candidate facilitiesLocation coordinatesCapacityTransportation costA(96.855 675, 33.507 342)122 00012.2B(97.213 021, 32.711 051)100 00010.0C(95.818 005, 32.604 522)140 00014.0D(96.328 710, 34.087 447)160 00016.0E(95.966 188, 33.411 089)185 00018.5F(96.548 062, 32.983 319)100 00010.0Table 2
Nominal demand for materials in disaster areas.
Disaster areaLocation coordinatesNominal demandLongbao(96.423 140, 33.268 340)23 950Xialaxiu(96.610 170, 32.663 040)3560Sahuteng(95.624 890, 33.043 210)4200Chengwen(97.095 760, 33.366 309)2280Qingshuihe(97.137 670, 33.803 491)4800Jiajibo(95.630 240, 33.844 551)29 050Xiangda(96.475 460, 32.203 150)3120Yuegai(95.811 340, 34.131 201)3880Jiegu(96.978 617, 33.001 271)4200Xiewu(97.355 080, 33.130 181)8700Zhaduo(96.743 690, 33.767 039)10 600Zhenqin(97.301 630, 33.407 661)13 700


## 4.2. Numerical Experiment Results

The following sections introduce the results of a series of experiments we conducted. The 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. The factors studied include the uncertainty of demand, the size of the customer’s size, and some key parameters of the relevant uncertainty set.


### 4.2.1. Determining Model Results

Firstly, we solve the deterministic model and give its optimal solution, as shown in Figure2.Figure 2
Facility location under nominal demand.It can be seen from Figure2 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.


### 4.2.2. 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 Figure2). Then, 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.Figure 3
Facility location under low uncertainty andΓ=3.Figure 4
Facility location under high uncertainty andΓ=3.Table 3
Comparison of results under uncertain demand and uncertain set.
UncertaintyCostOpen facility numberTime∗Gap∗Γ=3Γ=4Γ=5Γ=3Γ=4Γ=5072 737.573 061.473 528.3A, D, E, FA, D, E, FA, C, D, F34.520.28Low73 931.574 325.674 984.3A, C, D, FA, D, E, FA, D, E, F38.070.34High75 046.376 316.478 297.5A, B, D, E, FA, B, D, E, FA, B, D, E, F45.630.53∗Average results of 10 calculations.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 Figure3, 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.Figure4 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.Table3 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.Figure 5
Vehicle routing with high uncertainty level andΓ=5.Table3 reports average statistics for robust and deterministic solutions. The 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. This 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.


## 4.2.1. Determining Model Results

Firstly, we solve the deterministic model and give its optimal solution, as shown in Figure2.Figure 2
Facility location under nominal demand.It can be seen from Figure2 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.


## 4.2.2. 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 Figure2). Then, 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.Figure 3
Facility location under low uncertainty andΓ=3.Figure 4
Facility location under high uncertainty andΓ=3.Table 3
Comparison of results under uncertain demand and uncertain set.
UncertaintyCostOpen facility numberTime∗Gap∗Γ=3Γ=4Γ=5Γ=3Γ=4Γ=5072 737.573 061.473 528.3A, D, E, FA, D, E, FA, C, D, F34.520.28Low73 931.574 325.674 984.3A, C, D, FA, D, E, FA, D, E, F38.070.34High75 046.376 316.478 297.5A, B, D, E, FA, B, D, E, FA, B, D, E, F45.630.53∗Average results of 10 calculations.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 Figure3, 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.Figure4 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.Table3 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.Figure 5
Vehicle routing with high uncertainty level andΓ=5.Table3 reports average statistics for robust and deterministic solutions. The 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. This 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.


## 5. Conclusion

The 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-vehicle-routing problem with uncertain demand. The 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. Therefore, 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. Then, 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. The 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. The cost increase is also less than the high uncertainty level with Γ=3. This 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.The 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.



---

*Source: 9152989-2021-12-31.xml*