# Modeling of Emergency Supply Scheduling Problem Based on Reliability and Its Solution Algorithm under Variable Road Network after Sudden-Onset Disasters

Authors: Tinggui Chen; Shiwen Wu; Jianjun Yang; Guodong Cong; Gongfa Li
Journal: Complexity
 (2020)

---

## Abstract

It is common that many roads in disaster areas are damaged and obstructed after sudden-onset disasters. The phenomenon often comes with escalated traffic deterioration that raises the time and cost of emergency supply scheduling. Fortunately, repairing road network will shorten the time of in-transit distribution. In this paper, according to the characteristics of emergency supplies distribution, an emergency supply scheduling model based on multiple warehouses and stricken locations is constructed to deal with the failure of part of road networks in the early postdisaster phase. The detailed process is as follows. When part of the road networks fail, we firstly determine whether to repair the damaged road networks, and then a model of reliable emergency supply scheduling based on bi-level programming is proposed. Subsequently, an improved artificial bee colony algorithm is presented to solve the problem mentioned above. Finally, through a case study, the effectiveness and efficiency of the proposed model and algorithm are verified.


---

## Body

## 1. Introduction

The frequency of large-scale natural disasters is increasing in the world, which has caused significant disturbance to people’s daily lives and brought huge losses [1]. Especially, the extreme disasters such as earthquakes and debris flows even cause a large number of casualties in addition to enormous economic losses. Therefore, study on emergency supply scheduling after sudden-onset disasters has become a hot issue recently. Usually, when earthquake, debris flow, and other disasters occur, not only roads are damaged and obstructed but also a large number of emergency supplies are urgently needed in the affected areas in a short time. In addition, in the case of road obstruction, it is of great practical significance to study how to quickly dispatch emergency supplies, improve the reliability of the emergency plan, reduce the time of in-transit distribution, and lower costs of handling such incidents.Currently, most research focus on relatively simple problems such as transportation route programming and relief supply distribution under the condition of sufficient relief supplies [2, 3]. However, after sudden-onset disasters, adequate relief supplies are hardly guaranteed in most cases, so it is necessary to give a comprehensive consideration to the choices of routes, supplies distribution plans, relief center locations, and other issues and implement an overall optimization. In addition, most of the studies mentioned above are based on the single-level emergency rescue network with simple transport means and emergency supply kinds, lacking in-depth consideration of the actual situation. Actually, the postdisaster emergency supplies scheduling process usually has the following remarkable characteristics [4]: (1) part of road networks may be impassable after disasters; (2) a large demand for relief supplies in disaster areas may lead to the problem of insufficient supplies; (3) in order to reduce the cost and time of supply distribution, a certain number of distribution centers should be built to transport supplies from long distances to the vicinity of the disaster areas; (4) the rescue is urgent, while there exist numerous disaster areas [5]; (5) various means of transportation are required; (6) it is necessary to allocate a certain amount of supplies to repair the damaged road networks and form different transport networks so as to shorten the distribution time. All of the above characteristics increase the complexity of this problem.This paper aims to propose a scheduling model of reliable emergency supplies based on bi-level programming and then judge whether to repair the damaged road networks. In addition, a swarm intelligence algorithm named artificial bee colony (ABC) is adopted to solve the problem. Subsequently, a case study on optimizing rescue routing plans for a sudden-onset earthquake disaster is carried out to validate the effectiveness and efficiency of the proposed model and algorithm.The structure of this paper is organized as follows. Section2 is literature review. Section 3 proposes a multidepot scheduling model under bi-level programming. In Section 4, an improved ABC algorithm is used to solve the problem built in Section 3. Subsequently, a case study is given to testify the effectiveness and efficiency of the proposed model and algorithm in Section 5. Finally, Section 6 concludes the paper and points out future research directions.


## 2. Literature Review

Lots of literatures have studied the problem of postdisaster scheduling of emergency supplies. For example, Wang et al. [6] proposed a time-space network model to address the dynamic emergency logistics planning problem and decomposed the proposed model into two multiperiod multicommodity network flow problems. Qiu et al. [7] studied the emergency logistics problem with multicenter, multicommodity, and single-affected point, considering that the path near the disaster point might be damaged, the information about the paths’ state was not complete, and the travel time was uncertain. They further established the nonlinear programming model where objective function was the maximization of time-satisfaction degree. Cui et al. [8] considered supply chain disruption and path risk to achieve the minimum supply transportation time and supply transportation cost. The postdisaster supply transportation model under uncertainty conditions was established in their research as well. Chi et al. [9] put forward that different goals in emergency supply scheduling were correlated. Their study combined two scheduling objectives, i.e., time and resource satisfaction, into a timeliness function and designed a nonlinear evaluation model for emergency supply scheduling that incorporated multiple supply centers and a single-affected point to realize that function. Liu et al. [10] established a dynamic supply allocation model with the aim of minimizing response time based on time-varying emergency supply distribution framework of supply-and-demand constraints. Zhou et al. [11] presented a multiobjective optimization model and applied an improved genetic algorithm in response to the situation of multiperiod dynamic emergency resource scheduling. Within an integrated framework, Jeong et al. [12] developed emergency logistics networks (ELNs) based upon efficiency, risk, and robustness metrics. In brief, all the abovementioned studies attempted to solve the problem of multiple objectives existing in emergency supply scheduling on the premise of adequate relief supplies.Moreover, some literatures have adopted the method of multilevel programming to carry out studies. For instance, Rath and Gutjahr [13] put forward that, in the process of scheduling emergency supplies after disasters, a supply system with intermediate warehouses should be established. They used a three-objective optimization model and further proposed a “math-heuristic” algorithm. Wang and Hu [14] introduced a bi-level scheduling optimization model based on multiple supply subjects and proposed an improved nondominated sorted genetic algorithm-II (NSGAII) to solve the model. Wang et al. [15] proposed a multiperiod model for the allocation of emergency supplies and used a nonlinear utility function to illustrate the relationship between equity of allocations and the cost of emergency response. These literatures mostly considered multilevel distribution or multiperiod distribution, but gave insufficient attention to the transportation difficulties caused by road network blockage in the early postdisaster period.Furthermore, some scholars also studied the problem of damaged road networks after disasters. For instance, Ahmadi et al. [16] considered the failure of road networks, but they did not take the specific impact of road damage into account. Itabashi et al. [17] proposed a model to determine the optimal facility locations and routes for working vehicles when a road network was broken by a large-scale disaster. Rawls and Turnquist [18] established a two-stage stochastic mixed integer programming model for emergency rescue facility location according to the dual uncertainties of the amount of supply demand and the availability of transportation network after disaster. In addition, combining the elements of road damage after disaster, Liberatore et al. [19] designed an emergency supply distribution model which could be used when the road was under repair.Based on the abovementioned literature review, most studies are focused on targets such as scheduling time and cost, but few considered the reliable supply distribution path programming under the situation that the roads are damaged. Thus, this paper combines the bi-level programming model and the “soft time window” problem caused by road network failure and proposes a reliable emergency supply scheduling scheme, which is dispatched in possibly damaged roads. In this scheme, the distribution of relief supplies with different degrees of urgency and the situation of multiple relief points are fully considered. Also, the road network structure and supply distribution scheduling are jointly optimized to analyze and solve the reliable emergency supply scheduling model in this scheme. The analysis and solution work with the maximum satisfaction degree of multiple disaster points are at the upper level. On the other side, the minimal total cost of bi-level distribution as well as the shortest distribution time is at the lower level.


## 3. Multidepot Scheduling Model under Bi-Level Programming

In view of the characteristics of emergency rescue in the early postdisaster period, if the relief supplies are directly transported by road from the collector-distributor center, and it will be impossible to ensure the rapid and timely distribution of large quantities of relief supplies to the disaster points. In order to offer a more timely distribution of large amounts of supplies from peripheral collector-distributor centers to disaster-stricken areas, the paper considers different modes of transportation and diversified distribution conditions, thus establishing the middle-level distribution centers to distribute relief supplies. The middle-level distribution centers give full play to the advantages of various means of transportation, which enables a more efficient distribution of supplies. In order to simplify the process of modeling and apply the model to a wider context, this paper constructs a bi-level emergency supply scheduling model and describes the whole scheduling process from the upper level to lower level, respectively. At the upper level, the emergency supplies are transported from the distribution center to the middle-level distribution centers; while at the lower one, the emergency supplies are transported from the middle-level distribution centers to different disaster areas. Therefore, based on the actual situation, this paper establishes the mathematical model consisting of multiple supply depots and multiple rescue-needed points under the bi-level programming, in order to provide a reasonable supply scheduling plan for responding to sudden-onset disasters, especially to the damaged road networks caused by disasters. The research framework of this paper is presented in Figure1.Figure 1
Research framework.


### 3.1. Problem Description

After a disaster, the emergency rescue will be restrained by the factors of time and the amount of supplies. In order to gather the maximum quantities of emergency supplies and dispatch them to the disaster area quickly, the usual practice is to collect supplies from the collector-distributor center and transport them to the distribution center near the disaster area as soon as possible, and then distribute them to the disaster area in time. However, due to the impact of sudden-onset disasters, damages may occur in some segments of the road network. This situation requires reality-based consideration about whether to consume some supplies to repair the damaged segments, so as to shorten the transportation distance, reduce the distribution time, and prompt efficient distribution. Therefore, this paper mainly solves the problem of supply allocation and route programming in terms of limited resources and strives to ensure the stability and reliability of emergency logistics to the highest extent by comprehensively considering the time, cost, efficiency, and other objectives under the constraints of vehicles, load, and soft time windows.Considering the actual rescue situation, this paper divides up the distribution system of emergency supplies into the upper level and lower level. The former is responsible for dispatching the rescue supplies from the collector-distributor center to various distribution centers and allocating part of resources to repair the damaged road networks, while the latter is for delivering rescue supplies to the respective stricken areas. At the same time, the transportation route and supply distribution must be planned reasonably so as to achieve the optimal delivery. Moreover, the distribution of resources used for repairing the road network also needs to be taken into account while arranging upper-to-lower distribution. According to the different tasks that need to be completed, this paper adopts the method of decomposing the upper and lower levels to design the road network and plan supply distribution route, as shown in Figure2.Figure 2
Schematic diagram of double-layer emergency logistics after disaster.Located in different positions of the system, the upper and lower levels will make different decisions. The upper level mainly considers whether the damaged road network should be repaired when transporting supplies from the collector-distributor center to the middle-level distribution center, while the lower level mainly ponders how to quickly formulate the best transportation and distribution route and quickly distribute relief supplies to each disaster point. Because different tasks lead to different decisions and objectives, the upper level has priority to achieve the reasonable distribution of supplies to distribution centers or to the damaged road network with limited relief supplies, while the lower level aims at the fairness and rapidity of supply distribution. They are interrelated and together produce the optimal decision-making scheme. Supply distribution to different damaged segments will lead to the formation of different structures of road networks, which further causes different distribution time, as illustrated in Figure3. The remaining part of supplies is further distributed to other distribution center. According to the quantity of remaining relief supplies, the quantity of the supplies distributed to the affected areas is different, which will also affect the satisfaction degree of the affected areas. The upper and lower levels interact with each other. The ultimate goal is evaluated in terms of different distribution delay and satisfaction degree.Figure 3
Different road network status (partial). (a) Repair all. (b) Partial repair. (c) Not repaired all.
(a)
(b)
(c)


### 3.2. Model Assumptions

