Improved whale swarm algorithm for solving material emergency dispatching problem with changing road conditions

1.   Introduction

Over the past few years, the frequent occurrences of natural disasters in the world have caused serious damage to social construction and economic development, such as the flood disaster in India in 2013, the hail disaster in Yancheng in 2016, the 2019 hill fires in Australia and Hurricane Laura in the United States in 2020 ^[1,2,3]. In China, 69,227 people died in the "5.12" Wenchuan earthquake in 2008, with direct economic losses of $119.7 billion ^[4]. The "7.20" extremely heavy rain disaster in Zhengzhou in 2021 affected 14.786 million people, with direct economic losses reaching $17 billion ^[5]. Since most sudden natural disasters cannot be predicted in advance, strengthening the guarantee of emergency supplies after disasters has become the key to improving rescue efficiency. Therefore, in order to ensure the smooth development of rescue work, it is necessary to formulate a scientific and reasonable material emergency dispatching (MED) scheme, which can minimize the subsequent losses from the disaster.

With the in-depth study of MED, the combinations of MED and actual scenes have become the main research direction. Therefore, an array of variant problems have been derived, such as the problem with time windows ^[6,7], the problem with multiple distribution centers ^[8,9], the multi-objective optimization problem ^[10,11] and the dynamic routing problem (DRP) ^[12]. Among them, the DRP has attracted high attention from researchers ^{[13,14,15,16]}. Psaraftis et al. ^[14] analyzed the impact of the new generation of information technology on DRP and summarized the research literature according to 11 classification standards. Reference [16] classified the DRP models into dynamic demand and dynamic road conditions and introduced the research progress of its route update strategies and solving algorithms. DRP is proposed relative to the static routing problem (SRP). In the SRP, relevant information such as distribution points, demand and road conditions are known and fixed before planning the distribution path, and the vehicle route will not be changed once it is determined. However, some of the information in the DRP may be changed, so its dynamics bring a great challenge in optimizing the material scheduling path. More specifically, due to the continuous impact of disasters and secondary disasters, the emergency situations that road blockage or destruction make impassable vehicles occur frequently in the disaster scenario, which makes it necessary to dynamically adjust the planned path according to real-time information during material scheduling. Therefore, it is more realistic to consider the "changing road conditions" in the study of emergency material dispatching. In addition, the essence of MED is a vehicle routing problem (VRP). As a classical combinatorial optimization problem, the VRP has been proved to be an NP-hard problem. Therefore, while paying attention to the construction strategy of the MED model, scholars have also conducted a lot of research on its solution algorithm. At present, there are mainly the exact algorithm ^[17,18,19] and the heuristic algorithm ^[20,21,22]. Among them, the exact algorithm establishes the corresponding mathematical model for a specific problem and uses mathematical methods to find the optimal solution. Although the exact algorithm can achieve better results in solving the problem, its computational time explosively increases with an increase of the problem size, which makes it only suitable for solving small-scale problems ^[23]. The heuristic algorithm can give a feasible solution for the combinatorial optimization problem within an acceptable range, which has better robustness and feasibility in dealing with large-scale problems compared to the exact algorithm. However, its solution performance is mostly limited by the hand-designed algorithm rules. Especially in the cases of single distribution centers, no time constraints and changes in road conditions, the stability of the heuristic algorithm is deficient ^[24,25].

Although the research on "changing road conditions" has been widely carried out, finding high-performance algorithms to improve the execution efficiency, feasibility and reliability of the scheduling scheme is still an urgent practical demand problem in post-disaster material distribution. Therefore, considering the disadvantages of the traditional heuristic algorithm and the complexity of the material emergency dispatching with dynamic road conditions (MEDDRC) scheme, this paper proposes an improved whale swarm algorithm (IWSA) to solve the MEDDRC problem. The main contributions of this study are as follows:

1) The improved scan Clarke-Wright (ISCW) method is designed to plan the initial driving route of vehicles, which can obtain the initial scheduling scheme with minimum cost.

2) The group movement rule is proposed to update the individuals, which not only ensures that individuals in the population carry out a detailed local search of in different regions, but also enhances the global search capability of individuals.

3) To avoid individuals prematurely congregating in a certain area during the iterative process, a different weights strategy is constructed to expand the individual search space, which is helpful to maintain the population diversity.

4) The feasibility of the proposed algorithm in solving material emergency dispatch under dynamic changes in road conditions is verified through extensive experiments.

The remaining sections of this paper are organized as follows. Section 2 provides a review of the relevant literature on MED in recent years. Section 3 develops a mathematical model for the MEDDRC problem. A detailed description of the proposed IWSA and the main strategies are shown in Section 4. Section 5 analyzes the experimental results of the algorithm. Our conclusions and future work are given in Section 6.

2.   Literature review

As a key part of emergency management, the merits and demerits of the MED scheme will directly affect the smooth development of subsequent rescue work. For this reason, scholars have carried out a lot of research work around the MED. According to the different focuses of the research, the existing works mainly include two aspects of the model construction and solution.

As the specific situations of disasters are different, the objectives and influencing factors considered in material dispatching are also somewhat different. Therefore, in order to meet the special needs in different disaster situations, researchers have constructed corresponding mathematical models for specific application scenarios. For example, in solving the single-objective problem, Sung and Lee ^[26] achieved maximizing the number of rescuers with limited resources by modeling patient priority as an ambulance scheduling problem. Balcik et al. ^[27] proposed a mixed integer programming model for studying the last mile distribution problem, which can reduce the transportation cost to maximize the benefit of the recipients. However, in the actual disaster relief operations, there is a lack of comprehensive consideration of multiple objectives. For this reason, Vitoriano et al. ^[28] constructed a multi-criteria optimization model for the MED problem by considering objective criteria such as time, cost, security and equity. Zhang et al. ^[29] constructed a multi-objective optimization model for multiple reserve points and multiple relief supplies, which can solve the difficult problem of the allocation-scheduling of emergency relief supplies. In addition, Chen et al. ^[30] presented a multi-objective optimization model with a utility–priority–economy allocation strategy, which can effectively solve the MED problem. However, due to the complexity of the situation in the disaster area and the uneven distribution of the locations of the affected points, it is difficult to meet the urgent time requirements of MED for a single rescue point. For this reason, Wex et al. ^[31] constructed a rescue unit allocation and scheduling model for multiple rescue points with the objective of minimizing the total completion time and obtained the optimal solution by using the heuristics algorithm. Meanwhile, Song et al. ^[32] introduced the emergency point satisfaction and materials lack loss coefficient and built an emergency material dispatching model with a multi-emergency point, a multi-rescue point and multi-stage. However, most of the above models only consider the influence of static factors, without considering the time-variant of road conditions and traffic environment, which will lead to the vehicle not being able to smoothly complete the material distribution tasks according to the existing scheduling scheme. Therefore, Qu et al. ^[33] took the timeliness and fairness of emergency material distribution as the goal and established a multi-objective optimization model for distribution path planning under multi-period dynamic conditions. In addition, the update strategies of the scheduling scheme such as a fixed time interval ^[34] and a key point ^[35] are often adopted. For example, Okhrin et al. ^[36] adopted the method of updating the status of the current road section at regular intervals to solve the problem of real-time vehicle routing; however, it may cause an inconsistency between the time of information update and the time of actual route change. Therefore, Li et al. ^[37] proposed a new method of updating the route at key points to solve the problem of dynamic network vehicle distribution, which was more effective compared to the traditional updating mechanism.

