# Multiobjective Location Routing Problem considering Uncertain Data after Disasters

Authors: Keliang Chang; Hong Zhou; Guijing Chen; Huiqin Chen
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

The relief distributions after large disasters play an important role for rescue works. After disasters there is a high degree of uncertainty, such as the demands of disaster points and the damage of paths. The demands of affected points and the velocities between two points on the paths are uncertain in this article, and the robust optimization method is applied to deal with the uncertain parameters. This paper proposes a nonlinear location routing problem with half-time windows and with three objectives. The affected points can be visited more than one time. The goals are the total costs of the transportation, the satisfaction rates of disaster nodes, and the path transport capacities which are denoted by vehicle velocities. Finally, the genetic algorithm is applied to solve a number of numerical examples, and the results show that the genetic algorithm is very stable and effective for this problem.


---

## Body

## 1. Introduction

In recent years, man-made or natural disasters which caused huge casualties and economic losses occurred frequently in different regions and countries. Timely and effective rescue works are very important after disasters. The emergency logistics management mainly includes two aspects [1]: the facility location problem (FLP) and the vehicle routing problem (VRP). In fact, there is a close relationship between facility location problem and vehicle routing problem about emergency logistics, that is, location routing problem (LRP). Compared with the traditional transportation logistics, the facility location problem and vehicle routing problem about emergency logistics are more challenging and complex [2]. von Boventer [3] combined the facility location problem and the vehicle routing problem. Up to now, there have been abundant research results about LRP [4, 5].After disasters, there are a high degree of uncertainties in all aspects, such as the relief demands and the path transport capacities. Three major methods for dealing with the uncertain parameters are stochastic programming, fuzzy method, and robust optimization method. Ahmadi-Javid and Seddighi [6] studied a stochastic LRP problem. The materials which were provided by facilities are random variables. Under uncertain demands in the short time after disasters [7], large scale emergency scheduling problems were studied, including the selection of supply points, route selection, the decisions of the transport quantities, and the distribution methods, and the fuzzy set was applied to express the uncertain demands. Zare Mehrjerdi and Nadizadeh [8] studied the LRP problem with capacity restrictions of the vehicles and the distribution centers. Triangular fuzzy number represented the uncertainty demands. Lu [9] considered the robust weighted P-center model, and the travel time and demands of affected points were all uncertain. Wang et al. [10] proposed a dynamic time-space network model. It was a network flow model with multistages and multimaterials. Robust optimization method was applied to deal with the uncertain demand. Wang et al. [11] proposed a multiobjective location routing model with split delivery. It optimized the allocation of reliefs after earthquake. The objectives were to minimize the travel time, to minimize the total costs including the fixed costs of the distribution centers and the vehicles transportation costs and to maximize the path reliability, but the demands of disaster points and the probabilities of available arcs were certain. Koç et al. [12] studied the LRP problem with time windows. It was a certain problem with heterogeneous fleets and with spit delivery, and the objective was to minimize the total costs.It can be seen from the existing literatures that there are little researches on LRP problem with uncertain parameters which are solved by robust optimization. In this paper, the robust optimization method is used to deal with the uncertain demands and the uncertain velocities between two disaster points on the paths. In addition, the most scholars take the cost as the objective of LRP problem. Actually, the distances and the travel time between two points after disasters will be affected, and we will pay the price. Then, the distance and travel time can be understood as the travel cost. In this paper, the half-time windows constraints are quoted, so the time when the materials reach the demand points cannot be later than the specified time. Thus, the timeliness of emergency rescues is improved. After disasters, all kinds of materials are in short supplies, so it is very important to possess adequate relief supplies. Therefore, the minimization of the total distribution costs is one of the objectives, and the maximization of the worst path satisfaction rates is the second objective. After incidents, the transport network will be destroyed. In order to find a better path, the maximization of path transport capacities is the third objective in this paper.


## 2. The Problem Description