Combining the problem description with the reality, this paper proposes the following assumptions: (1) the coordinates of the collector-distributor center, distribution centers, and disaster areas are known; (2) when using vehicles to load supplies, this paper does not differentiate the types of supplies or consider special supplies, but mixes all kinds of them in the same vehicle so as to improve the loading capability of vehicles; (3) the total amount of relief supplies is limited and supply falls short of demand, so in view of different urgency degrees of various supplies, priority should be given to the distribution of supplies with higher urgency degrees; (4) the relief supplies needed by emergency rescue points do not change over time. If the demand information changes during this period, we can divide this period into several time intervals so as to guarantee constant demand by emergency rescue points in every time interval; (5) the time window factor of repair-needed segment is not considered. When emergency supplies are allocated to the damaged segments, the traffic of damaged segments is supposed to be restored. Specific model symbols are shown in Table1.Table 1
Model symbol description.
i
Represents the collector-distributor center of supplies at the upper level
I
Represents the set of collector-distributor center at the upper level
v
1
Represents the fixed cost of vehicles
v
2
Represents the convertible cost of vehicles on the way
C
Represents the set of all kinds of emergency supplies
c
Represents different types of emergency supplies
l
Represents damaged road segments
L
Represents the set of damaged road segments
n
Represents the total number of vehicles participating in supply scheduling
K
Represents the set of supply-needed disaster locations
K
Represents a certain supply-needed disaster location (k ∈ K)
M
Represents the set of transportation modes
m
Represents a certain mode of transportation (m ∈ M)
j
Represents a distribution center at the middle level
J
Represents the set of middle-level distribution centers (j ∈ J)
P
i
j
m
c
The quantity of emergency suppliesc transported from the collector-distributor center i to the emergency distribution center j by means of transportation m
S
i
c
The amount of supplyc provided by the collector-distributor center i (i ∈ I, c ∈ C)
T
h
j
The handling capacity/throughput of distribution centerc
T
i
j
m
The rescue time spent on transporting supplies from the collector-distributor centeri to the emergency distribution center j by means of transportation m
y
i
j
m
c
Transport emergency suppliesc from the collector-distributor center i to the distribution center j by means of transportation m
c
a
p
i
j
m
The single transporting capability by means of transportationm for dispatching supplies from the collector-distributor center i to the distribution center j
D
i
j
The distance a vehicle travels at (i, j) segment
t
i
j
m
The transporting time spent on transporting supplies from the collector-distributor centeri to the emergency distribution center j by means of transportation m
ω
c
The urgency degree of emergency supplyc ∑c∈Cωc=1
t
j
c
u
p
The time needed for distributing emergency supplyc from collector-distributor i at the upper level to the distribution center j
t
i
m
c
u
p
The time needed for distributing emergency supplyc from collector-distributor i to the middle level by means of transportation m
X
j
k
m
Variable 0-1. When vehiclem passes (j, k) segment, the variable is 1 or it becomes 0
a
m
j
The moment by vehiclem to drive out of the distribution center j
d
l
k
c
The latest demand-meeting time for supplyc required by supply-needed location k
δ
m
k
c
The time delayed when transporting supplyc to supply-needed location k by means of transportation m. If the transportation is completed before dlkc, the δmkc will be 0
t
j
k
m
The operation time of vehiclem in segment (j, k). If (j, k)∈L, then tjkm = ∞
U
k
c
The quantity unmet of demand for supplyc at supply-needed location k
c
max
Maximum load weight of vehicles at the lower level
Y
m
k
c
The quantity of supplyc transported to disaster-affected location k by means of transportation m
a
m
k
c
The time spent on transporting supplyc to location k by means of transportation m. If no supply c is transported to location k, amkc = 0
l
c
The quantity of supplyc needed for making urgent repair on road segment l
ξ
l
It is 1 when the decision is made to repair road segmentl, or it becomes 0
q
k
c
The quantity of demand for supplyc at location k


### 3.3. Model Establishment

In a certain period after disasters, part of road networks may be impassable and the relief supplies are limited, so it is reasonable to consider allocating limited supplies to the road networks and disaster areas. In that case, when supplies are allocated to the damaged segment, the traffic of this segment will be assumed as restored. In the model established in this paper, if no supply is allocated to the damaged segment, it is considered that this segment is not feasible, which may lead to the exceeding of required time limit and the consequent delay that will cause corresponding punishment. When supplies are allocated to the damaged segment, the road segment returns to the normal state, and the traveling time of this on road segment is calculated according to the normal standard.Based on the previous analysis, the multirescue model integrating bi-level programming is established to minimize the cost of distribution and maximize the satisfaction degree of stricken areas (i.e., to achieve the maximum quantity needed of supplies and the timeliness of distribution in disaster areas). The model is described as follows.


#### 3.3.1. The Model for the Upper Level

The objective of the upper level is as follows:(1)MinZ1=∑c∈Cωc∑k∈K∑c∈CTkcXjkmYmkc+∑c∈C∑l∈LTlclcξl+∑c∈C∑m∈M∑k∈KUkcδmkc,(2)s.t.∑j∈J∑m∈MPijmc=Sic,(3)∑i∈I∑n∈M∑c∈CPijmcωc≤XjThj,(4)timcup=∑j∈JtjkcupXjkc+Tlc,(5)Tkc=∑j∈JtjkcupXjkc+tkcdown,(6)Tlc=∑l∈Ltijcupξl,where the first two parts of Formula (1) at the upper level minimize the weighted time calculated by allocating emergency supplies to disaster areas and damaged segments. Due to the different types of supplies, receiving time, and receiving quantity required by different disaster areas and damaged segments and the different demand for emergency supplies at each supply-needed location, the total relative efficiency cost is affected by time, supply quantity, and urgency degree of relief supplies. In the last part of Formula (1), the quantity of the unmet supplies, the urgency degree, and the delay time are used to show that, in case of the same delay, the longer the delay time is and the higher the urgency degree of supplies is, the greater the impact is and consequently the higher the delay cost will be. Formula (2) denotes that the amount of supplies dispatched from the collector-distributor center is equal to the amount of resources supplied. Formula (3) represents the throughput limitation of emergency distribution centers. Formula (4) denotes the time spent on transporting supplies at the upper level from the collector-distributor center i to emergency distribution centers and damaged road segments by means of transportation m, which represents the total transit time of the transports with relief supplies to the disaster points. Formula (5) represents the time that the disaster-stricken location k needs for receiving emergency supply c. Formula (6) represents the time required for transporting rescue supplies from the collector-distributor center i to the segment in need of repair.


#### 3.3.2. The Model for the Lower Level

The objective of the lower level is as follows:(7)MinZ2=nv1+v2∑j∈J∑k∈K∑m∈MXjkmDij,(8)s.t.∑j∈J∑k∈K∑m∈MXjkm≤n,(9)∑j∈JXjkm≤1,(10)amj=0,(11)amj+tjkm−amk≤B1−Xjkm,(12)0≤amkc−δmkc≤dlkc,(13)δmkc≤∑j∈JMXjkm,(14)Ymkc≤cmax∑k∈KXjkm,(15)∑m∈MYmkc+Ukc−qkc=0,(16)∑j∈J∑k∈Kjc=∑c∈C∑k∈K∑m∈MYmkc+∑c∈C∑l∈Llcξl,(17)Ymkc≥0,Ukc≥0,δmkc≥0,(18)yijmc,Xj,ξl,Xjkm∈0,1.The lower-level objective function is Formula (7), which represents the sum of the fixed cost and the convertible cost of the transportation of supplies from the lower-level distribution center to the disaster location. Formula (8) denotes the limited number of transport vehicles dispatched from the distribution center. Formula (9) indicates that each vehicle goes to the affected area only once at most. Formulas (10)–(13) guarantee the feasibility of the orderly driving of vehicles, which means that once a vehicle arrives at the supply-needed point later than the required time limit, there will be time delay. Formulas (14) and (15) ensure the feasibility of emergency supply logistics, that is, the supplies transported by distribution vehicles are less than the total amount of relief supplies in distribution centers. Formula (16) demonstrates the conservation of the quantity of supplies, indicating that the amount of supplies used to repair damaged road segments plus the amount of supplies used to meet the demand of supply-needed location equals the sum of all emergency supplies. Formula (17) represents nonnegative constraints. Formula (18) indicates 0-1 variable constraints. The identifiers B and M in the abovementioned formulas are defined as integers large enough.


## 3.1. Problem Description

After a disaster, the emergency rescue will be restrained by the factors of time and the amount of supplies. In order to gather the maximum quantities of emergency supplies and dispatch them to the disaster area quickly, the usual practice is to collect supplies from the collector-distributor center and transport them to the distribution center near the disaster area as soon as possible, and then distribute them to the disaster area in time. However, due to the impact of sudden-onset disasters, damages may occur in some segments of the road network. This situation requires reality-based consideration about whether to consume some supplies to repair the damaged segments, so as to shorten the transportation distance, reduce the distribution time, and prompt efficient distribution. Therefore, this paper mainly solves the problem of supply allocation and route programming in terms of limited resources and strives to ensure the stability and reliability of emergency logistics to the highest extent by comprehensively considering the time, cost, efficiency, and other objectives under the constraints of vehicles, load, and soft time windows.Considering the actual rescue situation, this paper divides up the distribution system of emergency supplies into the upper level and lower level. The former is responsible for dispatching the rescue supplies from the collector-distributor center to various distribution centers and allocating part of resources to repair the damaged road networks, while the latter is for delivering rescue supplies to the respective stricken areas. At the same time, the transportation route and supply distribution must be planned reasonably so as to achieve the optimal delivery. Moreover, the distribution of resources used for repairing the road network also needs to be taken into account while arranging upper-to-lower distribution. According to the different tasks that need to be completed, this paper adopts the method of decomposing the upper and lower levels to design the road network and plan supply distribution route, as shown in Figure2.Figure 2
Schematic diagram of double-layer emergency logistics after disaster.Located in different positions of the system, the upper and lower levels will make different decisions. The upper level mainly considers whether the damaged road network should be repaired when transporting supplies from the collector-distributor center to the middle-level distribution center, while the lower level mainly ponders how to quickly formulate the best transportation and distribution route and quickly distribute relief supplies to each disaster point. Because different tasks lead to different decisions and objectives, the upper level has priority to achieve the reasonable distribution of supplies to distribution centers or to the damaged road network with limited relief supplies, while the lower level aims at the fairness and rapidity of supply distribution. They are interrelated and together produce the optimal decision-making scheme. Supply distribution to different damaged segments will lead to the formation of different structures of road networks, which further causes different distribution time, as illustrated in Figure3. The remaining part of supplies is further distributed to other distribution center. According to the quantity of remaining relief supplies, the quantity of the supplies distributed to the affected areas is different, which will also affect the satisfaction degree of the affected areas. The upper and lower levels interact with each other. The ultimate goal is evaluated in terms of different distribution delay and satisfaction degree.Figure 3
Different road network status (partial). (a) Repair all. (b) Partial repair. (c) Not repaired all.
(a)
(b)
(c)


## 3.2. Model Assumptions