In problem-solving, the quality of the scheduling scheme depends on the algorithm for solving the problem; a high-performance solving algorithm can significantly increase the possibility of obtaining an efficient and reasonable scheduling scheme. For this reason, researchers have conducted a lot of research on how to improve the performance of problem-solving algorithms. For example, Zhang et al. ^[38] used a local search technique and high-priority rule to design an algorithm based on linear programming to solve the optimal allocation of emergency resources, which can reduce the problem complexity by improving the linear programming solution after relaxation. Meanwhile, Queiroga et al. ^[39] proposed two branch cut and price (BCP) algorithms based on the backhaul vehicle scheduling problem, which incorporated advanced elements such as strong branching and route enumeration to reduce the computation time of the problem. Although the above work has achieved good results, the used algorithms have a high complexity in solving large-scale problems; therefore, most heuristic algorithms are used to solve such problems. For example, Xiang et al. ^[40] proposed the ant colony system (ACS) algorithm based on pairwise proximity learning to solve dynamic VRPs, which can predict the local visiting order of customers after the occurrence of changes according to the optimal routes found before the changes occur. Wang et al. ^[41] designed a multi-objective cellular genetic algorithm (GA) to optimize the post-disaster emergency resource allocation model, which achieved a balance between global and local optimization by introducing the auxiliary population and neighborhood structure in the cellular automata. In addition, Zhang et al. ^[42] proposed an improved particle swarm optimization (PSO) algorithm to solve the optimal scheduling problem with multiple objectives, which can increase the population diversity and expand the search range by adding perturbation operations to the update process of particle velocity. Although the heuristic algorithm has a high feasibility, it suffers from the disadvantage of falling into a local optimum prematurely, which makes it difficult to obtain an approximate optimal solution in the global scope. To this end, scholars adopted the whale swarm algorithm (WSA) ^[43] with an improved performance to solve the problem, and amended it on the basis of its superior iteration rules to enhance the global solution ability of the algorithm. For example, Dong and Ye ^[44] designed a new individual position movement strategy and a multi-neighborhood search strategy for WSA to solve the flow shop scheduling problem, which can balance the global and local search abilities of the algorithm to improve the solution performance. However, with the update and iteration of the population, the population diversity will rapidly decrease. For this reason, Zhang et al. ^[45] proposed a discrete WSA to solve the flow shop scheduling problem, and improved the individual movement rules of the WSA by constructing five mutation operators, which can maintain the diversity of the population and enhance the solving ability of the algorithm.

In addition, the relevant literatures on model construction and solution algorithms for MED problems are listed in Table 1. Among them, the superior iteration rules of WSA make its performance better than the traditional heuristic algorithm. However, due to the discrete characteristics of the MEDDRC problem, the algorithm rule in the traditional WSA cannot well solve this problem. Moreover, during the iteration of the algorithm, the convergence of individual whale movements in the WSA makes the population individuals still gather quickly in a certain region of the solution space, which in turn limits the ability of individual whales to search for solutions in the global range. Therefore, in this paper, the algorithm rules of WSA will be improved in terms of the initial route generation strategy, the individual movement method and the search space expansion, which can adapt to the MEDDRC problem and enhance the effectiveness of the algorithm.

Table 1.  Relevant literatures on the construction and solution for MED model.

Research Objectives	Solving Methods	Application Scenarios	Literature
Maximizing the number of people rescued	Column generation approach	Medical resource allocation	[26]
Minimizing transportation costs; Maximizing benefits to aid recipients	GAMS/Cplex	Humanitarian relief	[27]
Maximizing transportation safety and distribution equity; Minimizing transportation time and costs	GAMS/Cplex	Humanitarian relief	[28]
Minimizing distribution time; Maximizing material satisfaction and reliability	ACS、GA	Emergency relief supplies distribution and dispatch	[29]
Maximizing total allocation utility; Minimizing distribution costs	Goal programming、Fuzzy analytic hierarchy process	Multi-depot and multiple affected area emergency logistics system	[30]
Minimizing the sum of completion times	Heuristic algorithm	Natural disaster relief	[31]
Minimizing total dispatching cost; Maximizing emergency point satisfaction	PSO	Nonlinear continuous consumption MED problem	[32]
Minimizing transportation time; Maximizing distribution equity	GA	Post-earthquake MED problem	[33]
Minimizing travel time	GA	Dynamic VRP	[35]; [36]; [37]
Minimizing the rescue costs	Linear programming	Multiple-resource and multiple-depot MED problem	[38]
Minimizing distribution costs	BCP	VRP with backhauls	[39]
Minimizing total driving costs	ACS	Dynamic VRP	[40]
Minimize the disaster losses and the total risks of the selected routes in the rescue process	GA	MED for multi-period post-disaster	[41]
Minimizing distribution distance; Maximizing distribution equity	PSO	Emergency material allocation and vehicle routing plan generation	[42]
Minimizing completion costs	WSA	Flow shop scheduling problem	[44]; [45]

---

 | Show Table

DownLoad: CSV

3.   Problem description and model building

3.1.   Problem description

