Three-Dimensional Unmanned Aerial Vehicle Path Planning in Simulated Rugged Mountainous Terrain Using Improved Enhanced Snake Optimizer (IESO)

1. Introduction

The complex mountainous terrain, characterized by significant elevation changes and limited transportation options, often renders certain areas difficult to access using conventional observation, surveillance, and logistical support methods. Due to their agility and flexibility, unmanned aerial vehicles (UAVs) have become extensively utilized in these challenging environments. Consequently, high-quality and rapid flight path planning in such complex terrains is critically important for enhancing the efficiency of UAV missions [1,2].

Recent research in UAV path planning has introduced many different methodologies, including meta-heuristic optimization algorithms (MAs) [3], the quality-oriented hybrid A* and Q-learning algorithm [4], the artificial potential field [5], trajectory planning with matrix alignment Dijkstra’s algorithm [6,7]. Flores-Caballero et al. [8] approached path planning as a constrained global optimization problem, utilizing MAs in conjunction with the A* algorithm. However, their objective function was overly simplistic, concentrating exclusively on flight distance, and consequently, it did not achieve optimal practical outcomes. Shao et al. [9] proposed a hierarchical trajectory optimization framework that combines an adaptive particle swarm optimization algorithm with the Gaussian pseudospectral method. This methodology accounted for the UAV dynamics model and state constraints, thereby enabling the minimization of flight time. However, the simulated environment was overly simplistic and not applicable to complex terrains, such as mountainous regions or ravines. For large-scale collaborative searches, Chen et al. [10] utilized the ant colony optimization algorithm to efficiently plan paths for UAV swarms, reducing task completion times. Despite its rapid optimization, the algorithm exhibited limited convergence precision.

UAV path planning necessitates considering rugged terrain environments, such as mountainous regions, flight safety, and UAV-specific conditions, rendering it a quintessential multi-objective optimization problem. A literature review reveals that current research predominantly focuses on MAs due to their global search capabilities, which effectively address UAV motion constraints and path planning requirements, particularly in complex environments characterized by rugged mountainous terrain.

In this study, we examined the Enhanced Snake Optimizer (ESO) proposed by Yao et al. [11], which builds upon the original Snake Optimizer (SO) developed by Hashim and Hussien [12]. The ESO incorporates a randomized population initialization and mirroring imaging. This enhancement has been validated using classical test and engineering design functions, demonstrating its applicability to UAV path planning. However, in line with the No Free Lunch theorem [13], which posits that no single algorithm is universally optimal for all optimization problems, we acknowledge the limitations of ESO.

We propose three major ESO revisions at three stages that significantly improve optimization results, particularly in rugged terrains. Therefore, we name our method Improved Enhanced Snake Optimizer (IESO). First, the original SO employs a simple normal distribution random function to initialize the "population," which consists of the initial guess points. We replaced this with a Chebyshev chaotic map to improve the diversity and quality of the initial population. Second, we introduced a non-monotonic coefficient to adjust the "temperature," which controls the search merge and split processes. This modification addresses the SO’s persistent "cold temperature" state, balancing global exploration and local exploitation. Third, we added a boundary value to address occasional anomalies from the dynamic opposition learning mechanism, a key enhancement of ESO over SO.

With these three major revisions to ESO, we tested IESO against five state-of-the-art optimization algorithms, namely Enhanced Snake Optimizer (ESO) [11], Snake Optimizer (SO) [12], Golden Jackal Optimization (GJO) [14], Honey Badger Algorithm (HBA) [15], Grey Wolf Optimizer (GWO) [16], demonstrating superior performance in precision, convergence speed and planned trajectory quality.

2. Methodologies: Improved Enhanced Snake Optimizer (IESO) Algorithm and Simulation Design

This section is organized into four subsections: 2.1 briefly reviews the Snake Optimizer (SO). 2.2 presents a concise overview of the Enhanced Snake Optimizer (ESO). 2.3 offers a detailed explanation of our Improved Enhanced Snake Optimizer (IESO). 2.4 discusses the design of the simulation experiments and the evaluation criteria.

2.1. Basic Concepts of Snake Optimizer (SO)