Combining the problem description with the reality, this paper proposes the following assumptions: (1) the coordinates of the collector-distributor center, distribution centers, and disaster areas are known; (2) when using vehicles to load supplies, this paper does not differentiate the types of supplies or consider special supplies, but mixes all kinds of them in the same vehicle so as to improve the loading capability of vehicles; (3) the total amount of relief supplies is limited and supply falls short of demand, so in view of different urgency degrees of various supplies, priority should be given to the distribution of supplies with higher urgency degrees; (4) the relief supplies needed by emergency rescue points do not change over time. If the demand information changes during this period, we can divide this period into several time intervals so as to guarantee constant demand by emergency rescue points in every time interval; (5) the time window factor of repair-needed segment is not considered. When emergency supplies are allocated to the damaged segments, the traffic of damaged segments is supposed to be restored. Specific model symbols are shown in Table1.Table 1
Model symbol description.
i
Represents the collector-distributor center of supplies at the upper level
I
Represents the set of collector-distributor center at the upper level
v
1
Represents the fixed cost of vehicles
v
2
Represents the convertible cost of vehicles on the way
C
Represents the set of all kinds of emergency supplies
c
Represents different types of emergency supplies
l
Represents damaged road segments
L
Represents the set of damaged road segments
n
Represents the total number of vehicles participating in supply scheduling
K
Represents the set of supply-needed disaster locations
K
Represents a certain supply-needed disaster location (k ∈ K)
M
Represents the set of transportation modes
m
Represents a certain mode of transportation (m ∈ M)
j
Represents a distribution center at the middle level
J
Represents the set of middle-level distribution centers (j ∈ J)
P
i
j
m
c
The quantity of emergency suppliesc transported from the collector-distributor center i to the emergency distribution center j by means of transportation m
S
i
c
The amount of supplyc provided by the collector-distributor center i (i ∈ I, c ∈ C)
T
h
j
The handling capacity/throughput of distribution centerc
T
i
j
m
The rescue time spent on transporting supplies from the collector-distributor centeri to the emergency distribution center j by means of transportation m
y
i
j
m
c
Transport emergency suppliesc from the collector-distributor center i to the distribution center j by means of transportation m
c
a
p
i
j
m
The single transporting capability by means of transportationm for dispatching supplies from the collector-distributor center i to the distribution center j
D
i
j
The distance a vehicle travels at (i, j) segment
t
i
j
m
The transporting time spent on transporting supplies from the collector-distributor centeri to the emergency distribution center j by means of transportation m
ω
c
The urgency degree of emergency supplyc ∑c∈Cωc=1
t
j
c
u
p
The time needed for distributing emergency supplyc from collector-distributor i at the upper level to the distribution center j
t
i
m
c
u
p
The time needed for distributing emergency supplyc from collector-distributor i to the middle level by means of transportation m
X
j
k
m
Variable 0-1. When vehiclem passes (j, k) segment, the variable is 1 or it becomes 0
a
m
j
The moment by vehiclem to drive out of the distribution center j
d
l
k
c
The latest demand-meeting time for supplyc required by supply-needed location k
δ
m
k
c
The time delayed when transporting supplyc to supply-needed location k by means of transportation m. If the transportation is completed before dlkc, the δmkc will be 0
t
j
k
m
The operation time of vehiclem in segment (j, k). If (j, k)∈L, then tjkm = ∞
U
k
c
The quantity unmet of demand for supplyc at supply-needed location k
c
max
Maximum load weight of vehicles at the lower level
Y
m
k
c
The quantity of supplyc transported to disaster-affected location k by means of transportation m
a
m
k
c
The time spent on transporting supplyc to location k by means of transportation m. If no supply c is transported to location k, amkc = 0
l
c
The quantity of supplyc needed for making urgent repair on road segment l
ξ
l
It is 1 when the decision is made to repair road segmentl, or it becomes 0
q
k
c
The quantity of demand for supplyc at location k


## 3.3. Model Establishment

In a certain period after disasters, part of road networks may be impassable and the relief supplies are limited, so it is reasonable to consider allocating limited supplies to the road networks and disaster areas. In that case, when supplies are allocated to the damaged segment, the traffic of this segment will be assumed as restored. In the model established in this paper, if no supply is allocated to the damaged segment, it is considered that this segment is not feasible, which may lead to the exceeding of required time limit and the consequent delay that will cause corresponding punishment. When supplies are allocated to the damaged segment, the road segment returns to the normal state, and the traveling time of this on road segment is calculated according to the normal standard.Based on the previous analysis, the multirescue model integrating bi-level programming is established to minimize the cost of distribution and maximize the satisfaction degree of stricken areas (i.e., to achieve the maximum quantity needed of supplies and the timeliness of distribution in disaster areas). The model is described as follows.


### 3.3.1. The Model for the Upper Level

The objective of the upper level is as follows:(1)MinZ1=∑c∈Cωc∑k∈K∑c∈CTkcXjkmYmkc+∑c∈C∑l∈LTlclcξl+∑c∈C∑m∈M∑k∈KUkcδmkc,(2)s.t.∑j∈J∑m∈MPijmc=Sic,(3)∑i∈I∑n∈M∑c∈CPijmcωc≤XjThj,(4)timcup=∑j∈JtjkcupXjkc+Tlc,(5)Tkc=∑j∈JtjkcupXjkc+tkcdown,(6)Tlc=∑l∈Ltijcupξl,where the first two parts of Formula (1) at the upper level minimize the weighted time calculated by allocating emergency supplies to disaster areas and damaged segments. Due to the different types of supplies, receiving time, and receiving quantity required by different disaster areas and damaged segments and the different demand for emergency supplies at each supply-needed location, the total relative efficiency cost is affected by time, supply quantity, and urgency degree of relief supplies. In the last part of Formula (1), the quantity of the unmet supplies, the urgency degree, and the delay time are used to show that, in case of the same delay, the longer the delay time is and the higher the urgency degree of supplies is, the greater the impact is and consequently the higher the delay cost will be. Formula (2) denotes that the amount of supplies dispatched from the collector-distributor center is equal to the amount of resources supplied. Formula (3) represents the throughput limitation of emergency distribution centers. Formula (4) denotes the time spent on transporting supplies at the upper level from the collector-distributor center i to emergency distribution centers and damaged road segments by means of transportation m, which represents the total transit time of the transports with relief supplies to the disaster points. Formula (5) represents the time that the disaster-stricken location k needs for receiving emergency supply c. Formula (6) represents the time required for transporting rescue supplies from the collector-distributor center i to the segment in need of repair.


### 3.3.2. The Model for the Lower Level

The objective of the lower level is as follows:(7)MinZ2=nv1+v2∑j∈J∑k∈K∑m∈MXjkmDij,(8)s.t.∑j∈J∑k∈K∑m∈MXjkm≤n,(9)∑j∈JXjkm≤1,(10)amj=0,(11)amj+tjkm−amk≤B1−Xjkm,(12)0≤amkc−δmkc≤dlkc,(13)δmkc≤∑j∈JMXjkm,(14)Ymkc≤cmax∑k∈KXjkm,(15)∑m∈MYmkc+Ukc−qkc=0,(16)∑j∈J∑k∈Kjc=∑c∈C∑k∈K∑m∈MYmkc+∑c∈C∑l∈Llcξl,(17)Ymkc≥0,Ukc≥0,δmkc≥0,(18)yijmc,Xj,ξl,Xjkm∈0,1.The lower-level objective function is Formula (7), which represents the sum of the fixed cost and the convertible cost of the transportation of supplies from the lower-level distribution center to the disaster location. Formula (8) denotes the limited number of transport vehicles dispatched from the distribution center. Formula (9) indicates that each vehicle goes to the affected area only once at most. Formulas (10)–(13) guarantee the feasibility of the orderly driving of vehicles, which means that once a vehicle arrives at the supply-needed point later than the required time limit, there will be time delay. Formulas (14) and (15) ensure the feasibility of emergency supply logistics, that is, the supplies transported by distribution vehicles are less than the total amount of relief supplies in distribution centers. Formula (16) demonstrates the conservation of the quantity of supplies, indicating that the amount of supplies used to repair damaged road segments plus the amount of supplies used to meet the demand of supply-needed location equals the sum of all emergency supplies. Formula (17) represents nonnegative constraints. Formula (18) indicates 0-1 variable constraints. The identifiers B and M in the abovementioned formulas are defined as integers large enough.


## 3.3.1. The Model for the Upper Level

The objective of the upper level is as follows:(1)MinZ1=∑c∈Cωc∑k∈K∑c∈CTkcXjkmYmkc+∑c∈C∑l∈LTlclcξl+∑c∈C∑m∈M∑k∈KUkcδmkc,(2)s.t.∑j∈J∑m∈MPijmc=Sic,(3)∑i∈I∑n∈M∑c∈CPijmcωc≤XjThj,(4)timcup=∑j∈JtjkcupXjkc+Tlc,(5)Tkc=∑j∈JtjkcupXjkc+tkcdown,(6)Tlc=∑l∈Ltijcupξl,where the first two parts of Formula (1) at the upper level minimize the weighted time calculated by allocating emergency supplies to disaster areas and damaged segments. Due to the different types of supplies, receiving time, and receiving quantity required by different disaster areas and damaged segments and the different demand for emergency supplies at each supply-needed location, the total relative efficiency cost is affected by time, supply quantity, and urgency degree of relief supplies. In the last part of Formula (1), the quantity of the unmet supplies, the urgency degree, and the delay time are used to show that, in case of the same delay, the longer the delay time is and the higher the urgency degree of supplies is, the greater the impact is and consequently the higher the delay cost will be. Formula (2) denotes that the amount of supplies dispatched from the collector-distributor center is equal to the amount of resources supplied. Formula (3) represents the throughput limitation of emergency distribution centers. Formula (4) denotes the time spent on transporting supplies at the upper level from the collector-distributor center i to emergency distribution centers and damaged road segments by means of transportation m, which represents the total transit time of the transports with relief supplies to the disaster points. Formula (5) represents the time that the disaster-stricken location k needs for receiving emergency supply c. Formula (6) represents the time required for transporting rescue supplies from the collector-distributor center i to the segment in need of repair.


## 3.3.2. The Model for the Lower Level

The objective of the lower level is as follows:(7)MinZ2=nv1+v2∑j∈J∑k∈K∑m∈MXjkmDij,(8)s.t.∑j∈J∑k∈K∑m∈MXjkm≤n,(9)∑j∈JXjkm≤1,(10)amj=0,(11)amj+tjkm−amk≤B1−Xjkm,(12)0≤amkc−δmkc≤dlkc,(13)δmkc≤∑j∈JMXjkm,(14)Ymkc≤cmax∑k∈KXjkm,(15)∑m∈MYmkc+Ukc−qkc=0,(16)∑j∈J∑k∈Kjc=∑c∈C∑k∈K∑m∈MYmkc+∑c∈C∑l∈Llcξl,(17)Ymkc≥0,Ukc≥0,δmkc≥0,(18)yijmc,Xj,ξl,Xjkm∈0,1.The lower-level objective function is Formula (7), which represents the sum of the fixed cost and the convertible cost of the transportation of supplies from the lower-level distribution center to the disaster location. Formula (8) denotes the limited number of transport vehicles dispatched from the distribution center. Formula (9) indicates that each vehicle goes to the affected area only once at most. Formulas (10)–(13) guarantee the feasibility of the orderly driving of vehicles, which means that once a vehicle arrives at the supply-needed point later than the required time limit, there will be time delay. Formulas (14) and (15) ensure the feasibility of emergency supply logistics, that is, the supplies transported by distribution vehicles are less than the total amount of relief supplies in distribution centers. Formula (16) demonstrates the conservation of the quantity of supplies, indicating that the amount of supplies used to repair damaged road segments plus the amount of supplies used to meet the demand of supply-needed location equals the sum of all emergency supplies. Formula (17) represents nonnegative constraints. Formula (18) indicates 0-1 variable constraints. The identifiers B and M in the abovementioned formulas are defined as integers large enough.