For the MEDDRC problem studied in this paper, due to the existence of a dynamic factor such as the change of road conditions, it is difficult for the fixed vehicle scheduling scheme to complete the distribution task of all disaster-affected points in the shortest time. Therefore, in order to describe the dynamic characteristics of MEDDRC and to consider the complexity of information received by the vehicle while driving, the concept of a "key point" ^[46] is introduced and divided into two levels of key points. The first-level key point is the location of the disaster-affected point, and the second-level key point is the intersection of two roads (that is, the intersections of paths to different affected points are considered the secondary key points). The details are shown in Figure 1.

Figure 1.  Example of vehicle driving route.

DownLoad: Full-Size Img PowerPoint

When a vehicle arrives at the first-level key point, the road condition information of the disaster area is updated, and the relevant information is used to judge whether the vehicle scheduling scheme needs to be redesigned to visit the disaster-affected points that has not been visited currently. Due to the different locations of the vehicles sent from the distribution center at the same time, it is necessary to determine the position of each vehicle when reformulating the vehicle scheduling plan. Assuming that at time $t$, vehicles ${k_1}$ and ${k_2}$ are located on the first-level key point and the path from disaster-affected point $i$ to $j$, respectively, and there are n second-level key points between $i$ and $j$, the departure time of the vehicle at each key point are ${t_i}$, ${t_j}$, ${t_p}$($p$ is the $p$-th second-level key point, $1 \leqslant p \leqslant n$), respectively. If ${t_i} < t \leqslant {t_p}$, the key point $p$ is taken as the position of vehicle $k$, and if ${t_p} < t \leqslant {t_j}$, the key point $j$ is taken as the position of vehicle $k$. Figure 1 shows an example of two vehicles performing a distribution task. At time $t$, vehicle 1 arrives at key point 2, and vehicle 2 is located between key points ${p_3}$ and ${p_4}$. When updating the vehicle scheduling scheme, the starting position of vehicle 1 is the first-level key point 2, and the starting position of vehicle 2 is the second-level key point ${p_4}$.

3.2.   Mathematical model

There are many complex factors in the disaster area environment, which have brought many influences to the solution of the MEDDRC problem. Therefore, in order to reduce the complexity of the problem, the following constraints are given in this paper: 1) only one vehicle completes the material distribution task at each disaster-affected point; 2) the total demand of each affected point visited by each vehicle cannot exceed the maximum load of the vehicle; 3) the vehicles are fully loaded when they depart from the distribution center; and 4) the vehicle must leave after completing the material distribution task at the affected point. In addition, the objective function of the MEDDRC problem is to minimize the time cost spent by the vehicle, so the mathematical model is constructed as follows:

and subject to

The notations in the above formula are defined as shown in Table 2.

Table 2.  Notations for various parameters of the mathematical model.

Sets	Definition
${S_I}$	Set of key points
$R$	Set of disaster-affected points
$K$	Set of vehicles
Parameters
0	Number of distribution center
$C$	Maximum load of vehicle
$k$	The $k$-th vehicle
$i、j、h$	Node variable
${q_i}$	Material demand of node $i$
${c_d}$	The cost per unit time of vehicle travel
${d_{ij}}$	Distance between nodes $i$ and $j$
${v_{ij}}$	Running speed of the vehicle between nodes $i$ and $j$
${t_{ijk}}$	Load of vehicle $k$ when it reaches node $j$ from node $i$
Variables
${x_{ijk}}$	Path decision variable, ${x_{ijk}}$=1 when vehicle $k$ goes from $i$ to $j$, and ${x_{ijk}}$=0 otherwise
${y_{ik}}$	Vehicle decision variable, ${y_{ik}}$ = 1 when disaster-affected point $i$ is delivered by vehicle $k$, and ${y_{ik}}$ = 0 otherwise

---

 | Show Table

DownLoad: CSV

4.   Model solving

4.1.   Overview of the IWSA

In this paper, an IWSA is proposed to solve the MEDDRC problem. The framework of the IWSA is shown in Figure 2. Firstly, the ISCW algorithm is used to plan the initial driving route of the vehicle for obtaining a high-quality initial solution. Then, the group movement rule is used to update the individuals to ensure the detailed local search of population individuals in different regions. In addition, the individual search space is extended by constructing the different weights strategy to maintain the population diversity. The details of the algorithm are shown in Sections 4.2 and 4.3.

Figure 2.  Flow chart of IWSA for solving MEDDRC problem.

DownLoad: Full-Size Img PowerPoint

4.2.   Generation of the initial path

At the initial time, it is necessary to formulate the optimal path for the vehicle based on the current known location of the disaster-affected point, demand and road information. Assuming that the vehicle is running at a normal speed, the Clarke-Wright (C-W) algorithm ^[47] can generally be used to generate the initial path of the vehicle. Although the C-W algorithm has good results when used alone, its solution speed will decrease as the VRP scale increases. Therefore, in order to improve the performance of the algorithm to obtain a better initial path, this paper proposes an ISCW algorithm, which takes the angles between the lines connecting the disaster-affected points and the origin as the division basis, scans all disaster-affected points counterclockwise, and uses the improved C-W algorithm to plan the access path of the disaster-affected point. The basic idea is as follows:

1) Calculate the angle of the disaster-affected point in polar coordinates. Taking the distribution center as the origin, the coordinates of the disaster-affected point are converted from a rectangular coordinate system to a polar coordinate system, and rays are made from the origin to all the disaster-affected points, respectively. The line between the origin and a randomly selected disaster-affected point is taken as the polar axis. The remaining disaster-affected points are scanned counterclockwise starting from the polar axis, and the angle between the polar axis and the line connecting the disaster-affected point and the origin is calculated.

2) Determine the adjacent disaster-affected points with the largest included angle. Starting from the polar axis, the included angle between the lines formed by connecting the adjacent disaster-affected points to the origin is calculated in turn. The disaster-affected points corresponding to the two rays with the largest angle are selected, and their relative positions are judged (that is, the counterclockwise angle between the polar axis and the line connecting the origin and the disaster-affected point).

