Multi-Drone Optimal Mission Assignment and 3D Path Planning for Disaster Rescue

Abstract

In a three-dimensional (3D) disaster rescue mission environment, multi-drone mission assignments and path planning are challenging. Aiming at this problem, a mission assignment method based on adaptive genetic algorithms (AGA) and a path planning method using sine–cosine particle swarm optimization (SCPSO) are proposed. First, an original 3D digital terrain model is constructed. Second, common threat sources in disaster rescue environments are modeled, including mountains, transmission towers, and severe weather. Third, a cost–revenue function that considers factors such as drone performance, demand for mission points, elevation cost, and threat sources, is formulated to assign missions to multiple drones. Fourth, an AGA is employed to realize the multi-drone mission assignment. To enhance convergence speed and optimize performance in finding the optimal solution, an AGA using both the roulette method and the elite retention method is proposed. Additionally, the parameters of the AGA are adjusted according to the changes in the fitness function. Furthermore, the improved circle algorithm is also used to preprocess the mission sequence for AGA. Finally, based on the sine–cosine function model, a SCPSO is proposed for planning the optimal flight path between adjacent task points. In addition, the inertia and acceleration coefficients of linear weights are designed for SCPSO so as to enhance its performance to escape the local minimum, explore the search space more thoroughly, and achieve the purpose of global optimization. A multitude of simulation experiments have demonstrated the validity of our method.

1. Introduction

Disaster events such as fires, floods, and landslides are characterized by randomness, dynamics, and urgency. In most cases, these disaster events are inevitable. If urgent measures are not taken, it will cause significant economic losses and threaten human life. Therefore, how to reduce the losses caused by disaster events is the key issue in emergency rescue [1]. With the advancement of science and technology, drones can complete complex tasks such as three-dimensional (3D) map reconstruction [2], emergency mapping [3], and environmental assessment [4]. Compared with traditional disaster rescue applications [5], multi-drone-based disaster rescue schemes can use drones to achieve damage assessment [6], material delivery [7], etc. In the context of disaster rescue, it is challenging for rescuers to promptly reach the affected area due to the complex terrain, obstructions from mountains and rivers, road collapse, and other factors. At this time, due to the characteristics of drones, their mission planning is particularly important [8,9]. Therefore, drones have emerged as pioneers in disaster rescue missions, swiftly reaching disaster areas, and reducing the harm caused by disasters to the public [10]. There are many factors affecting drone disaster rescue, including the complexity of the mission itself, environmental factors [11] (such as terrain obstacles, weather conditions, or electromagnetic interference), and drone performance factors [12] (such as maximum load, energy consumption, or maximum flight distance). How to allocate disaster rescue missions to multiple drones reasonably and how to plan the flight path of drones to make them perform disaster rescue missions optimally and safely are the focus of current research.

Over the past few decades, various methods have been developed to solve mission assignments and path planning [13,14,15]. Existing mission assignment methods are divided into two categories: traditional mathematical programming and heuristic algorithms. Traditional mathematical programming methods consider constraints and resource availability to address optimal mission assignments, including enumeration algorithms [16] and dynamic programming [17]. Although these methods have the advantages of computational efficiency and are simple to implement, they oversimplify the drone model and can only be applied to simple mission scenarios. Heuristic algorithms iteratively optimize mission assignment schemes through strategies based on experience or rules, including genetic algorithm (GA) [18], evolutionary algorithm (EA) [19], tabu search (TS) algorithm [20], simulated annealing (SA) algorithm [21], etc. Path planning methods are divided into five categories: graph-based methods accomplish path planning by constructing a robust path graph, such as the Voronoi diagram method [22]. Random sampling search algorithms generate a path by iteratively sampling points between the start and end points and connecting them based on various constraints, such as the rapidly-exploring random tree (RRT) [23,24] and probabilistic roadmap (PRM) [25]. Node-based optimal search algorithms construct a node topology to represent the path and employ heuristic functions to facilitate efficient path searching, such as the Dijkstra algorithm [26], the A-star (A*) algorithm [27], and the harmony search (HS) algorithm [28]. The artificial potential field method [29,30] realizes path planning by constructing gravitational fields and repulsive fields to simulate the interaction forces between objects. Additionally, the bionic evolutionary algorithms optimize path planning solutions through the simulation of biological evolution processes, such as the particle swarm optimization (PSO) algorithm [31,32] and the ant colony optimization (ACO) algorithm [33,34].

Since mission planning is a NP-hard problem, it can be effectively solved using heuristic algorithms. However, these algorithms have the disadvantages of slow convergence speed and the tendency to fall into local optimums. So far, there is no complete solution to the problem of drone mission planning. Thus, the exploration of alternative and more efficacious remedies becomes indispensable. In [35], the GA was employed to achieve a multi-drone cooperative reconnaissance mission. In [36], the mission planning model was established according to the mission clustering of drones, and a mission assignment approach using the K-means clustering method of an improved simulated annealing algorithm was proposed. In [37], a hierarchical mission assignment method was developed that decomposed the primitive problem into multiple subproblems, then used mixed integer programming and ACO to solve these sub-problems. In [38], a hybrid algorithm based on PSO and the metropolis criterion was designed to reduce the local optimum of PSO. In [39], a differential evolution algorithm combined with a quantum particle swarm optimization algorithm was introduced for path planning in a drone marine environment. In [40], a pigeon swarm optimization considering path length, path curvature, and path risk was also created.

Clearly, extensive simulation results have demonstrated that when dealing with the problem of mission planning, the heuristic algorithm still suffers from the issue of premature convergence during the evolutionary process. Although the existing research has improved these algorithms and proposed many approaches with better performance, there are still some shortcomings among them: first, with the combination of algorithms, the calculation becomes more complicated. Second, global and local optimization abilities are two important factors to be considered in evaluating an algorithm [41], and it is difficult to maintain an effective balance between them. Third, most studies only consider two-dimensional environments, and the flight environment is too simple, which makes it difficult to apply the research results to actual flight. Therefore, this paper improves the GA and PSO and applies the improved algorithm to the mission planning of multiple drones in complex three-dimensional disaster rescue environments.

This paper presents a mission assignment and path planning system for multiple drones for disaster rescue applications. Additionally, it is mainly to coordinate and control multiple drones for efficient material distribution within a 3D disaster area. First, an original digital terrain and three common threat sources [42] for disaster rescue environments are designed, including the mountain, transmission tower, and severe weather. Then, a cost–revenue function considering drone performance, mission point requirements, elevation cost, and threat source factors is proposed. When using the AGA for mission assignment, the improved circle algorithm, adaptive crossover rate and mutation rate, and a strategy that uses both roulette and elite retention methods are used to increase the properties of the AGA. Finally, based on the sine–cosine function model, a SCPSO is proposed for searching the optimal flight path between adjacent task points. The inertia and acceleration coefficients of linear weights are designed to further maintain an effective balance of SCPSO between global exploration and local development.

