3D UAV Trajectory Planning Based on Improved Whale Optimization Algorithm

1. Introduction

With the rapid rise of low-altitude economy, the application of Unmanned Aerial Vehicle (UAV) [1], as one of its core technologies, in low-altitude airway planning has received wide attention. With its autonomous and controllable, highly flexible, as well as safe and reliable characteristics, UAV has shown great potential in many fields, especially in low-altitude economic application scenarios such as urban mapping, logistics and transport, urban monitoring, crop monitoring, environmental protection, and emergency rescue. With the development of smart cities and intelligent transport, the role of UAVs is not only limited to traditional surveying and mapping and monitoring tasks, but also extends to scenarios that require high timeliness and accuracy, such as urban logistics and distribution, emergency medical supplies transport, disaster monitoring and rescue.

The booming development of low-altitude economy has driven the accelerated progress of UAV technology, especially the innovation and breakthroughs in intelligence, autonomy and collaboration. Among them, trajectory planning, as one of the core technologies of UAS, has a crucial impact on the efficiency and safety of mission execution. In order to adapt to the complex low-altitude environment and diversified mission requirements, trajectory planning not only needs to quickly calculate collision-free and threat-free paths, but also must fully consider multiple factors such as the smoothness of the flight path, mission execution cost, fuel consumption, and flight time. This puts forward extremely high technical requirements for the application of UAVs in complex environments, such as logistics and distribution in dense urban building areas and emergency medical rescue. The path planning of UAVs must not only avoid fixed terrain obstacles, but must also have the ability to cope with dynamic environments, such as weather changes and moving obstacles in flight.

Currently, in order to solve the UAV low altitude path planning problem, many heuristic search algorithms are widely used, such as the ant colony algorithm [2], which improves the UAV’s obstacle avoidance ability and path optimization efficiency in complex environments by introducing a dynamic pheromone updating mechanism, and the particle swarm algorithm [3], which, by combining the particle swarm algorithm with simulated annealing algorithms, proposes a hybrid optimization method that significantly improves the path planning accuracy and convergence speed of UAVs in dynamic environments, Grey Wolf Algorithm [4]Some studies have adapted the predation and tracking strategies in the Grey Wolf Algorithm to make it more adapted to the complex 3D environments of UAVs, A* Algorithm[5,6], Whale Algorithm[7], and Dijkstra’s Algorithm [8], among others. Recently, many scholars have been making great efforts to solve this problem under special circumstances. Zhang etc. [9] incorporated an improved slap swarm algorithm to establish a rotary wing UAV trajectory plan model, which considered UAV patrol efficiency, trajectory cost and power consumption cost to overcome the shortcomings of low algorithmic search efficiency and unsmooth paths in open-air warehouse storage environments. Cherif etc. [10] formulate the cargo-UAV mission as a multi-objective problem to ensure a safe delivery mission, which is aiming to minimize energy consumption, reduce handoff events and guarantee cellular reliability along the trajectory. Yuan etc. [11] employed enhanced particle swarm optimization to optimize the UAV path planning with terrain constraints, which interpolated new constraint-altitude waypoints above the terrain based on altitude information and terrain resolution. While these algorithms try to find optimal or near-optimal paths for UAVs in complex environments by simulating group behaviour in nature, path search mechanisms, etc. However, these traditional algorithms often face challenges in practical applications, especially when dealing with complex 3D environments and dynamically changing obstacles, they are vulnerable to problems such as initial conditions, algorithmic parameter settings, and local optimal solutions, which affect the efficiency and accuracy of path planning.

To cope with this problem, this paper proposes an improved whale optimisation algorithm, which introduces a reverse learning mechanism and a nonlinear convergence factor to enhance the global search capability of the algorithm and avoid falling into local optimal solutions. Through experimental verification, the improved whale optimisation algorithm shows significant improvement in convergence accuracy and computational efficiency, with 22.1% improvement in convergence accuracy and 74.1% reduction in standard deviation compared with the traditional algorithm. The algorithm can better adapt to the complex three-dimensional environment and effectively improve the overall performance of UAVs in low-altitude route planning, providing more reliable and efficient technical support for the application of low-altitude economy. With the continuous innovation of technology, UAVs will play a more critical role in many low-altitude economic scenarios in the future, such as smart cities, precision agriculture, and emergency response. And efficient and intelligent route planning algorithms will be one of the key drivers of UAV technology breakthroughs.