The Snake Optimizer (SO) was initially proposed by Hashim and Hussien [12] in 2022. It is an optimization algorithm inspired by the snake mating process, which is influenced by "temperature" (mode switch) and "food" (change position) availability. Snakes mate only when the temperature is low, and food is abundant. When food is scarce, the SO algorithm prioritizes the food search, representing the exploration phase (global search). Once food is sufficient, the algorithm progresses to the exploitation phase (local search). During "hot" temperatures (T > 0.6), mating is deferred. When temperatures are low (T < 0.6), snakes engage in mating, based on a random probability $P_r$ (0, 1), whether to mate with the opposite sex or engage in same-sex aggression. When $P_r>0.6$, mating behaviour occurs between sexes; while $P_r$ < 0.6, males compete for mating rights in a fight mode, followed by a mating mode. Successful mating results in females laying eggs and hatching new snakes, thereby renewing the population (add more search points).

In SO, the "population" (N) represents the search points. The initial positions are generated using a random uniform distribution function, as defined in Equation 1. The population is then equally divided into two groups: males and females.

$$X_i=X_{min}+ran{d}_i\times (X_{max}-X_{min})$$

(1)

In SO, The "temperature" (T) determines the switch between global and local search modes. It is defined by Equation 2, where t is the current iteration and T is the maximum iteration and the same for the rest of this manuscript.

$$Temp=ex{p}^{-\mathit{t}/T}$$

(2)

In SO, a "food" parameter represents the feature position. The algorithm performs a broad, global search when no food is available. If food is present and the temperature is "hot", the algorithm refines the search point positions by eating food. When the temperature is "cold," the algorithm enters the exploitation phase, where males and females select optimal mates, resulting in convergence and repositioning new search points. The food quantity (Q) is defined in Eq. 3.

$$Q=c_{1}ex{p}^{(t-T)/T},c_1=0.5$$

(3)

When Q < 0.25, the snakes search for food by selecting any random position and updating their position for it, which is called the exploration phase. The male version is presented as Eq. 4, and $A_m$ is defined as Eq. 5. The female version is identical except for replacing the subscript _m to _f.

$$X_{i,m}(t+1)=X_{rand,m}\left(t\right)\pm c_2\times A_m\times ((X_{max}-X_{min})\times rand+X_{min}),c_2=0.05$$

(4)

$$A_m=ex{p}^{(-f_{rand,m}/f_{i,m})}$$

(5)

When Q > 0.6 and T > 0.6, the snake searches for food, keeping the global search. And it is presented in Eq. 6.

$$X_{i,j}(t+1)=X_{food}\pm c_3\times Temp\times rand\times ((X_{food}-X_{i,j}\left(t\right)),c_3=2,$$

(6)

where $X_{i,j}$ is the $i^{t}h$ snake position (male or female), $X_{food}$ is the position of the food.

When Q > 0.6 and T < 0.6, the snake will enter the flight mode (FM) or mating mode (MM), which means the search points are merging and splitting, which are presented in Eq. 7 and 8.

$$X_{i,m}(t+1)=X_{i,m}\left(t\right)+c_3\times FM\times rand\times Q\times X_{best,f}-X_{i,m}\left(t\right)),c_3=2$$

(7)

$$X_{i,f}(t+1)=X_{i,f}(t+1)+c_3\times FF\times rand\times Q\times X_{best,m}-X_{i,f}(t+1),c_3=2,$$

(8)

where $X_{i,m}$ is the $i^{t}h$ male position, $X_{best,f}$ refers to the best female position. FM and FF are the fighting ability of a snake (male and female), and it is defined in Eq. 9 and 10.

$$\mathrm{FM}={exp}^{-f_{best,f}/f_i}$$

(9)

$$\mathrm{FF}={exp}^{-f_{best,m}/f_i}$$

(10)

The mating mode will create new search points and replace the worst current ones. The equations are given in Eq. 11 and 12. The mating abilities are given in Eq. 13 and 14.

$$X_{i,m}(t+1)=X_{i,m}\left(t\right)+c_3\times M_m\times rand\times Q\times X_{i,f}\left(t\right)-X_{i,m}\left(t\right),c_3=2$$

(11)

$$X_{i,f}(t+1)=X_{i,f}\left(t\right)+c_3\times M_f\times rand\times Q\times X_{i,m}\left(t\right)-X_{i,f}\left(t\right),c_3=2$$

(12)

$$M_m={exp}^{-f_{i,f}/f_{i,m}}$$

(13)

$$M_f={exp}^{-f_{i,m}{f}_{i,f}}$$

(14)

The replace functions are given in Eq. 15 and 16.

$$X_{worst,m}=X_{min}+rand\times X_{max}-X_{min}$$

(15)

$$X_{worst,f}=X_{min}+rand\times X_{max}-X_{min}$$

(16)

2.2. Enhanced Snake Optimizer (ESO)