## 4. Algorithm Design for Problem Solving

The model established in Section3 not only considers the change of different road networks and the problem of supply distribution but also involves the optimization of vehicle routing from the distribution center to each disaster location. It is a typical combination optimization problem aiming at achieving multiple objectives. In addition, because of the large number of variables involved in the model and the high complexity of the solution, it is difficult to solve this kind of NP-hard problem by general exact algorithms. Since artificial bee colony algorithm simulates bee colony’s behavior of foraging in nature and boasts high solving efficiency, simple algorithm, fewer parameters, and strong robustness [20, 21], this paper designs an improved artificial bee colony algorithm to solve the problem.


### 4.1. Standard Artificial Bee Colony Algorithm

ABC algorithm is a novel global optimization algorithm based on swarm intelligence proposed by Karaboga [22] in 2005, whose basic principle is inspired by the honey collecting behavior of the bee colony. Bees carry out different activities according to their respective work and duties and realize the sharing and exchange of colony information, so as to find the optimal solution to problems.The standard ABC algorithm divides the artificial bee colony into three categories by simulating the actual mechanism of honey gathering, namely, employed bees, onlooker bees, and scout bees. The goal of the whole bee colony is to find the largest honey source.


#### 4.1.1. Employed Bees

Each employed bee is responsible for a certain honey source (solution vector) and searches the neighborhood of the honey source in iterations.


#### 4.1.2. Onlooker Bees

According to the abundance degree of honey source (adaptive value), onlooker bees are employed by the method of roulette to collect honey (search for new honey source).


#### 4.1.3. Scout Bees

If the honey source has not been improved after several updates, it will be abandoned and the employed bees will be transformed into scout bees to search for new honey source randomly.In the process of initialization,SN feasible solutions (equal to the number of employed bees) are generated randomly and the function value of fitness is calculated. The formula for randomly generating feasible solutions is as follows:(19)xpq=xmin,q+rand0,1xmax,q−xmin,q,where xpq,q=1,2,…,SN is a D-dimensional vector, representing the feasible solution generated randomly, and D is the number of optimized parameters (q∈1,2,…,D). Bees record their optimum values so far and launch searches in the neighborhood of current honey sources. The basic ABC formula for searching new honey sources nearby is as follows:(20)vpq=xpq+φpqxpq−kpq,where vpq denotes the solution of the new honey source (q∈1,2,…,D, k∈1,2,…,SN). k is generated randomly, and k≠p. φpq represents a random number in [−1, 1]. Onlooker bees select employed bees according to the probability formula as follows:(21)Qp=fitxp∑n=1SNfitxn,where fitxp refers to the abundance degree of honey source corresponding to the fitness value of the pth solution. The more abundant the honey in the source, the higher the probability to be selected by onlooker bees.In order to prevent the algorithm from falling into local optimum, when a honey source subjects to no more improvement after limited times of iteration, the honey source will be abandoned and recorded in the tabu list. At the same time, the corresponding employed bees will be transformed into scout bees and a new location is randomly generated to replace the original honey source according to Formula (19).


### 4.2. Improved Artificial Bee Colony Algorithm and Its Implementation

ABC algorithm was initially applied to continuous optimization problems, but not to discrete problems. Therefore, this paper improves the existing bee colony algorithm to adapt it to discrete problems and focus on the main idea of large-scale heuristic neighborhood search. When the employed bees and the onlooker bees in the ABC algorithm search for the honey source near the food source, the employed bees implement the local search for the best alternative and then rank their fitness values and retain the elite solutions. After that, by means of the crossover operator and mutation operator in genetic algorithm, employed bees will adjust the search range through self-adaption and improve the speed and accuracy of solving process, thus avoiding falling into local optimum and accelerating the rate of convergence.In addition, in order to enhance the search speed and convergence accuracy of the standard ABC algorithm, this paper improves the search strategy of the follower bees and the reconnaissance strategy of the scout bees, so as to strengthen the algorithm’s function of searching local areas, expand the search scope, and avoid the trap of local optimum. In the process of global searching, since the encoding method is integer number encoding, which is different from other methods such as the binary system, the neighborhood formula of standard ABC algorithm can be adopted to generate candidate solutions. In this paper, several common global search strategies are adopted. To be specific, firstly, the order of three locations is changed randomly, and then two search operators, namely, neighborhood inversion and neighborhood exchange, are used for global search. By combining these two search operators, the search range can be enlarged, the diversity of honey sources can be increased, and the global search ability can also be promoted.The emergency supply scheduling model with “soft time window” built in this paper considers multiple objective functions synthetically, making it difficult to achieve several optimal solutions at the same time, but only obtaining the feasible solution to the problem. Because emergency logistics is closely related to the importance of rescue itself, it features relatively weak economical efficiency. However, in the process of emergency rescue, the speed of rescue and the satisfaction rate of rescue supplies are of the most importance, so the method ofepsilon constraint [23] is used to solve the problem of multiobjective programming. Its basic idea is to regard the objective function Z1 as the main objective, then set an acceptable threshold value ε, and provide the prerequisite that Z2 < ε. If a proper threshold value ε is selected, the objective function Z2 will be transferred into the constraint and this optimization problem will become a single objective one, thus the efficient solution to the problem can be obtained.The concrete operating process of the algorithm is as follows.


#### 4.2.1. Encoding Strategy

According to the characteristics of the problem, first of all, the bee colonies are encoded by the integer-code method. The length corresponds to the total number of supply-needed disaster-stricken locations, and the solution range of the honey source is the set of integer serial numbers of all supply-needed locations. Concrete coding rules are as follows.Assume that there are two middle-level distribution centers and 15 supply-needed locations, which are labeled as number “1” and “2” and number “3, 4, 5, …, 16, 17,” respectively. First of all, a feasible solution with the length of 15 should be randomly constructed. After that, the distances between the first disaster location and the two distribution centers should be compared, and the number of the distribution center closer to the first disaster location should be inserted in the front of the code. Subsequently, accumulate the amount of the emergency supplies needed at the disaster-stricken locations on the route from the distribution center and compare the accumulated value with the maximum load of the vehicle. If the weight of relief supplies needed at several successive disaster locations along the route exceeds the maximum vehicle load, the number of distribution center will be inserted at the end of the route to indicate that the segment is coded as a vehicle driving route. Repeat the abovementioned steps until all the encoding information of honey source is processed. The encoding process is described in Figure4.Figure 4
An example of the coding process.An example of coding process shown in Figure4 contains three subpaths. The first subpath is “1-5-16-8-14-3-1,” which means that the vehicle starts from the distribution center No. 1 and returns to it after visiting “5-16-8-14-3.” Similarly, the second and third subpaths are “2-6-11-17-7-15-2” and “1-9-12-13-10-4-1.”


#### 4.2.2. Fitness Function

Fitness function is transformed from two objective functions contained in the model. Considering the characteristics and constraints of emergency logistics problems, the objective function is not directly set as the fitness function. The specific expression is as follows:(22)Fit=min1Z1,1Z1+B+λmaxdl−T,0,where the first part of Formula (22) denotes to apply penalty function to deal with constraints. If the objective function Z2 is larger than the set threshold value ε, a larger positive number is added to the objective function Z1 to punish the constraints, and then the reciprocal of the objective function is counted into the fitness function. This method can better highlight the difference and strengthen the selection function of the algorithm. The second part is designed to strengthen the impact of delay, in which λ is a large positive integer ensuring the diversity of chromosomes.


#### 4.2.3. Self-Adaptive Search Strategy

In order to reduce the blind search of bee colonies and improve the convergence speed, this paper adopts the crossover operator and mutation operator of genetic algorithm to improve the reconnaissance strategy of the scout bees and the selection mechanism of the follower bees. Subsequently, the order is arranged in accordance with the fitness degree, and the probability of the honey source being selected by the follower bees is calculated. The selection method of roulette is still used here, which is a method of random selection, incapable of guaranteeing that individuals with better adaptability in the contemporary population can be preserved. In order to improve this situation, the method ofelite protection individuals is added to the process of roulette selection. That means, the first two best individuals in each generation are directly selected and retained, and then the roulette is adopted to choose other individuals in the population, so as to ensure the convergence of the algorithm. When the follower bees select honey source according to the information provided by the employed bees, the self-adaptive probability is considered to carry out selective crossover and mutation. However, the crossover and mutation probability does not take the determined value in 0,1, but are dynamically adjusted according to the fitness value.The self-adaptive crossover probability can be obtained by(23)ηm=11+expkmin−f,where kk>0 refers to the initial crossover probability, min to the minimum of individual fitness, f to the average individual fitness value smaller than the average fitness value, and ηm to the crossover probability of adaptive change. In Formula (23), since min−f<0, the range of exponential function value expkmin−f belongs to 0,1. The larger the min−f is, the more dispersed the individual fitness value of the population is, which requires to generate new individuals through crossover so as to enlarge the search range. On the contrary, it indicates that the fitness value of individual is close to each other, so the crossover probability should be decreased to preserve excellent individuals.The self-adaptive mutation probability can be obtained by(24)ηn=1−11+expkmin−f,where kk>0 refers to the initial mutation probability, min to the minimum of individual fitness, f to the average individual fitness value smaller than the average fitness value, and ηn to the self-adaptive mutation probability. In Formula (24), since min−f<0, the range of exponential function value expkmin−f belongs to 0,1. When the fitness values are closer to each other, ηn is closer to 0.5 and it is easier to search for new individuals. When the fitness values differ greatly, the mutation probability is decreased in order to avoid damaging the good individual solutions and ensure smooth convergence.


#### 4.2.4. Algorithm Steps

To sum up, the detailed steps of the improved ABC algorithm are as follows:Step 1: initialize the population. Firstly, determine the population size and parameters. Set the colony size asSN, the number of employed bees as 50% of the total colony size, i.e., SN/2, onlooker bees as 50% of the total colony size, and the initial probability of crossover and mutation is k.Step 2: generateSN initial honey sources (i.e., feasible solutions) according to the encoding method set by ABC algorithm.Step 3: calculate the fitness values of all feasible solutions and rank them.Step 4: calculate the probability of each honey source being selected by the employed bees and use the method of roulette to determine which honey source the employed bees will choose.Step 5: calculate and rank the fitness of all honey sources, and then directly copy and preserve the first two individuals with the best fitness value.Step 6: calculate the adaptive crossover probability of each honey source and perform the crossover operation.Step 7: calculate the adaptive mutation probability of each honey source and perform the mutation operation.Step 8: evaluate the fitness value of the new honey source and apply the greed principle. If the new honey source is better than the original feasible solution, it will replace the original one and participate in the next round operation, otherwise the original honey source will be retained in the next round of operation.Step 9: if the honey source subjects to no improvement afterlimit (limit is the parameter used to prevent bees from falling into local optimization solutions) times of updates, it will be abandoned and the employed bees will be transformed into scout bees to search for new honey sources.Step 10: judge whether the maximum number of iterations has been reached. If so, terminate the process, otherwise move back to step 5.Based on the above analysis, the operational process of the improved artificial bee colony algorithm adopted in this paper is presented in Figure5.Figure 5
Improved artificial bee colony algorithm.


## 4.1. Standard Artificial Bee Colony Algorithm