Generally, the distribution network of reliefs after disasters is described as a graphG=V,E. V is a vertex set, and M={n+1,n+2,…,n+m} is a set of candidate distribution centers. We assume that there are no demands for a point in the set M. N={1,2,…,n} is disaster point set. E=i,j:i,j∈V,i≠j is effective arc set. K={1,2,…,k} is vehicle set. This paper considers the following three objectives:(1)
Minimize total costs including fixed costs and vehicle transportation costs.(2)
Maximize the minimum material satisfaction rates of demand points.(3)
Maximize the transport capacities of the worst path. (The transport capacities are represented by the velocities of the vehicle.)This paper has the following assumptions:(1)
The disaster points and the candidate distribution centers are known and the capacities of the candidate distribution centers are large enough.(2)
The available arcs and distances between two points on the transport network are known.(3)
Because the reliefs are calculated by volume, different types of relief supplies can be regarded as a kind of material.(4)
The relief demands of the demand points are greater than or equal to the amounts of supplies because of the materials shortage after disasters.(5)
In this paper, we can only consider the disaster points which can be serviced by vehicles, and the disaster points which can be serviced by the special transportation methods (e.g., helicopters) are ignored.The parameters and variables used in this paper are introduced in Notations.


## 3. The Multiobjectives LRP Model

### 3.1. Objective 1: Minimization of the Total Distribution Cost—min⁡f1

In the location routing problem, it is necessary to determine the number and location of the distribution centers and to arrange the disaster points to the distribution centers, and the corresponding vehicle routing will be decided. Therefore, the total costs of the relief distributions include the fixed costs of distribution centers and the vehicle transport costs.(1)f1=∑j∈Mfjxj+∑k∈K∑i,j∈Eckdijyijk.


### 3.2. Objective 2: Maximization of the Worst Path Satisfaction Rates—max⁡f2

We hope that reliefs can reach the demand point timely and effectively after disasters, so we should consider the relief satisfaction rates of the demand points. We hope to maximize the relief satisfaction rates in affected areas, and the fairness of the relief distributions is taken into account. Therefore, the second objective is maximization of the worst satisfaction rates.(2)f2=mini∈Nri,where ri=∑k∈Kqik/Di,∀i∈N, is called the relief satisfaction rate of demand point i. The demands of point i are uncertain: Di∈D_i,D¯i,∀i∈N,0≤D_i<D¯i.


### 3.3. Objective 3: Maximization of the Worst Path Transport Capacities—max⁡f3

The original paths after disasters are affected more or less. At this time, we take into account the path transport capacities, and the worst path transport capacities are maximized in this paper. The sum of velocities indicates the path transport capacities. The vehicle velocities are uncertain. Letvij∈v_ij,v¯ij,∀i,j∈E,0<v_ij<v¯ij.(3)f3=mink∈K∑i,j∈Evijyijk.The following formulas are the constraint conditions. In adition, the constraints conditions (4), (5), (6), (7), (8), (9), (16), (17), and (18) of the literature [11] are still used in our study. (4)tik-tjk+si+dij≤M1-yijk,∀i,j∈E,k∈K,(5)tik≤bi,∀i∈N,k∈K,(6)∑i∈N∑k∈Kqik≤Q,(7)∑i∈Nqik≤Lk,∀k∈K,(8)qik≥0,∀k∈K,i∈N,(9)yijk∈0,1,∀i,j∈E,k∈K.Formulas (4) and (5) are the time window constraints. M is a great positive number. This article requires that the sum of travel time and the service time is not bigger than the set time. Constraint (6) shows that the total amounts of reliefs delivered from the distribution centers to the disaster points shall not exceed the total amounts of available reliefs. Constraint (7) ensures that the relief volume transported to the affected areas by vehicle cannot exceed the load capacities of vehicles. Constraints (8)-(9) ensure that the decision variables are 0-1 nonnegative variables.There are subloop elimination constraints in order to avoid the subloop. The constraints proposed by Dror et al. [13] are applied in this paper. Let di indicate the outgoing degree of point i: di=∑k∈K∑j∈Vyijki∈N. (10)∑k∈K∑i,j∈Nyijk≤∑i∈Ndi-k,∑j∈M∑i∈Nyjik≤1,∀i,j∈E,k∈K,∑k∈K∑j∈Vyjik≥1,∀i∈N,∑j∈Vyjik≤1,∀i∈N,k∈K.


## 3.1. Objective 1: Minimization of the Total Distribution Cost—min⁡f1

In the location routing problem, it is necessary to determine the number and location of the distribution centers and to arrange the disaster points to the distribution centers, and the corresponding vehicle routing will be decided. Therefore, the total costs of the relief distributions include the fixed costs of distribution centers and the vehicle transport costs.(1)f1=∑j∈Mfjxj+∑k∈K∑i,j∈Eckdijyijk.