The Enhanced Snake Optimizer (ESO) was proposed by Yao et al. in 2023 [11]. The main modifications are the dynamic coefficients. First, the food quantity initialization function opens the fixed coefficient $c_1$ = 0.5 to a dynamic updated $c_1^{new}$, given as Eq. 17.

$$c_1^{new}=c_1+\frac{1}{10}\times cos\left(r_1^4\times \frac{\pi }{2}\right),$$

(17)

where $r_1$ is a random number between (0,1).

Then, the fixed temperature coefficient $c_2$ = 0.05 to a new dynamic $c_2^{new}$ as shown in Eq. .

$$c_2^{new}=c_2+\frac{1}{1000}\times cos\left(r_2^4\times \frac{\pi }{2}\right),$$

(18)

where $r_2$ is a random number between (0,1). And, the fixed $c_3$ = 2.

Another major enhancement is the adoption of mirror imaging strategy [17] and opposition-based learning [18]. The upper boundary (UB) and lower boundary (LB) were added and are presented in Eq. 19 and 20.

$$\frac{(UB+LB)/2-X}{X^{\ast }-(UB+LB)/2}=\frac{H}{H^{\ast }},$$

(19)

where $H^{\ast }$ and $X^{\ast }$ are the mirror positions of H and X.

$$X^{\ast }=\frac{UB+LB}{2}+\frac{UB+LB}{2\delta }-\frac{X}{\delta },$$

(20)

where $\delta =\frac{H}{H^{\ast }}$.

2.3. Improved Enhanced Snake Optimizer (IESO)

In our Improved Enhanced Snake Optimizer (IESO), we have implemented three significant enhancements based on the Enhanced Snake Optimizer (ESO). First, we adopted Chebyshev Population Initialization to initialize the snake’s initial population. Second, we introduced a non-monotonic factor to modulate the temperature. Third, we applied a Dynamic Boundary-based Opposition Learning (DBOL) approach to define the upper and lower boundaries. The general flowchart is illustrated in Figure 1.

2.3.1. Chebyshev Population Initialization

The SO and ESO randomly populate the initial population, leading to unstable target accuracy and varying convergence speeds across multiple searches. Achieving a more evenly distributed initial population within the solution space enhances the probability of locating the optimal value. Chaotic search, renowned for its randomness, ergodicity, and uniqueness, is preferred over random search strategies for generating initial populations due to its thorough coverage of the solution space. Therefore, this paper employs the Chebyshev chaotic map [19] for initializing the snake population to ensure a more comprehensive occupation of the entire solution space. The mathematical model of the Chebyshev chaotic map is as follows Eq. 21:

$$X_i(t+1)=cos\left(t{cos}^{-1}\left(X_i\left(t\right)\right)\right)$$

(21)

This comparison, illustrated in Figure 2a,b, demonstrates that the Chebyshev chaotic map produces a more uniformly distributed initial population. The Chebyshev chaotic map exhibits enhanced search capabilities, particularly when solutions to constrained problems are situated at the boundaries of the search space.

2.3.2. Non-Monotone Temperature Factor

In SO and ESO, the temperature factor influences the behaviour of the snake swarm in a monotonic manner. Figure 2 demonstrates that the temperature factor solely depends on the number of iterations, with a threshold value of approximately 220. Before this threshold, the temperature coefficient T consistently remains above 0.6, indicating a "warm" state and foraging activity. Conversely, as the optimization process progresses into its later stages, the temperature shifts to a "cold" state and remains "cold." This transition indicates that after a significant number of iterations, both SO and ESO are prone to local convergence failure.

To address this issue, we introduce a non-monotonic temperature factor by incorporating a random coefficient a into each iteration’s temperature coefficient T. This adjustment allows the modified temperature coefficient to balance more effectively between local development modes, as shown in Equation 22.

$$T=a\ast {\mathrm{e}}^{-\frac{t}{t_{max}}}$$

(22)

Figure 3. Temperature factor.

Figure 3. Temperature factor.

2.3.3. Dynamic Boundary-Based Opposition Learning (DBOL)

After carefully examining the mirror imaging strategy used in ESO [17] and the original opposition-based learning (OBL) [18], it is evident that the ESO mirror method effectively enhances the diversity and quality of the population by calculating the opposition solutions for the current position and then selecting the optimal solution.