ABC algorithm is a novel global optimization algorithm based on swarm intelligence proposed by Karaboga [22] in 2005, whose basic principle is inspired by the honey collecting behavior of the bee colony. Bees carry out different activities according to their respective work and duties and realize the sharing and exchange of colony information, so as to find the optimal solution to problems.The standard ABC algorithm divides the artificial bee colony into three categories by simulating the actual mechanism of honey gathering, namely, employed bees, onlooker bees, and scout bees. The goal of the whole bee colony is to find the largest honey source.


### 4.1.1. Employed Bees

Each employed bee is responsible for a certain honey source (solution vector) and searches the neighborhood of the honey source in iterations.


### 4.1.2. Onlooker Bees

According to the abundance degree of honey source (adaptive value), onlooker bees are employed by the method of roulette to collect honey (search for new honey source).


### 4.1.3. Scout Bees

If the honey source has not been improved after several updates, it will be abandoned and the employed bees will be transformed into scout bees to search for new honey source randomly.In the process of initialization,SN feasible solutions (equal to the number of employed bees) are generated randomly and the function value of fitness is calculated. The formula for randomly generating feasible solutions is as follows:(19)xpq=xmin,q+rand0,1xmax,q−xmin,q,where xpq,q=1,2,…,SN is a D-dimensional vector, representing the feasible solution generated randomly, and D is the number of optimized parameters (q∈1,2,…,D). Bees record their optimum values so far and launch searches in the neighborhood of current honey sources. The basic ABC formula for searching new honey sources nearby is as follows:(20)vpq=xpq+φpqxpq−kpq,where vpq denotes the solution of the new honey source (q∈1,2,…,D, k∈1,2,…,SN). k is generated randomly, and k≠p. φpq represents a random number in [−1, 1]. Onlooker bees select employed bees according to the probability formula as follows:(21)Qp=fitxp∑n=1SNfitxn,where fitxp refers to the abundance degree of honey source corresponding to the fitness value of the pth solution. The more abundant the honey in the source, the higher the probability to be selected by onlooker bees.In order to prevent the algorithm from falling into local optimum, when a honey source subjects to no more improvement after limited times of iteration, the honey source will be abandoned and recorded in the tabu list. At the same time, the corresponding employed bees will be transformed into scout bees and a new location is randomly generated to replace the original honey source according to Formula (19).


## 4.1.1. Employed Bees

Each employed bee is responsible for a certain honey source (solution vector) and searches the neighborhood of the honey source in iterations.


## 4.1.2. Onlooker Bees

According to the abundance degree of honey source (adaptive value), onlooker bees are employed by the method of roulette to collect honey (search for new honey source).


## 4.1.3. Scout Bees

If the honey source has not been improved after several updates, it will be abandoned and the employed bees will be transformed into scout bees to search for new honey source randomly.In the process of initialization,SN feasible solutions (equal to the number of employed bees) are generated randomly and the function value of fitness is calculated. The formula for randomly generating feasible solutions is as follows:(19)xpq=xmin,q+rand0,1xmax,q−xmin,q,where xpq,q=1,2,…,SN is a D-dimensional vector, representing the feasible solution generated randomly, and D is the number of optimized parameters (q∈1,2,…,D). Bees record their optimum values so far and launch searches in the neighborhood of current honey sources. The basic ABC formula for searching new honey sources nearby is as follows:(20)vpq=xpq+φpqxpq−kpq,where vpq denotes the solution of the new honey source (q∈1,2,…,D, k∈1,2,…,SN). k is generated randomly, and k≠p. φpq represents a random number in [−1, 1]. Onlooker bees select employed bees according to the probability formula as follows:(21)Qp=fitxp∑n=1SNfitxn,where fitxp refers to the abundance degree of honey source corresponding to the fitness value of the pth solution. The more abundant the honey in the source, the higher the probability to be selected by onlooker bees.In order to prevent the algorithm from falling into local optimum, when a honey source subjects to no more improvement after limited times of iteration, the honey source will be abandoned and recorded in the tabu list. At the same time, the corresponding employed bees will be transformed into scout bees and a new location is randomly generated to replace the original honey source according to Formula (19).


## 4.2. Improved Artificial Bee Colony Algorithm and Its Implementation

ABC algorithm was initially applied to continuous optimization problems, but not to discrete problems. Therefore, this paper improves the existing bee colony algorithm to adapt it to discrete problems and focus on the main idea of large-scale heuristic neighborhood search. When the employed bees and the onlooker bees in the ABC algorithm search for the honey source near the food source, the employed bees implement the local search for the best alternative and then rank their fitness values and retain the elite solutions. After that, by means of the crossover operator and mutation operator in genetic algorithm, employed bees will adjust the search range through self-adaption and improve the speed and accuracy of solving process, thus avoiding falling into local optimum and accelerating the rate of convergence.In addition, in order to enhance the search speed and convergence accuracy of the standard ABC algorithm, this paper improves the search strategy of the follower bees and the reconnaissance strategy of the scout bees, so as to strengthen the algorithm’s function of searching local areas, expand the search scope, and avoid the trap of local optimum. In the process of global searching, since the encoding method is integer number encoding, which is different from other methods such as the binary system, the neighborhood formula of standard ABC algorithm can be adopted to generate candidate solutions. In this paper, several common global search strategies are adopted. To be specific, firstly, the order of three locations is changed randomly, and then two search operators, namely, neighborhood inversion and neighborhood exchange, are used for global search. By combining these two search operators, the search range can be enlarged, the diversity of honey sources can be increased, and the global search ability can also be promoted.The emergency supply scheduling model with “soft time window” built in this paper considers multiple objective functions synthetically, making it difficult to achieve several optimal solutions at the same time, but only obtaining the feasible solution to the problem. Because emergency logistics is closely related to the importance of rescue itself, it features relatively weak economical efficiency. However, in the process of emergency rescue, the speed of rescue and the satisfaction rate of rescue supplies are of the most importance, so the method ofepsilon constraint [23] is used to solve the problem of multiobjective programming. Its basic idea is to regard the objective function Z1 as the main objective, then set an acceptable threshold value ε, and provide the prerequisite that Z2 < ε. If a proper threshold value ε is selected, the objective function Z2 will be transferred into the constraint and this optimization problem will become a single objective one, thus the efficient solution to the problem can be obtained.The concrete operating process of the algorithm is as follows.


### 4.2.1. Encoding Strategy

According to the characteristics of the problem, first of all, the bee colonies are encoded by the integer-code method. The length corresponds to the total number of supply-needed disaster-stricken locations, and the solution range of the honey source is the set of integer serial numbers of all supply-needed locations. Concrete coding rules are as follows.Assume that there are two middle-level distribution centers and 15 supply-needed locations, which are labeled as number “1” and “2” and number “3, 4, 5, …, 16, 17,” respectively. First of all, a feasible solution with the length of 15 should be randomly constructed. After that, the distances between the first disaster location and the two distribution centers should be compared, and the number of the distribution center closer to the first disaster location should be inserted in the front of the code. Subsequently, accumulate the amount of the emergency supplies needed at the disaster-stricken locations on the route from the distribution center and compare the accumulated value with the maximum load of the vehicle. If the weight of relief supplies needed at several successive disaster locations along the route exceeds the maximum vehicle load, the number of distribution center will be inserted at the end of the route to indicate that the segment is coded as a vehicle driving route. Repeat the abovementioned steps until all the encoding information of honey source is processed. The encoding process is described in Figure4.Figure 4
An example of the coding process.An example of coding process shown in Figure4 contains three subpaths. The first subpath is “1-5-16-8-14-3-1,” which means that the vehicle starts from the distribution center No. 1 and returns to it after visiting “5-16-8-14-3.” Similarly, the second and third subpaths are “2-6-11-17-7-15-2” and “1-9-12-13-10-4-1.”


### 4.2.2. Fitness Function

Fitness function is transformed from two objective functions contained in the model. Considering the characteristics and constraints of emergency logistics problems, the objective function is not directly set as the fitness function. The specific expression is as follows:(22)Fit=min1Z1,1Z1+B+λmaxdl−T,0,where the first part of Formula (22) denotes to apply penalty function to deal with constraints. If the objective function Z2 is larger than the set threshold value ε, a larger positive number is added to the objective function Z1 to punish the constraints, and then the reciprocal of the objective function is counted into the fitness function. This method can better highlight the difference and strengthen the selection function of the algorithm. The second part is designed to strengthen the impact of delay, in which λ is a large positive integer ensuring the diversity of chromosomes.


### 4.2.3. Self-Adaptive Search Strategy

In order to reduce the blind search of bee colonies and improve the convergence speed, this paper adopts the crossover operator and mutation operator of genetic algorithm to improve the reconnaissance strategy of the scout bees and the selection mechanism of the follower bees. Subsequently, the order is arranged in accordance with the fitness degree, and the probability of the honey source being selected by the follower bees is calculated. The selection method of roulette is still used here, which is a method of random selection, incapable of guaranteeing that individuals with better adaptability in the contemporary population can be preserved. In order to improve this situation, the method ofelite protection individuals is added to the process of roulette selection. That means, the first two best individuals in each generation are directly selected and retained, and then the roulette is adopted to choose other individuals in the population, so as to ensure the convergence of the algorithm. When the follower bees select honey source according to the information provided by the employed bees, the self-adaptive probability is considered to carry out selective crossover and mutation. However, the crossover and mutation probability does not take the determined value in 0,1, but are dynamically adjusted according to the fitness value.The self-adaptive crossover probability can be obtained by(23)ηm=11+expkmin−f,where kk>0 refers to the initial crossover probability, min to the minimum of individual fitness, f to the average individual fitness value smaller than the average fitness value, and ηm to the crossover probability of adaptive change. In Formula (23), since min−f<0, the range of exponential function value expkmin−f belongs to 0,1. The larger the min−f is, the more dispersed the individual fitness value of the population is, which requires to generate new individuals through crossover so as to enlarge the search range. On the contrary, it indicates that the fitness value of individual is close to each other, so the crossover probability should be decreased to preserve excellent individuals.The self-adaptive mutation probability can be obtained by(24)ηn=1−11+expkmin−f,where kk>0 refers to the initial mutation probability, min to the minimum of individual fitness, f to the average individual fitness value smaller than the average fitness value, and ηn to the self-adaptive mutation probability. In Formula (24), since min−f<0, the range of exponential function value expkmin−f belongs to 0,1. When the fitness values are closer to each other, ηn is closer to 0.5 and it is easier to search for new individuals. When the fitness values differ greatly, the mutation probability is decreased in order to avoid damaging the good individual solutions and ensure smooth convergence.


### 4.2.4. Algorithm Steps

To sum up, the detailed steps of the improved ABC algorithm are as follows:Step 1: initialize the population. Firstly, determine the population size and parameters. Set the colony size asSN, the number of employed bees as 50% of the total colony size, i.e., SN/2, onlooker bees as 50% of the total colony size, and the initial probability of crossover and mutation is k.Step 2: generateSN initial honey sources (i.e., feasible solutions) according to the encoding method set by ABC algorithm.Step 3: calculate the fitness values of all feasible solutions and rank them.Step 4: calculate the probability of each honey source being selected by the employed bees and use the method of roulette to determine which honey source the employed bees will choose.Step 5: calculate and rank the fitness of all honey sources, and then directly copy and preserve the first two individuals with the best fitness value.Step 6: calculate the adaptive crossover probability of each honey source and perform the crossover operation.Step 7: calculate the adaptive mutation probability of each honey source and perform the mutation operation.Step 8: evaluate the fitness value of the new honey source and apply the greed principle. If the new honey source is better than the original feasible solution, it will replace the original one and participate in the next round operation, otherwise the original honey source will be retained in the next round of operation.Step 9: if the honey source subjects to no improvement afterlimit (limit is the parameter used to prevent bees from falling into local optimization solutions) times of updates, it will be abandoned and the employed bees will be transformed into scout bees to search for new honey sources.Step 10: judge whether the maximum number of iterations has been reached. If so, terminate the process, otherwise move back to step 5.Based on the above analysis, the operational process of the improved artificial bee colony algorithm adopted in this paper is presented in Figure5.Figure 5
Improved artificial bee colony algorithm.