## 3.2. Objective 2: Maximization of the Worst Path Satisfaction Rates—max⁡f2

We hope that reliefs can reach the demand point timely and effectively after disasters, so we should consider the relief satisfaction rates of the demand points. We hope to maximize the relief satisfaction rates in affected areas, and the fairness of the relief distributions is taken into account. Therefore, the second objective is maximization of the worst satisfaction rates.(2)f2=mini∈Nri,where ri=∑k∈Kqik/Di,∀i∈N, is called the relief satisfaction rate of demand point i. The demands of point i are uncertain: Di∈D_i,D¯i,∀i∈N,0≤D_i<D¯i.


## 3.3. Objective 3: Maximization of the Worst Path Transport Capacities—max⁡f3

The original paths after disasters are affected more or less. At this time, we take into account the path transport capacities, and the worst path transport capacities are maximized in this paper. The sum of velocities indicates the path transport capacities. The vehicle velocities are uncertain. Letvij∈v_ij,v¯ij,∀i,j∈E,0<v_ij<v¯ij.(3)f3=mink∈K∑i,j∈Evijyijk.The following formulas are the constraint conditions. In adition, the constraints conditions (4), (5), (6), (7), (8), (9), (16), (17), and (18) of the literature [11] are still used in our study. (4)tik-tjk+si+dij≤M1-yijk,∀i,j∈E,k∈K,(5)tik≤bi,∀i∈N,k∈K,(6)∑i∈N∑k∈Kqik≤Q,(7)∑i∈Nqik≤Lk,∀k∈K,(8)qik≥0,∀k∈K,i∈N,(9)yijk∈0,1,∀i,j∈E,k∈K.Formulas (4) and (5) are the time window constraints. M is a great positive number. This article requires that the sum of travel time and the service time is not bigger than the set time. Constraint (6) shows that the total amounts of reliefs delivered from the distribution centers to the disaster points shall not exceed the total amounts of available reliefs. Constraint (7) ensures that the relief volume transported to the affected areas by vehicle cannot exceed the load capacities of vehicles. Constraints (8)-(9) ensure that the decision variables are 0-1 nonnegative variables.There are subloop elimination constraints in order to avoid the subloop. The constraints proposed by Dror et al. [13] are applied in this paper. Let di indicate the outgoing degree of point i: di=∑k∈K∑j∈Vyijki∈N. (10)∑k∈K∑i,j∈Nyijk≤∑i∈Ndi-k,∑j∈M∑i∈Nyjik≤1,∀i,j∈E,k∈K,∑k∈K∑j∈Vyjik≥1,∀i∈N,∑j∈Vyjik≤1,∀i∈N,k∈K.


## 4. The Solution Method of This Model

### 4.1. Data Uncertainty Description and Robust Solution

The travel velocities of vehicles on each arc are uncertain, so the interval of the vehicles travel velocity onarc(i,j) is vij∈v_ij,v¯ij,∀i,j∈E,0<v_ij<v¯ij. The relief demands on each point i are uncertain too. The relief demands on each point i are (11)Di∈D_i,D¯i,∀i∈N,0≤D_i<D¯i.Let z=(X,Y) be a feasible scheme for the above model, where X={xj,j∈M} and Y={yijk,(i,j)∈E}. Let Z be the feasible solutions set, so a scheme z∈Z corresponds to a group of qik(z) and a path r(z). S=[D_i,D¯i]×[v_ij,v¯ij],∀i∈N,(i,j)∈E indicates Cartesian product of two intervals. s∈S is called a scenario. In order to facilitate the description of the problem, f2, f3 are changed equally. (12)f2′=maxi∈NDi∑k∈Kqik,(13)f3′=max1∑i,j∈Evijyijk,∀k∈K.Given a scenario, problems (12) and (13) are certain. In certain problem the objective function 2 is(14)F2s=maxz∈Zf2′s,z.We need to define the maximum velocities sum of path r(z) under the scheme z and the scenario s. (15)vz,s,rz=max∑i,j∈rzvij,∀k∈K,(16)f3′′s,z=maxz∈Z1vz,s,rz.The objective function 3 is expressed by formula (17) under the scheme z and the scenario s. (17)F3s=maxz∈Zf3′′s,z.The robust deviations dev2(z,s) and dev3(z,s) about formulas (12) and (16) under the scheme z and the scenario s are indicated by formulas (18). Let z∗(s) be the optimal scheme of (14) and (17) under scenario s. (18)dev2z,s=f2′z,s-f2′z∗s,s,dev3z,s=f3′′z,s-f3′′z∗s,s.The robust costs of a scheme can be represented by the following problems: (19)rc2z=maxs∈Sdev2z,s,(20)rc3z=maxs∈Sdev3z,s.The second objective and the third objective of the above problem can be expressed as follows:(21)minz∈Zrc2z,minz∈Zrc3z.Formula (21) implies minimizing the maximum of the robust deviations.


