In the past decade, the proliferation of unmanned aerial vehicles (UAVs), commonly known as drones, has transformed a wide range of industries, including precision agriculture, search and rescue, infrastructure inspection, environmental monitoring, and military operations. One of the fundamental challenges that underpins the autonomy of unmanned aerial vehicles is path planning. The objective of path planning is to generate a feasible, safe, and efficient trajectory from a starting point to a target location while satisfying various operational constraints. In complex three-dimensional environments, the path planner must simultaneously consider terrain obstacles, radar threats, restricted flight zones, and the kinematic and dynamic limitations of the aircraft. As the scope of missions expands, the complexity of the environment and the number of objectives to be balanced increase significantly. Traditional algorithms such as Dijkstra, A*, and artificial potential field methods often fail to provide high-quality solutions in high-dimensional or cluttered environments due to premature convergence, poor scalability, or sensitivity to local minima. In this context, evolutionary algorithms have emerged as a powerful class of metaheuristics that can handle such challenging optimization problems. In this work, we present a comprehensive investigation into the design, implementation, and evaluation of evolutionary algorithms for path planning of unmanned aerial vehicles. We first formulate the path planning problem as a constrained multi-objective optimization problem. Then, we introduce two novel evolutionary algorithms: an enhanced teaching-learning particle swarm optimization algorithm (E-TLPSO) for single-objective path planning, and a multi-objective evolutionary algorithm based on decomposition with a hybrid penalty boundary intersection mechanism (MOEADHPBI). Through extensive simulations on multiple benchmark functions and realistic UAV scenarios, we demonstrate that the proposed methods achieve superior convergence, diversity, and feasibility compared to state-of-the-art algorithms. Our results highlight the potential of evolutionary algorithms as robust and scalable tools for autonomous navigation of unmanned aerial vehicles in complex environments.
The remainder of this paper is organized as follows. Section 1 provides the background and a review of related work. Section 2 formulates the UAV path planning problem mathematically. Section 3 details the proposed single-objective algorithm E-TLPSO. Section 4 presents the multi-objective algorithm MOEADHPBI. Section 5 reports and discusses the experimental results. Finally, Section 6 concludes the paper and outlines future research directions.
1. Introduction
Unmanned aerial vehicles (UAVs), also known as drones, have become indispensable tools in modern society. Their low cost, high flexibility, and ability to operate in dangerous environments make them attractive for a variety of tasks, including topographic mapping, object detection, crop management, and disaster rescue. For any autonomous UAV system, one of the most critical components is the path planning module. Path planning for unmanned aerial vehicles aims to identify a collision-free, flyable, and near-optimal trajectory from a start point to a target point, while respecting the dynamic constraints of the aircraft and the environmental limitations.
The path planning problem becomes increasingly challenging when the environment is three-dimensional and contains both static obstacles (e.g., mountains, buildings) and dynamic threats (e.g., radar zones, no-fly zones). In such real-world scenarios, the planner must balance multiple conflicting objectives, such as minimizing path length, avoiding threats, maintaining safe altitude, and satisfying smoothness constraints. Moreover, for multiple UAVs operating in the same airspace, the problem becomes even more complex because the paths must not only be individually optimal but also mutually collision-free. These requirements turn UAV path planning into a highly constrained multi-objective optimization problem (MOP).
Over the years, researchers have proposed numerous algorithms for UAV path planning. These can be broadly categorized into classical graph-based methods, artificial potential field methods, deep reinforcement learning methods, and evolutionary algorithms. Graph-based methods, such as Dijkstra and A*, are well-established and easy to implement, but their computational cost scales poorly with the dimensionality of the search space. Artificial potential field methods are efficient for online planning, yet they are notorious for getting trapped in local minima. Deep reinforcement learning offers great adaptivity but requires substantial training data and computational resources, and its decision-making process is often opaque. Evolutionary algorithms, which mimic natural selection and collective behavior, have gained considerable attention due to their strong global search capability, flexibility in handling constraints, and ability to solve multi-objective problems directly.
Among evolutionary algorithms, particle swarm optimization (PSO) has been widely applied to UAV path planning because of its simple mechanism and rapid convergence. However, standard PSO often suffers from premature convergence and poor diversity when dealing with complex, high-dimensional search spaces. Teaching-learning-based optimization (TLBO) is another evolutionary method that simulates the teaching and learning processes in a classroom. TLBO emphasizes exploration through its teaching and learning phases. By combining the exploitation strength of PSO with the exploration strength of TLBO, one can develop a hybrid algorithm that inherits the advantages of both. This observation motivated us to propose an enhanced teaching-learning particle swarm optimization algorithm, which we refer to as E-TLPSO.
In addition, because UAV path planning naturally involves multiple conflicting objectives, a more appropriate formulation is to treat it as a multi-objective optimization problem. Multi-objective evolutionary algorithms (MOEAs) can simultaneously optimize several objectives and provide a set of trade-off solutions, known as the Pareto optimal set. Among various MOEAs, decomposition-based approaches, such as MOEA/D, have attracted significant interest. MOEA/D decomposes the multi-objective problem into a number of scalar subproblems and solves them cooperatively. A crucial element in such frameworks is the scalarizing function, e.g., the penalty boundary intersection (PBI) function. While PBI is effective in preserving diversity, it may suffer from slow convergence in certain cases. The inverted penalty boundary intersection (IPBI) method, on the other hand, accelerates convergence by maximizing a scalarized function, but it may yield poorly distributed solutions. To overcome these limitations, we design a multi-objective evolutionary algorithm based on decomposition with a hybrid PBI mechanism, called MOEADHPBI. The proposed algorithm adaptively switches between PBI and IPBI according to the state of the population, and further incorporates a constraint repair technique to ensure the feasibility of the final solutions.
The main contributions of this work are threefold:
- We establish a comprehensive mathematical model for UAV path planning that includes a three-dimensional environment representation, obstacle and threat models, flight constraints, and multiple conflicting objectives.
- We propose E-TLPSO, a single-objective evolutionary algorithm that combines PSO and TLBO with chaotic initialization, nonlinear adaptive inertia weight adjustment, and dynamic teaching-learning mechanisms. E-TLPSO is capable of solving both single-UAV and multi-UAV path planning problems.
- We propose MOEADHPBI, a multi-objective evolutionary algorithm that integrates a hybrid penalty boundary intersection strategy, adaptive penalty factor adjustment, an evolutionary operation mechanism, a constraint repair technique, and a flight route selection mechanism. The algorithm shows significant performance gains on both benchmark problems and realistic UAV path planning scenarios.
2. Mathematical Modeling for UAV Path Planning
In this section, we describe the mathematical model used throughout this paper. The model encompasses the three-dimensional flight space, terrain and radar threats, the constraints imposed by the UAV dynamics, and the objective functions to be optimized.
2.1 Flight Space and Path Representation
We consider a cuboid flight space \(S\) defined by:
\[
S = \{(x,y,z)\ |\ x_{\min} \le x \le x_{\max},\ y_{\min} \le y \le y_{\max},\ z_{\min} \le z \le z_{\max}\}.
\tag{1}
\]
A flight path is represented as a sequence of \(M\) waypoints:
\[
\mathrm{Path} = \{(x_i, y_i, z_i)\ |\ (x_i, y_i, z_i) \in S,\ i=1,2,\ldots,M\}.
\tag{2}
\]
The first and last waypoints correspond to the start and target positions, respectively.
2.2 Obstacle Model
Mountains are modeled as three-dimensional Gaussian surfaces. For the \(k\)-th mountain with center \((x_k^0, y_k^0)\) and shape parameters \(\delta_k, \sigma_k, \tau_k\), the elevation at a horizontal position \((\bar{x}, \bar{y})\) is:
\[
\bar{z} = \tau_k \exp\left\{-\left(\frac{\bar{x}-x_k^0}{\delta_k}\right)^2 – \left(\frac{\bar{y}-y_k^0}{\sigma_k}\right)^2\right\}.
\tag{3}
\]
A path point \( (x_i,y_i,z_i) \) is considered to be in collision with the mountain if \(z_i\) is lower than the terrain elevation at that horizontal location.
Radar threats are modeled as spherical danger zones. For the \(r\)-th radar centered at \((\tilde{x}_r^0,\tilde{y}_r^0,\tilde{z}_r^0)\) with radius \(d_r\), the danger zone is:
\[
DZ_r = \{(\tilde{x},\tilde{y},\tilde{z})\ |\ \sqrt{(\tilde{x}-\tilde{x}_r^0)^2 + (\tilde{y}-\tilde{y}_r^0)^2 + (\tilde{z}-\tilde{z}_r^0)^2} \le d_r\}.
\tag{4}
\]
Any path point inside a danger zone is considered unsafe.
2.3 Flight Constraints
To ensure the flight path is physically realizable, the UAV’s pitch angle and yaw angle must be kept within acceptable ranges. For a path segment from \((x_i,y_i,z_i)\) to \((x_{i+1},y_{i+1},z_{i+1})\), the pitch angle \(\alpha_i\) and yaw angle \(\eta_i\) are computed as:
\[
\alpha_i = \arctan\left(\frac{z_{i+1}-z_i}{\sqrt{(x_{i+1}-x_i)^2+(y_{i+1}-y_i)^2}}\right),
\tag{5}
\]
\[
\eta_i = \arctan\left(\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\right).
\tag{6}
\]
We require \(\alpha_{\min} \le \alpha_i \le \alpha_{\max}\) and \(\eta_{\min} \le \eta_i \le \eta_{\max}\). Violations are counted as penalties.
2.4 Objective Functions
We consider two primary objectives: path length and radar threat. The total path length is:
\[
f_{\mathrm{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{7}
\]
To quantify radar threat, we sample \(m\) points along each path segment and check whether they fall inside any radar danger zone. The threat objective is:
\[
f_{\mathrm{threat}} = \sum_{i=1}^{M-1}\sum_{j=1}^{m} G_{i,j},
\tag{8}
\]
where \(G_{i,j}=1\) if the \(j\)-th sample of the \(i\)-th segment is inside a radar zone, and \(0\) otherwise.
Additionally, constraints related to flight height and angle consistency can be translated into penalty terms. For the single-objective formulation, we construct a weighted sum:
\[
F(\mathrm{Path}) = \rho_1 f_{\mathrm{height}} + \rho_2 f_{\mathrm{pitch}} + \rho_3 f_{\mathrm{yaw}} + \rho_4 f_{\mathrm{length}} + \rho_5 f_{\mathrm{threat}},
\tag{9}
\]
where \(\rho_i\) are user-defined weights. For the multi-objective formulation, we directly optimize the two objectives subject to a constraint on the total angle violations:
\[
\min F(\mathrm{Path}) = (f_{\mathrm{length}}(\mathrm{Path}), f_{\mathrm{threat}}(\mathrm{Path})) \quad \text{s.t.} \quad f_{\mathrm{pitch}} + f_{\mathrm{yaw}} \le \varepsilon.
\tag{10}
\]
This model provides a realistic trade-off between flight efficiency and operational safety, which we exploit in the design of our algorithms.
3. Enhanced Teaching-Learning Particle Swarm Optimization for Single-Objective UAV Path Planning
In this section, we present the E-TLPSO algorithm, which is specifically designed to solve the single-objective UAV path planning problem formulated in Equation (9). E-TLPSO integrates the exploitation capability of particle swarm optimization with the exploration capability of teaching-learning-based optimization. We describe the three main components: chaotic population initialization, dynamic teaching-learning strategy, and adaptive evolution with a nonlinear inertia weight.
3.1 Standard PSO and TLBO
In standard PSO, each particle \(i\) has a position vector \(X_i\) and a velocity vector \(V_i\). The velocity update at iteration \(t\) is:
\[
V_i(t+1) = \omega V_i(t) + c_1 r_1 (pbest_i(t) – X_i(t)) + c_2 r_2 (gbest(t) – X_i(t)),
\tag{11}
\]
and the position update is:
\[
X_i(t+1) = X_i(t) + V_i(t+1).
\tag{12}
\]
Here, \(\omega\) is the inertia weight, \(c_1\) and \(c_2\) are learning factors, \(r_1,r_2\) are random numbers in [0,1), and \(pbest_i, gbest\) are the personal and global best solutions, respectively.
In standard TLBO, the population evolves through a teaching phase and a learning phase. During the teaching phase, each individual learns from the best individual (teacher) and moves towards the mean of the population:
\[
X_{\mathrm{teach},i}(t) = X_i(t) + r_3 \left(X_{\mathrm{teacher}}(t) – T_f X_{\mathrm{center}}(t)\right),
\tag{13}
\]
where \(T_f\) is a teaching factor. During the learning phase, each individual learns from a randomly selected better individual:
\[
X_{\mathrm{learn},i}(t) = X_i(t) + r_4 \left(X_u(t) – X_v(t)\right) \quad \text{or} \quad X_i(t) + r_5 \left(X_v(t) – X_u(t)\right).
\tag{14}
\]
3.2 Proposed E-TLPSO Algorithm
We observed that standard PSO converges quickly but may lose diversity, while TLBO provides diversity but can be slow in refining solutions. The hybrid E-TLPSO combines the strengths of both. The algorithm proceeds as follows.
3.2.1 Chaotic Population Initialization
Instead of using random sequences, we adopt a chaotic map to generate initial positions. The chaotic factor is produced by:
\[
\gamma_i = \mu \gamma_{i-1} (1-\gamma_{i-1}), \quad \gamma_0 \in [0,1), \gamma_0 \notin \{0.25,0.5,0.75\},
\tag{15}
\]
with \(\mu = 4\). The \(j\)-th dimension of the initial position is then:
\[
X_{i,j}(0) = X_{\min,j} + \gamma_i (X_{\max,j} – X_{\min,j}).
\tag{16}
\]
Chaotic initialization enhances population diversity and can accelerate the convergence of the algorithm.
3.2.2 Dynamic Teaching and Learning
To further improve exploration, we replace the random numbers \(r_3, r_4, r_5\) in TLBO with chaotic factors generated by Equation (15). This yields the following dynamic teaching and learning rules:
\[
X_{\mathrm{teach},i}(t) = X_i(t) + \vartheta_i \left(X_{\mathrm{teacher}}(t) – T_f X_{\mathrm{center}}(t)\right),
\tag{17}
\]
\[
X_{\mathrm{learn},i}(t) = X_i(t) + \vartheta_i (X_u(t) – X_v(t)) \quad \text{or} \quad X_i(t) + \vartheta_i (X_v(t) – X_u(t)),
\tag{18}
\]
where \(\vartheta_i\) is a chaotic factor assigned to individual \(i\). This dynamic mechanism helps the algorithm escape local optima and promotes balanced exploration.
3.2.3 Adaptive Evolution with Nonlinear Inertia Weight
After the teaching-learning phase, we apply the PSO velocity and position updates. To adapt the inertia weight to the current search state, we propose a nonlinear adjustment strategy:
\[
\omega = (\omega_{\max} – \omega_{\min}) \times \frac{N_s}{N} + \zeta \omega_{\min},
\tag{19}
\]
where \(N_s\) denotes the number of particles that have updated their personal best in the current iteration, \(N\) is the population size, and \(\zeta\) is a chaotic factor. When many particles improve, the population is likely exploring promising regions; increasing \(\omega\) accelerates exploration. Conversely, when improvements are fewer, the population tends to exploit the best region, and a smaller \(\omega\) enhances fine local search.
The learning factors \(c_1\) and \(c_2\) are adjusted linearly over time:
\[
c_1 = (c_{\min} – c_{\max}) \frac{t}{t_{\max}} + c_{\max},
\tag{20}
\]
\[
c_2 = (c_{\max} – c_{\min}) \frac{t}{t_{\max}} + c_{\min}.
\tag{21}
\]
Such an adjustment balances global exploration and local exploitation throughout the optimization process.
The overall E-TLPSO procedure is summarized in Algorithm 1. The algorithm can handle both single-UAV and multi-UAV path planning by encoding multiple paths in the individual representation.
| Step | Description |
|---|---|
| 1 | Initialize population using chaotic map (Equation (16)). |
| 2 | Evaluate fitness of each individual using Equation (9). |
| 3 | Set pbest and gbest. |
| 4 | For \(t=1\) to \(t_{\max}\): |
| 5 | Update teacher and center of population. |
| 6 | Adjust \(c_1, c_2\) using Equations (20)-(21). |
| 7 | For each individual \(i\): |
| 8 | Perform dynamic teaching using Equation (17). |
| 9 | If improved, update \(X_i\); update gbest. |
| 10 | Perform dynamic learning using Equation (18). |
| 11 | If improved, update \(X_i\); update gbest. |
| 12 | Update velocity and position using PSO rules. |
| 13 | Update pbest and gbest; increment \(N_s\) if improved. |
| 14 | End for. |
| 15 | Update inertia weight using Equation (19). |
| 16 | End for. |
| 17 | Return gbest. |
For multi-UAV path planning, an individual is represented by a two-dimensional array where each row corresponds to the waypoint sequence of a particular UAV. The fitness of an individual equals the sum of costs over all UAV paths. Collision avoidance between UAVs is handled via a repair mechanism that adjusts conflicting segments until a minimum safety distance is satisfied.
4. Multi-Objective Evolutionary Algorithm Based on Decomposition with Hybrid Penalty Boundary Intersection
In many real-world missions, minimizing path length alone is insufficient because the unmanned aerial vehicle must also avoid threats. In this section, we consider the multi-objective formulation in Equation (10) and introduce a decomposition-based multi-objective evolutionary algorithm called MOEADHPBI.
4.1 Decomposition Framework
The algorithm starts by generating \(N\) uniformly spread weight vectors \(\lambda = \{\lambda_1,\ldots,\lambda_N\}\), each corresponding to a subproblem. For each weight vector \(\lambda_i\), its neighborhood set \(\text{Nes}(i)\) contains the indices of the \(T\) closest weight vectors in terms of Euclidean distance. The population is initialized randomly, and the ideal point \(Z^+\) and the nadir point \(Z^-\) are computed from the initial objective values.
At each iteration, for each subproblem \(i\), we generate an offspring solution using evolutionary operations (e.g., simulated binary crossover and polynomial mutation). The offspring is then used to update the ideal point and the population in the neighborhood or the whole population according to a selection probability. After a fixed number of iterations, a constraint repair technique is applied to ensure feasibility. Finally, a flight route selection mechanism is used to choose the most suitable path from the obtained Pareto set.
4.2 Hybrid Penalty Boundary Intersection
The traditional PBI scalarizing function is defined as:
\[
\min \quad g_{\text{pbi}}(x|\lambda, Z^+) = d_1 + \theta_1 d_2,
\]
where
\[
d_1 = \frac{\|(Z^+ – F(x))^T \lambda\|}{\|\lambda\|}, \quad d_2 = \|F(x) – (Z^+ – d_1 \lambda)\|.
\tag{22}
\]
PBI tends to produce well-distributed solutions but may slow down convergence.
The inverted PBI (IPBI) function is defined as:
\[
\max \quad g_{\text{ipbi}}(x|\lambda, Z^-) = d_1 – \theta_2 d_2,
\]
with
\[
d_1 = \frac{\|(Z^- – F(x))^T \lambda\|}{\|\lambda\|}, \quad d_2 = \|(Z^- – F(x)) – d_1 \lambda\|.
\tag{23}
\]
IPBI accelerates convergence by pulling solutions outward from the nadir point, but may sacrifice distribution uniformity.
We introduce a hybrid strategy that dynamically chooses between PBI and IPBI based on the current state of the population. Let \(\sigma_t\) be the standard deviation of the population objective vectors, and \(\delta_t\) the average distance of the population to \(Z^+\). We define the diversity change rate \(r_\sigma = \sigma_t/\sigma_{t-1}\) and the convergence change rate \(r_\delta = \delta_t/\delta_{t-1}\). The ratio \(r_\tau = r_\delta / r_\sigma\) indicates whether the population is improving more in convergence or diversity. A sigmoid function is used to compute the selection probability:
\[
r_\partial = \frac{1}{1 + e^{-(1-r_\tau)}}.
\tag{24}
\]
When \(r_\partial\) is larger than a threshold \(r_\partial^*\), we use PBI to enhance diversity; otherwise, we use IPBI to enhance convergence. The penalty factors \(\theta_1\) and \(\theta_2\) are adaptively tuned according to:
\[
\theta_1 = \frac{1}{1+\exp(-\varphi t/t_{\max})} \gamma_t + (\theta_{\max} – \theta_{\min}) \left(\frac{t}{t_{\max}}\right)^2,
\tag{25}
\]
\[
\theta_2 = \theta_{\max} \gamma_t + (\theta_{\min} – \theta_{\max})\frac{t}{t_{\max}},
\tag{26}
\]
where \(\gamma_t\) is a chaotic factor generated as in Equation (15). This adaptive adjustment helps balance convergence and diversity at different stages of the search.
4.3 Evolutionary Operation
For each subproblem, we generate an offspring solution using a local or global selection strategy. A random number is compared to a threshold \(\delta^*\). If the random number is less than or equal to \(\delta^*\), we select the neighborhood set \(\text{Nes}(i)\) as the mating pool; otherwise, we use the entire population. From the selected pool, two parents are chosen and simulated binary crossover (SBX) and polynomial mutation are applied to create an offspring. This operation encourages local refinement while preserving global exploration.
4.4 Constraint Repair Technique
Because the path planning problem contains angle constraints, many generated paths may violate the pitch or yaw limits. To transform infeasible solutions into feasible ones without severely degrading other objectives, we design a repair technique. For each violated angle, we clamp it to the nearest permissible bound:
\[
\alpha_{\text{new},i} = \max(\alpha_{\min}, \min(\alpha_i, \alpha_{\max})),
\quad
\eta_{\text{new},i} = \max(\eta_{\min}, \min(\eta_i, \eta_{\max})).
\tag{27}
\]
Then we recalculate the waypoint coordinates using the repaired angles:
\[
\begin{aligned}
x_i^* &= \max(x_{\min}, \min(x_{\max}, x_i + DS_i \cos(\eta_{\text{new},i}))) \\
y_i^* &= \max(y_{\min}, \min(y_{\max}, y_i + DS_i \sin(\eta_{\text{new},i}))) \\
z_i^* &= \max(z_{\min}, \min(z_{\max}, z_i + DS_i \tan(\alpha_{\text{new},i})))
\end{aligned}
\tag{28}
\]
where \(DS_i = \sqrt{(x_{i+1}-x_i)^2 + (y_{i+1}-y_i)^2}\). This repair technique is applied periodically to the current population to ensure that the final Pareto optimal set is feasible.
4.5 Flight Route Selection Mechanism
Once the algorithm terminates, we obtain a set of Pareto-optimal solutions. To select a single route for deployment, we consider the inherent priority of avoiding radar threats over path length. Our selection mechanism first sorts the solutions by \(f_{\text{threat}}\) in ascending order. If multiple solutions have the same threat value, we choose the one with the smallest \(f_{\text{length}}\). If both values are identical, we randomly select one. This conservative strategy guarantees a safe flight path, which is critical in many operational scenarios.
The complete MOEADHPBI procedure is outlined in Algorithm 2.
| Step | Description |
|---|---|
| 1 | Generate \(N\) weight vectors and compute neighborhoods. |
| 2 | Initialize population \(Pop\) and evaluate objectives. |
| 3 | Initialize ideal point \(Z^+\) and nadir point \(Z^-\). |
| 4 | For \(t=1\) to \(t_{\max}\): |
| 5 | Update penalty factors \(\theta_1,\theta_2\) using Equations (25)-(26). |
| 6 | Update \(Z^-\) and compute diversity/convergence rates. |
| 7 | For each subproblem \(i\): |
| 8 | Call EvolutionaryOperation to generate an offspring. |
| 9 | Update ideal point \(Z^+\). |
| 10 | Call AdaptiveHPBI to update the population. |
| 11 | End for. |
| 12 | If \(t \bmod t_{\text{mod}} == 0\), apply ConstraintRepair. |
| 13 | End for. |
| 14 | Return the best safe path selected from the final Pareto set. |
The computational complexity of MOEADHPBI is \(O(t_{\max} N^2)\), which is comparable to that of standard decomposition-based MOEAs and is acceptable for offline path planning applications.
5. Experimental Results and Analysis
We conducted comprehensive experiments to evaluate the performance of the proposed algorithms. All simulations were performed in MATLAB 2023a on a Windows 10 PC with an Intel i5 processor at 3.10 GHz and 32 GB of RAM. The experimental design includes (i) benchmark function tests to validate general optimization capability, (ii) single-UAV path planning tests, and (iii) multi-UAV path planning tests.
5.1 Benchmark Function Tests
We first evaluated E-TLPSO on 15 classical benchmark functions in 10 and 50 dimensions. The competing algorithms included several PSO variants, genetic algorithm, and recent metaheuristics. Table 1 summarizes the number of benchmark functions for which each algorithm obtained the best average value.
| Algorithm | Mean | Max | Min | Std |
|---|---|---|---|---|
| AIWPSO | 3 | 2 | 3 | 2 |
| TCSPSO | 2 | 2 | 3 | 3 |
| MPSO | 0 | 0 | 0 | 1 |
| AO | 0 | 0 | 0 | 0 |
| DE | 0 | 0 | 0 | 2 |
| GA | 0 | 0 | 0 | 0 |
| HBA | 4 | 5 | 1 | 2 |
| E-TLPSO | 6 | 6 | 8 | 5 |
In 50D, E-TLPSO achieved the best results on a larger number of functions than any other algorithm, especially in terms of minimum values and standard deviation, indicating excellent convergence and stability.
5.2 Evaluation of MOEADHPBI on Benchmarks
We compared MOEADHPBI with eight state-of-the-art MOEAs on ZDT, DTLZ, MaF, and WFG test suites. The population size was set to 100 and the maximum number of function evaluations was set to 50,000 for bi-objective and 100,000 for tri-objective problems. We report performance using the Inverted Generational Distance (IGD) and Hypervolume (HV) metrics. The IGD is defined as:
\[
\mathrm{IGD} = \frac{1}{c}\sum_{i=1}^{c} d_i,
\tag{29}
\]
where \(c\) is the number of true Pareto front points and \(d_i\) is the Euclidean distance from the \(i\)-th point to the nearest obtained solution. The HV metric measures the volume of the objective space dominated by the obtained solution set and is computed with respect to a reference point.
Table 2 shows the summarized Wilcoxon rank-sum test results on the ZDT and DTLZ suites. The symbol ‘+’ indicates that the corresponding comparison algorithm outperforms MOEADHPBI, ‘−’ indicates that MOEADHPBI is significantly better, and ‘≈’ indicates no significant difference.
| Metric | EAGMOEAD | FLEA | IMMOEAD | MOEADPaS | MOBCA | MaOEAIGD | CoMMEA | HREA |
|---|---|---|---|---|---|---|---|---|
| HV (win/loss) | 3/2 | 4/1 | 5/0 | 4/1 | 4/1 | 5/0 | 5/0 | 5/0 |
| IGD (win/loss) | 4/1 | 4/1 | 5/0 | 4/1 | 4/1 | 5/0 | 5/0 | 5/0 |
Our proposed MOEADHPBI showed significant improvements in both HV and IGD over most comparison algorithms. This is attributed to the hybrid scalarizing function that dynamically adjusts convergence and diversity, and to the adaptive penalty factor, which effectively guides the search across the entire Pareto front.
5.3 Single-UAV Path Planning Results
We evaluated E-TLPSO on three simulated terrains of increasing difficulty: Scenario 1 (simple terrain with low mountains and few radars), Scenario 2 (moderate complexity), and Scenario 3 (many high mountains and dense radar zones). The starting point and target point were fixed at (0,0,5) and (60,60,5), respectively. We compared E-TLPSO against seven algorithms: HCPSOA, CIPSO, AWPSO, HPSO, TVAPSO, GA, and HBA. Each algorithm ran 18 independent times with 100 iterations per run.
Table 3 reports the average path length (APL) and the final average optimal fitness value (FAOFV) of each algorithm in Scenarios 1-3 for single-UAV path planning.
| Scenario | Metric | HCPSOA | CIPSO | AWPSO | HPSO | TVAPSO | GA | HBA | E-TLPSO |
|---|---|---|---|---|---|---|---|---|---|
| Scenario 1 | APL (m) | 72.128 | 72.128 | 72.128 | 72.128 | 72.129 | 72.130 | 72.130 | 72.126 |
| FAOFV | 2163.85 | 2163.85 | 2163.83 | 2163.85 | 2163.87 | 2163.90 | 2163.93 | 2163.79 | |
| Scenario 2 | APL (m) | 74.47 | 74.33 | 74.00 | 74.18 | 74.28 | 74.54 | 74.43 | 72.95 |
| FAOFV | 2234.20 | 2229.76 | 2220.08 | 2225.38 | 2228.36 | 2236.22 | 2232.98 | 2188.46 | |
| Scenario 3 | APL (m) | 78.35 | 78.85 | 77.56 | 77.84 | 77.44 | 80.10 | 80.17 | 75.64 |
| FAOFV | 2347.94 | 2351.04 | 2333.69 | 2342.77 | 2340.31 | 2402.86 | 2405.09 | 2269.30 |
From Table 3, we observe that in all three scenarios, E-TLPSO produced the shortest average path and the lowest FAOFV among all compared algorithms. In Scenario 3, which is the most complex, E-TLPSO reduced the FAOFV by about 5.6% compared to CIPSO, demonstrating its robustness in difficult environments. The convergence curves (not shown for brevity) also indicated that E-TLPSO converged faster and reached a lower global optimum than the other algorithms.
5.4 Multi-UAV Path Planning Results
We further tested E-TLPSO on multi-UAV path planning scenarios with a formation of five UAVs. Two scenarios were used: Scenario 4 (terrain with mountains but no radar threats) and Scenario 5 (both mountains and radar threats). The maximum iteration number was set to 300. The total cost was the sum of the path costs of all five UAVs, and a collision repair mechanism was applied for all algorithms.
Table 4 lists the average path length (AL) over the five UAVs and the FAOFV obtained by each algorithm in Scenarios 4 and 5.
| Scenario | Metric | AWPSO | CIPSO | HCPSOA | HPSO | TVAPSO | GA | HBA | E-TLPSO |
|---|---|---|---|---|---|---|---|---|---|
| Scenario 4 | AL (m) | 83.17 | 81.01 | 80.67 | 80.01 | 79.13 | 81.79 | 82.41 | 76.51 |
| FAOFV | 17260.37 | 13630.49 | 13802.22 | 10549.68 | 14217.42 | 15102.15 | 16525.80 | 10272.03 | |
| Scenario 5 | AL (m) | 81.26 | 81.54 | 81.16 | 78.93 | 78.68 | 81.99 | 82.11 | 74.47 |
| FAOFV | 15271.70 | 16151.08 | 13125.33 | 11327.14 | 20790.51 | 15044.57 | 17343.07 | 11165.30 |
E-TLPSO achieved the shortest AL in both scenarios. In Scenario 4, its FAOFV was 37.82% lower than that of HBA, and in Scenario 5, it was 46.15% lower than that of TVAPSO. These results demonstrate the scalability of E-TLPSO for multi-UAV cooperative path planning.
5.5 Performance of MOEADHPBI on UAV Path Planning
We assessed MOEADHPBI on four UAV scenarios, denoted as Scenarios 6-9, with diverse terrain and radar threats. We compared it with eight constrained multi-objective evolutionary algorithms: BiCo, CCMO, cDPEA, CMOPSO, MCCMO, MSCMO, DPPPS, and MOEA2DE. Each algorithm ran 30 independent times, and we report the averaged HV and IGD values.
Table 5 shows the results for Scenarios 6 and 7, and Table 6 shows the results for Scenarios 8 and 9. The best result in each column is highlighted in bold.
| Algorithm | Scenario 6 | Scenario 7 | ||
|---|---|---|---|---|
| HV | IGD | HV | IGD | |
| BiCo | 1.67e-1 (±1.31e-1) | 2.46e+1 (±1.72e+1) | 2.11e-1 (±1.23e-1) | 1.91e+1 (±1.76e+1) |
| CCMO | 2.68e-1 (±5.08e-2) | 1.05e+1 (±5.59e+0) | 2.68e-1 (±5.25e-2) | 1.05e+1 (±5.76e+0) |
| cDPEA | 3.02e-1 (±3.65e-2) | 6.81e+0 (±3.99e+0) | 3.15e-1 (±3.09e-2) | 5.51e+0 (±3.35e+0) |
| CMOPSO | 2.60e-1 (±1.67e-2) | 1.15e+1 (±1.84e+0) | 2.56e-1 (±1.93e-2) | 1.18e+1 (±2.12e+0) |
| MCCMO | 2.70e-1 (±5.48e-2) | 1.03e+1 (±6.00e+0) | 2.67e-1 (±5.38e-2) | 1.07e+1 (±5.91e+0) |
| MSCMO | 2.67e-1 (±5.58e-2) | 1.07e+1 (±6.11e+0) | 2.62e-1 (±5.20e-2) | 1.12e+1 (±5.69e+0) |
| DPPPS | 3.01e-1 (±3.33e-2) | 7.08e+0 (±3.53e+0) | 2.98e-1 (±2.83e-2) | 7.51e+0 (±3.05e+0) |
| MOEA2DE | 2.03e-1 (±4.31e-3) | 1.02e+1 (±7.07e-1) | 2.01e-1 (±2.94e-3) | 1.03e+1 (±7.40e-1) |
| MOEADHPBI | 6.12e+2 (±8.32e+0) | 1.23e-1 (±6.43e-2) | 6.59e+2 (±6.60e+0) | 1.35e-1 (±6.70e-2) |
| Algorithm | Scenario 8 | Scenario 9 | ||
|---|---|---|---|---|
| HV | IGD | HV | IGD | |
| BiCo | 1.75e-1 (±1.08e-1) | 2.61e+1 (±1.19e+1) | 1.70e-1 (±1.00e-1) | 2.62e+1 (±1.20e+1) |
| CCMO | 2.11e-1 (±4.98e-2) | 1.77e+1 (±5.05e+0) | 1.84e-1 (±6.00e-2) | 2.03e+1 (±6.22e+0) |
| cDPEA | 2.38e-1 (±4.52e-2) | 1.59e+1 (±4.47e+0) | 2.23e-1 (±4.70e-2) | 1.72e+1 (±4.19e+0) |
| CMOPSO | 2.27e-1 (±1.84e-2) | 1.59e+1 (±1.93e+0) | 2.16e-1 (±1.51e-2) | 1.71e+1 (±1.54e+0) |
| MCCMO | 2.00e-1 (±4.66e-2) | 1.87e+1 (±4.72e+0) | 1.98e-1 (±6.92e-2) | 2.00e+1 (±6.50e+0) |
| MSCMO | 1.99e-1 (±5.45e-2) | 1.91e+1 (±5.44e+0) | 1.58e-1 (±7.52e-2) | 2.59e+1 (±1.49e+1) |
| DPPPS | 2.39e-1 (±5.29e-2) | 1.31e+1 (±6.14e+0) | 2.23e-1 (±2.71e-2) | 1.68e+1 (±2.29e+0) |
| MOEA2DE | 9.24e-2 (±4.10e-3) | 9.46e+0 (±3.20e-14) | 9.23e-2 (±4.21e-3) | 1.03e+1 (±4.62e-14) |
| MOEADHPBI | 7.79e+2 (±1.71e-1) | 1.17e-1 (±6.07e-2) | 8.16e+2 (±5.12e-12) | 1.11e-1 (±6.24e-2) |
It is evident that MOEADHPBI dramatically outperforms all comparison algorithms in terms of both HV and IGD. For instance, in Scenario 8, the HV value of MOEADHPBI is more than 3000 times larger than that of the best comparison algorithm (DPPPS), and its IGD value is about 1.2% of the best comparison algorithm’s value. This outstanding performance stems from the hybrid penalty boundary intersection strategy, which effectively balances convergence and diversity, and from the constraint repair technique that ensures all returned paths are feasible.
5.6 Ablation Study
To investigate the contribution of each component in E-TLPSO, we conducted an ablation study. Three variants were introduced:
- E-TLPSO-C: replaces chaotic initialization with random initialization.
- E-TLPSO-L: replaces the nonlinear inertia weight with a linear time-varying one.
- E-TLPSO-D: removes the dynamic teaching-learning strategy.
Table 7 reports the FAOFV and average path length/AL for all variants in five UAV scenarios.
| 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 (m) | 72.158 | 72.128 | 72.134 | 72.126 | |
| Scenario 2 | FAOFV | 2191.33 | 2200.06 | 2190.30 | 2188.46 |
| APL (m) | 73.04 | 73.34 | 73.01 | 72.95 | |
| Scenario 3 | FAOFV | 2325.39 | 2349.28 | 2312.54 | 2269.30 |
| APL (m) | 77.51 | 78.31 | 77.09 | 75.64 | |
| Scenario 4 | FAOFV | 10531.48 | 14415.20 | 11348.25 | 10272.03 |
| AL (m) | 76.22 | 81.24 | 76.35 | 76.51 | |
| Scenario 5 | FAOFV | 11246.80 | 12656.98 | 11245.43 | 11165.30 |
| AL (m) | 74.93 | 80.86 | 74.87 | 74.47 |
In all scenarios, the complete E-TLPSO achieves the lowest FAOFV. In complex scenarios (e.g., Scenario 3 and Scenario 5), the differences are pronounced. The dynamic teaching-learning strategy and nonlinear inertia weight play critical roles in maintaining search diversity and preventing premature convergence.
6. Conclusion and Future Work
In this paper, we have presented a comprehensive study on evolutionary algorithms for path planning of unmanned aerial vehicles. We formulated the UAV path planning problem as a constrained multi-objective optimization problem and developed two novel algorithms: E-TLPSO for single-objective path planning and MOEADHPBI for multi-objective path planning.
E-TLPSO integrates particle swarm optimization and teaching-learning-based optimization with chaotic initialization, dynamic teaching-learning, and adaptive nonlinear inertia weight adjustment. It demonstrated superior convergence and solution quality for both single-UAV and multi-UAV path planning tasks in various simulated environments. The algorithm’s ability to escape local optima and maintain diversity makes it particularly effective in complex terrain with dense obstacles and threats.
MOEADHPBI employs a hybrid penalty boundary intersection strategy that dynamically switches between PBI and IPBI based on the population’s state. The adaptive penalty factor and constraint repair technique ensure that the final Pareto set is both well-distributed and feasible. Extensive experiments on benchmark functions and realistic UAV scenarios showed that MOEADHPBI significantly outperforms eight state-of-the-art constrained multi-objective evolutionary algorithms in terms of HV and IGD. In particular, the algorithm achieves orders-of-magnitude improvement in HV and IGD metrics in the tested UAV path planning scenarios.
There are several promising directions for future work. First, we intend to validate the proposed algorithms on real physical UAV platforms using the DJI Mobile SDK and Onboard SDK. This will require developing middleware to bridge the gap between simulation and real-world deployment, including handling sensor noise and communication delays. Second, we plan to investigate multi-scale search strategies that can better balance global exploration and local refinement in high-dimensional path planning problems. Third, we aim to extend the current framework to cooperative path planning of heterogeneous UAV swarms with communication constraints. Moreover, the computational efficiency of the proposed algorithms can be improved by leveraging parallel computing techniques, such as GPU acceleration, to enable real-time path replanning in dynamic environments. We believe that the work presented here provides a solid foundation for these future advancements in autonomous navigation of unmanned aerial vehicles.