## 4.2.1. Encoding Strategy

According to the characteristics of the problem, first of all, the bee colonies are encoded by the integer-code method. The length corresponds to the total number of supply-needed disaster-stricken locations, and the solution range of the honey source is the set of integer serial numbers of all supply-needed locations. Concrete coding rules are as follows.Assume that there are two middle-level distribution centers and 15 supply-needed locations, which are labeled as number “1” and “2” and number “3, 4, 5, …, 16, 17,” respectively. First of all, a feasible solution with the length of 15 should be randomly constructed. After that, the distances between the first disaster location and the two distribution centers should be compared, and the number of the distribution center closer to the first disaster location should be inserted in the front of the code. Subsequently, accumulate the amount of the emergency supplies needed at the disaster-stricken locations on the route from the distribution center and compare the accumulated value with the maximum load of the vehicle. If the weight of relief supplies needed at several successive disaster locations along the route exceeds the maximum vehicle load, the number of distribution center will be inserted at the end of the route to indicate that the segment is coded as a vehicle driving route. Repeat the abovementioned steps until all the encoding information of honey source is processed. The encoding process is described in Figure4.Figure 4
An example of the coding process.An example of coding process shown in Figure4 contains three subpaths. The first subpath is “1-5-16-8-14-3-1,” which means that the vehicle starts from the distribution center No. 1 and returns to it after visiting “5-16-8-14-3.” Similarly, the second and third subpaths are “2-6-11-17-7-15-2” and “1-9-12-13-10-4-1.”


## 4.2.2. Fitness Function

Fitness function is transformed from two objective functions contained in the model. Considering the characteristics and constraints of emergency logistics problems, the objective function is not directly set as the fitness function. The specific expression is as follows:(22)Fit=min1Z1,1Z1+B+λmaxdl−T,0,where the first part of Formula (22) denotes to apply penalty function to deal with constraints. If the objective function Z2 is larger than the set threshold value ε, a larger positive number is added to the objective function Z1 to punish the constraints, and then the reciprocal of the objective function is counted into the fitness function. This method can better highlight the difference and strengthen the selection function of the algorithm. The second part is designed to strengthen the impact of delay, in which λ is a large positive integer ensuring the diversity of chromosomes.


## 4.2.3. Self-Adaptive Search Strategy

In order to reduce the blind search of bee colonies and improve the convergence speed, this paper adopts the crossover operator and mutation operator of genetic algorithm to improve the reconnaissance strategy of the scout bees and the selection mechanism of the follower bees. Subsequently, the order is arranged in accordance with the fitness degree, and the probability of the honey source being selected by the follower bees is calculated. The selection method of roulette is still used here, which is a method of random selection, incapable of guaranteeing that individuals with better adaptability in the contemporary population can be preserved. In order to improve this situation, the method ofelite protection individuals is added to the process of roulette selection. That means, the first two best individuals in each generation are directly selected and retained, and then the roulette is adopted to choose other individuals in the population, so as to ensure the convergence of the algorithm. When the follower bees select honey source according to the information provided by the employed bees, the self-adaptive probability is considered to carry out selective crossover and mutation. However, the crossover and mutation probability does not take the determined value in 0,1, but are dynamically adjusted according to the fitness value.The self-adaptive crossover probability can be obtained by(23)ηm=11+expkmin−f,where kk>0 refers to the initial crossover probability, min to the minimum of individual fitness, f to the average individual fitness value smaller than the average fitness value, and ηm to the crossover probability of adaptive change. In Formula (23), since min−f<0, the range of exponential function value expkmin−f belongs to 0,1. The larger the min−f is, the more dispersed the individual fitness value of the population is, which requires to generate new individuals through crossover so as to enlarge the search range. On the contrary, it indicates that the fitness value of individual is close to each other, so the crossover probability should be decreased to preserve excellent individuals.The self-adaptive mutation probability can be obtained by(24)ηn=1−11+expkmin−f,where kk>0 refers to the initial mutation probability, min to the minimum of individual fitness, f to the average individual fitness value smaller than the average fitness value, and ηn to the self-adaptive mutation probability. In Formula (24), since min−f<0, the range of exponential function value expkmin−f belongs to 0,1. When the fitness values are closer to each other, ηn is closer to 0.5 and it is easier to search for new individuals. When the fitness values differ greatly, the mutation probability is decreased in order to avoid damaging the good individual solutions and ensure smooth convergence.


## 4.2.4. Algorithm Steps

To sum up, the detailed steps of the improved ABC algorithm are as follows:Step 1: initialize the population. Firstly, determine the population size and parameters. Set the colony size asSN, the number of employed bees as 50% of the total colony size, i.e., SN/2, onlooker bees as 50% of the total colony size, and the initial probability of crossover and mutation is k.Step 2: generateSN initial honey sources (i.e., feasible solutions) according to the encoding method set by ABC algorithm.Step 3: calculate the fitness values of all feasible solutions and rank them.Step 4: calculate the probability of each honey source being selected by the employed bees and use the method of roulette to determine which honey source the employed bees will choose.Step 5: calculate and rank the fitness of all honey sources, and then directly copy and preserve the first two individuals with the best fitness value.Step 6: calculate the adaptive crossover probability of each honey source and perform the crossover operation.Step 7: calculate the adaptive mutation probability of each honey source and perform the mutation operation.Step 8: evaluate the fitness value of the new honey source and apply the greed principle. If the new honey source is better than the original feasible solution, it will replace the original one and participate in the next round operation, otherwise the original honey source will be retained in the next round of operation.Step 9: if the honey source subjects to no improvement afterlimit (limit is the parameter used to prevent bees from falling into local optimization solutions) times of updates, it will be abandoned and the employed bees will be transformed into scout bees to search for new honey sources.Step 10: judge whether the maximum number of iterations has been reached. If so, terminate the process, otherwise move back to step 5.Based on the above analysis, the operational process of the improved artificial bee colony algorithm adopted in this paper is presented in Figure5.Figure 5
Improved artificial bee colony algorithm.


## 5. A Case Study

In this section, in order to verify the efficiency and effectiveness of proposed model and algorithm, an actual case is analyzed by means of MATLABR2016b.


### 5.1. Background

After a sudden-onset earthquake disaster in a certain area, 20 disaster-stricken locations are waiting for emergency rescue supplies and two distribution centers in this area can play the role middle-level distribution centers. In addition, there is a collector-distributor center to allocate supplies. Two distribution centers in the middle level are numbered “1” and “2” in this case. The disaster-stricken locations are numbered “3” to “22.” Moreover, it is necessary to take into account the possibility of road obstruction due to collapse or cracks. The descriptions above are shown in Figure6.The locations of all distribution centers, supply-needed locations, and upper-level supply collector-distributor center are known, including the number, location coordinate of all distribution centers and stricken locations, time window requirements for relief supplies at all disaster locations, and the amount of emergency supplies each disaster-stricken location needs. Disaster-stricken people and rescue workers in disaster locations are in urgent need of various foods, drinking water, and other daily necessities, as well as the tools and supplies needed to repair the road network. A variety of transportation modes can be adopted to transport supplies from the collector-distributor center to distribution centers and damaged road networks. But only vehicles can be used to transport supplies from distribution centers to various disaster locations due to condition limitation. The maximum load capacity of vehicles is known. Due to the inaccessibility of part of the road networks, it is necessary to allocate a certain amount of supplies for emergency repair in order to significantly reduce the distribution delay and to pass the road within the normal time. In view of the scarcity of relief supplies in the early stage of the earthquake, the amount of relief supplies should be as many as the living supplies needed in the disaster area. The starting time of transportation task is 0, and the earliest time window and the latest time window for relief supplies at the disaster locations depend on the actual situation. The concrete information is shown in Table2.Table 2
Transportation network information table.
No.
Coordinate position
Earliest time window
Latest time window
Required relief supplies
1
(0, 0)
0
2
(30, 30)
0
3
(20, 24)
0.2
1.6
15
4
(4, 1)
0.55
2.4
50
5
(2, 15)
0.05
1.8
30
6
(10, 7)
0.25
1.4
30
7
(10, 27)
0.35
3
40
8
(13, 8)
0.5
3.6
40
9
(6, 26)
0.95
2.8
10
10
(14, 16)
0.3
2
20
11
(12, 10)
0.45
3.4
25
12
(23, 27)
0.1
2.8
60
13
(14, 27)
0.8
4
40
14
(5, 24)
0.65
2.8
10
15
(12, 2)
0.75
2.2
10
16
(15, 12)
0.9
3.8
30
17
(11, 30)
0.4
3.2
35
18
(18, 8)
1
2.2
45
19
(5, 16)
0.6
1.8
50
20
(15, 15)
0.15
1.6
10
21
(25, 26)
0.85
2.6
30
22
(19, 26)
0.7
2.4
30


### 5.2. Parameter Setting

In this case, the parameters are set as follows. 600 units of emergency supplies are prepared in the collector-distributor center to ensure the needs of disaster areas. The fixed cost of a vehicle is 300 yuan and the driving cost is 5 yuan per hour. The maximum load of a vehicle is 150 units and the average speed is 35 km/h. The penalty cost is 200 yuan per hour when the time window requirement of the disaster point is not met. The two damaged road segments need 12 units and 18 units to get restored, respectively. If the damaged segment is not repaired and the road is impassable, the distance of the segment is set to be a large enough integerM, which is 1000 in this case. On the basis of the model built in Section 3 and the actual situation, the parameters of the adaptive artificial bee colony algorithm are set as follows. The population size is set as 200, the initial probability rate of adaptive crossover and mutation as 0.01, the maximum number of iterations as 100, the maximum number of searches as 100, and the program runs 20 times to get the optimal results.


### 5.3. Result Analysis

#### 5.3.1. No Repairs of the Road Network

First of all, the situation that all damaged segments are not repaired should be considered, that is, all relief supplies are distributed to the disaster-stricken areas and no supplies are allocated to the damaged segment. It is assumed that Segment 1 between Supply-Needed Location 6 and Location 11 and Segment 2 between Distribution Center 2 and Supply-Needed Location 21 in the transportation network are blocked (shown as the dotted line in Figure6), and a penalty function is set as 1000. In the case of road congestion, the optimal path of the emergency scheduling problem is [15-18-16-21-8-11-6-12-13-7-3-10-5-9-22-4-14-19-17]. A total of five vehicles are needed to complete the whole transportation plan, as shown in Table 3.Figure 6
Road network structure.Table 3
Subpath transport route without road repair.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
1-15-18-16-21-1
54.29
0.97
78.92
2
1-8-11-6-1
541.25
27.73
1029.707
3
2-12-13-7-2
40.84
0.33
40.84
4
2-3-10-5-9-22-20-2
77.60
0.40
91.33
5
1-4-14-19-17-1
65.17
0.43
83.32


#### 5.3.2. Repairs of the Partial Road Network