### 4.2. Robust Cost Analysis

Because the uncertain demands and uncertain transport velocities are represented by continuous intervals, the scenario set is an infinite set. Therefore, it is very difficult to get the estimate values ofrc2(z), rc3(z). To find the worst scenario under scheme z is an urgent problem. Let si(z) be a scenario introduced by scheme z on the point i, and r(z) is the corresponding path for scheme z. The following two assumptions are satisfied:(1)
The demands of pointi equal the corresponding upper bound Di=Di¯, and the demands of other points equal the corresponding lower bound Dj=Dj¯,j≠i.(2)
The transport velocities of arc on the pathr(z) equal the corresponding lower bounds vij=v_ij,(i,j)∈r(z), and the transport velocities of arc on the other paths equal the corresponding upper bounds vij=v¯ij,(i,j)∉r(z). si(z) and si can be used to make the following changes.Lemma 1.
With the above introduction ofsi and s¯, for a scheme z,(22)rc2z=f2′s¯,z-f2′s¯,z∗s¯=f2′si,z-f2′si,z∗si,rc3z=f3′′s¯,z-f3′′s¯,z∗s¯=f3′′si,z-f3′′si,z∗si,where z∗(s¯) and z∗(si) are optimal schemes under the scenarios s¯ and si.Proof.
Taking the second objective function as an example, the third objective function can be proved by the same method. Letz∗(s¯) be the optimal solution of problem (14), and (s¯,z∗(s¯)) is the optimal solution of problem (19).
Main Claim. (s¯,z∗(s¯)) is the optimal solution of problem (19), so si is the worst scenario under scheme z, and z∗(s¯) is the optimal scheme of F2(si). Then the following formulas can be obtained:(23)f2′si,z-f2′si,z∗si=f2′si,z-f2′si,z∗s¯=f2′s¯,z-f2′s¯,z∗s¯,f3′si,z-f3′si,z∗si=f3′si,z-f3′si,z∗s¯=f3′s¯,z-f3′s¯,z∗s¯.
It can be seen that the main claim of Lemma1 can be proved converting s¯ to si. We can finish this conversion by the following steps.
Step 1. The upper bound D¯i of point i substitutes Di(s¯).
Before converting, since(s¯,z∗(s¯)) is the optimal solution of (19), f2′(s¯,z)-f2′(s¯,z∗(s¯))≥0. Therefore, ∑k∈Kqik(z∗(s¯))≥∑k∈Kqik(z(s¯)). In the first conversion step, f2′(s¯,z∗(s¯)) cannot increase by more than (D¯i-Di(s¯))(1/∑k∈Kqikz∗s¯)≤(D¯i-Di(s¯))(1/∑k∈Kqikzs¯). In addition, in Step 1, the value of f2′(s¯,z)-f2′(s¯,z∗(s¯)) cannot decrease because D¯i≥D¯i(s¯) and cannot increase because (s¯,z∗(s¯)) is the optimal solution of problem (19). Therefore, the value of f2′(s¯,z)-f2′(s¯,z∗(s¯)) does not change in Step 1.
Step 2. The lower bound D_l substitutes Dl(s¯), where ∀l∈N,l≠i, and the values of f2′(s¯,z) and f2′(s¯,z∗(s¯)) do not change. Therefore, (s¯,z∗(s¯)) remains the optimal solution of (19). Thus, the main claim is proven, as in Lemma 1.The following theorem is obtained by Lemma1, which greatly simplifies formulas (19) and (20).Theorem 2.
For anyz∈Z, (24)rc2z=maxi∈ND¯i∑k∈Kqik-f2′si,z∗si,rc3z=minz∈Z1∑i,j∈rzv_ij-f3′′si,z∗si.