3) Determine the starting point of scanning and group the disaster-affected points. According to the two disaster-affected points with the largest angle determined in 2), the disaster-affected point with the largest angle with the polar axis is selected as the scanning starting point of counterclockwise, and the remaining disaster-affected points are grouped according to the maximum load constraint of the vehicle. That is, if adding a new scanned disaster-affected point makes the total demand in the current group more than the maximum load of the vehicle, the disaster-affected point is added to a new group, and the remaining disaster-affected points are continue scanned.

4) Determine the visiting order of the disaster-affected points. According to the grouping of the disaster-affected points, the improved C-W algorithm is used to determine the service order of the disaster-affected points in each group, that is, the visiting order is determined sequentially by the saving value. Specifically, for a group of disaster-affected points that do not violate the vehicle load constraint, when the two selected disaster-affected points with the largest savings value are not related to the planned disaster-affected points, the improved C-W algorithm directly adds them to the set of planned disaster-affected points instead of selecting new vehicles to visit.

4.3.   Update strategy of the path

When the road conditions change, the main works of IWSA proposed in this paper include: 1) the initial population is generated based on the shortest path mutation to adapt to the dynamic characteristics of the problem; 2) when the population is iteratively updated, the guide individual and the specific movement mode of each individual in the population are determined by the group movement strategy; and 3) in the process of population iteration, when the population diversity is lower than a certain threshold, the method of moving individuals based on different weights is used to screen the individuals and generate new individuals.

4.3.1.   Population initialization method based on hybrid operator

Due to the update of road condition information, the start node is not the distribution center when re-planning the unvisited disaster-affected points, and the vehicles are not necessarily fully loaded at this time. Therefore, in order to adapt to the specificity of the problem, this paper uses the hybrid operator based on the shortest path mutation to initialize the population. The shortest path algorithm is used to generate the initial solution of the subproblem, and the initial solution is expanded by using the 2-opt operator ^[48] and the split mutation operator ^[49] to obtain the initial population. The specific method is shown in Figure 3.

Figure 3.  Principle of 2-opt operator and split mutation operator.

DownLoad: Full-Size Img PowerPoint

In Figure 3, the 2-opt operator flips the order of nodes between two randomly selected positions $l1$ and $l2$ to generate new individuals. The split mutation operator randomly selects a split point $l$, keeps the node order before and after the split point $l$ unchanged, and uses the improved C-W algorithm to replan the remaining nodes for generating two new individuals.

4.3.2.   Individual update method based on group movement strategy

The individual movement method in WSA is mainly suitable for solving continuous problems; however, MEDDRC is a discrete problem. Therefore, in order to adapt to its discreteness and improve the quality of the problems solution, this paper constructs a group movement strategy for updating individuals. The main ideas are as follows:

1) Grouping of individuals in the population. Hierarchical clustering is used to divide the individuals in the current population equally according to their differences, and to minimize the total individual differences within the group. The clustering objective is as in Eq (7), and the individual with the best quality in each group is recorded as ${G_b}$. Eq (8) is used to calculate the difference between individuals $i$ and $j$, and Eq (9) is used to measure the quality of individuals in the population $P$.

In Eq (7), $c$ is the class number, $x$ is the sample vector, ${\mu _i}$ is the mean vector of the $i$-th class samples, and ${S_i}$ is the set of the $i$-th class samples. In Eq (8), ${r_{ij}}$ represents the same path length traveled by vehicles in individual $i$ and $j$, n is the total number of disaster-affected points, and ${m_i}$ and ${m_j}$ are the number of vehicles used in individual $i$ and $j$, respectively. In Eq (9), $dis{t_i}$ is the total distance traveled by the vehicles to finish the material distribution task for all disaster-affected points represented by individual $i$, and ${f_q}\left(i \right)$ represents the quality of individual $i$.

2) Update the optimal individual in the long term memory unit (LTMU). If the current population is the initial population, the highest quality individual ${W_b}$ and the higher quality individual with the largest difference from ${W_b}$ are selected from it to join the LTMU; otherwise, the highest quality individual ${p_b}$ in the current population is found and compared with ${W_b}$. If individual ${p_b}$ is better than ${W_b}$, the ${W_b}$ is replaced by ${p_b}$.

3) Update other individuals in the LTMU. The high-quality individual with the largest difference from ${W_b}$ is selected from the current population to be added to LTMU, and the remaining individuals except in the ${W_b}$ original LTMU are deleted.

4) Individual movement. Within the group, all individuals in the group except the optimal individual ${G_b}$ are moved with ${G_b}$ as the guide. Between groups, the individual from the LTMU, which is closest to the optimal individual ${G_b}$ of each group and better quality, is selected as its guide individual and moved. For the optimal individual ${p_b}$ in the population, if it is the same as the optimal individual ${W_b}$ in LTMU, the 2-opt operator is used to generate its offspring individuals.

In order to better illustrate the idea of the group movement strategy, Figure 4 shows an algorithm diagram with 12 disaster-affected points and LTMU size of 4.

Figure 4.  Diagram of the group movement strategy.

DownLoad: Full-Size Img PowerPoint

According to their differences and similarities, when the individual moves to its guide individual, this paper adopts the following movement rules. First, the same path segment in individual $i$ and its guide individual $j$ is found, such as $a \to b \to c$, and the path is regarded as a whole node $A$, whose initial visit position is node $a$, and its end visit position is node $c$. Then, the improved C-W algorithm is used to plan the path between the node $A$ and other nodes to obtain the offspring individual $k$. Figure 5 shows a specific example.

Figure 5.  Diagram of individual movement.

DownLoad: Full-Size Img PowerPoint

In Figure 5, the same path segments of individual $i$ and guide individual $j$ are 1→2 and 6→7, so they are regarded as two whole nodes, and the path planning is carried out together with the remaining nodes, so as to obtain the new offspring individual $k$.

4.3.3.   Expand the search space based on different weights

In the MEDDRC problem, the solution space of the problem grows explosively with the increase of the number of nodes. The movement of individual whale is a local small-range search based on the original individual and guided by the guide individual. Therefore, by iteratively updating of the population, the individuals gradually gather in a certain area, which makes it difficult for the superposition of the movable space of all individuals in the population to cover the entire solution space of the problem when the individual move toward the guide individual. If the individuals in the population gather in a certain area of the solution space prematurely, the population diversity will be lost, which may cause the algorithm to fall into a local optimum, making it difficult to find the global optimal solution of the problem. Therefore, in order to maintain a high population diversity and to ensure that individuals can be distributed in a large area of the solution space, this paper uses the strategy of different weight to expand the individual search space. When the population diversity is below a certain threshold, half of the individuals with a poor quality and a smaller total difference value in the current population are discarded. Then the search space of the remaining individuals in the population is expanded using different weight strategies to obtain new offspring individuals.