2. Modelling

2.1. Environmental Modelling

In order to visualise the results and the realism of the simulation environment, this paper selects the study to compare and verify the algorithm in a 3D environment. The 3D space is gridded and numerically coded to simulate the mountainous terrain environment, while the simulation object also includes the no-fly zone. In this paper, the 3D coordinates are planned as $n_x$*$n_y$*$n_z$ parts, and all the height information is represented by matrix H, with $z_{ij}$ represents the coordinates $\left(i,j\right)$ of the three-dimensional height, in this study $n=100$, then the simulated terrain map is shown below:

Figure 1. The simulation generates mountainous terrain maps, and each simulation generates new maps to verify that the algorithm works for different environments.

Figure 1. The simulation generates mountainous terrain maps, and each simulation generates new maps to verify that the algorithm works for different environments.

The specific expression for the matrix H is given below:

$$H=\left[\begin{array}{cc}\begin{array}{cc}{z}_{11}& z_{12}\end{array}& \begin{array}{cc}\dots & z_{1 n}\end{array}\\ \begin{array}{cc}{z}_{21}& z_{22}\end{array}& \begin{array}{cc}\dots & z_{2 n}\end{array}\\ \begin{array}{c}\begin{array}{cc}⋮ & ⋮\end{array}\\ \begin{array}{cc}{z}_{n 1}& z_{n 2}\end{array}\end{array}& \begin{array}{c}\begin{array}{cc}\ddots & ⋮\end{array}\\ \begin{array}{cc}\dots & z_{n n}\end{array}\end{array}\end{array}\right]$$

(1)

In the real environment, UAV flights are subject to many restrictions, such as crowded areas, military control areas and radar areas, etc. In this paper, these restricted areas are collectively called no-fly zones. In this paper, these restricted areas are collectively referred to as no-fly zones, which need to be avoided during UAV trajectory planning, and are set as cylindrical models in contact with the terrain surface model, and the simulation of no-fly zones is shown below.

Figure 2. Add a simulated no-fly zone for cylinders to the simulated terrain map, where cylinders are not allowed to fly over at infinite heights.

Figure 2. Add a simulated no-fly zone for cylinders to the simulated terrain map, where cylinders are not allowed to fly over at infinite heights.

Simulation of all cylinder model information is represented by a matrix as shown in Equation 2.

$$O=\left[\begin{array}{ccc}{x}_1& y_1& r_1\\ \begin{array}{c}{x}_2\\ \begin{array}{c}⋮\\ x_k\end{array}\end{array}& \begin{array}{c}{y}_2\\ \begin{array}{c}⋮\\ y_k\end{array}\end{array}& \begin{array}{c}{r}_2\\ \begin{array}{c}⋮\\ r_k\end{array}\end{array}\end{array}\right]$$

(2)

In the formula: $\left(x_k,y_k\right)$, $r_k$ are the coordinates of the centre and radius of the$k$ and the radius of the centre of the threat zone, respectively. From the security point of view, the altitude of the threat zone is set to be infinite, i.e., no flight is allowed over the threat zone.

2.2. Consideration Functions

In order to find the optimal trajectory, it is also necessary to take into account the influence of different factors, and this paper uses a cost function to evaluate the influence of different elements.

In this paper, the total cost of the research design model is the sum of terrain constraint cost, boundary constraint cost, no-fly zone constraint cost, flight path cost, flight altitude cost, flight altitude change cost, and flight path smoothing cost. The objective is to minimize the cost of the solution

2.2.1. Terrain Constraints