### 4.3. The Solutions of This Model

Many practical problems need to optimize multiple objectives simultaneously. Sometimes these goals often compete with each other or contradict each other, so the definition of Pareto optimal solution is introduced.(1) Pareto Dominance. Considering all objectives, if solution x1 is at least as equal as x2, and better than x2 with at least one objective value, solution x1 dominates x2 (denoted as x1≻x2). For minimizing (f1,…,fw),x1≻x2 if(25)∀w∈1,2,…,W:fwx1≤fwx2∧∀w′∈1,2,…,W:fw′x1≤fw′x2.(2) Pareto Optimum. A solution x1 is called Pareto optimal or nondominated solution if and only if there is no any solution x2 that satisfies x2≻x1.(3) Pareto Front. Furthermore, if x1 is Pareto optimal (nondominated), then f(x1)={f1(x1),…,fw(x1)} is said to be the nondominated vector. The set of all nondominated vectors is called Pareto front (or nondominated frontier).The papers [14–17] are the common methods to solve multiobjective optimization problems. In this paper, genetic algorithm [17] is applied to solve the uncertain problem with half-time windows. The specific processes are as follows.(1) Initial Population. According to the characteristics of the LRP problem, each chromosome includes three substrings. (26)Xgt=xg11t,xg12t,…,xg1kt︸xg1t,xg21t,xg22t,…,xg2kt︸xg2t,xg31t,xg32t,…,xg3nt︸xg3t,g=1,2,…,NP.We can calculate three objective function values through decomposing Xgt [11]. Choose the worst scenario according to Sections 4.1 and 4.2. Transport time is calculated based on the distances between the points and the velocities of the vehicle (Ignoring the distribution time of distribution centers). If the sum of travel time and service time are greater than bi, the transport capacities of path are punished, and the transport capacities are updated to bi (transport capacity)/(transport time + service time).(2) Mutation Operation. Variation vectors can be obtained through mutating. In order to avoid being trapped in local optimum, this paper introduces the inversion sequence variation method for each substring of chromosome [18]. It randomly chooses two notes within a chromosome and then reverses their contents. The following two formulas are the examples. Parent is [5 6∣2 3 7 8 9∣1 4], and child is [5 6∣9 8 7 3 2∣1 4] after reversing.(3) Crossover Operation. Crossover operation is applied to obtain the trail vector Ugt after mutation operator. According to the literature [19], the two points crossover method is used in ug3, and single point crossover method is applied to produce ug1, ug2. ug1, ug2, and ug3 which constitute the trail vector Ugt. (27)Ugt=ug11t,ug12t,…,ug1kt︸ug1t,ug21t,ug22t,…,ug2kt︸ug2t,ug31t,ug32t,…,ug3nt︸ug3t.(4) Selection Operation. The weighted method [20] is applied to sort individuals in population. λ1 is selected randomly on interval (0.1,0.3], and λ3 is obtained randomly on interval (0.3,0.5]. Because the satisfaction rate is less than 1, it is easy to lose with the total costs and the total path velocities, so we select λ′ on interval [500,1000] randomly. Then, we let λ2=λ′×N (N is the number of disaster points).(5) Steps of Genetic Algorithm
(S1)
Generate initial populationX0 of size NP randomly.
(S2)
The variation population is obtained by the variation process, and trail population is obtained by the crossover operation.
(S3)
Combine the parent populationXt and trial population Ut together to form population Rt.
(S4)
Compute the objective values for each chromosome inRt.
(S5)
Determine the values ofλ1, λ2, and λ3 and compute the value of F=λ1(-f1)+λ2f2+λ3f3 and sort individuals based on F.
(S6)
Select the first NP individuals as the next generation population.
(S7)
Stop the procedure if the generationt is bigger than maxgen (maximum of iteration times); else turn to (S2).


## 4.1. Data Uncertainty Description and Robust Solution