The main contributions of this paper include: (1) A modeling method for multi-drone 3D disaster rescue is presented. The original digital terrain is defined. Three threat sources are proposed. A new cost–revenue function is established and formulated as a constrained, multi-objective optimization problem. (2) The AGA and SCPSO algorithms are proposed. A strategy of using both roulette and the elite retention method is proposed, and the capacity of AGA is improved by combining the improved circle algorithm. The inertia and acceleration coefficients of linear weights are designed for SCPSO to increase optimization efficiency. By integrating the mission assignment method with the path planning algorithm, multi-drone mission allocation and path planning in a complex 3D environment can be implemented.

In the following sections, first, the description of the multi-drone disaster rescue problem will be presented in Section 2. Second, the digital terrain, threat sources, cost–revenue function, AGA, and SCPSO modeling methods are introduced in Section 3. Third, a series of experiments and simulations under complex 3D terrain are carried out and discussed in Section 4. Finally, the conclusion is given in Section 5.

2. Problem Descriptions

The method proposed in this paper is utilized to deploy multiple drones that originate from a base and perform disaster rescue missions at different target locations within the disaster area. The primary aim of these missions is to efficiently distribute vital supplies to the disaster-stricken region. Figure 1 shows a sketch of multi-drone disaster relief. Multiple mission points are situated in different locations within the disaster area, each requiring a specific quantity of rescue materials. Commencing from a designated starting point, multiple drones initiate the mission, transporting the materials to each mission point in the assigned order, and subsequently returning. The flight of the drone is affected by many factors, such as the mountains, transmitting towers, and severe weather. The mountain may affect the flight safety of drones. The transmitting tower may affect the communication system of the drone. Moreover, the severe weather may change the original flight path of the drone. In addition, flight distance and load will also affect the assignment of the entire mission.

In the aforementioned application scenario, multiple drones traverse N mission nodes in a specific sequence, starting from the base. Once the mission is accomplished, all drones return to the base. Assuming that N mission points are {x₁, x₂, …, x_N}, the number of drones is M. Since the drone needs to return to the base after the mission is completed, we need to copy the base node as an M virtual base node. This enables the decomposition of the multi-path generation problem into one-way travel for a group of drones. This approach ensures that the flight path from the first drone to the last drone is connected from beginning to end, forming a single loop, as shown in Figure 2. Among them, O is the base node, and O₁, O₂, and O₃ are the virtual base nodes.

If the mission of drone u includes the flight track from the mission point x_i to the mission point x_j, the decision variable is defined as $x_{ij}^u$ = 1, otherwise $x_{ij}^u$ = 0. The total flight distance of multiple drones is expressed as Equation (1). Equations (2) and (3) ensure that each mission point is visited only once during a mission cycle. Equations (4) and (5) specify that each drone can only take off and land at the base. Formula (6) ensures that the total flight displacement of each drone does not exceed its maximum flight displacement constraint. Equation (7) means that the load of each drone does not exceed its maximum load constraint.

$$distance={\displaystyle \sum _{u=1}^M{\displaystyle \sum _{i=0}^N{\displaystyle \sum _{j=0, j\ne i}^N{x}_{ij}^u·s_{ij}}}}$$

(1)

$${\displaystyle \sum _{j=0,j\ne i}^N{\displaystyle \sum _{u=1}^M{x}_{ij}^u=1}}, i=1, \dots , N$$

(2)

$${\displaystyle \sum _{i=0,i\ne j}^N{\displaystyle \sum _{u=1}^M{x}_{ij}^u=1}}, j=1, \dots , N$$

(3)

$${\displaystyle \sum _{i=1}^N{x}_{0i}^u=1}, \forall u$$

(4)

$${\displaystyle \sum _{i=1}^N{x}_{i0}^u=1}, \forall u$$

(5)

$${\displaystyle \sum _{i=0}^N{\displaystyle \sum _{j=0,j\ne i}^N{x}_{ij}^u·s_{ij}\le s_{max}}}, \forall u$$

(6)

$${\displaystyle \sum _{i=0}^N{\displaystyle \sum _{j=0,j\ne i}^N{x}_{ij}^u·q_i\le q_{max}}}, \forall u$$

(7)

$$x_{ij}=\left\{\begin{array}{l}1, if drone u travel from i to j;\\ 0, otherwise.\end{array}\right. i\ne j, i, j\in \left\{0,1,2,\dots ,N\right\}, u\in \left\{1,2,\dots ,M\right\}$$

(8)

$$s_{ij}=\sqrt{{(X_i-X_j)}^2+{(Y_i-Y_j)}^2+{(Z_i-Z_j)}^2}$$

(9)

where s_ij is the displacement between any pairs of target points (x_i, x_j); (X_i, Y_i, Z_i) and (X_j, Y_j, Z_j) are the 3D coordinates of mission points x_i and x_j, respectively; s_max means the maximum displacement of drone driving; q_i represents the material demand of mission point x_i; and q_max means the maximum load of the drone.

3. Environment Modeling and Optimization Method

3.1. Proposed Flowchart

For multi-drone disaster rescue, it is crucial to design the mission sequence required by multiple drones for the whole rescue mission—that is, to plan the optimal mission assignment. In addition, to ensure a safe and effective drone flight, it is necessary to consider various threat sources when planning the drone’s flight path. Figure 3 shows the main process of multi-drone disaster rescue. First, before the mission planning of drones, it is essential to model the drone disaster relief environment. Then, a cost–revenue function is proposed based on drone performance and rescue scenarios. Finally, the AGA is used to realize the optimal mission assignment of multi-drones, and the SCPSO is employed for planning the optimal flight path between adjacent target points.

3.2. Construction of Original Digital Terrain

In the mission area, the safety constraints of the terrain should be considered while carrying out the material allocation. The original digital terrain refers to the basic surface appearance and terrain height in the 3D environment. Therefore, before carrying out multi-drone mission planning in a complex 3D disaster rescue environment, the original digital terrain must be constructed first. Without loss of generality, one of the reference models of the original digital terrain can be computed by the method in Formula (10). Other 3D terrains can also be obtained by satellite surveying and mapping.

$$z_1(x,y)=sin(y+a)+b\hspace{0.17em}sin(x)+c\hspace{0.17em}cos(d\sqrt{y^2+x^2})+e\hspace{0.17em}sin(e\sqrt{y^2+x^2})+f\hspace{0.17em}cos(y)$$

(10)

where (x, y) is the two-dimensional coordinate of the horizontal projection plane; z₁(x, y) represents the height information corresponding to the horizontal projection plane coordinates; a, b, c, d, e, and f mean the undetermined terrain coefficients. These coefficients control the amplitude of the ups and downs of the map and can simulate the actual terrain of various complex landforms. Currently, our hypothetical disaster rescue mission is conducted in the mountainous region of southwestern China, which is characterized by a complex climate and an extensive seismic belt, making it prone to natural disasters such as landslides and earthquakes. To describe the terrain of southwest China, we employ the subjective evaluation method for continuous parameter adjustments. Ultimately, the original terrain coefficients are set as a = 1.5 π, b = 5, c = 3, d = 5, e = 3, and f = 10, aiming to simulate the authentic topography of the region. Figure 4 is the corresponding simulated digital terrain.

3.3. The Modeling Method of Threat Source

It is assumed that there are three external threat sources for drone flight in our system: the mountains, transmission towers, and severe weather.

3.3.1. The Mountain Threat Source

The flight safety of the drone will be affected by the mountains. If the flight trajectory of the drone is not properly controlled, it may collide with the mountains and cause a crash. Therefore, it is important to plan the optimal scheme and path so that the drone can safely pass through the mountain area. The expression of the mountain mathematical model in this paper is shown in Equation (11). Additionally, a bi-Gaussian mixture model (bi-GMM) [43] (see Equations (12) and (13)) is used as the cost function of the mountain threat sources. This is due to its capacity to offer a smooth and continuous approximation of their shapes, rendering it better suited for real-world terrain in contrast to the traditional cone model that simulates mountain distribution.

$$z_2(x,y)={\displaystyle \sum _{N\_m=1}^{n_1}{h}_{N\_m}exp\left[-{\left(\frac{x-x_{O,N\_m}}{x_{sl,N\_m}}\right)}^2-{\left(\frac{y-y_{O,N\_m}}{y_{sl,N\_m}}\right)}^2\right]}$$

(11)

$$f_m(X_i,Y_i)=\frac{1}{\pi ·x_{sl,N\_m}·y_{sl,N\_m}}·exp\left[-{\left(\frac{X_i-x_{O,N\_m}}{x_{sl,N\_m}}\right)}^2-{\left(\frac{Y_i-y_{O,N\_m}}{y_{sl,N\_m}}\right)}^2\right]$$

(12)

$$C_m={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{N\_m=1}^{n_1}{k}_1·f_m(X_i,Y_i)}}$$