However, the performance exhibits a certain degree of uncertainty. While opposition solutions may have greater potential for global optimization, the fixed search upper and lower boundaries limit the ability to estimate the direction of new solutions. We propose a Dynamic Boundary-based Opposition Learning (DBOL) to address this issue. DBOL leverages the fact that elite individuals carry more effective information than ordinary ones to refine the opposition process. It fully utilizes the valuable experiences of historical individuals while maintaining a global perspective for new individuals. Through greedy selection, a superior solution is chosen between the current and opposition solutions to participate in the next iteration. The DBOL defines the boundary as follows: Eq. 23:

$$X_i^{j^{\ast }}=r_1\left(U{B}_i^j+L{B}_i^j\right)-X_i^j,$$

(23)

where $r_1$ is a random value in the interval (0, 1), $U{B}^j=max\left(X_i^j\right)$, $LBj=min\left(X_i^j\right)$, and $U{B}^j$ and $L{B}^j$ are the upper and lower dynamic boundary in the $j^{th}$ dimension, respectively.

2.4. Simulation Experiments Design and the Evaluation Criteria

2.4.1. Simulated Experiment

This study’s simulation experiments are divided into two parts: The first part involves selecting the benchmark test functions with a search dimension (D=1000). The objective is to compare different algorithms’ convergence speeds and optimization precision to ascertain whether IESO demonstrates a superior capability to evade local optima when dealing with large-scale optimization problems. The second part entails creating a simulated rugged mountainous terrain to test the path planning. In this experience, the performance of various algorithms is assessed when solving three-dimensional path planning problems based on the three optimization criteria (discussed in 2.4.2). It aims to further validate the feasibility of IESO in determining optimal flight paths within intricate, complex landscapes.

We selected thirteen testing functions shown in Table 1.

We selected five algorithms, namely GJO, GWO, HBA, ESO and SO and their parameters are listed in Table 2.

We randomly generated the digital elevation model of a simulated rugged mountainous terrain, with map dimensions of 100 km by 150 km by 3 km. The maximum elevation difference is at least 2 km, representing a terrain with significant undulations and numerous ravines. Arbitrary no-fly zones were also added to increase the complexity. Flight waypoints were set to evaluate the effectiveness of the flight paths generated by different algorithms.

The experimental setup utilized a Windows 10, 64-bit operating system with an AMD Ryzen 7 4800H CPU running at 2.9 GHz and 16GB of RAM. The algorithms were implemented on the Matlab2020b platform.

2.4.2. Evaluation Criteria

Considering the intricate nature of UAV three-dimensional path planning, this study formulates the path planning objective function by incorporating factors such as flight distance, path smoothness, and flight safety. The shortest flight path is defined as Eq. 24

$$F_f=\sum _{k=1}^D\sqrt{{\left(x_{k+1}-x_k\right)}^2+{\left(y_{k+1}-y_k\right)}^2+{\left(z_{k+1}-z_k\right)}^2}$$

(24)

, where D represents the number of waypoints and $(x_k,y_k,z_k)$ denotes the cooridnates of the $k^{t}h$ waypoints. Ensuring smoothness involves minimizing abrupt changes in yaw angle and altitude. Given mountainous areas’ rugged and diverse terrain, flight paths must adhere to maximum climb angles and rates. The path smoothness is given as Eq. 25, 26, and 27.

$${\phi }_k=arctan\left(\frac{∥l_k{l}_{k+1}∥}{l_k{l}_{k+1}}\right)$$

(25)

$${\varphi }_k=arctan\left(\frac{z_{k+1}-z_k}{\sqrt{{\left(x_{k+1}-x_k\right)}^2+{\left(y_{k+1}-y_k\right)}^2}}\right)$$

(26)

$$F_e=\sum _{k=1}^{D-2}{\phi }_k+\sum _{k=1}^{D-1}\left({\varphi }_k-{\varphi }_{k-1}\right),$$

(27)

where $l_k$ represents the distance between two waypoints. Eq 25 and 26 represent the yaw angle ${\phi }_k$ and pitch angle ${\varphi }_k$, respectively.

A meticulously designed path should prioritize the safe operation of UAVs. The objective function design must incorporate a safety margin, enabling UAVs to avoid obstacles or restricted no-fly zones. As shown in Figure 4 the obstacle within the airspace has a center coordinate of $O_p$ and radius $R_p$. The perpendicular distance $d_p$ from the flight node of the UAV to the obstacle should exceed the safety distance $S_f$, ensuring the UAV operates outside designated shaded regions for safety. The flight safety formula is provided as Eq. 28