The travel velocities of vehicles on each arc are uncertain, so the interval of the vehicles travel velocity onarc(i,j) is vij∈v_ij,v¯ij,∀i,j∈E,0<v_ij<v¯ij. The relief demands on each point i are uncertain too. The relief demands on each point i are (11)Di∈D_i,D¯i,∀i∈N,0≤D_i<D¯i.Let z=(X,Y) be a feasible scheme for the above model, where X={xj,j∈M} and Y={yijk,(i,j)∈E}. Let Z be the feasible solutions set, so a scheme z∈Z corresponds to a group of qik(z) and a path r(z). S=[D_i,D¯i]×[v_ij,v¯ij],∀i∈N,(i,j)∈E indicates Cartesian product of two intervals. s∈S is called a scenario. In order to facilitate the description of the problem, f2, f3 are changed equally. (12)f2′=maxi∈NDi∑k∈Kqik,(13)f3′=max1∑i,j∈Evijyijk,∀k∈K.Given a scenario, problems (12) and (13) are certain. In certain problem the objective function 2 is(14)F2s=maxz∈Zf2′s,z.We need to define the maximum velocities sum of path r(z) under the scheme z and the scenario s. (15)vz,s,rz=max∑i,j∈rzvij,∀k∈K,(16)f3′′s,z=maxz∈Z1vz,s,rz.The objective function 3 is expressed by formula (17) under the scheme z and the scenario s. (17)F3s=maxz∈Zf3′′s,z.The robust deviations dev2(z,s) and dev3(z,s) about formulas (12) and (16) under the scheme z and the scenario s are indicated by formulas (18). Let z∗(s) be the optimal scheme of (14) and (17) under scenario s. (18)dev2z,s=f2′z,s-f2′z∗s,s,dev3z,s=f3′′z,s-f3′′z∗s,s.The robust costs of a scheme can be represented by the following problems: (19)rc2z=maxs∈Sdev2z,s,(20)rc3z=maxs∈Sdev3z,s.The second objective and the third objective of the above problem can be expressed as follows:(21)minz∈Zrc2z,minz∈Zrc3z.Formula (21) implies minimizing the maximum of the robust deviations.


## 4.2. Robust Cost Analysis

Because the uncertain demands and uncertain transport velocities are represented by continuous intervals, the scenario set is an infinite set. Therefore, it is very difficult to get the estimate values ofrc2(z), rc3(z). To find the worst scenario under scheme z is an urgent problem. Let si(z) be a scenario introduced by scheme z on the point i, and r(z) is the corresponding path for scheme z. The following two assumptions are satisfied:(1)
The demands of pointi equal the corresponding upper bound Di=Di¯, and the demands of other points equal the corresponding lower bound Dj=Dj¯,j≠i.(2)
The transport velocities of arc on the pathr(z) equal the corresponding lower bounds vij=v_ij,(i,j)∈r(z), and the transport velocities of arc on the other paths equal the corresponding upper bounds vij=v¯ij,(i,j)∉r(z). si(z) and si can be used to make the following changes.Lemma 1.
With the above introduction ofsi and s¯, for a scheme z,(22)rc2z=f2′s¯,z-f2′s¯,z∗s¯=f2′si,z-f2′si,z∗si,rc3z=f3′′s¯,z-f3′′s¯,z∗s¯=f3′′si,z-f3′′si,z∗si,where z∗(s¯) and z∗(si) are optimal schemes under the scenarios s¯ and si.Proof.
Taking the second objective function as an example, the third objective function can be proved by the same method. Letz∗(s¯) be the optimal solution of problem (14), and (s¯,z∗(s¯)) is the optimal solution of problem (19).
Main Claim. (s¯,z∗(s¯)) is the optimal solution of problem (19), so si is the worst scenario under scheme z, and z∗(s¯) is the optimal scheme of F2(si). Then the following formulas can be obtained:(23)f2′si,z-f2′si,z∗si=f2′si,z-f2′si,z∗s¯=f2′s¯,z-f2′s¯,z∗s¯,f3′si,z-f3′si,z∗si=f3′si,z-f3′si,z∗s¯=f3′s¯,z-f3′s¯,z∗s¯.
It can be seen that the main claim of Lemma1 can be proved converting s¯ to si. We can finish this conversion by the following steps.
Step 1. The upper bound D¯i of point i substitutes Di(s¯).
Before converting, since(s¯,z∗(s¯)) is the optimal solution of (19), f2′(s¯,z)-f2′(s¯,z∗(s¯))≥0. Therefore, ∑k∈Kqik(z∗(s¯))≥∑k∈Kqik(z(s¯)). In the first conversion step, f2′(s¯,z∗(s¯)) cannot increase by more than (D¯i-Di(s¯))(1/∑k∈Kqikz∗s¯)≤(D¯i-Di(s¯))(1/∑k∈Kqikzs¯). In addition, in Step 1, the value of f2′(s¯,z)-f2′(s¯,z∗(s¯)) cannot decrease because D¯i≥D¯i(s¯) and cannot increase because (s¯,z∗(s¯)) is the optimal solution of problem (19). Therefore, the value of f2′(s¯,z)-f2′(s¯,z∗(s¯)) does not change in Step 1.
Step 2. The lower bound D_l substitutes Dl(s¯), where ∀l∈N,l≠i, and the values of f2′(s¯,z) and f2′(s¯,z∗(s¯)) do not change. Therefore, (s¯,z∗(s¯)) remains the optimal solution of (19). Thus, the main claim is proven, as in Lemma 1.The following theorem is obtained by Lemma1, which greatly simplifies formulas (19) and (20).Theorem 2.
For anyz∈Z, (24)rc2z=maxi∈ND¯i∑k∈Kqik-f2′si,z∗si,rc3z=minz∈Z1∑i,j∈rzv_ij-f3′′si,z∗si.