(13)

where z₁(x, y) means the peak elevation corresponding to point (x, y); n₁ is the number of mountain threat sources; i is the ith mission point; N_m means the N_mth mountain threat source; h_{N_m} is the height of the N_mth mountain threat source; (x_{O,N_m}, y_{O,N_m}) represents the central coordinate of the N_mth mountain threat source; x_{sl,N_m} and y_{sl,N_m} mean the slope parameters of the N_mth mountain threat source along the X-axis and Y-axis, respectively. Greater values of x_{sl,N_m} and y_{sl,N_m} indicate a flatter profile for the corresponding mountain threat source, while smaller values indicate a steeper profile. By manipulating the parameters h_{N_m}, x_{sl,N_m}, and y_{sl,N_m}, it is possible to simulate mountain threat sources with varying heights and contours. The symbol f_m(X_i, Y_i) represents the probability density function of the bi-GMM model; (X_i, Y_i) is the horizontal projection coordinate of the ith mission point; C_m is the cost function of the mountain threat source; N is the number of mission points; and k₁ is the corresponding weight coefficient.

Finally, Equation (14) merges the mountain threat source with the original digital terrain to form a 3D environment-equivalent digital terrain, as shown in Figure 5.

$$z(x,y)=max\left[z_1(x,y),z_2(x,y)\right]$$

(14)

where z(x, y) is the terrain height corresponding to the point (x, y); max means the function of the maximum value.

3.3.2. The Transmission Tower Threat Source

The electromagnetic interference emitted by the transmission tower can disrupt the communication and navigation systems of drones [44], potentially leading to the failure of accurate mission execution. Therefore, the position and influence range of the transmission tower needs to be considered when the drone is flying so as to reduce the risk of flight. The electromagnetic interference emitted by the tower can be regarded as the diffusion of a spherical model into space. In addition, the emitted electromagnetic interference should not have a specific limit, beyond which the risk of damage is zero. Therefore, a probability density model of the transmission tower threat source is proposed in this paper, as shown in Equations (15)–(17).

$$f_t(d_{i,N\_t}^t,R_{N\_t}^t)=\left\{\begin{array}{ll}exp\left(-\frac{d_{i,N\_t}^t}{{\sigma }_{i,N\_t}}\right),& d_{i,N\_t}^t>R_{N\_t}^t\\ \infty ,& d_{i,N\_t}^t\le R_{N\_t}^t\end{array}\right.$$

(15)

$${\sigma }_{i,N\_t}=\frac{R_{N\_t}^t}{lg\left(\frac{L_i}{2}\right)}$$

(16)

$$C_t={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{N\_t=1}^{n_2}{k}_2·f_{tower}(d_{i,N\_t}^t,R_{N\_t}^t)}}$$

(17)

where f_t($d_{i,N\_t}^t$, $R_{i,N\_t}^t$) represents the probability density model of the signal transmission tower threat source; N_t means the N_tth transmission tower threat source; i is the ith mission point; $d_{i,N\_t}^t$ is the distance from the ith mission point to the threat source center of the N_tth transmitting tower; $R_{i,N\_t}^t$ is the minimum safe distance of the N_tth transmission tower; σ_{i,N_t} means the control parameter of the probability density model of the transmitting tower; L_i represents the distance from the current mission point to the ith mission point; C_t is the transmission tower threat source cost function; N is the number of mission points; n₂ is the number of transmission tower threat sources; and k₂ is the corresponding weight coefficient.

3.3.3. The Severe Weather Threat Source

Usually, when the drone is performing its mission, it may encounter local severe weather (such as storms, heavy rain, lightning, etc.). If the drone is compelled to operate in that environment, it may damage the motor or sensor of the drone or even crash. Therefore, drones need to change lanes to avoid severe weather. In this paper, severe weather events such as storms, rainstorms, and lightning are abstracted and simplified into cylindrical threat areas. A probability density model of the source of the threat of severe weather is shown in Equations (18) and (19).

$$f_w(k_{N\_w}^w,d_{i,N\_w}^w,r_{N\_w}^w,R_{N\_w}^w)=\left\{\begin{array}{l}\infty ,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{d}_{i,N\_w}^w\le r_{N\_w}^w\\ \frac{k_{N\_w}^w}{{(R_{N\_w}^w-r_{N\_w}^w)}^2}{(d_{i,N\_w}^w-R_{N\_w}^w)}^2,r_{N\_w}^w<d_{i,N\_w}^w\le R_{N\_w}^w\\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{d}_{i,N\_w}^w>R_{N\_w}^w \end{array}\right.$$

(18)

$$C_w={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{N\_w=1}^{n_3}{k}_3·f_w(k_{N\_w}^w,d_{i,N\_w}^w,r_{N\_w}^w,R_{N\_w}^w)}}$$

(19)