In order to avoid collision with the terrain during UAV flight, the height of each point on each trajectory should be greater than the terrain height. In the experiment if it is lower than the terrain height, a positive value is generated indicating the degree of violation of the terrain constraint, otherwise it is 0.

$$R_1={\sum }_1^{N}max\left(H-N+P,0\right)$$

(3)

where $R_1$ represents the degree of violation of terrain constraints, H represents terrain altitude information, N represents trajectory altitude information, P represents a safety threshold to ensure that the altitude of the path point will not be lower than the altitude of the terrain, the actual value is adjusted according to the implementation of the task and the flight environment, usually set a few metres to several tens of metres, in this paper, the value of 1.

2.2.2. Boundary Constraints

In order to avoid the planned UAV trajectory from going beyond the flight region, the path algorithm is investigated with the requirement that each point on the trajectory does not exceed the bounds of the specified spatial coordinates. In the experiments if the boundaries are exceeded, a positive value is generated indicating the degree of violation of the boundary constraints, otherwise it is 0.

$$R_2={\sum }_{i=1}^M{\sum }_{j=1}^{2}max\left(2-C\left(i,j\right),0\right)$$

(4)

$$R_3={\sum }_{i=1}^M{\sum }_{j=1}^{2}max\left(C\left(i,j\right)-T\left(j\right)+2,0\right)$$

(5)

where $R_2$ denotes the total extent to which the path point exceeds the left and lower boundaries of the map; $R_3$ denotes the total extent to which the path points exceed the right and upper boundaries of the map;$C$ denotes the matrix of the 3D coordinates of the total path points after interpolation; $M$ denotes the total number of path points after interpolation;$i$ Table S shows the first $i$ path point, where$i$ has two coordinate values of x and y; when $j=1$ when $i$ x coordinate value is taken, when $j=2$ when $i$ takes the y coordinate value; $T$ Is a vector containing the extent of the map.

2.2.3. No-Fly Zone Constraint Costs

In order to avoid UAVs from colliding with buildings and staying away from restricted areas during flight, all no-fly zones are set as an infinite-height impenetrable cylinder model during the research process in this paper. The distance from each point on the trajectory to the centre of the no-fly zone circle is calculated and then this distance is subtracted from the radius. If this difference is less than zero, it is replaced by zero to ensure that there are no negative values.

$$R_4={\sum }_{i=1}^{N_f}{\sum }_{j=1}^{N_p}\max \left(R_i-L_{ij},0\right)$$

(6)

where $R_4$ denotes the degree to which the indicated path point has entered the no-fly zone, and $N_f$ denotes the number of no-fly zones, and $N_p$ denotes the track point, and$R_i$ denotes the radius of the it no-fly zone, and $L_{ij}$ denotes the distance from the point on the track to the centre of the no-fly zone circle.

2.2.4. Cost of Flight Trajectory

In order to reduce the distance travelled, the shortest paths are sought during trajectory planning in order to save energy. This is done by calculating the square root of the sum of the squares of the differences in the coordinates of two neighboring trajectory points.

$$R_5={\sum }_{i=1}^{N-1}{D}_i$$

(7)

$$D_i=\sqrt{{\left(x_{i+1}-x_i\right)}^2+{\left(y_{i+1}-y_i\right)}^2+{\left(z_{i+1}-z_i\right)}^2}$$

(8)

where $R_5$ denotes the degree of flight distance, the$D_i$ is the first$i$ and$i+1$ distance between the path points.

2.2.5. Cost of Flight Altitude Change

In order to ensure the safety, stability and energy efficiency of the UAV, the trajectory should be planned with as little altitude fluctuation as possible. This is achieved by calculating the variance of the altitude of the trajectory points during design.