## 4.3. The Solutions of This Model

Many practical problems need to optimize multiple objectives simultaneously. Sometimes these goals often compete with each other or contradict each other, so the definition of Pareto optimal solution is introduced.(1) Pareto Dominance. Considering all objectives, if solution x1 is at least as equal as x2, and better than x2 with at least one objective value, solution x1 dominates x2 (denoted as x1≻x2). For minimizing (f1,…,fw),x1≻x2 if(25)∀w∈1,2,…,W:fwx1≤fwx2∧∀w′∈1,2,…,W:fw′x1≤fw′x2.(2) Pareto Optimum. A solution x1 is called Pareto optimal or nondominated solution if and only if there is no any solution x2 that satisfies x2≻x1.(3) Pareto Front. Furthermore, if x1 is Pareto optimal (nondominated), then f(x1)={f1(x1),…,fw(x1)} is said to be the nondominated vector. The set of all nondominated vectors is called Pareto front (or nondominated frontier).The papers [14–17] are the common methods to solve multiobjective optimization problems. In this paper, genetic algorithm [17] is applied to solve the uncertain problem with half-time windows. The specific processes are as follows.(1) Initial Population. According to the characteristics of the LRP problem, each chromosome includes three substrings. (26)Xgt=xg11t,xg12t,…,xg1kt︸xg1t,xg21t,xg22t,…,xg2kt︸xg2t,xg31t,xg32t,…,xg3nt︸xg3t,g=1,2,…,NP.We can calculate three objective function values through decomposing Xgt [11]. Choose the worst scenario according to Sections 4.1 and 4.2. Transport time is calculated based on the distances between the points and the velocities of the vehicle (Ignoring the distribution time of distribution centers). If the sum of travel time and service time are greater than bi, the transport capacities of path are punished, and the transport capacities are updated to bi (transport capacity)/(transport time + service time).(2) Mutation Operation. Variation vectors can be obtained through mutating. In order to avoid being trapped in local optimum, this paper introduces the inversion sequence variation method for each substring of chromosome [18]. It randomly chooses two notes within a chromosome and then reverses their contents. The following two formulas are the examples. Parent is [5 6∣2 3 7 8 9∣1 4], and child is [5 6∣9 8 7 3 2∣1 4] after reversing.(3) Crossover Operation. Crossover operation is applied to obtain the trail vector Ugt after mutation operator. According to the literature [19], the two points crossover method is used in ug3, and single point crossover method is applied to produce ug1, ug2. ug1, ug2, and ug3 which constitute the trail vector Ugt. (27)Ugt=ug11t,ug12t,…,ug1kt︸ug1t,ug21t,ug22t,…,ug2kt︸ug2t,ug31t,ug32t,…,ug3nt︸ug3t.(4) Selection Operation. The weighted method [20] is applied to sort individuals in population. λ1 is selected randomly on interval (0.1,0.3], and λ3 is obtained randomly on interval (0.3,0.5]. Because the satisfaction rate is less than 1, it is easy to lose with the total costs and the total path velocities, so we select λ′ on interval [500,1000] randomly. Then, we let λ2=λ′×N (N is the number of disaster points).(5) Steps of Genetic Algorithm
(S1)
Generate initial populationX0 of size NP randomly.
(S2)
The variation population is obtained by the variation process, and trail population is obtained by the crossover operation.
(S3)
Combine the parent populationXt and trial population Ut together to form population Rt.
(S4)
Compute the objective values for each chromosome inRt.
(S5)
Determine the values ofλ1, λ2, and λ3 and compute the value of F=λ1(-f1)+λ2f2+λ3f3 and sort individuals based on F.
(S6)
Select the first NP individuals as the next generation population.
(S7)
Stop the procedure if the generationt is bigger than maxgen (maximum of iteration times); else turn to (S2).