where f_w($k_{N\_w}^w, d_{i,N\_w}^w$, $r_{N\_w}^w$, $R_{N\_w}^w$) means the probability density model of the threat source of severe weather; N_w is the N_wth severe weather threat source; i is the ith mission point; symbol $k_{N\_w}^w$ means the risk coefficient of the N_wth severe weather; $d_{i,N\_w}^w$ represents the distance from the ith mission point to the center of the N_wth severe weather threat source; $r_{N\_w}^w$ is the minimum safe distance of the N_wth severe weather threat source; $R_{N\_w}^w$ is the maximum impact range of the N_wth severe weather threat source; C_w means the cost function of the severe weather threat source; N is the number of mission points; n₃ is the number of severe weather threat sources; and k₃ is the corresponding weight coefficient.

3.4. The Cost–Revenue Function

In this paper, a cost–benefit function for multi-drone mission assignments is designed. When constructing the cost function, the influence of threat sources is considered, and the calculation methods are given in Equations (13), (17) and (19). The drone’s flight distance and flight elevation are also considered. Equations (20) and (21) give their calculation methods. Equation (22) represents the total cost function. The revenue function is expressed as the revenue from completing a mission. The greater the demand for materials at the mission point, the higher the revenue. Therefore, the revenue function can be expressed by Equation (23). Finally, the total cost–revenue function is expressed as Equation (24).

$$C_{distance}=k_4·distance$$

(20)

$$C_{high}=k_5·{\displaystyle \sum _{u=1}^M{\displaystyle \sum _{i=0}^N{\displaystyle \sum _{j=0,j\ne i}^N{x}_{ij}^u·m_{ij}^u·\left|h_{ij}\right|}}}$$

(21)

$$C_{total}=C_m+C_t+C_w+C_{distance}+C_{high}$$

(22)

$$R_{total}=k_r·q_i$$

(23)

$$T={\omega }_1{C}_{total}-{\omega }_2{R}_{total}$$

(24)

where C_distance is the distance cost; k₄ is the weight; C_high is the elevation cost; k₅ is the weight; the symbol $m_{ij}^u$ means the quality of the drone u flying from the ith mission point to the jth mission point; h_ij represents the height difference between the ith mission point and the jth mission point; C_total means the total cost; R_total is the total revenue; k_r is a weight, which is set to k_r = 2 in this paper; q_i means the material demand of the ith mission point; T represents a cost–revenue function; and ω₁ and ω₂ are weights, which are set to ω₁ = 1 and ω₂ = 1 in this paper.

3.5. Mission Planning Algorithms: AGA and SCPSO

An AGA method is designed for a multi-drone mission assignment, and a SCPSO algorithm is designed for 3D path planning.

3.5.1. The Mission Assignment of Multi-Drone Based on AGA

GA is a computational model that simulates the biological evolution processes of natural selection and genetics. In GA, each problem solution is represented as a set of chromosomes. Through selection, crossover, and mutation operations in each generation, GA gradually improves the quality of its solutions. The following are the calculation steps of the AGA algorithm.

• The population initialization

The initial solution to the mission assignment is generated by random numbers. Let us define the population size N_AGA, the maximum iteration time T_AGA, the maximum crossover rate P_C,max, the minimum crossover rate P_C,min, the maximum mutation rate P_M,max, and the minimum mutation rate P_M,min.

The initial chromosome is a Hamilton cycle [45]. In order to accelerate the convergence speed of the algorithm and obtain a better initial solution for each generation, the improved circle algorithm is used in AGA. That is, to judge whether two groups of adjacent gene points v_p−1v_p(p = 2, …, N − 2) and v_rv_r+1(r = p + 1, …, N − 1) in the chromosome meet the Equation (25). If it is satisfied, then v_p and v_r are arranged in the reverse order; otherwise, the next judgment is made until all gene points are traversed.

$$d_{ch}\left(v_{p-1}{v}_r)+d_{ch}\right(v_p{v}_{r+1})<d_{ch}(v_{p-1}{v}_p)+d_{ch}\left(v_r{v}_{r+1}\right)$$

(25)

where d_ch means the distance between two gene points (mission points); v_p is the pth gene point; and v_r is the rth gene point.

2.

The calculation of chromosome fitness

The chromosome fitness is the mission assignment fitness function in this paper, and its calculation is shown in Equation (26).

$$f_{tAGA,chr}^{AGA}=T$$

(26)

where $f_{tAGA,chr}^{AGA}$ represents the fitness function of the chrth chromosome in the tAGAth generation; T is defined in Equation (24).

3.

The evolution of AGA

This step implements three operations: selection, crossover, and mutation. First, in the selection stage, the better the fitness, the greater the probability of producing excellent individuals. The method of roulette to choose (see Equation (27)) is used in this paper. Second, the two-point crossover is used in the cross-point phase; that is, two crossover points are randomly set for the chromosome, and then the chromosome gene sequence is cross-changed. Finally, the mutation phase is realized by randomly exchanging two genes in an individual.

$$P_{S,tAGA,chr}=\frac{f_{tAGA,chr}^{AGA}}{{\displaystyle \sum _{chr=1}^{N_{AGA}}{f}_{tAGA,chr}^{AGA}}}$$

(27)

where P_S,tAGA,chr means the selection rate for the chrth chromosome of the tAGAth generation.

Without losing generality, we take seven gene points as an example. Figure 6 displays an example of AGA evolution. In Figure 6, two parent chromosomes are {7, 5, 3, 6, 2, 1, 4} and {4, 3, 1, 5, 2, 7, 6}, respectively, where the number represents the mission point. These chromosomes produce new fetus {7, 3, 1, 5, 6, 2, 4} through the two-point crossover. Finally, two gene points are randomly switched to obtain child {7, 6, 1, 5, 3, 2, 4}.

In order to improve the convergence speed and optimization effect of the algorithm, some computational strategies are designed in this paper. The crossover rate and mutation rate of AGA do not use setting values but are adaptively selected according to fitness changes, as shown in Equations (28) and (29).

