As the application domains of unmanned aerial vehicle (UAV) technology continue to expand, the complexity of flight environments and the scale of obstacle distributions have increased dramatically, making path planning an ever more challenging optimization problem. Efficient path planning is fundamentally critical for the successful deployment of unmanned aerial vehicle systems in real-world missions, particularly in scenarios such as terrain mapping, object detection, precision agriculture, and disaster rescue. However, a significant number of existing path planning algorithms are primarily designed for single-UAV scenarios and frequently suffer from several inherent limitations, including susceptibility to premature convergence toward local optima, slow convergence rates, and difficulties in scaling effectively to multi-UAV cooperative path planning tasks. Furthermore, when confronted with complex missions involving multiple conflicting flight objectives, current algorithms often struggle to strike a satisfactory balance between minimizing path length and ensuring adequate obstacle avoidance, which consequently leads to suboptimal solutions that fail to meet the stringent requirements of optimal path planning in complex environments.
This research addresses these critical challenges from two complementary perspectives. First, to overcome the limitations of existing algorithms that are mostly tailored to single-scenario applications, we propose an Enhanced Teaching-Learning Particle Swarm Optimization algorithm (E-TLPSO). This novel algorithm is capable of effectively handling both single-UAV and multi-UAV path planning tasks. It strategically leverages the exploitation capability of Particle Swarm Optimization (PSO) and the exploration capability of Teaching-Learning-Based Optimization (TLBO). Through three synergistic improvements—chaotic population initialization, nonlinear adaptive parameter adjustment, and dynamic teaching-learning strategies—the algorithm achieves robust optimization performance. E-TLPSO follows a single-objective optimization paradigm, whose core step involves aggregating multiple objectives and constraints in unmanned aerial vehicle path planning into a single unconstrained fitness function using a set of weight coefficients. Nevertheless, the setting of these weight coefficients is intrinsically subjective and empirical, and when dealing with complex path planning problems, this approach can easily lead to entrapment in local optima, making it difficult to obtain high-quality solutions within limited computational time.
Second, to address the issue of suboptimal solution quality in single-objective optimization algorithms, we propose a Multi-Objective Evolutionary Algorithm with Decomposition-based Hybrid Penalty Boundary Intersection (MOEADHPBI). Recognizing that unmanned aerial vehicle path planning is inherently a multi-objective optimization problem where multiple objectives often conflict with each other and require simultaneous trade-off optimization, this algorithm integrates the Penalty Boundary Intersection (PBI) and Inverted Penalty Boundary Intersection (IPBI) methods. To dynamically balance convergence and diversity, an adaptive penalty factor adjustment strategy is introduced. Concurrently, an evolutionary operation mechanism is employed to enhance both global exploration (diversity) and local exploitation (convergence) capabilities. Considering that solutions violating constraints may emerge during the evolutionary process, a constraint repair technique is designed to ensure that the final Pareto-optimal solution set is completely feasible, effectively guiding the search process toward the feasible region.
Through systematic experiments conducted in MATLAB, incorporating nine distinct terrain and threat complexity scenarios, we validate the superiority of both proposed algorithms in terms of convergence precision, solution quality, and comprehensive multi-objective performance metrics.
1. Mathematical Modeling of UAV Path Planning
The practical application of unmanned aerial vehicle path planning requires comprehensive consideration of complex flight environments, flight constraints, and flight objectives. In a three-dimensional space, the path planning algorithm must not only avoid terrain obstacles and dangerous zones but also satisfy the physical performance limitations of the UAV. To ensure the feasibility and practicality of the planned path, a complete mathematical model is established, encompassing flight environment representation, constraint condition definitions, and optimization objective design.
1.1 Three-Dimensional Environment Model
The core objective of unmanned aerial vehicle path planning (UPP) is to find a globally optimal path through optimization algorithms within a complex three-dimensional flight environment, while ensuring that the UAV can safely and efficiently reach the target location. This research considers a scenario involving UAVs, mountains, and radars. The allowable flight area can be modeled as the flight space \(S\), as defined in Equation (1):
$$ S = \{(x, y, z) | x_{min} \leq x \leq x_{max}, y_{min} \leq y \leq y_{max}, z_{min} \leq z \leq z_{max}\} \tag{1} $$
A complete UAV flight path can be represented as a sequence of waypoints:
$$ Path = \{(x_i, y_i, z_i) | (x_i, y_i, z_i) \in S, i = 1, 2, \ldots, M\} \tag{2} $$
where \(M\) denotes the total number of waypoints, including both the start and end points, and \((x_i, y_i, z_i)\) represents the coordinates of the \(i\)-th waypoint.
For the \(k\)-th mountain obstacle (\(k = 1, 2, \ldots, n_{mountain}\)), the terrain elevation at point \((\bar{x}, \bar{y})\) is modeled using a Gaussian-like function:
$$ \bar{z} = h_k(\bar{x}, \bar{y}) = \tau_k \exp\left\{\frac{(\bar{x}-\bar{x}_k^0)^2}{\delta_k} + \frac{(\bar{y}-\bar{y}_k^0)^2}{\sigma_k}\right\} \tag{3} $$
where \((\bar{x}_k^0, \bar{y}_k^0)\) determines the peak center, \(\delta_k\) and \(\sigma_k\) control the spread in the x and y directions respectively, and \(\tau_k\) denotes the maximum height of the \(k\)-th mountain.
For radar threats, each radar \(r\) \((r = 1, 2, \ldots, n_{radar})\) is characterized by a hemispherical detection zone:
$$ DZ_r = \{(\tilde{x}_r, \tilde{y}_r, \tilde{z}_r) | \sqrt{(\tilde{x}_r – \tilde{x}_r^0)^2 + (\tilde{y}_r – \tilde{y}_r^0)^2 + (\tilde{z}_r – \tilde{z}_r^0)^2} \leq d_r\} \tag{4} $$
where \((\tilde{x}_r^0, \tilde{y}_r^0, \tilde{z}_r^0)\) represents the center of the radar detection region and \(d_r\) represents its radius.
1.2 Constraint Functions
UAV path planning requires the exploration of globally optimal solutions under multiple constraints. These constraints primarily encompass flight altitude, pitch angle, and yaw angle considerations. By scientifically designing these constraint functions, the feasibility and optimality of the planned path can be ensured.
For a given flight path \(Path = \{(x_i, y_i, z_i) | i = 1, 2, \ldots, M\}\), the UAV flight altitude cost \(f_{height}\) is designed to ensure that the UAV flies above all mountain peaks. For \(m\) equally spaced points \(\{(x_{ij}, y_{ij}, z_{ij}) | j = 1, 2, \ldots, m\}\) extracted from the flight segment between \((x_i, y_i, z_i)\) and \((x_{i+1}, y_{i+1}, z_{i+1})\), each point is associated with a penalty value \(A_{i,j}\):
$$ A_{i,j} = \begin{cases} 1 & \text{if } \exists k \text{ such that } z_{ij} \leq h_k(x_{ij}, y_{ij}) \\ 0 & \text{otherwise} \end{cases} \tag{5} $$
The total altitude cost is then:
$$ f_{height} = \sum_{i=1}^{M-1} \sum_{j=1}^{m} A_{i,j} \tag{6} $$
For pitch angle \(\alpha_i\) and yaw angle \(\eta_i\) constraints, the costs are defined as:
$$ f_{pitch} = \sum_{i=1}^{M-1} B_i \quad \text{where} \quad B_i = \begin{cases} 0 & \text{if } \alpha_{min} \leq \alpha_i \leq \alpha_{max} \\ 1 & \text{otherwise} \end{cases} \tag{7} $$
$$ f_{yaw} = \sum_{i=1}^{M-1} C_i \quad \text{where} \quad C_i = \begin{cases} 0 & \text{if } \eta_{min} \leq \eta_i \leq \eta_{max} \\ 1 & \text{otherwise} \end{cases} \tag{8} $$
The pitch angle \(\alpha_i\) and yaw angle \(\eta_i\) can be computed as:
$$ \alpha_i = \arctan\frac{z_{i+1}-z_i}{\sqrt{(x_{i+1}-x_i)^2 + (y_{i+1}-y_i)^2}} \tag{9} $$
$$ \eta_i = \arctan\frac{y_{i+1}-y_i}{x_{i+1}-x_i} \tag{10} $$
1.3 Objective Functions
The optimization objectives in UPP directly influence the performance and effectiveness of the path planning algorithm. This research considers two key objectives: path length and radar threat.
The path length cost \(f_{length}\) is minimized to enhance operational efficiency:
$$ f_{length} = \sum_{i=1}^{M-1} \sqrt{(x_{i+1} – x_i)^2 + (y_{i+1} – y_i)^2 + (z_{i+1} – z_i)^2} \tag{11} $$
The radar threat cost \(f_{threat}\) quantifies exposure risks. For each of the \(m\) points sampled from the flight segment, a cost value \(G_{i,j}\) is defined:
$$ f_{threat} = \sum_{i=1}^{M-1} \sum_{j=1}^{m} G_{i,j} \quad \text{where} \quad G_{i,j} = \begin{cases} 1 & \text{if } \exists r \text{ such that } (x_{ij}, y_{ij}, z_{ij}) \in DZ_r \\ 0 & \text{otherwise} \end{cases} \tag{12} $$
This mathematical modeling framework transforms complex flight problems into computable optimization problems, providing a solid theoretical foundation for the subsequent algorithm design.
2. Enhanced Teaching-Learning Particle Swarm Optimization for Single-Objective UAV Path Planning
Particle Swarm Optimization (PSO) achieves optimization search by simulating the foraging behavior of bird flocks and fish schools. With advantages such as simple implementation and fast convergence speed, PSO has demonstrated good performance in unmanned aerial vehicle path planning. PSO combines two important mechanisms—the personal best mechanism and the global best mechanism—to find optimal or near-optimal solutions. Through the personal best mechanism, each particle in PSO continuously improves the quality of its solution by moving toward its historical best position. Meanwhile, through the global best mechanism, the entire population converges toward the current best solution found by any particle. These two mechanisms enable PSO to find good solutions quickly and efficiently, but they also make it susceptible to falling into local optima when solving UPP problems.
To address this limitation, we introduce Teaching-Learning-Based Optimization (TLBO). TLBO finds optimal or near-optimal solutions through its teaching phase and learning phase. In the teaching phase, the entire population is guided toward the current average level position. In the learning phase, the algorithm promotes information exchange among individuals through a random learning mechanism, thereby providing better global exploration capability.
2.1 Problem Definition for Single-Objective UPP
To find the optimal flight path, we integrate the flight constraints and flight objectives described in the previous section by constructing a fitness function \(F(\cdot)\) that evaluates the overall quality of any flight path:
$$ F(Path) = \rho_1 f_{height} + \rho_2 f_{pitch} + \rho_3 f_{yaw} + \rho_4 f_{length} + \rho_5 f_{threat} \tag{13} $$
where \(\rho_1, \ldots, \rho_5\) represent weight coefficients that need to be set in advance based on specific task or environmental requirements. Through this transformation, the UPP problem is converted into an unconstrained minimization problem, with the goal of finding a path (for single-UAV scenarios) or multiple paths (for multi-UAV scenarios) that minimizes the fitness value.
2.2 Standard Particle Swarm Optimization
In PSO, for the \(i\)-th particle at the \(t\)-th iteration (\(t = 1, 2, \ldots, t_{max}\)), the position \(X_i(t)\) and velocity \(V_i(t)\) are updated according to:
$$ V_i(t+1) = \omega \times V_i(t) + c_1 \times r_1 \times (p_{best_i}(t) – X_i(t)) + c_2 \times r_2 \times (g_{best}(t) – X_i(t)) \tag{14} $$
$$ X_i(t+1) = X_i(t) + V_i(t+1) \tag{15} $$
where \(\omega\) is the inertia weight, \(c_1\) and \(c_2\) are learning factors, \(r_1\) and \(r_2\) are random values in the interval [0,1), and \(p_{best_i}(t)\) and \(g_{best}(t)\) represent the best solutions found so far by the \(i\)-th particle and the entire population, respectively.
2.3 Standard Teaching-Learning-Based Optimization
TLBO simulates the teaching and learning process in a classroom, comprising two phases. In the teaching phase, each individual learns from the teacher:
$$ X_{teach_i}(t) = X_i(t) + r_3 \times (X_{teacher}(t) – T_f \times X_{center}(t)) \tag{16} $$
where \(r_3\) is a random value in [0,1), \(T_f \in [1,2)\) is a coefficient affecting teaching capability, and \(X_{center}(t) = \frac{1}{N} \sum_{i=1}^{N} X_i(t)\) is the average of all individuals.
In the learning phase, each individual learns from two randomly selected individuals \(X_u(t)\) and \(X_v(t)\):
$$ X_{learn_i}(t) = \begin{cases} X_i(t) + r_4 \times (X_u(t) – X_v(t)) & \text{if } F(X_u(t)) < F(X_v(t)) \\ X_i(t) + r_5 \times (X_v(t) – X_u(t)) & \text{otherwise} \end{cases} \tag{17} $$
2.4 The Proposed E-TLPSO Algorithm
The proposed E-TLPSO algorithm synergistically integrates PSO and TLBO, leveraging TLBO’s advantages in exploring promising regions of the search space to enhance PSO’s exploration capability. The framework consists of three key components: population initialization, dynamic teaching-learning (derived from TLBO), and adaptive evolution (derived from PSO).
2.4.1 Chaotic Population Initialization
To increase the diversity and uniform distribution of the initial population, a chaotic method is adopted. The chaotic factor \(\gamma_i\) is generated through the logistic map:
$$ \gamma_i = \mu \times \gamma_{i-1} \times (1 – \gamma_{i-1}) \tag{18} $$
where \(\mu\) is the chaotic coefficient (typically set to 4), and the random seed \(\gamma_0\) is randomly selected from the interval [0,1). The \(j\)-th dimension of the \(i\)-th individual is then determined by:
$$ X_{i,j}(0) = X_{min,j} + \gamma_i \times (X_{max,j} – X_{min,j}) \tag{19} $$
2.4.2 Dynamic Teaching and Learning
To address the issue of absolute randomness in pseudo-random number generators, a chaotic mechanism is introduced. The chaotic factor \(\vartheta_i\) replaces the random numbers \(r_3\), \(r_4\), and \(r_5\) in the TLBO update equations:
$$ X_{teach_i}(t) = X_i(t) + \vartheta_i \times (X_{teacher}(t) – T_f \times X_{center}(t)) \tag{20} $$
$$ X_{learn_i}(t) = X_i(t) + \vartheta_i \times (X_u(t) – X_v(t)) \tag{21} $$
2.4.3 Adaptive Evolution
A novel inertia weight adjustment method is proposed, which considers the population state:
$$ \omega = (\omega_{max} – \omega_{min}) \times N_s/N + \zeta \times \omega_{min} \tag{22} $$
where \(\omega_{min}\) and \(\omega_{max}\) denote the minimum and maximum inertia weights, \(N_s\) is the number of particles whose personal best solutions were updated in the current iteration, and \(\zeta\) is a chaotic factor. When \(N_s/N\) is large, it suggests that the current evolutionary direction is favorable, and the inertia weight is increased to encourage further exploration. Conversely, when \(N_s/N\) is small, the inertia weight is decreased to enhance local exploitation.
The learning factors are adjusted using a time-varying linear adjustment method:
$$ c_1 = (c_{min} – c_{max}) \times t/t_{max} + c_{max} \tag{23} $$
$$ c_2 = (c_{max} – c_{min}) \times t/t_{max} + c_{min} \tag{24} $$
2.5 Computational Complexity Analysis
The total time complexity of E-TLPSO is \(O(t_1DN + t_2DN + t_3DN)\). Since the population size \(N\) is typically much larger than the optimization dimension \(D\) in single-objective optimization problems, the total time complexity can be simplified to \(O(t_{max}N)\), which is consistent with most existing single-objective evolutionary algorithms.
3. Multi-Objective UAV Path Planning with MOEADHPBI
Single-objective evolutionary algorithms aggregate multiple objectives and constraints into a single unconstrained optimization problem using weight factors. However, when handling multi-objective and multi-constraint problems, these algorithms have limitations in weight coefficient setting, requiring extensive experiments to determine approximately optimal weights.
To overcome these limitations, this chapter defines the unmanned aerial vehicle path planning problem as a constrained multi-objective optimization problem:
$$ \min F(Path) = (f_{length}(Path), f_{threat}(Path)) \tag{25} $$
$$ s.t. \quad f_{pitch}(Path) + f_{yaw}(Path) \leq \varepsilon $$
where \(f_{length}(Path)\) evaluates the flight distance, \(f_{threat}(Path)\) quantifies the radar detection risk, and \(\varepsilon\) is a preset threshold for constraint violation.
3.1 Individual Encoding and Population Initialization
For UPP, a one-dimensional array with elements \((x_{ij}, y_{ij}, z_{ij})\) \((j = 1, 2, \ldots, M)\) represents the \(i\)-th individual, where \(M\) is the number of path points. During initialization, \(N\) weight vectors \(\lambda = \{\lambda_1, \lambda_2, \ldots, \lambda_N\}\) and population individuals are randomly generated. For the \(i\)-th weight vector \(\lambda_i \in \lambda\), its neighborhood index set \(Nes_i\) consists of indices corresponding to the \(T\) closest weight vectors, where closeness is determined by the Euclidean distance \(din_{ij} = \|\lambda_i – \lambda_j\|\).
3.2 Evolutionary Operation Mechanism
To help population individuals obtain better planned paths, an evolutionary operation mechanism is employed after initialization. A random number \(r\) is generated and compared with the threshold \(\delta^*\). If \(r \leq \delta^*\), the \(i\)-th neighborhood index set \(Nes_i\) is selected as the evolutionary operation index set \(Ne\), meaning local exploitation is performed. Otherwise, the entire population indices are selected as \(Ne\), performing global exploration.
The new individual \(Off_r\) is generated using simulated binary crossover and polynomial mutation based on the selected parents from \(Pop(Ne)\).
3.3 Adaptive Hybrid Penalty Boundary Intersection Strategy
The traditional PBI method is defined as:
$$ \min_x \{g_{pbi}(x|\lambda, Z^+) = d_1 + \theta_1 d_2\} \tag{26} $$
where:
$$ d_1 = \frac{\|(Z^+ – F(x))^T \lambda\|}{\|\lambda\|}, \quad d_2 = \|F(x) – (Z^+ – d_1\lambda)\| \tag{27} $$
The IPBI method, in contrast, maximizes the scalarizing function value to push solutions outward from the reference point \(Z^-\):
$$ \max_x \{g_{ipbi}(x|\lambda, Z^-) = d_1 – \theta_2 d_2\} \tag{28} $$
where:
$$ d_1 = \frac{\|(Z^- – F(x))^T \lambda\|}{\|\lambda\|}, \quad d_2 = \|(Z^- – F(x)) – d_1\lambda\| \tag{29} $$
To dynamically balance convergence and diversity, the selection probability \(r_\partial\) is computed:
$$ r_\partial = \frac{1}{1 + e^{-(1-r_\tau)}} \tag{30} $$
where \(r_\tau = r_\delta / r_\sigma\) is the population distribution change rate, \(r_\delta\) represents the convergence change rate, and \(r_\sigma\) represents the diversity change rate.
The adaptive penalty factors are adjusted as:
$$ \theta_1 = \frac{1}{1 + \exp(-\varphi * t/t_{max})} * \gamma_t + (\theta_{max} – \theta_{min}) * (t/t_{max})^2 \tag{31} $$
$$ \theta_2 = \theta_{max} * \gamma_t + (\theta_{min} – \theta_{max}) * t/t_{max} \tag{32} $$
When \(r_\partial > r_\partial^*\), the adaptive PBI is adopted to enhance diversity. Otherwise, the IPBI mechanism is executed to improve convergence.
3.4 Constraint Repair Technique
Path planning is a complex multi-objective optimization problem, and the planned flight paths may violate constraints. To address this, a constraint repair technique is proposed to convert infeasible solutions into feasible ones. For individuals with constraint violation values exceeding \(\varepsilon\), the pitch and yaw angles are computed, and when they exceed the constraint range, they are repaired:
$$ \begin{cases} \alpha_{new,i} = \max(\alpha_{min}, \min(\alpha_i, \alpha_{max})) \\ \eta_{new,i} = \max(\eta_{min}, \min(\eta_i, \eta_{max})) \end{cases} \tag{33} $$
The corrected flight path point coordinates are then calculated:
$$ \begin{cases} x_i^* = \max(x_{min}, \min(x_{max}, x_i + DS_i \cdot \cos(\eta_{new,i}))) \\ y_i^* = \max(y_{min}, \min(y_{max}, y_i + DS_i \cdot \sin(\eta_{new,i}))) \\ z_i^* = \max(z_{min}, \min(z_{max}, z_i + DS_i \cdot \tan(\alpha_{new,i}))) \end{cases} \tag{34} $$
where \(DS_i = \sqrt{(x_{i+1} – x_i)^2 + (y_{i+1} – y_i)^2}\). To reduce computational cost, the constraint repair technique is executed at regular intervals (\(t_{mod}\)).
3.5 Flight Path Selection Mechanism
After the algorithm reaches the termination condition, a set of Pareto solutions is generated. A flight path selection mechanism is designed to select the most suitable solution from these solutions, prioritizing those without collisions. The mechanism selects the solution with the smallest collision risk value first, and if equal, the one with the smallest path length is selected.
3.6 Computational Complexity Analysis
The total time complexity of MOEADHPBI is \(O(t_1lN^2 + t_2lN^2 + t_3lN + t_4lN)\). Since the population size \(N\) is typically much larger than the number of optimization objectives \(l\), the total time complexity can be simplified to \(O(t_{max}lN^2)\), consistent with most existing multi-objective evolutionary algorithms.
4. Experimental Results and Analysis
4.1 Simulation Setup
The simulation environment is based on MATLAB 2023a platform on Windows 10 operating system, equipped with an Intel i5 processor (3.10 GHz) and 32GB memory. The system model parameters are listed in Table 1.
| Parameter | Value | Parameter | Value |
|---|---|---|---|
| \(\alpha_{min}\) | -60° | \(\alpha_{max}\) | 60° |
| \(\eta_{min}\) | -60° | \(\eta_{max}\) | 60° |
| \(x_{min}\) | 0m | \(x_{max}\) | 60m |
| \(y_{min}\) | 0m | \(y_{max}\) | 60m |
| \(z_{min}\) | 0m | \(z_{max}\) | 10m |
For a comprehensive evaluation, nine scenarios with varying terrain and threat complexities are considered. The obstacle information for these scenarios includes mountain peaks modeled by Gaussian functions and radar threats modeled by hemispherical detection zones.
4.2 Benchmark Function Optimization Performance
Before applying E-TLPSO to UPP, we first verified its basic optimization capability on 15 classical benchmark functions in 10 dimensions (10D) and 50 dimensions (50D). We compared E-TLPSO with AIWPSO, TCSPSO, MPSO, AO, DE, GA, and HBA algorithms. Each algorithm performed 25 independent runs on each benchmark function.
The statistical results demonstrate that E-TLPSO performs excellently in handling unimodal, multimodal, and hybrid multimodal optimization problems. Table 2 summarizes the performance comparison in 10D and 50D:
| Dimension | Metric | AIWPSO | TCSPSO | MPSO | AO | DE | GA | HBA | E-TLPSO |
|---|---|---|---|---|---|---|---|---|---|
| 10D | Mean | 3 | 1 | 0 | 0 | 0 | 0 | 9 | 2 |
| Max | 5 | 1 | 0 | 0 | 0 | 0 | 7 | 2 | |
| Min | 1 | 1 | 0 | 0 | 0 | 0 | 6 | 7 | |
| Std | 5 | 2 | 0 | 0 | 3 | 0 | 4 | 1 | |
| 50D | Mean | 3 | 2 | 0 | 0 | 0 | 0 | 4 | 6 |
| Max | 2 | 2 | 0 | 0 | 0 | 0 | 5 | 6 | |
| Min | 3 | 3 | 0 | 0 | 0 | 0 | 1 | 8 | |
| Std | 2 | 3 | 1 | 0 | 2 | 0 | 2 | 5 |
The results verify the effectiveness and superiority of E-TLPSO in solving high-dimensional optimization problems.
4.3 Single-UAV Path Planning Performance
We evaluated E-TLPSO against representative algorithms including CIPSO, TVAPSO, HCPSOA (PSO variants for UPP), HPSO, AWPSO (general optimization strategies), and GA, HBA (other evolutionary algorithms for UPP). Each algorithm was independently run 18 times in each scenario.
Table 3 presents the performance results across three different scenarios:
| Scenario | Metric | HCPSOA | CIPSO | AWPSO | HPSO | TVAPSO | GA | HBA | E-TLPSO |
|---|---|---|---|---|---|---|---|---|---|
| Scenario 1 | APL | 72.128 | 72.128 | 72.128 | 72.128 | 72.129 | 72.130 | 72.130 | 72.126 |
| ANI | 100 | 100 | 100 | 80 | 30 | 20 | 10 | 50 | |
| FAOFV | 2163.85 | 2163.85 | 2163.83 | 2163.85 | 2163.87 | 2163.90 | 2163.93 | 2163.79 | |
| Scenario 2 | APL | 74.47 | 74.33 | 74.00 | 74.18 | 74.28 | 74.54 | 74.43 | 72.95 |
| ANI | 100 | 100 | 100 | 100 | 20 | 50 | 70 | 70 | |
| FAOFV | 2234.20 | 2229.76 | 2220.08 | 2225.38 | 2228.36 | 2236.22 | 2232.98 | 2188.46 | |
| Scenario 3 | APL | 78.35 | 78.85 | 77.56 | 77.84 | 77.44 | 80.10 | 80.17 | 75.64 |
| ANI | 100 | 100 | 100 | 100 | 30 | 100 | 100 | 40 | |
| FAOFV | 2347.94 | 2351.04 | 2333.69 | 2342.77 | 2340.31 | 2402.86 | 2405.09 | 2269.30 |
In Scenario 3 (the most complex single-UAV environment), E-TLPSO achieved an APL of 75.64m and FAOFV of 2269.30, significantly outperforming all comparison algorithms.
4.4 Multi-UAV Path Planning Performance
For multi-UAV path planning, a formation of five UAVs is considered. The algorithms are evaluated based on the total cost of the five generated paths. Paths must be collision-free, with a minimum inter-path distance threshold of 2m.
| UAV ID | AWPSO | CIPSO | HCPSOA | HPSO | TVAPSO | GA | HBA | E-TLPSO |
|---|---|---|---|---|---|---|---|---|
| UAV1 | 82.94 | 79.32 | 79.63 | 77.86 | 79.60 | 81.02 | 82.47 | 76.26 |
| UAV2 | 83.19 | 81.65 | 81.11 | 80.88 | 78.50 | 83.83 | 84.39 | 76.32 |
| UAV3 | 82.57 | 82.53 | 79.61 | 80.62 | 78.58 | 80.97 | 80.79 | 76.79 |
| UAV4 | 84.19 | 79.83 | 81.96 | 80.39 | 79.93 | 82.38 | 83.53 | 76.55 |
| UAV5 | 82.94 | 81.70 | 81.02 | 80.29 | 79.05 | 80.74 | 80.87 | 76.61 |
| AL | 83.17 | 81.01 | 80.67 | 80.01 | 79.13 | 81.79 | 82.41 | 76.51 |
E-TLPSO not only achieves the shortest average length (AL=76.51m) for all five paths but also produces the shortest individual path for each UAV. In Scenario 4, E-TLPSO reduced FAOFV by 37.82% compared to HBA.
| Scenario | Metric | HCPSOA | CIPSO | AWPSO | HPSO | TVAPSO | GA | HBA | E-TLPSO |
|---|---|---|---|---|---|---|---|---|---|
| Scenario 4 | ANI | 120 | 120 | 210 | 120 | 90 | 210 | 210 | 30 |
| FAOFV | 13802.22 | 13630.49 | 17260.37 | 10549.68 | 14217.42 | 15102.15 | 16525.80 | 10272.03 | |
| Scenario 5 | ANI | 210 | 210 | 210 | 90 | 60 | 150 | 180 | 60 |
| FAOFV | 13125.33 | 16151.08 | 15271.7 | 11327.14 | 20790.51 | 15044.57 | 17343.07 | 11165.30 |
4.5 Ablation Study of E-TLPSO
To validate the effectiveness of the three important components of E-TLPSO, we designed an ablation study with three variant algorithms: E-TLPSO-C (random initialization instead of chaotic), E-TLPSO-L (linear inertia weight adjustment instead of nonlinear), and E-TLPSO-D (without dynamic teaching-learning strategies).
| Scenario | Metric | E-TLPSO-C | E-TLPSO-L | E-TLPSO-D | E-TLPSO |
|---|---|---|---|---|---|
| Scenario 1 | FAOFV | 2164.75 | 2163.83 | 2164.01 | 2163.77 |
| APL | 72.158 | 72.128 | 72.134 | 72.126 | |
| Scenario 2 | FAOFV | 2191.33 | 2200.06 | 2190.30 | 2188.46 |
| APL | 73.04 | 73.34 | 73.01 | 72.95 | |
| Scenario 3 | FAOFV | 2325.39 | 2349.28 | 2312.54 | 2269.30 |
| APL | 77.51 | 78.31 | 77.09 | 75.64 | |
| Scenario 4 | FAOFV | 10531.48 | 14415.20 | 11348.25 | 10272.03 |
| AL | 76.22 | 81.24 | 76.35 | 76.51 | |
| Scenario 5 | FAOFV | 11246.80 | 12656.98 | 11245.43 | 11165.30 |
| AL | 74.93 | 80.86 | 74.87 | 74.47 |
The results demonstrate that E-TLPSO outperforms all three ablation variants in almost all scenarios, particularly in complex environments. The dynamic teaching-learning strategy proves essential for maintaining solution quality, while chaotic initialization and nonlinear parameter adjustment contribute to escaping local optima.
4.6 Multi-Objective Benchmark Function Performance
We validated MOEADHPBI on four benchmark function sets: ZDT (5 problems), DTLZ (7 problems), MaF (13 problems), and WFG (9 problems). The algorithm was compared with eight unconstrained MOEAs: EAGMOEAD, FLEA, IMMOEAD, MOEADPaS, MOBCA, MaOEAIGD, CoMMEA, and HREA.
Tables 7-10 summarize the statistical results using two performance metrics: Inverted Generational Distance (IGD) and HyperVolume (HV).
$$ IGD = \frac{\sum_{i=1}^{c} d_i}{c} \tag{35} $$
where \(c\) represents the number of sampling points on the true PF, and \(d_i\) represents the minimum Euclidean distance from the \(i\)-th sampling point to the approximate PF obtained by the algorithm.
$$ HV(U) = VOL\left(\bigcup_{x \in U} [f_1(x), z_1^+] \times \ldots \times [f_l(x), z_l^+]\right) \tag{36} $$
| Problem | Metric | Algorithm | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| EAGMOEAD | FLEA | IMMOEAD | MOEADPaS | MOBCA | MaOEAIGD | CoMMEA | HREA | MOEADHPBI | ||
| ZDT1 | HV | 0.9864 | 0.9851 | 0.9912 | 0.9897 | 0.9873 | 0.9786 | 0.9752 | 0.9801 | 0.9923 |
| IGD | 0.0241 | 0.0263 | 0.0182 | 0.0201 | 0.0235 | 0.0312 | 0.0348 | 0.0296 | 0.0157 | |
| ZDT2 | HV | 0.9812 | 0.9823 | 0.9887 | 0.9851 | 0.9826 | 0.9743 | 0.9718 | 0.9765 | 0.9908 |
| IGD | 0.0287 | 0.0276 | 0.0198 | 0.0224 | 0.0271 | 0.0345 | 0.0378 | 0.0325 | 0.0162 | |
| ZDT3 | HV | 0.9897 | 0.9962 | 0.9931 | 0.9918 | 0.9925 | 0.9821 | 0.9795 | 0.9837 | 0.9973 |
| IGD | 0.0203 | 0.0145 | 0.0172 | 0.0191 | 0.0186 | 0.0278 | 0.0312 | 0.0265 | 0.0134 | |
| ZDT4 | HV | 0.9953 | 0.9941 | 0.9862 | 0.9835 | 0.9821 | 0.9612 | 0.9587 | 0.9652 | 0.9968 |
| IGD | 0.0152 | 0.0178 | 0.0256 | 0.0283 | 0.0297 | 0.0389 | 0.0423 | 0.0351 | 0.0141 | |
| ZDT6 | HV | 0.9902 | 0.9915 | 0.9878 | 0.9842 | 0.9838 | 0.9725 | 0.9691 | 0.9756 | 0.9938 |
| IGD | 0.0189 | 0.0167 | 0.0215 | 0.0248 | 0.0253 | 0.0321 | 0.0356 | 0.0302 | 0.0148 | |
According to the Wilcoxon rank-sum test results, MOEADHPBI demonstrates statistically significant superiority over all comparison algorithms in both HV and IGD metrics across the ZDT test suite.
| Problem | Metric | Algorithm | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| EAGMOEAD | FLEA | IMMOEAD | MOEADPaS | MOBCA | MaOEAIGD | CoMMEA | HREA | MOEADHPBI | ||
| DTLZ1 | HV | 0.9821 | 0.9835 | 0.9956 | 0.9912 | 0.9878 | 0.9652 | 0.9618 | 0.9687 | 0.9972 |
| IGD | 0.0253 | 0.0238 | 0.0124 | 0.0162 | 0.0198 | 0.0345 | 0.0382 | 0.0316 | 0.0112 | |
| DTLZ2 | HV | 0.9842 | 0.9856 | 0.9918 | 0.9885 | 0.9849 | 0.9721 | 0.9962 | 0.9754 | 0.9938 |
| IGD | 0.0228 | 0.0215 | 0.0163 | 0.0192 | 0.0217 | 0.0286 | 0.0354 | 0.0271 | 0.0135 | |
| DTLZ3 | HV | 0.9786 | 0.9802 | 0.9991 | 0.9945 | 0.9923 | 0.9412 | 0.9385 | 0.9458 | 0.9996 |
| IGD | 0.0295 | 0.0278 | 0.0082 | 0.0126 | 0.0154 | 0.0428 | 0.0462 | 0.0395 | 0.0078 | |
| DTLZ4 | HV | 0.9815 | 0.9829 | 0.9896 | 0.9862 | 0.9837 | 0.9685 | 0.9652 | 0.9712 | 0.9931 |
| IGD | 0.0268 | 0.0251 | 0.0187 | 0.0213 | 0.0238 | 0.0312 | 0.0348 | 0.0296 | 0.0142 | |
| DTLZ5 | HV | 0.9752 | 0.9786 | 0.9854 | 0.9821 | 0.9798 | 0.9623 | 0.9598 | 0.9667 | 0.9912 |
| IGD | 0.0315 | 0.0298 | 0.0221 | 0.0245 | 0.0268 | 0.0352 | 0.0387 | 0.0335 | 0.0168 | |
| DTLZ6 | HV | 0.9738 | 0.9762 | 0.9835 | 0.9801 | 0.9772 | 0.9587 | 0.9552 | 0.9615 | 0.9896 |
| IGD | 0.0332 | 0.0315 | 0.0248 | 0.0273 | 0.0295 | 0.0378 | 0.0412 | 0.0362 | 0.0193 | |
| DTLZ7 | HV | 0.9235 | 0.9652 | 0.9786 | 0.9712 | 0.9684 | 0.9325 | 0.9287 | 0.9356 | 0.9892 |
| IGD | 0.0452 | 0.0286 | 0.0215 | 0.0248 | 0.0262 | 0.0358 | 0.0393 | 0.0336 | 0.0182 | |
4.7 Multi-Objective UAV Path Planning Performance
We compared MOEADHPBI with eight state-of-the-art constrained MOEAs: BiCo, CCMO, cDPEA, CMOPSO, DPPPS, MCCMO, MSCMO, and MOEA2DE across four scenarios (Scenarios 6-9).
| Algorithm | Scenario 6 | Scenario 7 | ||
|---|---|---|---|---|
| HV | IGD | HV | IGD | |
| BiCo | 1.6678e-1 (±1.31e-1) | 2.4603e+1 (±1.72e+1) | 2.1087e-1 (±1.23e-1) | 1.9091e+1 (±1.76e+1) |
| CCMO | 2.6793e-1 (±5.08e-2) | 1.0528e+1 (±5.59e+0) | 2.6818e-1 (±5.25e-2) | 1.0508e+1 (±5.76e+0) |
| cDPEA | 3.0243e-1 (±3.65e-2) | 6.8088e+0 (±3.99e+0) | 3.1502e-1 (±3.09e-2) | 5.5089e+0 (±3.35e+0) |
| CMOPSO | 2.5952e-1 (±1.67e-2) | 1.1452e+1 (±1.84e+0) | 2.5623e-1 (±1.93e-2) | 1.1830e+1 (±2.12e+0) |
| MCCMO | 2.6986e-1 (±5.48e-2) | 1.0346e+1 (±6.00e+0) | 2.6698e-1 (±5.38e-2) | 1.0650e+1 (±5.91e+0) |
| MSCMO | 2.6699e-1 (±5.58e-2) | 1.0660e+1 (±6.11e+0) | 2.6223e-1 (±5.20e-2) | 1.1197e+1 (±5.69e+0) |
| DPPPS | 3.0118e-1 (±3.33e-2) | 7.0774e+0 (±3.53e+0) | 2.9796e-1 (±2.83e-2) | 7.5069e+0 (±3.05e+0) |
| MOEA2DE | 2.0266e-1 (±4.31e-3) | 1.0222e+1 (±7.07e-1) | 2.0062e-1 (±2.94e-3) | 1.0271e+1 (±7.40e-1) |
| MOEADHPBI | 6.1201e+2 (±8.32e+0) | 1.2290e-1 (±6.43e-2) | 6.5912e+2 (±6.60e+0) | 1.3453e-1 (±6.70e-2) |
In Scenarios 6 and 7, MOEADHPBI achieves HV values exceeding 612 and 659 respectively, while the best-performing comparison algorithm (cDPEA) achieves only approximately 0.30. This represents a more than 2000-fold improvement in HV, demonstrating the quantity-level advantage of the Pareto solution set obtained by MOEADHPBI in terms of coverage breadth and distribution uniformity.
| Algorithm | Scenario 8 | Scenario 9 | ||
|---|---|---|---|---|
| HV | IGD | HV | IGD | |
| BiCo | 1.7516e-1 (±1.08e-1) | 2.6143e+1 (±1.19e+1) | 1.6954e-1 (±1.00e-1) | 2.6188e+1 (±1.20e+1) |
| CCMO | 2.1136e-1 (±4.98e-2) | 1.7704e+1 (±5.05e+0) | 1.8448e-1 (±6.00e-2) | 2.0347e+1 (±6.22e+0) |
| cDPEA | 2.3814e-1 (±4.52e-2) | 1.5871e+1 (±4.47e+0) | 2.2338e-1 (±4.70e-2) | 1.7167e+1 (±4.19e+0) |
| CMOPSO | 2.2710e-1 (±1.84e-2) | 1.5931e+1 (±1.93e+0) | 2.1625e-1 (±1.51e-2) | 1.7132e+1 (±1.54e+0) |
| MCCMO | 1.9990e-1 (±4.66e-2) | 1.8749e+1 (±4.72e+0) | 1.9803e-1 (±6.92e-2) | 2.0020e+1 (±6.50e+0) |
| MSCMO | 1.9902e-1 (±5.45e-2) | 1.9076e+1 (±5.44e+0) | 1.5844e-1 (±7.52e-2) | 2.5932e+1 (±1.49e+1) |
| DPPPS | 2.3943e-1 (±5.29e-2) | 1.3093e+1 (±6.14e+0) | 2.2274e-1 (±2.71e-2) | 1.6766e+1 (±2.29e+0) |
| MOEA2DE | 9.2396e-2 (±4.10e-3) | 9.4595e+0 (±3.20e-14) | 9.2333e-2 (±4.21e-3) | 1.0251e+1 (±4.62e-14) |
| MOEADHPBI | 7.7942e+2 (±1.71e-1) | 1.1730e-1 (±6.07e-2) | 8.1556e+2 (±5.12e-12) | 1.1070e-1 (±6.24e-2) |
In Scenario 8, the IGD value of MOEADHPBI is 0.1173, which is only approximately 1.2% of the best comparison algorithm result (MOEA2DE, 9.46), representing a relative performance improvement exceeding 98%. Similarly, in Scenario 9, MOEADHPBI achieves an even more impressive HV value of 815.56, with the comparison algorithms achieving at best around 0.24, representing an improvement of over 3000-fold.
5. Conclusion and Future Work
This research has made significant contributions to the field of unmanned aerial vehicle path planning through the development of two complementary evolutionary algorithms. The enhanced teaching-learning particle swarm optimization algorithm (E-TLPSO) demonstrates superior performance in both single-UAV and multi-UAV path planning scenarios. By seamlessly integrating the exploitation capability of PSO with the exploration capability of TLBO, and augmenting it with chaotic initialization, nonlinear adaptive parameter adjustment, and dynamic teaching-learning strategies, E-TLPSO achieves excellent convergence speed and solution quality. The experimental results confirm that E-TLPSO can plan high-quality, collision-free paths with the shortest path lengths and best average fitness values compared to state-of-the-art algorithms, particularly in complex scenarios with dense radar threats and high-altitude mountains.
The multi-objective evolutionary algorithm with decomposition-based hybrid penalty boundary intersection (MOEADHPBI) addresses the fundamental limitations of single-objective approaches by treating unmanned aerial vehicle path planning as a true multi-objective optimization problem. Through the intelligent fusion of PBI and IPBI methods, coupled with adaptive penalty factor adjustment, evolutionary operation mechanisms, and constraint repair techniques, MOEADHPBI achieves an excellent balance between convergence and diversity. The extensive validation across benchmark functions and complex path planning scenarios demonstrates that MOEADHPBI consistently outperforms state-of-the-art constrained multi-objective evolutionary algorithms, achieving orders-of-magnitude improvements in both HyperVolume and Inverted Generational Distance metrics.
These findings collectively underscore the significant potential of evolutionary algorithms in addressing the increasingly complex path planning challenges faced by unmanned aerial vehicle systems. The proposed algorithms provide robust and efficient solutions for both single-UAV and multi-UAV scenarios, contributing to the advancement of autonomous flight capabilities.
Future research directions will focus on several key areas: (1) constructing a physical validation platform to test the algorithms on real unmanned aerial vehicle systems, integrating Mobile SDK and Onboard SDK interfaces; (2) introducing multi-scale search strategies to achieve better balance between global search and local exploitation in high-dimensional path planning problems, particularly through coarse global search focused on critical waypoints and fine local search for waypoint refinement; (3) developing heterogeneous multi-UAV cooperative planning methods for obstacles and communication constraints in contested environments; and (4) optimizing computational efficiency through GPU-parallel evolutionary computing, exploring lightweight network models for path evaluation to reduce computational complexity while maintaining solution accuracy.