Population diversity can be used to represent the differences among individuals in a population. When the individuals in the population gather in a certain area, the differences among individuals are small. The population diversity is currently low, which makes the total search space of the individual during moving is also smaller. On the contrary, when the individuals in the population are distributed in a large area of the solution space, the differences among individuals are larger. Currently, the population diversity is higher, which makes the total search space of the individual during moving is also larger. Therefore, the population diversity can be used to reflect the size of the total search space of the individuals in the population during moving, and the population diversity can be represented by inter-individual differences. In this paper, Eq (10) is constructed based on Eq (9) to measure population diversity.

In Eq (10), ${d_{ij}}$ is the difference between individual $i$ and $j$ calculated by Eq (9), where $s$ is the population size.

Everyone in the population can reflect the number of vehicles required for material distribution to all the disaster-affected points and the order of vehicle visits to the disaster-affected sites. In addition, it can be known from Eq (8) that the individual quality is inversely proportional to the total distance traveled by the vehicle represented by the individual, that is, the longer the total distance traveled by the vehicle represented by the individual, the worse its quality. In the same individual, each vehicle departs from the distribution center and returns after completing the distribution task. Due to the different locations of the disaster-affected points, the distance traveled by each vehicle is different, so the distance traveled by each vehicle in the individual has a different degree of influence on its quality. In the same vehicle depot (the vehicle depot refers to the complete travel path of a vehicle, including the disaster location and distribution sequence of the vehicle), the relative position of the disaster-affected point influences the total distance traveled by the vehicle. This paper considers the influence of different parts of the individual when expanding the individual search space, and proposes a method for moving individuals based on different weights. The specific idea is as follows:

1) For individuals in the population, Eq (11) is used to calculate the proportion of the running distance $dist$ of each vehicle in the individual in the total running distance $Dist$, thus the weight of each vehicle depot is obtained, and the vehicle depot ${l_c}$ with the largest weight is selected.

In Eq (11), if $dist$ is the driving distance of a certain vehicle depot, then $Dist$ is the total driving distance. If $dist$ is the running distance required by vehicles to visit a certain node, then $Dist$ is the running distance of the vehicle depot where the node is located.

2) Equation (11) is used to calculate the proportion of the total distance $dist$ between the remaining nodes in ${l_c}$ except the distribution center node and its anteroposterior nodes in the total running distance $Dist$ of ${l_c}$, so as to obtain the weight of each node in ${l_c}$. Then, the node $a$ with the largest weight is deleted, and the available load $q$ of the vehicles in the current vehicle depot is calculated.

3) All the nodes except the distribution center in the remaining vehicle depots are traversed in turn, and node $b$, which can minimize the total distance of ${l_c}$ by inserting it at node $a$, is selected. Meanwhile, it is satisfied that the demand of node $b$ does not exceed the available load of vehicles in ${l_c}$, and the demand of node $a$ does not exceed the available load of the vehicle depot where node $b$ is located. The positions of nodes $a$ and $b$ are exchanged to generate new individuals.

5.   Experimental analysis

The suddenness and destructiveness of natural disasters has brought severe challenges to governments around the world. Moreover, for a populous country with a complex geographical environment like China, a timely and effective post-disaster response is particularly important. Therefore, in order to provide decision support for managers to formulate effective disaster relief strategies, this paper constructs a MEDDRC model according to the rescue characteristics of natural disaster events, and designs the ISWA solving algorithm. In the following sections, the effectiveness of IWSA in solving the MEDDRC model will be verified through experiments.

5.1.   Experimental environment and parameter setting

Our experiments were coded in Python 3.8, compiled on PyCharm, and run on an Intel (R) Core (TM) i5-8500T CPU @ 2.10GHz, Windows 10 operating system. In this paper, the experiments are conducted using the randomly generated dataset, and the change of road condition is represented by randomly modifying the distance between the affected points. The number of nodes in the data set is 10, 20, 50 and 100 respectively, and the size of the data set is 50. The node positions are randomly selected in [0,100]², and the node demand is randomly generated in ^[1,9], and the vehicle loads corresponding to different node numbers are 20, 30, 40 and 50, respectively.

The parameters in this paper mainly include the population size, the size of the long-term memory unit, the grouping number of individuals in the population and the threshold of population diversity, in which the relationship among the grouping number $g\_num$, the population size $pop\_num$ and the number of hierarchical clustering $n$ can be shown in Eq (12).

A larger number of groups means that the division of individuals in the population is more detailed, and the search space is relatively small when the individuals in the group move. A smaller number of groups means that the division of individuals in the population is wider, and the search space is relatively larger when the individuals in the group move. In order to achieve better performance of the algorithm, this paper used a pre-experimental method to determine some parameters, that is, combining the existing works and numerous tests to obtain a set of parameters with relatively good results. The pre-experiments were carried out with a node size of 20, a data set size of 50 and a vehicle load of 30. The final parameters are set as follows: the population size is 30, the long-term memory unit size is 6, the number of groups is 2, and the population diversity threshold is 0.4.

5.2.   ISCW validity verification

In order to verify the performance of the ISCW in the model, the ISCW, scan algorithm and C-W algorithm are used to generate the initial vehicle distribution routes of the MEDDRC problem respectively, and the results are shown in Table 3.

Table 3.  Performance comparison of initial paths generated by different algorithms.

Scales	Batches	ISCW			scan algorithm			Average cost deviation	C-W algorithm			Average cost deviation
		Average distance (km)	Average number of vehicles	Average time (s)	Average distance (km)	Average number of vehicles	Average time (s)		Average distance (km)	Average number of vehicles	Average time (s)
10	1	540.810	3.33	3.758 × 10−3	561.839	3.3	2.299 × 10−4	3.743%	540.948	3.77	6.386 × 10−3	0.026%
	2	522.060	3.26	6.899 × 10−4	536.310	3.25	9.995 × 10−5	2.657%	525.022	3.74	2.098 × 10−3	0.564%
	3	534.432	3.09	7.496 × 10−4	545.404	3.09	5.995 × 10−5	2.012%	544.666	3.7	2.059 × 10−3	1.879%