$$P_{C,tAGA,chr}=\left\{\begin{array}{ll}{k}_{AGA1}·P_{C,max}-\frac{(P_{C,max}-P_{C,min})(f_{tAGA,max}^{AGA}-f_{tAGA,bigger}^{AGA})}{(f_{tAGA,max}^{AGA}-f_{tAGA,mean}^{AGA})},& f_{tAGA,bigger}^{AGA}\ge f_{tAGA,mean}^{AGA}\\ k_{AGA2}·P_{C ,max},& f_{tAGA,bigger}^{AGA}<f_{tAGA,mean}^{AGA}\end{array}\right.$$

(28)

$$P_{M,tAGA,chr}=\left\{\begin{array}{ll}{k}_{AGA3}·P_{M ,max}-\frac{(P_{M,max}-P_{M,min})(f_{tAGA,max}^{AGA}-f_{tAGA,thr}^{AGA})}{(f_{tAGA,max}^{AGA}-f_{tAGA,mean}^{AGA})},& f_{tAGA,thr}^{AGA}\ge f_{tAGA,mean}^{AGA}\\ k_{AGA4}·P_{M ,max},& f_{tAGA,thr}^{AGA}<f_{tAGA,mean}^{AGA}\end{array}\right.$$

(29)

where P_C,tAGA,chr means the crossover rate of the chrth chromosome in the tAGAth generation; P_M,tAGA,chr represents the mutation rate of the chrth chromosome in the tAGAth generation; k_AGA1, k_AGA2, k_AGA3; and k_AGA4 are the control coefficients, which are set to k_AGA1 = 1, k_AGA2 = 1.05, k_AGA3 = 1, and k_AGA4 = 1.25 in this paper. P_C,max and P_M,max represent the maximum values of crossover rate and mutation rate; P_C,min and P_M,min represent the minimum values of crossover rate and mutation rate; $f_{tAGA,max}^{AGA}$ means the individual with the largest fitness in the tAGAth generation; $f_{tAGA,bigger}^{AGA}$ represents the bigger fitness in the two chromosomes to be crossed in the tAGAth generation and $f_{tAGA,mean}^{AGA}$ is the average fitness of all chromosomes in the tAGAth generation.

In order to prevent the degradation of the genetic algorithm, the elite retention method for each iteration is adopted in this paper. That is, we retain the optimal individual of the generation for the next generation, do not participate in the evolution process, and compare it with the optimal individual of the next generation. If it is better than the next generation of optimal individuals, the algorithm continues to retain them; otherwise, the next-generation optimal individual is retained as an elite.

3.5.2. The Path Planning between Mission Points Based on SCPSO

PSO is derived from the study of bird predation behavior. The basic idea of the algorithm is to initialize a set of particles within the solution space that are explored in the search space. The characteristics of particles are expressed by velocity, position, and fitness. The velocity of each particle represents the moving direction of the particle; the position indicates the candidate solution of the optimization problem; and the fitness represents the quality of a solution. Each particle continuously adjusts its position through individual experience and group collaboration to find the optimal solution. In the path planning problem of this paper, each particle represents the coordinate of the path-control point. The following are the calculation steps of the SCPSO algorithm.

• The population initialization

In this step, the number of particles N_SC, the number of path-control points N_C, and the maximum number of iterations T_SC are defined. The initial particle position (initial solution) is randomly generated in space.

2.

The fitness function calculation

SCPSO calculates the fitness function of each particle at each iteration. The fitness function of path planning in this paper is expressed as the total distance of drone flight, as shown in Equations (30) and (31).

$$f_{tSC,par}^{SCPSO}={\displaystyle \sum _{u=1}^M{\displaystyle \sum _{i=0}^N{\displaystyle \sum _{j=0,j\ne i}^N{x}_{ij}^u·l_{ij,tSC,par}}}}$$

(30)

$$\begin{gathered} l_{ij,tSC,par}=PT·{\displaystyle \sum _{s=1}^{np-1}‖B_{tSC,par,s+1}^{ij}-B_{tSC,par,s}^{ij}‖} \\ PT=\left\{\begin{array}{ll}1000,& if the path point is inside the threat source\\ 1,& otherwise\end{array}\right. \end{gathered}$$

(31)

where $f_{tSC,par}^{SCPSO}$ represents the fitness function of the parth particle in the tSCth generation; PT is a penalty factor; l_ij,tSC,par means the length of the three-dimensional path planned by the parth particle of the tSCth generation of SCPSO between the ith mission point and the jth mission point; np is the number of path points of the three-dimensional path curve generated by the path-control point N_C; $B_{tSC,par,s}^{ij}$ means the sth path point of the three-dimensional path planned by the parth particle of the tSCth generation of SCPSO between the ith mission point and the jth mission point.

3.

The iteration of SCPSO

Particle velocity is affected by inertia, individual optimal values, and global optimal values. Each particle adjusts the speed according to the current position, individual optimal value, and global optimal value, and then updates its position. The SCPSO algorithm proposed in this paper uses the properties of sine and cosine functions to improve the optimization ability of the algorithm by adaptively changing the amplitude of sine and cosine functions. In addition, in order to further balance the global exploration and local development capabilities of the algorithm, the inertia coefficient and acceleration coefficient of linear weights are also designed. Equations (32)–(37) show the speed and position update methods of the SCPSO algorithm.

$$v_{par,dim}^{tSC+1}=w·v_{par,dim}^{tSC}+c_1·rand\left(\right)·(p_{best,par,dim}-x_{par,dim}^{tSC})+x_{sc}$$

(32)

$$x_{sc}=\left\{\begin{array}{l}{R}_1·sin(R_2)·(g_{best,dim}-x_{par,dim}^{tSC}), R_3\le 0.5\\ R_1·cos(R_2)·(g_{best,dim}-x_{par,dim}^{tSC}), R_3>0.5\end{array}\right.$$

(33)

$$x_{par,dim}^{tSC+1}=x_{par,dim}^{tSC}+v_{par,dim}^{tSC}$$

(34)

$$w=w_{max}-(w_{max}-w_{min})·\frac{tSC}{T_{SC}}$$

(35)

$$c_1=c_{1max}-(c_{1max}-c_{1min})·\frac{tSC}{T_{SC}}$$

(36)

$$R_1=R_{1min}+(R_{1max}-R_{1min})·\frac{tSC}{T_{SC}}$$

(37)

where $v_{par,dim}^{tSC}$ is the velocity of the parth particle of the tSCth generation on the dimth dimension; w represents the inertia weight coefficient; c₁ means the individual acceleration coefficient, which is adaptively adjusted by Equations (35) and (36), respectively; rand() is the random number of [0, 1]; p_best,par,dim means the optimal position of the parth particle on the dimth dimension; $x_{par,dim}^{tSC}$ represents the position of the parth particle of the tSCth generation on the dimth dimension; x_sc means the sine and cosine components; R₁ is the control parameter, which mainly controls the amplitude of the sine and cosine functions and adjusts adaptively through Equation (37); R₂ and R₃ are random numbers obeying a uniform distribution, which are set to R₂ ∈ [0, 0.5π] and R₃ ∈ [0, 1] in this paper; g_best,dim is the global optimal position on the dimth dimension; w_max and w_min are the maximum and minimum values of inertia weight; c_1max and c_1min are the maximum and minimum values of individual acceleration coefficient; and R_1max and R_1min are the maximum and minimum values of sine and cosine amplitude.

Finally, the individual optimal value and the global optimal value are updated according to the calculation results of SCPSO, and then the algorithm enters the next iteration.

4. Results and Discussions

Based on the above model research and algorithm designs, Python programming is used for simulation experiments on our PC (Intel(R) Core(TM) i7-10700K CPU @ 3.80 GHz, 32 GB RAM, NVIDIA GeForce RTX 3090). For the purpose of proving the effectiveness of the presented system and method, first the GA and the AGA are used to evaluate the mission assignment. Then, the evaluation experiments of three-dimensional path planning are carried out, in which PSO and SCPSO are compared. Third, the simulation experiment of a multi-drone disaster rescue system is carried out by combining mission assignment and path planning.

4.1. Comparison of GA and AGA

In this section, GA and AGA algorithms are compared. Mission points are used for gene coding on chromosomes. The initial population size N_AGA of the algorithm is 200, the maximum crossover rate P_C,max, and minimum mutation rate P_C,min are 0.8 and 0.4, and the maximum crossover rate P_M,max, and minimum mutation rate P_M,min are 0.01 and 0.001. The experiment of multi-drone mission assignment is carried out in a 100 × 100 × 100 area. The drone base station coordinate is (2, 35, 1).

Table 1. Parameters of each threat source.

No.	Center Coordinate, Slope, k1, and Height of Mountain Threat Source	Center Coordinate, k2, and Radius of Transmission Tower Threat Source	Center Coordinate, k3, $k_{N\_w}^w$, rN_w, and $R_{N\_w}^w$ of Severe Weather Threat Source
1	(85, 85), (8, 8), 10, 40	(40, 40), 20, 18	(70, 50), 10, 0.6, 15, 25
2	(20, 80), (10, 12), 10, 60	(70, 20), 20, 12	(15, 45), 10, 0.4, 10, 20
3	(30, 20), (7, 7), 10, 50	(55, 80), 20, 15	(90, 20), 10, 0.2, 8, 18

shows the parameters of each threat source. The experiment of three drones performing 10 mission points is evaluated and tested. The minimum fitness function and the mean fitness function of each generation with the maximum iteration time T_AGA of 100, 200, and 300 are studied by performing 50 mission assignment experiments. Mission execution order, fitness function, and program running time are the main indexes of the mission assignment evaluation experiment. The fitness function is defined in Equation (24), where it is evident that a lower fitness function value indicates a better calculation outcome for the algorithm.

The maximum displacement of the drone is 300, the body weight is 1, and the maximum load is 4. Mission point coordinates and material requirements are shown in

Table 2. Coordinates and requirements of 10 mission points.

No.	Coordinate Mission Point	Requirement of Mission Point	No.	Coordinate Mission Point	Requirement of Mission Point
1	(5, 74, 10)	0.9	6	(94, 9, 7)	0.7
2	(21, 69, 35)	0.7	7	(58, 30, 1)	0.9
3	(39, 83, 4)	1.3	8	(29, 14, 37)	0.5
4	(90, 53, 1)	0.8	9	(95, 80, 10)	1.1
5	(4, 91, 3)	0.4	10	(80, 82, 34)	0.6

.

Table 3. Comparison of performance evaluation indexes between GA and AGA (10 mission points).

Fitness Function (TAGA = 100)
Method	OV	WV	MV	SDV
GA	711.7506	763.4414	734.2795	13.9377
AGA	677.0911	703.6567	692.9914	7.1411
Program running time (s) (TAGA = 100)
Method	OV	WV	MV	SDV
GA	10.8673	10.9425	10.8894	0.0155
AGA	10.9270	10.9460	10.9330	0.0038
Fitness function (TAGA = 200)
Method	OV	WV	MV	SDV
GA	699.2446	752.1145	725.0800	11.8666
AGA	677.0911	703.0757	691.9700	7.1607
Program running time (s) (TAGA = 200)
Method	OV	WV	MV	SDV
GA	21.6633	21.7874	21.7073	0.0225
AGA	21.7839	21.8891	21.8036	0.0147
Fitness function (TAGA = 300)
Method	OV	WV	MV	SDV
GA	694.5612	745.4266	719.5153	13.8968
AGA	677.0911	698.4715	689.9750	6.1316
Program running time (s) (TAGA = 300)
Method	OV	WV	MV	SDV
GA	32.4697	32.6227	32.5286	0.0322
AGA	32.6415	32.6807	32.6671	0.0090

presents the comparison results of performance evaluation indexes between GA and AGA. The results indicate that the optimal value (OV), the worst value (WV), and the mean value (MV) of AGA are lower than those of GA except for the running time. Additionally, the lower standard deviation value (SDV) proves that AGA can obtain the optimal assignment stably. Figure 7 shows the statistical results of the minimum fitness per generation and the mean fitness per generation of GA and AGA with 100, 200, and 300 iterations, respectively. It can be seen from Figure 7 that the convergence speed and optimal value of AGA are much better than those of GA.

Table 4. Multi-drone optimal mission assignment of GA and AGA (10 mission points).

Method	No. of Drone	TAGA = 100	TAGA = 200	TAGA = 300
GA	1	5, 1	5, 1	3, 5, 1
	2	3, 4, 9, 10	3, 7, 6, 8	7, 6, 8
	3	7, 6, 8, 2	4, 9, 10, 2	2, 10, 9, 4
AGA	1	3, 5, 1	3, 5, 1	3, 5, 1
	2	7, 6, 8	7, 6, 8	7, 6, 8
	3	4, 9, 10, 2	4, 9, 10, 2	4, 9, 10, 2

gives out the optimal mission execution order comparison result: when the iteration time is small, there are significant differences in the optimal allocation results between GA and AGA. However, when the iteration time is large, the optimal allocation results of GA and AGA become similar. AGA can achieve better results with fewer iterations.

4.2. The Evaluation Experiment of SCPSO

In this section, the evaluation experiment of the 3D path planning algorithm is carried out, and the PSO and SCPSO algorithms are compared. The initial population size (N_SC) of the algorithm is 50. The number of path-control points (N_C) is 5. The maximum inertia weights w_max and minimum inertia weights w_min are 0.9 and 0.4, respectively. The maximum individual acceleration coefficients c_1max and minimum individual acceleration coefficients c_1min are 2.5 and 1, respectively. Additionally, the maximum sine and cosine amplitudes R_1max and minimum sine and cosine amplitudes R_1min are 2.25 and 1. The experiment of multi-drone path planning is carried out in a 100 × 100 × 100 area, and the simulation time is 50. The starting point coordinate of the drone is (2, 35, 1), and the end point coordinate is (95, 80, 10). There are three threat sources: mountains, transmission towers, and severe weather. The parameters of each threat source are shown in Table 1.

Table 5. Comparison of performance evaluation indexes between PSO and SCPSO.

Fitness Function (TSC = 100)
Method	OV	WV	MV	SDV
PSO	117.6105	146.7281	129.5862	5.0713
SCPSO	108.7102	132.6111	115.7923	6.7508
Program running time (s) (TSC = 100)
Method	OV	WV	MV	SDV
PSO	51.1639	64.1155	57.9294	3.2033
SCPSO	61.5123	66.3144	64.5302	1.2511
Fitness function (TSC = 200)
Method	OV	WV	MV	SDV
PSO	116.1025	137.6843	125.6526	4.5177
SCPSO	108.6614	124.6490	114.0725	5.0003
Program running time (s) (TSC = 200)
Method	OV	WV	MV	SDV
PSO	106.0202	124.9831	116.0116	5.0245
SCPSO	122.7706	134.3059	129.5099	2.1687
Fitness function (TSC = 300)
Method	OV	WV	MV	SDV
PSO	114.6301	132.5005	122.0011	4.6042
SCPSO	108.6507	121.9915	113.1423	4.8832
Program running time (s) (TSC = 300)
Method	OV	WV	MV	SDV
PSO	147.9369	186.1021	170.9336	8.1782
SCPSO	192.3480	204.1927	198.4601	2.7036

presents the comparison results of performance evaluation indexes between PSO and SCPSO. Figure 8, Figure 9 and Figure 10 show the 3D path planning examples and statistical results of the minimum fitness and the mean fitness of each generation of PSO and SCPSO when the maximum iteration times T_SC are 100, 200, and 300. The fitness function and program running time are the main results of this experiment. In 3D path planning, the fitness function refers to the distance of drone flight. Therefore, the smaller the fitness function value, the better the calculation effect of the algorithm. All the aforementioned evaluation experiment results demonstrate that SCPSO is capable of achieving a shorter flight path.

4.3. The Disaster Relief Simulation Experiment

This section conducts a multi-drone disaster rescue simulation evaluation experiment. Four algorithms, GA + PSO, GA + SCPSO, AGA + PSO, and AGA + SCPSO, are used for testing. The comprehensive evaluation of the multi-drone disaster rescue simulation evaluation experiment is carried out in a 100 × 100 × 100 area. The experiment of three drones performing 10 missions is evaluated and tested. The drone base station coordinate is (2, 35, 1), the number of path-control points N_C is 2, the iteration times of GA and AGA are 300, the iteration times of PSO and SCPSO are 100, and the time of simulations is 20. The mission point, threat source, and other parameters of the algorithm are the same as those in Section 4.1 and Section 4.2.

Table 6. Multi-drone optimal mission results (the number of mission points is 10 and T_AGA = 300).

Method	No. of Drone	Optimal Mission Assignment Result
GA	1	3, 5, 1
	2	7, 6, 8
	3	2, 10, 9, 4
AGA	1	3, 5, 1
	2	7, 6, 8
	3	4, 9, 10, 2

shows the multi-drone optimal mission results of GA and AGA when the maximum iteration time is 300.

Table 7. The fitness function and path length statistics of four methods.

Method	Optimal Fitness Function Value	Path Length
		OV	WV	MV	SDV
GA + PSO	694.5612	694.7434	739.5507	716.5673	11.2946
GA + SCPSO	694.5612	681.5401	737.6858	711.0448	14.5077
AGA + PSO	677.0911	691.8288	739.5353	716.0973	12.7313
AGA + SCPSO	677.0911	680.2483	737.5755	710.2630	14.7289

shows the evaluation results of four methods. It can be seen from Table 7 that the performance of AGA + SCPSO is the best, i.e., its fitness function is optimal and the path is the shortest. Figure 11 shows the visualization results of the drone disaster rescue experiment. It can be seen that the path planned by SCPSO is better, that is, its path is shorter and the route is smoother.

4.4. Discussions

After the disaster, due to the barrier of mountains and rivers, it is a challenging task for relief forces to reach the disaster area for the first time. Due to the urgency of disaster rescue, prompt execution of search and rescue operations within a short timeframe is often critical. The drone has the advantages of search and rescue [46], efficient monitoring [47], material distribution, and recovery of communication [48], which can be used for a fast disaster response [49]. Compared with a single drone, a multi-drone can undertake different rescue missions simultaneously and reduce rescue time. Therefore, it is necessary to study how to reasonably allocate different mission points for drones and plan their flight paths before implementing rescue. To this end, a multi-drone mission planning system and method that comprehensively consider factors such as flight distance, drone performance, threat source factors, and disaster site requirements are proposed in this paper. The research results can provide theoretical support for future disaster relief.

In this paper, three common threat sources in disaster rescue environments are modeled: mountains, transmission towers, and severe weather. Due to the characteristics of mountains, the mountain threat source is simplified by the bi-GMM. In future research, neural networks based on the truncated sign distance function (TSDF) can be used to model real mountain areas [50], making the model closer to real 3D mountain scenes. Although many anti-interference technologies have been developed to reduce the electromagnetic interference of drones [51], the signal is susceptible to weather, terrain, and other factors; so it is best to avoid entering the corresponding influence ranges. In the future, the interference of transmission towers and mining areas can also be studied. Modeling severe weather is a complex task due to its complexity and uncertainty and the vulnerability of drones to changes in airflow and pressure. In this paper, to enable fast computation, the range of influence of the severe weather threat source is abstracted using a cylinder. In the future, the difference in the response of drones’ own characteristics to weather phenomena and the temperature model [52] can be considered.

A cost–revenue function is formulated to facilitate multi-drone disaster relief mission allocation. The cost function represents the cost and risk of the multi-drone rescue process. In addition, due to the influence of weight on overcoming the work conducted by gravity, the order of mission execution is different, and the energy consumption of the drone is also varied. Therefore, the elevation cost is taken into account in the cost function. The revenue function represents the relief value of the task completed by the drone. The cost–revenue function represents the difference between the cost and revenue functions. Therefore, a smaller cost–revenue function indicates a better calculation effect for our model. In this study, the flight distance, flight height, threat source factors, drone performance constraints, and mission point revenue are all considered in the cost–revenue function. In the future, additional aspects can be taken into consideration when designing the cost–revenue function, such as the time required for a drone to perform different missions or the smoothness of the drone’s flight trajectory.

In this paper, AGA and SCPSO are used to assign missions and plan paths for drones, respectively. In order to further verify the effectiveness of AGA, we further supplement the comparative experiments of GA and AGA. In the 100 × 100 × 100 area, the experiment of 5 drones performing 20 missions is evaluated. The maximum displacement of the drone is 400, the body weight is 1, and the maximum load is 5. Mission point coordinates and material requirements are shown in

Table 8. The coordinates and requirements of 20 mission points.

No.	Coordinate Mission Point	Requirement of Mission Point	No.	Coordinate Mission Point	Requirement of Mission Point
1	(5, 74, 10)	0.9	11	(10, 11, 11)	0.6
2	(21, 69, 35)	0.7	12	(40, 65, 2)	1.0
3	(39, 83, 4)	1.3	13	(56, 7, 6)	0.7
4	(90, 53, 1)	0.8	14	(64, 99, 8)	0.3
5	(4, 91, 3)	0.4	15	(28, 10, 11)	0.2
6	(94, 9, 7)	0.7	16	(18, 32, 4)	0.5
7	(58, 30, 1)	0.9	17	(5, 59, 1)	0.4
8	(29, 14, 37)	0.5	18	(19, 80, 60)	0.2
9	(95, 80, 10)	1.1	19	(84, 31, 1)	1.1
10	(80, 82, 34)	0.6	20	(41, 5, 18)	1.2

. The remaining parameters are the same as those in Section 4.1. A total of 50 mission assignment experiments are performed.

Table 9. Comparison of performance evaluation indexes between GA and AGA (20 mission points).

Fitness Function (TAGA = 100)
Method	OV	WV	MV	SDV
GA	1300.4399	1436.9930	1367.9740	31.1004
AGA	1005.1311	1063.0131	1038.3372	13.4497
Program running time (s) (TAGA = 100)
Method	OV	WV	MV	SDV
GA	11.4349	11.4740	11.4482	0.0084
AGA	11.4968	11.5775	11.5298	0.0116
Fitness function (TAGA = 200)
Method	OV	WV	MV	SDV
GA	1293.7675	1414.6581	1347.7407	29.4725
AGA	1001.0166	1056.1460	1035.4154	13.5783
Program running time (s) (TAGA = 200)
Method	OV	WV	MV	SDV
GA	22.8013	22.8360	22.8165	0.0084
AGA	22.9302	23.0175	22.9803	0.0210
Fitness function (TAGA = 300)
Method	OV	WV	MV	SDV
GA	1227.6617	1383.7537	1335.1681	29.8367
AGA	987.6911	1049.3092	1027.6919	14.2750
Program running time (s) (TAGA = 300)
Method	OV	WV	MV	SDV
GA	34.1791	34.2579	34.1992	0.0152
AGA	34.3697	34.4926	34.4398	0.0322

presents a comparison of the performance evaluation indexes between GA and AGA for the corresponding disaster relief mission. It can be seen that in the case of increased task complexity, except for the program running time, the performance of AGA is significantly better than that of GA. Figure 12 shows the statistical results of the minimum fitness of each generation and the mean fitness of each generation for the corresponding GA and AGA when iteration times are 100, 200, and 300. It can be seen that the convergence speed, optimal value, and average value of AGA are much better than those of GA. Compared with the previous results (see Section 4.1), it can be seen that the higher the mission complexity, the better the performance of AGA. This means that the performance gap between AGA and GA is larger, and AGA is more suitable for mission assignments with high mission complexity.

Table 10. Multi-drone optimal mission assignment of GA and AGA (20 mission points).

Method	No. of Drone	TAGA = 100	TAGA = 200	TAGA = 300
GA	1	12, 2, 3, 14, 6, 13, 15	1, 12	7, 15, 11, 8
	2	11, 17	17, 6, 19, 4, 7	13, 20
	3	19, 7, 4, 20, 8	10, 9, 14, 20, 11, 8, 15	17, 3, 12
	4	16, 18, 10, 9	16	16, 6, 4, 19, 10, 9
	5	5, 1	13, 5, 3, 18, 2	5, 1, 14, 18, 2
AGA	1	11	9, 4, 19, 6, 7	7, 6, 19, 4, 12
	2	15, 8, 20, 13, 6, 19, 7	11, 15, 20, 13, 8	8, 20, 13, 15, 11
	3	16	16	16
	4	4, 9, 10, 14, 18, 2	12, 3, 14, 10, 18, 2	3, 14, 9, 10, 18, 2
	5	17, 1, 5, 3, 12	17, 5, 1	5, 1, 17

presents the corresponding optimal mission execution orders for GA and AGA: as task complexity increases, the results show significant differences.

After repeated experiments, it can be found that GA is difficult to obtain optimal mission assignment with high complexity; differently, the AGA proposed in this paper can obtain optimal allocation only with fewer iterations. In order to further verify the effectiveness of SCPSO, we further supplement the comparative experiments of AGA + PSO and AGA + SCPSO under the same simulation experiment of 5 drones performing 20 mission points. The algorithm parameters are the same as those in Section 4.3. A total of 20 simulation experiments are performed.

Table 11. The path length statistics of different methods.

Path Length
Method	OV	WV	MV	SDV
AGA + PSO	936.6954	1016.0094	972.2957	19.4389
AGA + SCPSO	920.5319	1012.6603	967.2762	27.2980

shows the evaluation indexes of AGA + PSO and AGA + SCPSO. It can be seen that the performance of AGA + SCPSO is better in environments with higher mission complexity, i.e., the path is shorter and the route is smoother. Compared with the previous results in Section 4.3, it can be seen that SCPSO performs better when the mission complexity increases, and SCPSO is more suitable for path planning problems with high mission complexity. Figure 13 shows the visualization results of the corresponding drone disaster rescue experiment.

The effectiveness of the proposed model is demonstrated in this paper through a series of experiments. The fitness function and processing time are used in the experiment. When the number of missions is small, AGA generally reaches the optimal value within 10 generations, while GA needs hundreds of generations to reach the optimal value, and its optimal value is not as good as the optimal value of AGA. When the number of missions increases, AGA generally reaches the optimal value within 100 generations, while GA will use even 300 generations to reach the optimal value. Although AGA uses more computing resources and takes more computational time, it obviously takes less time if the time to reach the optimal value is discussed instead of the same number of iterations. For PSO and SCPSO, since SCPSO balances the ability of global and local optimizations, its convergence speed and optimal value are better than PSO. Clearly, all of the aforementioned experimental results demonstrate that our proposed method exhibits strong performance.

The method proposed in this paper has at least three advantages. First, the proposed system is designed to address the real-world scenarios of multi-drone disaster relief. The typical threat sources of disaster rescue scenarios, the requirements of disaster areas, and the performance of drones are considered when designing the cost–revenue function. This study can provide a basis for future drone disaster rescue. Second, for the problem of multi-drone disaster rescue, an optimization method with higher calculation accuracy is proposed, and this method is more suitable for optimization problems with high complexity. Third, the model has good scalability. Our model can simulate common disaster rescue scenarios and can be applied to various disaster relief environments with minimal modifications. Despite its advantages, our method also has limitations. For instance, our approach is specifically tailored for scenarios where complete environmental information is available. The presence of partially unknown information in the environment can potentially lead to mission planning failures. Moreover, an excessive amount of information regarding drones and mission points may have a detrimental impact on the real-time performance of the algorithm. These could be addressed in future work.

5. Conclusions

This paper realizes the problem of multi-drone mission planning in a complex 3D environment. The original digital terrain of drone flight and three common threat sources are constructed in this paper, including mountains, transmission towers, and severe weather. When constructing the cost–revenue function, factors such as the performance of the drone, the requirements of each mission point, the elevation cost, and the threat source are considered. AGA is designed to solve the problem of multi-drone mission assignments. The improved circle algorithm, adaptive crossover rate and mutation rate, and a strategy that uses both roulette and elite retention methods are used to improve the efficiency of our method. The experimental results demonstrate that AGA exhibits stronger optimization abilities than GA. SCPSO is designed to plan the optimal path between adjacent mission points. The inertia and acceleration coefficients of linear weights are designed to maintain an effective balance between global exploration and local development, further enhancing the performance of SCPSO. The experimental results show that the SCPSO algorithm can plan an effective and safe path for drones. This study can provide a basis for future disaster rescue decisions.