## 5. Numerical Experiment

Parameter setting isNP=5×(m+n); mutation probability is pm=0.7; crossover probability is pc=0.7.The combination of the number of candidate distribution centersm and the number of disaster points n is(28)m,n=3-10,3-100,3-1000,6-100,6-1000,10-100,10-1000.The coordinates of each point are randomly generated in the plane, and the distances from pointi to point j equal the distances from point j to point i. The fixed costs of the distribution centers (yuan) are randomly generated on interval (0,3000]. The demands of the disaster points and the velocities of the available arcs are uncertain, and uncertain parameters are processed in accordance with the method of Sections 4.1 and 4.2. The required numbers of large and medium and small vehicles for different number of disaster points are shown in Table 1. Let the time when the disaster happened be the zero time. The last time (hour) for each demand point to be served is randomly generated on the interval (0,100]. Table 2 gives the parameters of the vehicles. Table 3 shows the results of the calculation, and the approximate Pareto front is shown in Figure 1 (3-10).Table 1
Required number of vehicles.
Vehicle type
10 points
100 points
1000 points
Large vehiclek1
10
100
1000
Medium vehiclek2
6
60
600
Small vehiclek3
9
90
900Table 2
Vehicle parameters.
Vehicle type
Vehicle capacityLk (cm3)
General velocityvk (km/h)
Unit transportation costck  (yuan/km)
Large vehiclek1
600
×
250
×
175
50
10.0
Medium vehiclek2
280
×
200
×
145
30
3.1
Small vehiclek3
231
×
150
×
130
20
1.7Table 3
Calculation results.
Examples
Total costs
Satisfaction rate (%)
Path velocity
Mean value of approximate Pareto fronts
Approximate Pareto solution
Mean value of approximate Pareto fronts
Approximate Pareto solution
Mean value of approximate Pareto fronts
Approximate Pareto solution
3-10
8150.1
5811
55.48
55.44
643.8333
1098
3-100
34625
25263
64.28
62.77
8644.1
8761
3-1000
158190
139560
59.35
59.80
89865
88913
6-100
42807
38976
60.57
64.49
8554
8575
6-1000
129560
116650
54.89
60.01
11145
12270
10-100
43869
48567
58.99
59.62
8097.4
8329
10-1000
137980
118960
58.67
61.13
87660
89437Figure 1
Obtained by genetic algorithm (3-10).


## 6. Conclusion and Prospect

In this paper, a multiobjective nonlinear location routing model with half-time windows is proposed. The arrival time of the reliefs to the demand points cannot be later than the specified time, so the timeliness of emergency reliefs is enhanced. The affected points can be visited more than one time in this article. After the disasters, all kinds of materials are in short supplies, so it is very important to have adequate relief supplies. Therefore, the minimization of the total distribution costs is one of the objectives, and the maximization of the worst path satisfaction rates is the second objective. After incidents, the transport network will be influenced. In order to find a better path, the maximization of path transport capacities is the third objective in this paper. After disasters situations are very complex with a high degree of uncertainties, such as demands, transportation time, and path through velocities. Therefore, this paper assumes that the demands of the disaster nodes and the transportation velocities of the available path are uncertain, and the robust optimization is applied to deal with the uncertainty. The genetic algorithm is applied to solve a number of numerical examples; the results show that the algorithm is very stable and effective for this problem. Finally, the method of solving the problem can also apply simulated annealing algorithm or nondominated sorting genetic algorithm.



---

*Source: 1703608-2017-03-23.xml*