20	1	776.420	3.97	1.669 × 10−3	840.303	4.01	9.993 × 10−5	7.602%	865.543	5.88	8.405 × 10−3	10.297%
	2	780.678	4.02	1.529 × 10−3	829.346	4.04	2.299 × 10−4	5.868%	859.095	5.83	8.454 × 10−3	9.128%
	3	778.583	4.12	1.739 × 10−3	826.095	4.07	1.199 × 10−4	5.751%	855.665	5.83	9.054 × 10−3	9.008%
50	1	1469.121	7.11	5.427 × 10−3	1725.696	7.11	3.899 × 10−4	14.868%	1707.334	12.25	7.152 × 10−2	13.952%
	2	1490.557	7.05	4.837 × 10−3	1743.931	7.1	3.000 × 10−4	14.529%	1750.497	12.46	6.466 × 10−2	14.849%
	3	1508.282	7.05	4.967 × 10−3	1780.993	7.04	3.600 × 10−4	15.312%	1760.988	12.22	6.835 × 10−2	14.350%
100	1	2465.197	11	1.166 × 10−2	3097.789	10.95	6.296 × 10−4	20.421%	3021.683	22.14	3.222 × 10−1	18.416%
	2	2479.745	11.12	1.195 × 10−2	3156.978	11.1	6.696 × 10−4	21.452%	3054.719	22.04	3.132 × 10−1	18.822%
	3	2465.101	10.99	1.151 × 10−2	3082.964	10.97	6.296 × 10−4	20.041%	3035.305	22.39	3.079 × 10−1	18.786%

---

 | Show Table

DownLoad: CSV

As can be seen from Table 3, under different problem scales, the average distance solved by ISCW is smaller than that of both the scan algorithm and the C-W algorithm for the same set of samples. When the number of nodes is 10, 20, 50 and 100, the average distance calculated by ISCW is reduced by 2.804 and 0.823%, 6.407 and 9.478%, 14.903 and 14.384%, 20.638 and 18.675%, respectively, on average compared with scan algorithm and C-W algorithm. The reduction degree of the average distance solved by the ISCW increases greatly with the increase of problem scale, which is attributed to the improvement of the scan algorithm. In the scan algorithm, a disaster-affected point is randomly selected, and other disaster-affected points are scanned in a certain direction according to the angle between the lines joining the disaster-affected points and the distribution center. By default, the loop formed between the disaster-affected points in similar directions is short; however, the influence of the size of the angle between the lines connecting the disaster-affected points and the distribution center is ignored. In the ISCW, one of the two disaster points with the largest included angle is selected as the starting point to scan, which can avoid dividing the two nodes with the largest included angle in the same vehicle and reduce the distance of the solution to a certain extent.

In addition, Table 3 also shows that the average number of vehicles calculated by the ISCW is roughly the same as that of the scan algorithm for the same set of samples, but less than that of the C-W algorithm. This is because in the C-W algorithm, when the two disaster-affected points to be selected are not related to the disaster-affected points in the current planning set, a vehicle will be reassigned to access, which may cause unnecessary use of the vehicle or waste of used vehicle capacity. However, in the ISCW, since each group of disaster-affected points is obtained according to the improved scan algorithm, the total demand of the disaster-affected points must be less than or equal to the maximum load of the vehicle. Therefore, if the C-W algorithm is used to optimize the path for each group of disaster-affected points, the remaining disaster-affected points can be directly added to the planned disaster-affected points set when there are new vehicles to be added, which can avoid the waste of vehicle capacity. Additionally, in terms of solution time, the calculation time of ISCW is between that of scan algorithm and C-W algorithm. The ISCW takes more time to calculate than that of the scan algorithm. This is because the ISCW also uses the improved C-W algorithm to optimize the path of the disaster-affected point except for the scan algorithm, so it takes more time than just using the scan algorithm. Compared with the C-W algorithm, the ISCW takes less time because the C-W algorithm needs to select the appropriate nodes from all the disaster-affected points; in turn, according to the saving value, the set of disaster-affected points is relatively large. The ISCW firstly uses the scanning algorithm to group the disaster-affected points, which reduces the set of disaster-affected points compared with directly using C-W algorithm, and can shorten the time spent in selecting nodes to a certain extent. Therefore, compared with using the C-W algorithm, the solution time of the ISCW is shorter.

5.3.   Comparative analysis of algorithms

To verify the performance of the IWSA in solving the MEDDRC problem, the IWSA, the ACS ^[50] and the adaptive chaos genetic algorithm (ACGA) ^[51] are used to solve the problem in this paper respectively, and the calculation results are shown in Table 4 and Figures 6 and 7.

Table 4.  Calculation results of three algorithms on five groups of data with different numbers of nodes.

Number of nodes	Dataset	Best solution (km)			Deviation
		IWSA	ACS	ACGA	IWSA/ACS	IWSA/ACGA
10	1	1425	1682	1537	15.28%	7.29%
	2	1974	2020	1706	2.27%	−15.70%
	3	1613	1791	1684	9.97%	4.26%
	4	1690	1404	1755	−20.35%	3.71%
	5	1311	1724	2140	23.96%	38.72%
20	1	2828	3222	3515	12.22%	19.54%
	2	2439	2816	2127	13.39%	−14.67%
	3	1991	2978	3132	33.13%	36.41%
	4	1862	2305	2820	19.19%	33.96%
	5	2793	2751	3239	−1.53%	13.77%
50	1	4171	7457	8764	44.07%	52.41%
	2	3982	8799	7861	54.74%	49.34%
	3	3540	7145	6663	50.46%	46.88%
	4	4079	8299	7623	50.85%	46.49%
	5	4178	6743	7165	38.04%	41.69%
100	1	6279	17650	15137	64.42%	58.52%
	2	6880	15287	15973	55.00%	56.93%
	3	6505	17772	15919	63.40%	59.13%
	4	6606	16373	15208	59.65%	56.56%
	5	5587	15432	16756	63.80%	66.66%

---

 | Show Table

DownLoad: CSV

Figure 6.  Calculation results of IWSA and ACS with different number of nodes. Note:mean represents the mean of the results for the same set of sample data.