Under the situation that Segment 1 is repaired urgently and Segment 2 is not repaired, the optimal path is [6-11-13-9-4-19-20-21-12-22-7-3-15-18-10-16-8-5-17-14]. Five vehicles are needed to complete the whole transportation plan, as shown in Table4.Table 4
Subpath transport route on the repairing number 1.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
1-6-11-13-9-1
42.38
0
67.67
2
1-4-19-20-21-1
45.33
0.48
80.14
3
2-12-22-3-2
35.27
0.36
42.90
4
1-15-18-10-16-1
46.24
0.81
52.93
5
1-8-5-17-14-1
61.39
0
322.44Under the situation that Segment 2 is repaired urgently and Segment 1 is not repaired, the optimal path is [7-13-12-5-19-10-20-16-6-11-8-4-3-14-9-15-18-22-17-21]. Five vehicles are also needed to complete the whole transportation plan, as shown in Table5.Table 5
Subpath transport route on the repairing number 2.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
2-3-22-10-16-14-2
41.03
0.37
70.53
2
1-4-18-11-6-1
551.76
28.78
1038.31
3
2-21-13-17-9-2
32.76
1.12
52.43
4
2-12-20-5-7-2
62.26
0
69.68
5
1-15-19-8-1
50.90
0.40
54.396


#### 5.3.3. Repairs of the Whole Road Network

Under the situation that the whole road network is repaired urgently, the optimal path for emergency supply scheduling is [21-10-19-5-14-6-11-8-16-12-13-22-4-15-7-17-9-3-18-20]. Five vehicles are needed to complete the whole transportation plan, as shown in Table6.Table 6
Subpath transport route with all roads repaired.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
2-21-10-19-5-14-2
57.51
0.67
68.63
2
1-6-11-8-16-1
36.39
0.26
41.73
3
2-12-13-22-2
29.25
0.40
33.42
4
1-4-15-7-17-9-1
73.51
0.83
73.51
5
2-3-18-20-2
30
0.21
56.62Compared the value of objective functionZ1 and the transportation cost Z2 in different situations, the results are shown in Table 7.Table 7
Comparison of target values under four road networks.
Objective function
Repair road number 1
Repair road number 2
Repair all
Not repaired all
Z
1
1243.03
2069.97
266.76
1800.9
Z
2
3112.18
7926.72
2869.53
8115.63In Table7, the ascending order of the objective function values Z1 is as follows: Z1 in the situation that all road networks are repaired, Z1 in the situation that Segment 1 is repaired, Z1 in the situation that no road segment is repaired, and Z1 in the situation that Segment 2 is repaired. The data obtained in the situation of repairing all road networks are obviously smaller than those of repairing no road network or only repairing partial road networks. When all road networks are restored, the target value Z1 is smaller than that under other situations. This is because the fitness function has adjusted the target value and enlarged it when it does not meet the threshold value, which is conducive to enhance the difference among algorithm selections.Similarly, in Table7, the ascending order of the transportation cost Z2 in different road networks is as follows: Z2 in the situation that all road networks are repaired, Z2 in the situation that Segment 1 is repaired, Z2 in the situation that Segment 2 is repaired, and Z2 in the situation that no road segment is repaired. The transportation cost of the road networks that undergoes complete repair is lower than that of the road networks that are not repaired. It is because that the two damaged segments are both the shortest key route between distribution centers to the disaster-stricken areas. Once damaged, transportation costs will rise dramatically and transportation efficiency will decline. The benefit of repairing the roads is greater than the cost, so all the results obtained in the situation of repairing all road networks are optimal. If the relief supplies are only enough to restore part of road segments, it is advisable to repair the road segments which can produce greater benefits. That means the supplies can be sent to more disaster-stricken locations in the shortest time through the chosen road segments. Therefore, it is clear that in the case of road-varying network, it is of great practical value and significance to analyze and compare the whole road network and explore the path cost and benefit under different road networks.


### 5.4. Performance Comparison of Different Algorithms

Few similar examples and algorithms are available at present to handle the problem of scheduling emergency supplies under conditions of multiple distribution centers and road-varying network, so the algorithms cannot be compared with others directly. In order to analyze the performance of this algorithm, in the path optimization stage, the genetic algorithm and the standard ABC algorithm are used to compare with the improved ABC algorithm under the same condition that the whole road network is not restored. In consideration of the difference between genetic algorithm and the standard ABC algorithm, different population numbers after 500 times of iteration are compared, and the results of comparative tests are shown in Table8 and Figure 7.Table 8
Performance comparison between improved ABC algorithm and others.
Population
Improved ABC algorithm
Genetic algorithm
ABC algorithm
100
796
763
814
200
1496
1618
1589
500
3569
4015
3863Figure 7
Comparison of different optimization algorithms.Through comparison, it is obvious that the genetic algorithm has a slight advantage over the improved ABC algorithm when the population size is small, but when the population size increases gradually, the algorithm proposed in this paper significantly overcomes the genetic algorithm. Comparing the improved ABC algorithm with the standard one, we find that the former has advantages in terms of any population sizes. In addition, in Figure7, comparing the convergence speeds and optimal target values using these three algorithms, we can find that the improved ABC algorithm has faster convergence speed and better target value, which proves that the strategy of adaption in this paper is prominently effective.


## 5.1. Background

After a sudden-onset earthquake disaster in a certain area, 20 disaster-stricken locations are waiting for emergency rescue supplies and two distribution centers in this area can play the role middle-level distribution centers. In addition, there is a collector-distributor center to allocate supplies. Two distribution centers in the middle level are numbered “1” and “2” in this case. The disaster-stricken locations are numbered “3” to “22.” Moreover, it is necessary to take into account the possibility of road obstruction due to collapse or cracks. The descriptions above are shown in Figure6.The locations of all distribution centers, supply-needed locations, and upper-level supply collector-distributor center are known, including the number, location coordinate of all distribution centers and stricken locations, time window requirements for relief supplies at all disaster locations, and the amount of emergency supplies each disaster-stricken location needs. Disaster-stricken people and rescue workers in disaster locations are in urgent need of various foods, drinking water, and other daily necessities, as well as the tools and supplies needed to repair the road network. A variety of transportation modes can be adopted to transport supplies from the collector-distributor center to distribution centers and damaged road networks. But only vehicles can be used to transport supplies from distribution centers to various disaster locations due to condition limitation. The maximum load capacity of vehicles is known. Due to the inaccessibility of part of the road networks, it is necessary to allocate a certain amount of supplies for emergency repair in order to significantly reduce the distribution delay and to pass the road within the normal time. In view of the scarcity of relief supplies in the early stage of the earthquake, the amount of relief supplies should be as many as the living supplies needed in the disaster area. The starting time of transportation task is 0, and the earliest time window and the latest time window for relief supplies at the disaster locations depend on the actual situation. The concrete information is shown in Table2.Table 2
Transportation network information table.
No.
Coordinate position
Earliest time window
Latest time window
Required relief supplies
1
(0, 0)
0
2
(30, 30)
0
3
(20, 24)
0.2
1.6
15
4
(4, 1)
0.55
2.4
50
5
(2, 15)
0.05
1.8
30
6
(10, 7)
0.25
1.4
30
7
(10, 27)
0.35
3
40
8
(13, 8)
0.5
3.6
40
9
(6, 26)
0.95
2.8
10
10
(14, 16)
0.3
2
20
11
(12, 10)
0.45
3.4
25
12
(23, 27)
0.1
2.8
60
13
(14, 27)
0.8
4
40
14
(5, 24)
0.65
2.8
10
15
(12, 2)
0.75
2.2
10
16
(15, 12)
0.9
3.8
30
17
(11, 30)
0.4
3.2
35
18
(18, 8)
1
2.2
45
19
(5, 16)
0.6
1.8
50
20
(15, 15)
0.15
1.6
10
21
(25, 26)
0.85
2.6
30
22
(19, 26)
0.7
2.4
30


## 5.2. Parameter Setting

In this case, the parameters are set as follows. 600 units of emergency supplies are prepared in the collector-distributor center to ensure the needs of disaster areas. The fixed cost of a vehicle is 300 yuan and the driving cost is 5 yuan per hour. The maximum load of a vehicle is 150 units and the average speed is 35 km/h. The penalty cost is 200 yuan per hour when the time window requirement of the disaster point is not met. The two damaged road segments need 12 units and 18 units to get restored, respectively. If the damaged segment is not repaired and the road is impassable, the distance of the segment is set to be a large enough integerM, which is 1000 in this case. On the basis of the model built in Section 3 and the actual situation, the parameters of the adaptive artificial bee colony algorithm are set as follows. The population size is set as 200, the initial probability rate of adaptive crossover and mutation as 0.01, the maximum number of iterations as 100, the maximum number of searches as 100, and the program runs 20 times to get the optimal results.


## 5.3. Result Analysis

### 5.3.1. No Repairs of the Road Network

First of all, the situation that all damaged segments are not repaired should be considered, that is, all relief supplies are distributed to the disaster-stricken areas and no supplies are allocated to the damaged segment. It is assumed that Segment 1 between Supply-Needed Location 6 and Location 11 and Segment 2 between Distribution Center 2 and Supply-Needed Location 21 in the transportation network are blocked (shown as the dotted line in Figure6), and a penalty function is set as 1000. In the case of road congestion, the optimal path of the emergency scheduling problem is [15-18-16-21-8-11-6-12-13-7-3-10-5-9-22-4-14-19-17]. A total of five vehicles are needed to complete the whole transportation plan, as shown in Table 3.Figure 6
Road network structure.Table 3
Subpath transport route without road repair.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
1-15-18-16-21-1
54.29
0.97
78.92
2
1-8-11-6-1
541.25
27.73
1029.707
3
2-12-13-7-2
40.84
0.33
40.84
4
2-3-10-5-9-22-20-2
77.60
0.40
91.33
5
1-4-14-19-17-1
65.17
0.43
83.32


### 5.3.2. Repairs of the Partial Road Network

Under the situation that Segment 1 is repaired urgently and Segment 2 is not repaired, the optimal path is [6-11-13-9-4-19-20-21-12-22-7-3-15-18-10-16-8-5-17-14]. Five vehicles are needed to complete the whole transportation plan, as shown in Table4.Table 4
Subpath transport route on the repairing number 1.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
1-6-11-13-9-1
42.38
0
67.67
2
1-4-19-20-21-1
45.33
0.48
80.14
3
2-12-22-3-2
35.27
0.36
42.90
4
1-15-18-10-16-1
46.24
0.81
52.93
5
1-8-5-17-14-1
61.39
0
322.44Under the situation that Segment 2 is repaired urgently and Segment 1 is not repaired, the optimal path is [7-13-12-5-19-10-20-16-6-11-8-4-3-14-9-15-18-22-17-21]. Five vehicles are also needed to complete the whole transportation plan, as shown in Table5.Table 5
Subpath transport route on the repairing number 2.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
2-3-22-10-16-14-2
41.03
0.37
70.53
2
1-4-18-11-6-1
551.76
28.78
1038.31
3
2-21-13-17-9-2
32.76
1.12
52.43
4
2-12-20-5-7-2
62.26
0
69.68
5
1-15-19-8-1
50.90
0.40
54.396


### 5.3.3. Repairs of the Whole Road Network