$$F_s=\left\{\begin{array}{cc}0,& d_p\ge S\hfill \\ S-d_p,& R_p<d_p<S\hfill \\ \infty ,& d_p<R_p\hfill \end{array}\right.$$

(28)

By integrating these cost functions with weights, an objective function F for the multi-objective path planning problem is formed in Eq. 29

$$F={\omega }_1{F}_f+{\omega }_2{F}_e+{\omega }_3{F}_s$$

(29)

Figure 4. Flight safety to avoid obstacles or restricted no-fly zones.

Figure 4. Flight safety to avoid obstacles or restricted no-fly zones.

3. Results and Discussion

The test results across six different algorithms are summarized, with the algorithm achieving the lowest average optimization result for each test function highlighted in bold, which is shown in Table 3. In the rank-sum test, the symbols “+”, “-”, and “=” denote that IESO optimization is superior, inferior, or equivalent, respectively, compared to other algorithms. Moreover, the statistical analysis indicates that, overall, the optimization performance of IESO surpasses that of most tested algorithms. The integration of multiple strategies significantly enhances the convergence precision and stability of IESO beyond that of the standard SO, preliminarily confirming the effectiveness of these strategies. For any given test function, none of the algorithms significantly outperforms IESO. The optimization curves for various algorithms on the test function with D=1000 are illustrated in Figure 4.

Among the thirteen test functions, our IESO algorithm significantly outperforms other algorithms. This advantage is attributed to the non-monotonic decrease in the temperature factor, which fosters interaction among the foraging, combat, and mating modes. Specifically, the foraging mode utilizes the current global best position to generate high-quality solutions throughout the mid-to-high iterations, accelerating convergence. The only exception is $f_{12}$, where ESO slightly outperforms our IESO. Its numerous local minima characterize this function, posing a considerable challenge for global optimization algorithms. The analysis suggests that introducing the mirroring mechanism has somewhat expanded the exploration horizon of the SO. However, as stated in Wang et al., [20], due to limited executions and dynamically narrowing boundaries over iterations, the development mode, which relies on information exchange between the male and female snake populations, might still result in the optimizer becoming trapped in local optima. Wang et al. suggested incorporating Lévy flights into the ESO for mutating optimal individuals. As a distinctive heavy-tailed distribution, Lévy flights, with their long and short step lengths, help the algorithm escape local optima. In our IESO, we employed a dynamic boundary to the opposite learning mechanism to enhance performance further, which has better results escaping the local optima than ESO. However, in some cases, it may not be fully utilized.

Figure 5. Iteration fitness curves for different test functions.

Figure 5. Iteration fitness curves for different test functions.

A simulated digital elevation model of rugged mountainous terrain was obtained, with the map dimensions measuring 100 km by 150 km by 3 km. The maximum elevation difference in this area exceeds 2 km, featuring a terrain with significant undulations and numerous ravines. Randomly generated "no-fly zones" within this terrain are added and are illustrated in Figure 6 (purple areas). The launch point is set at coordinates (10, 90, 1.115), and the destination is at (130, 10, 1.373), with all units in kilometres (km). The fundamental challenge of path planning lies in identifying suitable waypoints during the flight to achieve an optimal flight route. To this end, this study sets 15 flight waypoints within the three-dimensional map and then employs interpolation to obtain detailed flight paths between adjacent waypoints. The parameters to be optimized are the spatial coordinates of these 15 waypoints.

The six algorithms—GJO, HBA, GWO, SO, ESO, and our IESO—were selected to plan the optimized routing on this simulated terrain. As shown in Figure 7, IESO exhibits a significantly faster convergence time than the other algorithms, at least 20% faster than GJO and over 35% faster than ESO. Additionally, Figure 8 and Figure 9 demonstrate that the route proposed by IESO is the shortest and successfully avoids all no-fly zones, as seen from both 2D and 3D perspectives. It is at least 10% shorter than the competitors. Therefore, we can conclusively state that IESO significantly improves flight planning in complex environments, such as rugged mountain terrains with substantial undulations and numerous ravines, while effectively avoiding arbitrary no-fly zones.

4. Conclusions

In this study, we propose an Improved Enhanced Snake Optimizer (IESO) to optimize 3D UAV path planning in complex conditions such as rugged mountainous terrains. The IESO refines the Enhanced Snake Optimizer (ESO) by incorporating Chebyshev chaotic mapping, a non-monotonically decreasing temperature factor, and dynamic boundary-based opposition learning (DBOL).

We validated the results in simulated experiments using thirteen different test functions under three criteria against five other innovative optimization algorithms. The IESO consistently outperformed these algorithms, demonstrating superior convergence accuracy and speed, which is over 30% shorter in distance and 10% faster in computing speed.