DownLoad: Full-Size Img PowerPoint

Figure 7.  Calculation results of IWSA and ACGA with different number of nodes.

DownLoad: Full-Size Img PowerPoint

It can be seen from Table 4 that when the three algorithms are used to solve the five groups of sample data under different node numbers, the IWSA has better solution performance than the ACS and the ACGA with the increase of number of nodes. Among the five groups of samples, when the number of nodes is 10, the best solutions of groups 1, 2, 3 and 5 samples solved by the IWSA are better than that of the ACS, and the best solutions of groups 1, 3, 4 and 5 samples solved by the IWSA are better than that of the ACGA. When the number of nodes is 20, the best solutions of groups 1, 2, 3 and 4 samples solved by the IWSA are better than that of the ACS, and the best solutions of groups 1, 3, 4 and 5 samples solved by the IWSA are better than that of the ACGA. This is because, ideally, the search area of individuals in the population is positively correlated with the size of the problem solution space. When the number of nodes is small, the solution space of the problem is relatively small, and the area that the individuals in the population can explore also becomes smaller when solving the problem. With the iterative update of the population, the probability of individuals in the population exploring the best solution with different algorithm rules is higher, so there will be cases in some sample results where the best solutions obtained by both the ACS and the ACGA are better than that of the IWSA. When the number of nodes is 50 and 100, the best solutions for the five groups of samples obtained by the IWSA are better than that of the ACS, and as the number of nodes increases the deviation of the best solutions between the IWSA and the ACS gradually increases, with the average deviation of 47.63 and 61.25%, respectively. Similarly, for the ACGA, with the increasing number of nodes, the deviation of the best solutions between the IWSA and the ACGA also increases gradually, and the average deviation reaches 47.36 and 59.56%, respectively.

Figure 6 shows that the best solutions obtained by the IWSA and the ACS have different degrees of superiority and inferiority under different number of nodes. Among them, when the number of nodes is 10, there are 71 groups of results calculated by the IWSA better than that of the ACS, that is, relative to ACS, the solving rate of the best solution of the IWSA (refers to the percentage of sample groups of one algorithm better than other algorithm in the experimental results) is 71%. When the number of nodes is 20, the solving rate of the IWSA is 98%. When the number of nodes is 50 and 100, the results solved by the IWSA are better than that of the ACS in 100 groups of samples, respectively, that is, the solving rate of the best solution of the IWSA is 100%. In addition, Figure 7 shows that when the number of nodes is 10, 87% of the best solution obtained by the IWSA is better than that of the ACGA in 100 groups of samples. When the number of nodes is 20, the best solution obtained by the IWSA is better than that of the ACGA in 93% of samples. When the number of nodes is 50 and 100, the best solution curves calculated by the IWSA are all below the average curve, while the best solution curves obtained by the ACGA are all above the average curve, which indicates that the best solutions obtained by the IWSA are better than that of the ACGA.

Combining Table 4 and Figures 6 and 7, the IWSA, the ACS and the ACGA have different performance in solving the same problem, which is related to the scale of the problem. With the increase of the number of nodes, the solving rate of the IWSA gradually increases, and it is better than that of the ACS and the ACGA. As we know, the optimization algorithm is prone to fall into the local extreme value prematurely, which makes individuals quickly cluster in a small area during the search process. Moreover, we also found that the population diversity indeed decreased rapidly through the calculation of Eq (10) during the experiments. Therefore, by abandoning some individuals in our improvement, the reduction speed of population diversity is effectively delayed, which allows the algorithm has enough time to find the global optimal solution in the search space. In addition, due to the irregular distribution of initial positions of individuals in the solution space caused by the randomness of the heuristic algorithm itself, the moving positions of individuals are also unstable in the process of iterative updating. Therefore, the individual can only move in the direction that may be better under the corresponding rules, and its final position is also related to the specific movement rules and parameters of the algorithm. For small-scale problems, the solution space is relatively small, and the area that the algorithm needs to explore in the optimization process is also smaller. At this time, when three algorithms are used to solve the problem, the best solutions obtained by the ACS and the ACGA are better than that of the IWSA. However, due to the advantages of individual movement rules designed in the IWSA, individuals can search more carefully in the solution space, so that the solving rate of the best solution of the IWSA is still higher than that of the ACS and the ACGA. The performance of the IWSA is much better than that of the ACS and the ACGA in solving large-scale problems. This may be because the grouping mechanism of IWSA allows individuals to conduct more detailed searches within the smaller area of the solution space. Moreover, when the individual search space is extended based on different weights, it can search for optimization multiple times in a larger search area based on the better individual. When the diversity of the population is below a certain threshold, the expansion of the search space can not only prevent individuals from congregating in a small area to avoid the algorithm falling into the local optimum prematurely, but also improve the probability of individuals exploring the best solution to a certain extent.

5.4.   Time complexity analysis of IWSA

Although the advanced individual optimization method enables the IWSA to obtain the best solution with higher quality in solving, its grouping mechanism and the detailed exploration of individuals in the population among the solution space make its time cost higher than that of the ACS and the ACGA. In order to better analyze the time cost of the IWSA, 100 groups of data under different node numbers are divided into 10 groups, and the IWSA, the ACS and the ACGA are used to calculate. The average time cost of each group of samples is shown in Table 5.

Table 5.  Calculation time of the three algorithms.

Groups	Calculation time of IWSA (s)				Calculation time of ACS (s)				Calculation time of ACGA (s)
	10	20	50	100	10	20	50	100	10	20	50	100
1	2.788	7.824	33.550	115.154	0.240	0.968	10.456	56.289	0.081	0.114	0.252	0.599
2	2.613	7.365	31.796	123.056	0.243	1.286	9.755	59.531	0.085	0.119	0.278	0.555
3	2.472	8.083	32.471	117.133	0.214	0.997	10.300	58.496	0.085	0.128	0.269	0.624
4	1.973	8.137	31.961	114.464	0.242	0.977	10.636	55.421	0.082	0.115	0.269	0.624
5	2.581	7.061	34.702	113.447	0.214	0.962	10.539	56.041	0.107	0.123	0.271	0.609
6	2.862	7.213	32.138	114.730	0.247	0.995	9.707	55.431	0.080	0.127	0.283	0.553
7	2.102	8.021	32.105	114.785	0.230	0.936	9.811	56.834	0.089	0.119	0.249	0.607
8	2.121	7.749	32.648	114.843	0.238	1.027	9.566	61.343	0.090	0.116	0.267	0.610
9	2.201	7.377	31.975	118.408	0.220	0.954	8.880	60.475	0.081	0.122	0.273	0.589
10	1.974	7.760	31.797	114.001	0.257	0.974	9.464	53.117	0.087	0.117	0.261	0.576
Average calculation time (s)	2.369	7.659	32.514	116.002	0.234	1.008	9.911	57.298	0.087	0.120	0.267	0.595