Under the situation that the whole road network is repaired urgently, the optimal path for emergency supply scheduling is [21-10-19-5-14-6-11-8-16-12-13-22-4-15-7-17-9-3-18-20]. Five vehicles are needed to complete the whole transportation plan, as shown in Table6.Table 6
Subpath transport route with all roads repaired.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
2-21-10-19-5-14-2
57.51
0.67
68.63
2
1-6-11-8-16-1
36.39
0.26
41.73
3
2-12-13-22-2
29.25
0.40
33.42
4
1-4-15-7-17-9-1
73.51
0.83
73.51
5
2-3-18-20-2
30
0.21
56.62Compared the value of objective functionZ1 and the transportation cost Z2 in different situations, the results are shown in Table 7.Table 7
Comparison of target values under four road networks.
Objective function
Repair road number 1
Repair road number 2
Repair all
Not repaired all
Z
1
1243.03
2069.97
266.76
1800.9
Z
2
3112.18
7926.72
2869.53
8115.63In Table7, the ascending order of the objective function values Z1 is as follows: Z1 in the situation that all road networks are repaired, Z1 in the situation that Segment 1 is repaired, Z1 in the situation that no road segment is repaired, and Z1 in the situation that Segment 2 is repaired. The data obtained in the situation of repairing all road networks are obviously smaller than those of repairing no road network or only repairing partial road networks. When all road networks are restored, the target value Z1 is smaller than that under other situations. This is because the fitness function has adjusted the target value and enlarged it when it does not meet the threshold value, which is conducive to enhance the difference among algorithm selections.Similarly, in Table7, the ascending order of the transportation cost Z2 in different road networks is as follows: Z2 in the situation that all road networks are repaired, Z2 in the situation that Segment 1 is repaired, Z2 in the situation that Segment 2 is repaired, and Z2 in the situation that no road segment is repaired. The transportation cost of the road networks that undergoes complete repair is lower than that of the road networks that are not repaired. It is because that the two damaged segments are both the shortest key route between distribution centers to the disaster-stricken areas. Once damaged, transportation costs will rise dramatically and transportation efficiency will decline. The benefit of repairing the roads is greater than the cost, so all the results obtained in the situation of repairing all road networks are optimal. If the relief supplies are only enough to restore part of road segments, it is advisable to repair the road segments which can produce greater benefits. That means the supplies can be sent to more disaster-stricken locations in the shortest time through the chosen road segments. Therefore, it is clear that in the case of road-varying network, it is of great practical value and significance to analyze and compare the whole road network and explore the path cost and benefit under different road networks.


## 5.3.1. No Repairs of the Road Network

First of all, the situation that all damaged segments are not repaired should be considered, that is, all relief supplies are distributed to the disaster-stricken areas and no supplies are allocated to the damaged segment. It is assumed that Segment 1 between Supply-Needed Location 6 and Location 11 and Segment 2 between Distribution Center 2 and Supply-Needed Location 21 in the transportation network are blocked (shown as the dotted line in Figure6), and a penalty function is set as 1000. In the case of road congestion, the optimal path of the emergency scheduling problem is [15-18-16-21-8-11-6-12-13-7-3-10-5-9-22-4-14-19-17]. A total of five vehicles are needed to complete the whole transportation plan, as shown in Table 3.Figure 6
Road network structure.Table 3
Subpath transport route without road repair.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
1-15-18-16-21-1
54.29
0.97
78.92
2
1-8-11-6-1
541.25
27.73
1029.707
3
2-12-13-7-2
40.84
0.33
40.84
4
2-3-10-5-9-22-20-2
77.60
0.40
91.33
5
1-4-14-19-17-1
65.17
0.43
83.32


## 5.3.2. Repairs of the Partial Road Network

Under the situation that Segment 1 is repaired urgently and Segment 2 is not repaired, the optimal path is [6-11-13-9-4-19-20-21-12-22-7-3-15-18-10-16-8-5-17-14]. Five vehicles are needed to complete the whole transportation plan, as shown in Table4.Table 4
Subpath transport route on the repairing number 1.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
1-6-11-13-9-1
42.38
0
67.67
2
1-4-19-20-21-1
45.33
0.48
80.14
3
2-12-22-3-2
35.27
0.36
42.90
4
1-15-18-10-16-1
46.24
0.81
52.93
5
1-8-5-17-14-1
61.39
0
322.44Under the situation that Segment 2 is repaired urgently and Segment 1 is not repaired, the optimal path is [7-13-12-5-19-10-20-16-6-11-8-4-3-14-9-15-18-22-17-21]. Five vehicles are also needed to complete the whole transportation plan, as shown in Table5.Table 5
Subpath transport route on the repairing number 2.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
2-3-22-10-16-14-2
41.03
0.37
70.53
2
1-4-18-11-6-1
551.76
28.78
1038.31
3
2-21-13-17-9-2
32.76
1.12
52.43
4
2-12-20-5-7-2
62.26
0
69.68
5
1-15-19-8-1
50.90
0.40
54.396


## 5.3.3. Repairs of the Whole Road Network

Under the situation that the whole road network is repaired urgently, the optimal path for emergency supply scheduling is [21-10-19-5-14-6-11-8-16-12-13-22-4-15-7-17-9-3-18-20]. Five vehicles are needed to complete the whole transportation plan, as shown in Table6.Table 6
Subpath transport route with all roads repaired.
Route number
Driving route
Efficiency cost
Delay cost
Transportation distance
1
2-21-10-19-5-14-2
57.51
0.67
68.63
2
1-6-11-8-16-1
36.39
0.26
41.73
3
2-12-13-22-2
29.25
0.40
33.42
4
1-4-15-7-17-9-1
73.51
0.83
73.51
5
2-3-18-20-2
30
0.21
56.62Compared the value of objective functionZ1 and the transportation cost Z2 in different situations, the results are shown in Table 7.Table 7
Comparison of target values under four road networks.
Objective function
Repair road number 1
Repair road number 2
Repair all
Not repaired all
Z
1
1243.03
2069.97
266.76
1800.9
Z
2
3112.18
7926.72
2869.53
8115.63In Table7, the ascending order of the objective function values Z1 is as follows: Z1 in the situation that all road networks are repaired, Z1 in the situation that Segment 1 is repaired, Z1 in the situation that no road segment is repaired, and Z1 in the situation that Segment 2 is repaired. The data obtained in the situation of repairing all road networks are obviously smaller than those of repairing no road network or only repairing partial road networks. When all road networks are restored, the target value Z1 is smaller than that under other situations. This is because the fitness function has adjusted the target value and enlarged it when it does not meet the threshold value, which is conducive to enhance the difference among algorithm selections.Similarly, in Table7, the ascending order of the transportation cost Z2 in different road networks is as follows: Z2 in the situation that all road networks are repaired, Z2 in the situation that Segment 1 is repaired, Z2 in the situation that Segment 2 is repaired, and Z2 in the situation that no road segment is repaired. The transportation cost of the road networks that undergoes complete repair is lower than that of the road networks that are not repaired. It is because that the two damaged segments are both the shortest key route between distribution centers to the disaster-stricken areas. Once damaged, transportation costs will rise dramatically and transportation efficiency will decline. The benefit of repairing the roads is greater than the cost, so all the results obtained in the situation of repairing all road networks are optimal. If the relief supplies are only enough to restore part of road segments, it is advisable to repair the road segments which can produce greater benefits. That means the supplies can be sent to more disaster-stricken locations in the shortest time through the chosen road segments. Therefore, it is clear that in the case of road-varying network, it is of great practical value and significance to analyze and compare the whole road network and explore the path cost and benefit under different road networks.


## 5.4. Performance Comparison of Different Algorithms

Few similar examples and algorithms are available at present to handle the problem of scheduling emergency supplies under conditions of multiple distribution centers and road-varying network, so the algorithms cannot be compared with others directly. In order to analyze the performance of this algorithm, in the path optimization stage, the genetic algorithm and the standard ABC algorithm are used to compare with the improved ABC algorithm under the same condition that the whole road network is not restored. In consideration of the difference between genetic algorithm and the standard ABC algorithm, different population numbers after 500 times of iteration are compared, and the results of comparative tests are shown in Table8 and Figure 7.Table 8
Performance comparison between improved ABC algorithm and others.
Population
Improved ABC algorithm
Genetic algorithm
ABC algorithm
100
796
763
814
200
1496
1618
1589
500
3569
4015
3863Figure 7
Comparison of different optimization algorithms.Through comparison, it is obvious that the genetic algorithm has a slight advantage over the improved ABC algorithm when the population size is small, but when the population size increases gradually, the algorithm proposed in this paper significantly overcomes the genetic algorithm. Comparing the improved ABC algorithm with the standard one, we find that the former has advantages in terms of any population sizes. In addition, in Figure7, comparing the convergence speeds and optimal target values using these three algorithms, we can find that the improved ABC algorithm has faster convergence speed and better target value, which proves that the strategy of adaption in this paper is prominently effective.


## 6. Conclusions

Aiming at the problem of reliable emergency supply scheduling in the rescue operation of sudden-onset natural disaster rescue, this paper combines the actual situation of multiple distribution centers and multiple supply-needed locations and designs a bi-level model. In the upper model, we take into account supply distribution and the urgency of time; while in the lower one, we take into account route programming and the fairness of distribution delay. Based on the conceptual requirements of network variability, path optimization, and time window, a multiobjective scheduling model is constructed, and a specific problem-solving strategy is proposed. Finally, the following conclusions are drawn. ① The scheduling schemes under different road network states can be obtained by virtue of the bi-level programming model, which gives better consideration to both supply distribution and route programming. ② When relief supplies are limited, different distribution schemes should be worked out to produce different benefits. Nevertheless, efficiency and fairness are contrary to each other, which require comprehensive consideration by decision makers. ③ In the early stage after the disaster, relief supplies cannot fully meet the needs of the disaster-stricken areas and damaged segments, so rational programming of supply distribution or the increase of resource input should play its part in overcoming the difficulty. ④ Different strategies for repairing the road network will lead to different revenue results, which are related to the overall structure of emergency supply network. Especially when the damaged segment is exactly the shortest route from the distribution center to the disaster area, the benefit of repairing the segment will be greater than the cost of supplies, which makes it more reasonable to repair the damaged segment. ⑤ In view of the model’s characteristics of multiple objectives and multiple parameters, the improved artificial bee colony algorithm is adopted to solve the problem. Finally, an actual case is given to verify the feasibility and validity of the proposed model and algorithm, which provides a reference for scheduling emergency supplies based on bi-level programming under conditions of multiple distribution centers and variable road networks.However, this paper only considers the supply scheduling problem when the supply demand is fixed. There are some other problems involved that need to be further explored as follows:(1)
In actual sudden-onset disasters, a lot of information is uncertain. In the initial stage after the disaster, the demand for relief supplies is ambiguous. In that case, some assumptions need to be more lenient, new models need to be constructed under other related constraints, and new effective parameter values need to be determined. Moreover, under the circumstance of uncertain information, fairness as well as effective route programming plan should be discussed and explored.(2)
Study the problem of the vulnerability of road network nodes [24, 25]. To solve the problem, complex network theory should be introduced to explore the model for measuring the vulnerability of road network and disaster site in emergency logistics, especially the risk impact on cities of different importance. By reducing the number of the damaged road network or increasing that of disaster points, the impact of attack strategies on the emergency logistics network after different disasters should be analyzed, which is of great significance for further prevention of sudden-onset natural disasters.(3)
Study the dynamic evolution characteristics of the emergency logistics network. Because of the derivative and secondary of nature disasters, the number of disaster-stricken areas is increasing, demand will also change dynamically after a sudden-onset disaster, and the corresponding scheduling network is also dynamic in real time. Therefore, establishing the corresponding network evolution model and revealing the change rules of emergency logistics network will lay a foundation for the later reliable scheduling.



---

*Source: 7501891-2020-01-20.xml*