$$R_6={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{\left(z_i-\overline{z}\right)}^2$$

(9)

where $R_6$ denotes the degree of UAV altitude fluctuation.

2.2.6. Flight Altitude Cost

In order to maintain the advantages of UAVs flying at low altitude, such as good manoeuvrability, stealth, clarity of shots, low energy consumption, etc., the trajectory planning should try to ensure a low altitude, which is achieved by calculating the average value of the altitude of the trajectory points during the design.

$$R_7={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{z}_i$$

(10)

where $R_7$ denotes the degree of UAV flight altitude, and$i$ is the number of path points.

2.2.7. Cost of Flight Path Smoothing

In order to ensure the safety of the UAV and its surroundings during flight, and the feasibility of the operation and its own structure, it is required that the path should be kept smooth during trajectory planning.

$${v1}_i=P_{i+1}-P_i$$

(11)

$${v2}_i=P_{i+2}-P_{i+1}$$

(12)

$${\delta }_i={\cos }^{-1}\left({\displaystyle \frac{{v1}_i·{v2}_i}{‖{v1}_i‖‖{v2}_i‖}}\right)$$

(13)

$$R_8={\sum }_{i=2}^n{\delta }_i$$

(14)

where $R_7$ denotes the degree of angle of change in all directions, which is used to measure the smoothness of the path.${v1}_i$ is the same as ${v2}_i$ represents the vector difference between two neighbouring points in the same direction as the UAV flight direction, and ${\delta }_i$ represents the difference between ${v1}_i$ is the inverse cosine of the cosine of the angle between ${v2}_i$ the inverse cosine of the cosine of the angle.

3. Whale Optimization Algorithm

Whale Optimization Algorithm (WOA) is a meta-heuristic intelligent optimization algorithm based on the foraging behavior of whales, which was proposed by Seyedali Mirjalili and Andrew Lewis in 2016 [12]. The mathematical model is based on the three foraging behaviors of humpback whales: encircling predation, bubble attack and food seeking. The flow chart of the calculation is shown below:

Figure 3. Whale optimization algorithm flowchart, the process reflects the stages of initialisation, encircling prey, bubble net predation method, searching for prey.

Figure 3. Whale optimization algorithm flowchart, the process reflects the stages of initialisation, encircling prey, bubble net predation method, searching for prey.

3.1. Surrounding the Prey

The locally optimal solution for encircling the prey is the one closest to the food, and the other individuals gradually approach the locally optimal individual to encircle the food. The mathematical model for this stage is as follows:

$$D=\left|C·x_{\ast }-x_i^{t-1}\right|$$

(15)

$$x_i^t=x_{\ast }-A·D$$

(16)

Eq:$x_{\ast }$ is the current optimal solution and $x_i^t$ denotes the position of the ith whale in generation t. A and C are coefficient vectors with the expression:

$$A=2a·r-a$$

(17)

$$C=2·r$$

(18)

where: r is a random number uniformly distributed on [0, 1] and $\alpha$ is the convergence factor, whose expression is:

$$a=2-{\displaystyle \frac{2t}{T_{max}}}$$

(19)

where: t is the current number of iterations, and ${\mathrm{T}}_{max}$ is the maximum number of iterations.

3.2. Subsection Bubble Attacks

There are two attack mechanisms that occur during bubble net predation, the constrictive encirclement mechanism and the spiral updating position. The constriction-enclosure mechanism is encircling predation, which is achieved by lowering the value of the$\alpha$ value. Now assume that the probability of both is 0.5. The mathematical model for this stage is as follows:

$$x_i^t=x_{\ast }-A·\left|C·x_{\ast }-x_i^{t-1}\right|\left(p<0.5\right)$$

(20)

$$x_i^t=x_{\ast }-A·\left|C·x_{\ast }-x_i^{t-1}\right|·e^{bl}·\cos \left(2\pi l\right)\left(p\ge 0.5\right)$$

(21)

Eq: $p$ is a random number on [0, 1] obeying a uniform distribution, and $b$ is a constant coefficient, and$l$ is a random number between [-1, 1].

3.3. Random Search

Referring to the whale at this point $A$ has a value greater than 1 or less than -1, the search agent’s position in the exploration phase will be updated here based on the randomly selected search agent, not on the optimal search agent so far [13].

$$X_{t+1}=x_{rand}-A·D_{rand}$$

(22)

$$D_{rand}=\left|C·x_{rand}-x_t\right|$$

(23)

where $x_{rand}$ denotes the coordinates of any whale.

4. Improving the Whale Optimization Algorithm

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

The WOA algorithm has a simple structure with only a few parameters to be adjusted and is suitable for continuous and discrete optimization problems. However, on some complex high-dimensional problems, especially in solving spatially complex and multi-peaked problems, the convergence of WOA [14] may be slower, and it still falls into the problem of locally optimal solutions even through a larger number of iterations. In this paper, an Improved Whale Optimization Algorithm (IWOA) is proposed. The IWOA algorithm uses reverse learning to increase the diversity of the initial population of the algorithm, and at the same time, the convergence factor is changed from linearly decreasing to nonlinearly decreasing, and a random number generation mechanism is integrated into the iterative process, which coordinates the global search capability and local development capability of the algorithm. The global search ability and local development ability of the algorithm are coordinated [15]. The flowchart of the IWOA algorithm is shown below:

Figure 4. Improved whale optimization algorithm flowchart, the process compared with the original algorithm in the initial stage with the reverse learning mechanism to expand the search range, the late addition of non-linear convergence factors to improve the ability to solve the local search.

Figure 4. Improved whale optimization algorithm flowchart, the process compared with the original algorithm in the initial stage with the reverse learning mechanism to expand the search range, the late addition of non-linear convergence factors to improve the ability to solve the local search.

4.1. Reverse Learning Initialization

Inverse learning is a commonly used global optimization strategy that is currently proven and applied in many swarm intelligence algorithms [16]. The idea of reverse learning is to use the current solution and its inverse solution to improve the search efficiency. In the study, it is assumed that the dimension of the search space is D, the population size is N, the upper and lower bounds of the search space are respectively $a$, the $b$, the initial population $X$ of each individual is $x_{ij}$ and the inverse solution of the population $X\mathrm{\text{'}}$ Each individual of the initial population is $x_{ij}^{\mathrm{\text{'}}}$ and the formula is as follows:

$$x_{ij}^{\mathrm{\text{'}}}=l\left(a_j+b_j-x_{ij}\right)$$

(24)

$$l=\left(a_2-1\right)·r_{rand}+1\left(-2<r_{rand}<-1\right)$$

(25)

where $l$ is the backward learning factor, is $\left[\mathrm{0,2}\right]$ the random number on the.

$$X=SelectBest\left(X,X^{\prime }\right)$$

(26)

At this point $X$ is the population after calculating the fitness update.

4.2. Nonlinear Convergence Factors

In WOA $a$ is the convergence factor [17], which decreases linearly from 2 to 0, as formulated in equation (19). At $a$ approaching 2 due to the complexity of the space the step size of its searching agent decreases rapidly as a result, which makes the exploration range of the global searching phase shrink rapidly, resulting in missing some potential globally optimal regions. The convergence of the $a$ the convergence of the algorithm from global search to local search can also be a premature shift from global search to local search resulting in local optimums. The algorithm converges to a local optimum at$a$ approaching 0, the linearly decreasing convergence property slows down the change of its own value, which makes the step size of the local search not obvious enough, leading to the weakening of the ability to search in the vicinity of the local optimal solution. Aiming at the above problems this paper proposes a nonlinear convergence factor with the following specific formula:

$$a_{nonlinear}=\alpha +{\displaystyle \frac{1-\beta ·\left(t/T_{max}\right)}{1-\gamma ·\left(t/T_{max}\right)}}$$

(27)

In order to ensure that the $t=0$ time $a_{nonlinear}$ is close to 2, at $t=T_{max}$ time $a_{nonlinear}$ is close to 0, the design $\alpha =1$ The $\beta =0.5$.$\gamma =50$. The resulting smaller value of t in the early stages when the need to decline is slow, with little change in the numerator denominator, makes the $a_{nonlinear}$ less variation in the early stage and the global search capability is improved. In the middle stage the algorithm’s ability to balance between global and local search is improved because the convergence factor varies nonlinearly. In the late stage when fast descent is required, a larger value of t makes the $a_{nonlinear}$ rapid descent in the later stages and enhances the local exploitation capability.

4.3. Random Number Generation Mechanism

In population intelligence algorithms, it is common to generate random numbers to improve the algorithm’s global and local search capabilities [18]. In this paper, the random number generation mechanism is used to optimize the algorithm in individual initialization, position update. During individual initialization, the random number$P$ is used to determine the prey action to be selected when the individual is updated. For position update, the random numbers$r_1$ is used to determine the predator action to be selected during the individual update,$r_2$ parameter$A$, $C$ are shown in Eqs. (17) and (18), and the index of the individual is selected by the random number. For the calculation of the nonlinear convergence factor, the random number $r_{rand}$ is guaranteed to $l$ has randomness as shown in equation (24), so that it can produce different spiral trajectories and avoid the algorithm to fall into local optimum.

5. Simulation Verification and Result Analysis

In this study, the 3D space is planned as a 100*100*100 grid and 30 sittings of peaks are randomly generated as shown below:

Figure 5. Terrain model.

Figure 5. Terrain model.

No-fly zones are also planned in space, which are non-flyable cylinders with base coordinates and radii as shown in Table 1:

5.1. Simulation Results

The whale optimization algorithm was modified as described above to study the UAV trajectory. Both algorithms assume a starting point at (35, 2) and an end point at (75, 80), with both heights varying according to randomly generated terrain. The calculation is set three times, the population size is set to 80, and each calculation undergoes 500 iterations. The performance information of the two algorithms is shown in Table 2, the simulation results are shown in Figure 6, and the convergence results are shown in Figure 7.

5.2. Analysis of Results

Firstly, as shown in Figure 4, the simulation results of the two algorithms in a unified environment, when solving the 3D trajectory planning problem, both algorithms are able to plan a smooth trajectory between the start and the end points and avoid collision with the mountain and the threatening area. Secondly, from Figure 6(a), we can intuitively see that the IWOA algorithm chooses a shorter route to bypass the mountain and plan a shorter trajectory. And from Figur. 6(b), we can intuitively see that the IWOA algorithm plans a smoother trajectory when avoiding the mountain and the no-fly zone, which can effectively improve the safety of the UAV’s flight operations. This is due to the fact that the WOA algorithm enters the local search too early when solving the trajectory and falls into the local optimum. Finally, Table 2 shows that the optimal value of the IWOA algorithm increases by about 0.882%, the worst value decreases by about 54.8%, the mean value decreases by about 22.1%, and the standard deviation decreases by about 74.1% compared with the WOA algorithm through five calculations.

Comprehensive analysis of the above shows that the algorithm can be effectively optimized by incorporating the improved methods of inverse learning, non-linear convergence factor, random number generation mechanism, which achieves the enhancement of global and local search capabilities IWOA algorithm can plan a smaller total cost, shorter total path and smoother trajectory compared to the WOA algorithm when solving such problems as three-dimensional trajectory planning.

Four common optimization algorithms were selected LO for Lemur optimization algorithm [19] SWO is the spider wasp optimization algorithm COA is the coyote optimization algorithm [20] KOA is the Keplerian Optimization Algorithm [21] to verify the performance of IWOA algorithm. In the comparison test, the populations of the IWO algorithm and WOA algorithm proposed in this paper and the remaining six algorithms are set to 50, and the maximum number of iterations are set to 100, and the test results are shown in Figure 8.

6. Results Conclusion

In this paper, based on the WOA algorithm, IWOA is proposed by incorporating inverse learning, adding nonlinear convergence factors, and random number generation mechanism. IWOA is able to plan UAV trajectories better than WOA in simulation experiments in compliance with the constraints, which indicates that the improved method solves the local optimum problem, improves the ability to deal with individual boundaries, and ensures the applicability of the algorithm on different problems. In future research, the introduction of more influencing factors to enrich the model will be considered on the one hand, and UAV trajectory planning in dynamic environment will also be considered on the other.