---

 | Show Table

DownLoad: CSV

Table 5 shows that the solution time cost of the three algorithms increases as the number of nodes increases, and the time cost required by the IWSA to solve the problem is the highest. The average time cost required by both the IWSA and the ACGA increases with the number of nodes. The relationship between them increases by multiple, reaching 195 times when the number of nodes is 100. On the contrary, the relationship between the average time cost required by both the IWSA and the ACS decreases by multiple. When the number of nodes is 10, 20, 50 and 100, the relationship between them is roughly a multiple of 10, 7, 3 and 2.

The multiplicative relationships between the time cost consumed by different algorithms show different trends with the change of the number of nodes. This is because the specific solution rules designed by different algorithms are different during the solution process, such as population initialization method, individual update operations and so on. As a result, the solution times required by the three algorithms are quite different for the same number of iterations and population size, that is, different algorithms have different time complexity due to different specific solution rules. Compared with the ACS and the ACGA, the update operation of population individuals in the IWSA is more complicated. After the hierarchical clustering algorithm is used to group individuals in the population, the individuals move to their guide individuals according to the movement rules of the group. When the population diversity is lower than the threshold, the search area needs to be expanded based on the current individuals, which makes everyone carry out a more detailed search in its neighboring area. Although it improves the probability of individuals searching for the best solution, it also increases the time cost.

Moreover, in the IWSA, the number of execution steps of the individual update operation is determined by the dimension of the individual (the number of nodes that are not currently accessed). Therefore, in the best case, the maximum number of execution steps of its update operation is $number = COUNT \times MAXInter \times pop\_num \times n$ ($COUNT$ and $MAXInter$ are the maximum number of iterations, $n$ is the number of nodes (problem size)), and the best time complexity is ${\rm O}(n)$. In the worst case, the maximum number of execution steps of its update operation is $number = COUNT \times MAXInter \times pop\_num \times n \times pop\_num \times n$, and the worst time complexity is ${\rm O}({n^2})$. However, due to the different movement rules when updating population individuals, the time complexities of ACS and ACGA are ${\rm O}(n)$ and ${\rm O}(1)$, respectively. The time costs of the three algorithms are different in solving problem of the same size due to the difference in the time complexity caused by the different solution rules within the algorithms, which makes the time difference more apparent in numerical terms with the increase of the problem size.

6.   Conclusions

Based on the dynamic change characteristics of road conditions in disaster areas, this paper proposes an IWSA for updating the material scheduling scheme to solve the shortcomings of the current MEDDRC problem such as the poor timeliness of updating the scheduling scheme and the existing algorithms easy to fall into the local extreme value prematurely in solving. The ISCW is constructed to formulate the initial vehicle routing, and the group movement strategy and different weights strategy are proposed to optimize the WSA from two aspects of individual movement mode and expanding the search space of individuals in the population respectively. The IWSA realizes the purpose of updating the vehicle scheduling scheme in real-time under the change of road conditions, which can provide support for decision-makers to formulate effective material emergency scheduling strategies in major natural disaster events.

At the same time, when the number of nodes is 10, 20, 50 and 100, the simulation calculation is carried out on four randomly generated data sets containing 100 groups of samples. Through the analysis of the experimental data, the following conclusions are obtained:

1) When generating the initial driving path of the vehicle, the best solution obtained by the ISCW is better than that of both the scan algorithm and the C-W algorithm. Moreover, with the increase of the problem scale, the reduction degree of the average distance solved by the ISCW increases greatly.

2) In the case of road conditions changing, the IWSA is used to plan the vehicle scheduling scheme. When the number of nodes is 10 and 20, the proportion of the best solution obtained by the IWSA is better than that of both the ACS and the ACGA is 71, 87% and 98, 93%, respectively. When the number of nodes is 50 and 100, the best solutions obtained by IWSA are better than that of both the ACS and the ACGA. And with the increase of the number of nodes, the average distance of the best solution obtained by the IWSA is much better than that of both the ACS and the ACGA.

3) In terms of algorithm complexity, as the individuals in the IWSA move within the solution space, everyone performs detailed searches in their neighboring area, so a large amount of time is spent on searching for the optimal solution. Therefore, compared with the ACS and the ACGA, the IWSA spends a higher time cost.

4) Although the proposed algorithm can solve the MEDDRC problem and has good results, it still has the disadvantage of high time complexity, which is related to its special solution rules. Besides, the values of the parameters can also have a large impact on the performance of the algorithm. Therefore, a suitable set of parameters is also particularly important for optimization algorithms. It is our next priority that how to find an optimal set of algorithm parameters to improve the model solving performance while reducing the computation time.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was supported by Natural Science Innovative Fund of Henan University of Technology (2021ZKCJ18), Joint Fund of Science and Technology R & D Program of Henan Province (222103810078) and Science and Technology Research Project of Henan Province (232102211010).

Conflict of interest

The authors declare there is no conflict of interest.

Supplementary

Table .  The abbreviations used in this article and their full names.

Abbreviations	Full Name
MED	Material Emergency Dispatching
DRP	Dynamic Routing Problem
SRP	Static Routing Problem
VRP	Vehicle Routing Problem
MEDDRC	Material Emergency Dispatching with Dynamic Road Conditions
WSA	Whale Swarm Algorithm
IWSA	Improved Whale Swarm Algorithm
C-W	Clarke-Wright
ISCW	Improved Scan Clarke-Wright
BCP	Branch Cut and Price
GA	Genetic Algorithm
PSO	Particle Swarm Optimization
LTMU	Long Term Memory Unit
ACS	Ant Colony System
ACGA	Adaptive Chaos Genetic Algorithm

---

 | Show Table

DownLoad: CSV