Multi-UAV Formation Path Planning Based on Compensation Look-Ahead Algorithm

Abstract

This study primarily studies the shortest-path planning problem for unmanned aerial vehicle (UAV) formations under uncertain target sequences. In order to enhance the efficiency of collaborative search in drone clusters, a compensation look-ahead algorithm based on optimizing the four-point heading angles is proposed. Building upon the receding-horizon algorithm, this method introduces the heading angles of adjacent points to approximately compensate and decouple the triangular equations of the optimal trajectory, and a general formula for calculating the heading angles is proposed. The simulation data indicate that the model using the compensatory look forward algorithm exhibits a maximum improvement of 12.9% compared to other algorithms. Furthermore, to solve the computational complexity and sample size requirements for optimal solutions in the Dubins multiple traveling salesman model, a path-planning model for multiple UAV formations is introduced based on the Euclidean traveling salesman problem (ETSP) pre-allocation. By pre-allocating sub-goals, the model reduces the computational scale of individual samples while maintaining a constant sample size. The simulation results show an 8.4% and 17.5% improvement in sparse regions for the proposed Euclidean Dubins traveling salesman problem (EDTSP) model for takeoff from different points.

1. Introduction

Unmanned aerial vehicles are extensively used in collaborative reconnaissance, joint strike operations, aerial displays, express logistics, and agricultural spraying, among other domains. This has consequently elevated the control systems for multi-UAV (unmanned aerial vehicle) clusters, denoted as unmanned aerial systems (UAS), to the focal point of current international and domestic research [1,2]. Within the field of UAV research, path planning assumes a position of paramount importance. A recurrent inquiry pertains to the determination of target sequences and the acquisition of the most expeditious flight path or time [3]. This problem bears resemblance to the variant of the traveling salesman problem (TSP) known as the multiple traveling salesman problem (MTSP), which is conventionally studied but harbors certain distinctions [4]. In the research of MTSP, the paths between target points adhere to a geometrically straight distribution, whereas the motion of UAVs adheres to curved trajectories constrained by maximum curvature [5]. This is particularly pronounced in the case of fixed-wing UAVs, where trajectory curvature is bounded, demanding a curvature configuration in the model that aligns more closely with practical requirements [6,7]. Consequently, researchers frequently substitute the path planning for UAVs with Dubins curves [8]. Despite the approximate similarity between UAV path planning and MTSP, the former poses a more formidable challenge due to kinematic conditions such as curvature constraints, necessitating the formulation of novel solution methodologies by researchers [9].

We combined the Dubins problem with the multiple traveling salesman problem (MTSP), coining it the Dubins traveling salesman problem (DTSP). Solving the DTSP entails not only determining the target sequence but also considering that distinct heading angles yield divergent path-planning outcomes. Consequently, the determination of heading angles must precede the establishment of the target sequence [10]. Using a greedy algorithm for the brute-force solution of the DTSP is imprudent, as the computational burden exponentially escalates with the augmentation of target points [11]. Therefore, resolving a pre-defined target sequence remains a common strategy for addressing curvature-constrained path planning. Ma et al. (2007) proffered a receding-horizon algorithm grounded in the retrograde principle, computing a smooth trajectory through the sequence of path points; however, the authors have not yet developed a complete general formula that can characterize the three path points and initial directions of all possible geometric configurations, so as to uniquely determine the optimal trajectory for any problem [12]. Pantelis Isaiah et al. (2015) synergized k-step anticipation with the retrograde principle and disrupted the target sequence by introducing a k-step foresight algorithm; however, this algorithm is only applicable to situations where there are fewer target points [13]. Meanwhile, Izack Cohen et al. (2017) introduced a novel class of discretized look-ahead algorithms (DLAAs), discretizing the heading orientation angles (HOAs) at various waypoints and formulating an integer optimization problem. The resolution is then ascertained using an integer programming solution, but this algorithm suffers from high computational complexity and low computational efficiency in large-scale formation calculations [14].

This article mainly studies the problem of path planning for fixed-wing unmanned aerial vehicle formation. Firstly, we propose a heading angle calculation algorithm based on four-point optimization and called the heading angle compensation look-ahead algorithm. Building upon the receding-horizon algorithm, we introduce the influence of the heading angle of any adjacent point, be it antecedent or subsequent, upon the selection of the intermediary point and propose a general formula for calculating the heading angle. This algorithm is compared with other algorithms, such as the 2-opt algorithm [15], the alternating algorithm (AA) [16], and the Euclidean distance traveling salesman problem (ETSP). Secondly, we have adopted a pre-allocation mechanism to address the computational problem of path planning in large-scale drone formations. Within the context of multiple UAV formations, we use the ETSP algorithm for the pre-allocation of target points, thereby creating several subsets of designated points. Subsequently, we apply the compensatory foresight algorithm introduced in the second section to each subset, striking a balance between computational scale and the optimization of pathways.

2. Problem Description and Necessary Theories

The typical scenario of UAV swarm mission path planning depicted in Figure 1 in this paper involves multiple fixed-wing UAVs conducting point-based search and reconnaissance within a designated area (${\mathrm{N}}_{\mathrm{u}\mathrm{a}\mathrm{v}}\le {\mathrm{N}}_{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}$). Initially, after obtaining coordinates of various target points scattered within the area via satellite reconnaissance or other means, the command center assigns these task points to the leader UAV based on the distribution of the target points. The leader UAV then plans paths for all target points, considering the maneuverability of the UAVs. Subsequently, the planned task paths are distributed to each follower UAV [17]. Following the assigned task paths, multiple UAVs will conduct a search of all target points.

However, constrained by the size of the target area and the maneuverability of fixed-wing UAVs themselves, in practical applications, UAV formations often encounter scenarios where they must fly around dense distributions of target points. This requires UAVs to perform extensive curved maneuvers to reach the target point coordinates, significantly reducing the efficiency of the UAV formation’s search and increasing the search time. Therefore, we recommend determining the sequence of target points and the heading angle for each UAV towards every target point before path planning to optimize the UAV formation’s path planning.

This remains challenging for various reasons. Academically, Dubins curves are commonly used to simulate the shortest paths between target points [18]. However, the inherent uncertainty of Dubins curves results in different paths between target points based on different angles, meaning there are countless possible paths between these points. Therefore, there is a need to design a new algorithm that can provide a general formula for calculating heading angles to determine the optimal heading angle for each target point. Additionally, in large-scale UAV formations, as the number of target points increases, the computational complexity of the model increases exponentially. Therefore, efficient path planning while maintaining manageable computational complexity is also a critical consideration.

2.1. Dubins Path

The trajectory of fixed-wing unmanned aerial vehicles is a linear variation, delineating its trajectory as a curve [19].

$$\mathrm{X}\left(\mathrm{t}\right)={\left[\mathrm{x}\left(\mathrm{t}\right),\mathrm{y}\left(\mathrm{t}\right),\mathsf{\theta }\left(\mathrm{t}\right)\right]}^{\mathrm{T}}\in {\mathrm{R}}^2\times {\mathrm{S}}^1$$

(1)

This paper focuses solely on the motion of unmanned aerial vehicles in a two-dimensional environment, where x and y represent the Cartesian coordinates of the drone’s position, and θ represents the angle of the UAV relative to the X-axis, which is the heading angle of the drone.

Assuming a model where the UAV’s velocity is set at 1 and the maximum curvature is 1/γ, the kinematic equations can be expressed as follows:

$$\dot{\mathrm{x}}=\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }$$

(2)

$$\dot{\mathrm{y}}=\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }$$

(3)

$$\dot{\mathsf{\theta }}=\mathrm{u}/\mathsf{\gamma }$$

(4)

In this context, $\mathsf{\theta }\in \left[-1/\mathsf{\gamma }1/\mathsf{\gamma }\right]$, and the constant u is the control variable that varies with time for the drone, with an absolute value of $\left|\mathrm{u}\right|\le 1$.

In mathematical modeling, the simulation of an unmanned aerial vehicle’s curved motion with curvature 1/γ is typically achieved through the utilization of Dubins curves [20]. J-D. Boissonnat et al. had a detailed account of the specific derivations for various scenarios involving Dubins curves [21]. Dubins curves are commonly divided into two major families: the CCC family, which includes RLR and LRL, and the CSC family, which includes RSR, LSL, RSL, and LSR [22]. C represents an arc with curvature not exceeding 1/γ, S represents a straight line, L represents a counterclockwise curve, and R represents a clockwise curve. The Dubins curve is depicted by drawing the maximum curvature radius circle tangent to the initial state, followed by the delineation of two tangent lines connecting the circles. Through defined coordinate positions and heading angles, the shortest path can be determined [23].

The Dubins shortest-path problem belongs to the realm of solving minimization problems in traditional optimal control theory. In this problem, the initial and final states of the UAV are uncertain, and the control vector is a continuous function, making it amenable to solution using variational methods. In the UAV motion model presented in this paper, the UAV’s velocity is set to 1. The performance functional for the shortest path of Dubins curves can be represented by an integral, as shown in Equation (5).

$$\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{J}={\int }_{{\mathrm{t}}_0}^{{\mathrm{t}}_{\mathrm{f}}}\mathrm{L}(\mathrm{x},\mathrm{y},\mathsf{\theta })\mathrm{d}\mathrm{t}={\int }_{{\mathrm{t}}_0}^{{\mathrm{t}}_{\mathrm{f}}}1\mathrm{d}\mathrm{t}$$

(5)

In Equation (5), ${\mathrm{t}}_0$ and ${\mathrm{t}}_{\mathrm{f}}$ represent the initial and terminal times of the UAV. We set ${\mathsf{\phi }}_{\mathrm{x}}$, ${\mathsf{\phi }}_{\mathrm{y}}$, ${\mathsf{\phi }}_{\mathsf{\theta }}$ as parameters for the shortest path. Combining Equations (2)–(4), which represent the UAV’s kinematic equations, with Equation (5), which represents the performance functional, we construct the Hamiltonian function based on the theorems of optimal control theory. Therefore, the Hamiltonian function for the UAV’s shortest path is given as

$$\mathrm{H}\left(\mathrm{x},\mathrm{y},\mathsf{\theta }\right)=1+{\mathsf{\phi }}_{\mathrm{x}}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }+{\mathsf{\phi }}_{\mathrm{y}}\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }+{\mathsf{\phi }}_{\mathsf{\theta }}\mathrm{u}/\mathsf{\gamma }$$

(6)

According to the Pontryagin maximum principle, we can obtain the necessary conditions for the shortest path of the Dubins curve [24].

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}\left[\begin{array}{c}{\mathsf{\phi }}_{\mathrm{x}}\\ {\mathsf{\phi }}_{\mathrm{y}}\\ {\mathsf{\phi }}_{\mathsf{\theta }}\end{array}\right]=\left(\begin{array}{c}0\\ 0\\ {\mathsf{\phi }}_{\mathrm{x}}\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }-{\mathsf{\phi }}_{\mathrm{y}}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }\end{array}\right)$$

(7)

The above equation indicates that ${\mathsf{\phi }}_{\mathrm{x}}$ and ${\mathsf{\phi }}_{\mathrm{y}}$ are piecewise constant within the UAV’s motion cycle [${\mathrm{t}}_0$, ${\mathrm{t}}_{\mathrm{f}}$]. Due to the linearity of the Hamiltonian in the control function u, there are two distinct types of optimal control [25]. In the situation that ${\mathsf{\phi }}_{\mathsf{\theta }}\ne 0$, u takes on the value of $\mathrm{u}=-\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({\mathsf{\phi }}_{\mathsf{\theta }}\right)$. Conversely, when ${\mathsf{\phi }}_{\mathsf{\theta }}=0$, u ∈ [−1, 1]. The former condition corresponds to the maximum curvature arc with a radius of r, while the latter aligns with the trajectory of a straight line where the direction is given by $\mathrm{t}\mathrm{a}\mathrm{n}\mathsf{\theta }={\mathsf{\phi }}_{\mathrm{x}}/{\mathsf{\phi }}_{\mathrm{y}}$.

Sussmann and his colleagues, in their seminal paper in 1991, established the existence of a solution to the Dubins path problem, wherein the optimal trajectories manifest as either C_αC_βC_δ or C_αS_dC_β, where 0 ≤ α, δ < 2π, π < β < 2π, and d ≥ 0 [26]. Given the inherent non-uniqueness of solutions to the Dubins path, generally, multiple trajectories can converge to the minimal value [27]. As discerned from Equation (7), the shortest path of the Dubins curve is intricately linked to the heading angle. The choice of distinct heading angles results in divergent selections within the path details. Consequently, determining the heading angles in the Dubins curve emerges as a pivotal factor in ascertaining the unmanned aerial vehicle’s shortest path [28]. In Section 2.2, we introduce a heading angle model that exhibits a notable enhancement when compared to existing algorithms. The specific data comparisons are presented in Section 4.

2.2. Compensatory Look-Ahead Algorithm

In Section 2.1, the necessary conditions for the shortest path and the existence of the optimal solution are demonstrated. It is also shown that the shortest path of Dubins curves is closely related to the heading angle. Determining the heading angle in Dubins curves is crucial for determining the shortest path for unmanned aerial vehicles. Compared to path planning between two waypoints, smooth curve optimization between three waypoints better meets practical requirements and yields superior planning results. Ma et al. introduced a midpoint $\left({\mathrm{x}}_1,{\mathrm{y}}_1\right)$ between the initial point $\left({\mathrm{x}}_0,{\mathrm{y}}_0\right)$ and the terminal point $\left({\mathrm{x}}_2,{\mathrm{y}}_2\right)$ of the Dubins curve. Referring to the analysis in [26], they obtained the boundary relationships of the midpoint, as shown below.

$$\begin{array}{c}{\mathsf{\phi }}_{\mathrm{x}}\left({\mathrm{t}}_1^{+}\right)={\mathsf{\phi }}_{\mathrm{x}}\left({\mathrm{t}}_1^{-}\right)+{\mathrm{k}}_1\\ {\mathsf{\phi }}_{\mathrm{y}}\left({\mathrm{t}}_1^{+}\right)={\mathsf{\phi }}_{\mathrm{y}}\left({\mathrm{t}}_1^{-}\right)+{\mathrm{k}}_2\\ \begin{array}{c}{\mathsf{\phi }}_{\mathsf{\theta }}\left({\mathrm{t}}_1^{+}\right)={\mathsf{\phi }}_{\mathsf{\theta }}\left({\mathrm{t}}_1^{-}\right)\\ \mathrm{H}\left({\mathrm{t}}_1^{+}\right)=\mathrm{H}\left({\mathrm{t}}_1^{-}\right)\end{array}\end{array}$$

(8)

In Equation (8), ${\mathrm{k}}_1$ and ${\mathrm{k}}_2$ are constants. Combined with Equation (7), it becomes evident that ${\mathsf{\phi }}_{\mathrm{x}}$ and ${\mathsf{\phi }}_{\mathrm{y}}$ exhibit a discontinuous jump at $\left({\mathrm{x}}_1,{\mathrm{y}}_1\right)$ and are not continuous.

Figure 2 presents a schematic diagram of the optimal path. The red trajectory represents the Dubins optimal path, which passes through the initial point, midpoint, and terminal point. ${\mathsf{\theta }}_0$ is the heading angle (relative to the X-axis) of the initial point $\left({\mathrm{x}}_0,{\mathrm{y}}_0\right)$. ${\mathsf{\theta }}_1$ and ${\mathsf{\theta }}_2$ are the angles of the line where they are located (relative to the X-axis). ${\mathsf{\theta }}_{\mathrm{m}}$ represents the heading angle of the midpoint $\left({\mathrm{x}}_1,{\mathrm{y}}_1\right)$. Ma et al. reached the conclusion that a line drawn from the midpoint $\left({\mathrm{x}}_1,{\mathrm{y}}_1\right)$ to the point of intersection A $\left({\mathrm{x}}_3,{\mathrm{y}}_3\right)$ of the two line segments bisects the angle of intersection, which is shown by the blue line segment in Figure 2. On the optimal path, ${\mathsf{\phi }}_{\mathsf{\theta }}=0$ and ${\mathrm{d}\mathsf{\phi }}_{\mathsf{\theta }}/\mathrm{d}\mathrm{t}=0$. Substituting these values into Equation (8) allows the system of simultaneous equations to be solved. The results can then be utilized in Equation (9).

$$\begin{array}{c}{\mathsf{\phi }}_{\mathrm{x}}=-\cos {\mathsf{\theta }}_1\\ {\mathsf{\phi }}_{\mathrm{y}}=-\sin {\mathsf{\theta }}_1\\ \begin{array}{c}{\mathsf{\phi }}_{\mathrm{x}}+{\mathrm{k}}_1=-\cos {\mathsf{\theta }}_2\\ {\mathsf{\phi }}_{\mathrm{y}}+{\mathrm{k}}_2=-\sin {\mathsf{\theta }}_2\end{array}\end{array}$$

(9)

From Equation (9), we can derive that ${\mathrm{k}}_1/{\mathrm{k}}_2=-\tan {{(\mathsf{\theta }}_1+\mathsf{\theta }}_2)/2$. The Hamiltonian function is continuous along the optimal trajectory at point (${\mathrm{x}}_1,{\mathrm{y}}_1$). Therefore, it can be represented as follows:

$$\begin{gathered} \underset{‖\mathrm{u}‖\le 1}{\min }\left[1+{\mathsf{\phi }}_{\mathrm{x}}\cos {\mathsf{\theta }}_{\mathrm{m}}+{\mathsf{\phi }}_{\mathrm{y}}\sin {\mathsf{\theta }}_{\mathrm{m}}+{\mathsf{\phi }}_{\mathsf{\theta }}\mathrm{u}\right]= \\ \underset{‖\mathrm{u}‖\le 1}{\min }\left[1+({\mathsf{\phi }}_{\mathrm{x}}+{\mathrm{k}}_1)\cos {\mathsf{\theta }}_{\mathrm{m}}+{\mathsf{\phi }}_{\mathrm{y}}+{\mathrm{k}}_2\sin {\mathsf{\theta }}_{\mathrm{m}}+{\mathsf{\phi }}_{\mathsf{\theta }}\mathrm{u}\right] \end{gathered}$$

(10)

From the above equation, it can be deduced that

$${\mathrm{k}}_1\cos {\mathsf{\theta }}_{\mathrm{m}}+{\mathrm{k}}_2\sin {\mathsf{\theta }}_{\mathrm{m}}=0$$

(11)

From Equations (7) and (9), Equation (10) can be obtained.

$${\mathsf{\theta }}_{\mathrm{m}}=\frac{{\mathsf{\theta }}_1+{\mathsf{\theta }}_2}{2}$$

(12)

By Equation (10), the midpoint heading angle of the shortest path can be obtained, and the solution to the coupled trigonometric equations can be found. While this can be accomplished numerically, it is challenging to describe using a general formula. In their paper, Ma et al. did not provide a complete equation that could represent all possible geometric configurations of three path points and initial direction, thus uniquely determining the optimal trajectory for any problem. Solving the coupled trigonometric equations can be performed numerically but is difficult to describe with a general formula. Therefore, based on the horizon backtracking principle, we propose a four-point optimization-based algorithm for calculating the heading angle (Algorithm 1) and provide a general formula for computing the heading angle size at the target point.

First, let us consider the scenario where the distances between the initial, intermediate, and terminal waypoints are relatively long. Figure 3 illustrates the path planning between the initial, intermediate, and terminal waypoints, where the solid red lines represent Dubins trajectories. Point F (${\mathrm{x}}_4,{\mathrm{y}}_4$) and point D (${\mathrm{x}}_5,{\mathrm{y}}_5$) denote the tangent points of two circular trajectories. While the midpoint heading angle ${\mathsf{\theta }}_{\mathrm{m}}$ can be solved using Equation (10), when the distances between the initial and intermediate points are considerable, the error between the angle ${\mathsf{\theta }}_1$ (relative to the X-axis) of the line segment connecting point F (${\mathrm{x}}_4,{\mathrm{y}}_4$) and point D (${\mathrm{x}}_5,{\mathrm{y}}_5$) and the angle ${\mathsf{\theta }}_3$ (relative to the X-axis) of the line segment connecting the initial point B (${\mathrm{x}}_0,{\mathrm{y}}_0$) and the tangent point D (${\mathrm{x}}_5,{\mathrm{y}}_5$) (red dashed line segment) can be neglected. Solving for the coordinates of tangent point F (${\mathrm{x}}_4,{\mathrm{y}}_4$) is difficult, so the direction angle ${\mathsf{\theta }}_{\mathrm{m}}$ (relative to the X-axis) of the intermediate point can be directly solved using the coordinates of the initial point B (${\mathrm{x}}_0,{\mathrm{y}}_0$) substituted into Equation (12).

From Figure 3, the total range L between the three points after task planning can be obtained, as shown below:

$$\begin{array}{c}\mathrm{L}={\mathrm{L}}_{\mathrm{B}\mathrm{D}}+{\mathrm{L}}_{\mathrm{E}\mathrm{C}}+{\mathrm{S}}_{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{D}\mathrm{E}}\\ >\mathrm{a}+\mathrm{b}+\frac{\mathsf{\alpha }+\mathsf{\beta }}{180}\mathsf{\pi }\mathrm{R}-2\mathrm{R}\sin \frac{\mathsf{\alpha }}{2}-2\mathrm{R}\sin \frac{\mathsf{\beta }}{2}\\ >\mathrm{a}+\mathrm{b}+\frac{\mathsf{\alpha }+\mathsf{\beta }}{180}\mathsf{\pi }\mathrm{R}-4\mathrm{R}\sin \frac{\mathsf{\alpha }+\mathsf{\beta }}{4}\cos \frac{\mathsf{\alpha }-\mathsf{\beta }}{4}\end{array}$$

(13)

From Equation (13), it can be seen that the value of L is related to α, β, as shown in Equation (14):

$$\mathrm{S}=\frac{\mathsf{\alpha }+\mathsf{\beta }}{180}\mathsf{\pi }\mathrm{R}-4\mathrm{R}\sin \frac{\mathsf{\alpha }+\mathsf{\beta }}{4}\cos \frac{\mathsf{\alpha }-\mathsf{\beta }}{4}$$

(14)

It is discernible from the equation that as the sum of α and β decreases, reaching equality (α = β), the value of Equation (14) attains its minimum. Illustrated graphically, when the center O of the circle lies on the angle bisector of angle A, the total course length L achieves its minimum value, denoted as ${\mathrm{L}}_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{a}+\mathrm{b}+\frac{\mathsf{\alpha }+\mathsf{\beta }}{180}\mathsf{\pi }\mathrm{R}-4\mathrm{R}\sin \frac{\mathsf{\alpha }+\mathsf{\beta }}{4}$. The heading angle ${\mathsf{\theta }}_{\mathrm{m}}$ at the midpoint $\left({\mathrm{x}}_1,{\mathrm{y}}_1\right)$ is given by ${\mathsf{\theta }}_{\mathrm{m}}=\frac{{\mathsf{\theta }}_3+{\mathsf{\theta }}_2}{2}$.

Next, let us consider the scenario where the distances between the initial, intermediate, and terminal waypoints are relatively close. Figure 3 depicts the path planning between the initial, intermediate, and terminal waypoints, where the solid red lines represent Dubins trajectories. However, when the distances between the initial and intermediate points are close, the error between the angle ${\mathsf{\theta }}_1$ (relative to the X-axis) of the line segment connecting point F (${\mathrm{x}}_4,{\mathrm{y}}_4$) and point D (${\mathrm{x}}_5,{\mathrm{y}}_5$) and the angle ${\mathsf{\theta }}_3$ (relative to the X-axis) of the line segment connecting the initial point B (${\mathrm{x}}_0,{\mathrm{y}}_0$) and the tangent point D (${\mathrm{x}}_5,{\mathrm{y}}_5$) (red dashed line segment) cannot be neglected. Therefore, we cannot directly substitute the coordinates of the initial point B (${\mathrm{x}}_0,{\mathrm{y}}_0$) into Equation (12) to solve for the direction angle ${\mathsf{\theta }}_{\mathrm{m}}$ (relative to the X-axis) of the intermediate point. In this scenario, the heading angles of the initial and terminal points will affect the values of ${\mathsf{\theta }}_1$ and ${\mathsf{\theta }}_2$.

From Figure 3, it can be observed that when the waypoints are close to each other, the actual ${\mathsf{\theta }}_3$ can significantly differ from the original ${\mathsf{\theta }}_1$. In such cases, continuing to use ${\mathsf{\theta }}_1$ for calculations would lead to significant errors. It is challenging to obtain a solution for the fully coupled trigonometric equations. To address this issue, we propose a four-point optimization algorithm for calculating the heading angles, aiming to minimize the deviation caused by error terms.

In the scenario where the initial point (${\mathrm{x}}_0,{\mathrm{y}}_0$) is close to the midpoint, the heading angle at the initial point will affect the heading angle at the midpoint. In this paper, we introduce an adjacent point (${\mathrm{x}}_6,{\mathrm{y}}_6$) to the initial point (${\mathrm{x}}_0,{\mathrm{y}}_0$). We preprocess the target point (${\mathrm{x}}_6,{\mathrm{y}}_6$) using the optimization algorithm for the scenario where the initial point is far from the midpoint, as described in Figure 3, to obtain a tangent point Q₁ (${\mathrm{x}}_4,{\mathrm{y}}_4$). Then, we use this tangent point Q₁ (${\mathrm{x}}_4,{\mathrm{y}}_4$) instead of the initial point (${\mathrm{x}}_0,{\mathrm{y}}_0$) as the new target point. We repeat the previous step for the midpoint (${\mathrm{x}}_1,{\mathrm{y}}_1$). This process allows us to obtain a compensated heading angle ${\mathsf{\theta }}_{\mathrm{m}}$ for the midpoint. As shown in Figure 4, ${\mathsf{\theta }}_1$ represents the angle of the line segment connecting tangent points Q₁ (${\mathrm{x}}_4,{\mathrm{y}}_4$) and Q₂ (${\mathrm{x}}_5,{\mathrm{y}}_5$), and ${\mathsf{\theta }}_2$ represents the angle of the line segment connecting tangent point Q₃ (${\mathrm{x}}_7,{\mathrm{y}}_7$) and the terminal point C. O₁ and O₂ are the centers of the two circular trajectories. The solid black line represents the process of solving for the substitute point Q₁, and the dashed black line represents the process of using the substitute point Q₁ instead of the initial point B to solve for the direction angle ${\mathsf{\theta }}_{\mathrm{m}}$ of the midpoint.

$${\mathsf{\theta }}_{\mathrm{m}}=\frac{{\mathsf{\theta }}_1+{\mathsf{\theta }}_2}{2}$$

(15)

Algorithm 1. Compensation look-ahead algorithm

Input: B$\left({\mathbf{x}}_0,{\mathbf{y}}_0\right),\mathbf{C}\left({\mathbf{x}}_1,{\mathbf{y}}_1\right),\mathbf{D}\left({\mathbf{x}}_6,{\mathbf{y}}_6\right),\mathbf{E}\left({\mathbf{x}}_2,{\mathbf{y}}_2\right),\mathbf{R}$ Output: compensation heading angle ${\mathsf{\theta }}_{\mathbf{m}}$ 1:   Center of circle ${\mathbf{O}}_1=\overrightarrow{\mathbf{B}}+\left(\right(\overrightarrow{\mathbf{B}\mathbf{C}}+\overrightarrow{\mathbf{B}\mathbf{D}})/\sqrt{{\left|\mathbf{B}\mathbf{C}\right|}^2+{\left|\mathbf{B}\mathbf{D}\right|}^2})/2$ 2:   Distance from point to tangent point $\left|\mathbf{C}{\mathbf{Q}}_1\right|=\sqrt{{\left|\mathbf{C}{\mathbf{O}}_1\right|}^2-{\mathbf{R}}^2}$ 3:   Acquiring tangent point ${\mathbf{Q}}_1=\left(\frac{\overrightarrow{\mathbf{C}{\mathbf{O}}_1}}{\left|\mathbf{C}{\mathbf{O}}_1\right|}\right)\ast \mathbf{C}\mathbf{o}\mathbf{s}\left(\mathbf{S}\mathbf{i}\mathbf{n}\right)\left(\mathbf{a}\mathbf{s}\mathbf{i}\mathbf{n}\left(\frac{\mathbf{R}}{\left|\mathbf{C}{\mathbf{Q}}_1\right|}\right)\right)$ 4:   Center of circle ${\mathbf{O}}_2=\overrightarrow{\mathbf{c}}+\left(\right(\overrightarrow{\mathbf{C}{\mathbf{Q}}_1}+\overrightarrow{\mathbf{C}\mathbf{E}})/\sqrt{{\left|\mathbf{C}{\mathbf{Q}}_1\right|}^2+{\left|\mathbf{C}\mathbf{E}\right|}^2})/2$ 5:   Distance from point to tangent point $\left|\mathbf{B}{\mathbf{Q}}_2\right|=\sqrt{{\left|\mathbf{B}{\mathbf{O}}_2\right|}^2-{\mathbf{R}}^2}$ 6:   Distance from point to tangent point $\left|\mathbf{B}{\mathbf{Q}}_3\right|=\sqrt{{\left|\mathbf{B}{\mathbf{O}}_2\right|}^2-{\mathbf{R}}^2}$ 7:   Acquiring tangent point ${\mathbf{Q}}_2=\left(\frac{\overrightarrow{\mathbf{B}{\mathbf{O}}_2}}{\left|\mathbf{B}{\mathbf{O}}_2\right|}\right)\ast \mathbf{C}\mathbf{o}\mathbf{s}\left(\mathbf{S}\mathbf{i}\mathbf{n}\right)\left(\mathbf{a}\mathbf{s}\mathbf{i}\mathbf{n}\left(\frac{\mathbf{R}}{\left|\mathbf{B}{\mathbf{Q}}_2\right|}\right)\right)$ 8:   Acquiring tangent point ${\mathbf{Q}}_3=\left(\frac{\overrightarrow{\mathbf{E}{\mathbf{O}}_2}}{\left|\mathbf{E}{\mathbf{O}}_2\right|}\right)\ast \mathbf{C}\mathbf{o}\mathbf{s}\left(\mathbf{S}\mathbf{i}\mathbf{n}\right)\left(\mathbf{a}\mathbf{s}\mathbf{i}\mathbf{n}\left(\frac{\mathbf{R}}{\left|\mathbf{E}{\mathbf{Q}}_3\right|}\right)\right)$ 9:   Obtain compensation heading angle ${\mathsf{\theta }}_{\mathbf{m}}=\left|\mathbf{a}\mathbf{c}\mathbf{o}\mathbf{s}\left(\frac{({\mathbf{x}}_2-{\mathbf{x}}_7)}{\left|{\mathbf{Q}}_3\mathrm{E}\right|}\right)\right|+\left|\mathbf{a}\mathbf{c}\mathbf{o}\mathbf{s}\left(\frac{({\mathbf{x}}_5-{\mathbf{x}}_6)}{\left|{\mathbf{Q}}_1{\mathbf{Q}}_2\right|}\right)\right|\mp 90°$ 10:   end

3. Model Decomposition

3.1. MTSP Model Based on Shortest Distance

The multiple traveling salesman problem (MTSP) is an extension of the traveling salesman problem (TSP) [29,30]. In the MTSP problem, there are n couriers assigned to visit m destinations, where m exceeds n. Each destination is assigned to only one courier, with the objective of minimizing the combined travel distance of all couriers [31]. Task planning of a drone cluster can be seen as an instantiation of the MTSP. It can be succinctly described as follows: provide an undirected graph G = (V, A), where V is the set of nodes, and A is the set of arcs representing road segments. The arcs of the undirected graph G have non-negative costs, and the objective function is to divide V into g non-empty subsets ${\left\{{\mathrm{S}\mathrm{u}\mathrm{b}}_{\mathrm{i}}\right\}}_{\mathrm{g}}^1$ and find the minimum cost through each target point in each subset ${\mathrm{S}\mathrm{u}\mathrm{b}}_{\mathrm{i}}$ [32].

The total range from the starting point of the drone cluster to the return of all drones is used as the optimization objective of the MTSP model. There are n task target points, m drones, and the coordinates of all target points are represented by $\mathrm{T}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}=\left[({\mathrm{x}}_1,{\mathrm{y}}_1),({\mathrm{x}}_2,{\mathrm{y}}_2),\cdots ,({\mathrm{x}}_{\mathrm{n}},{\mathrm{y}}_{\mathrm{n}})\right]$. Therefore, the distance between the two points after the drone cluster task planning can be expressed as

$${\mathrm{d}}_{\mathrm{i},\mathrm{j}}=\left|{\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}\right|+\left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{j}}\right|$$

(16)

The drones in this paper are all fixed-wing drones, which are different from rotary-wing drones. The trajectory change in fixed-wing drones belongs to a linear process. Therefore, the distance equation between the two points needs to add constraints that meet the trajectory of the drone. Here, we simulate the actual trajectory-change ability of the drone using a Dubins curve with a turning radius of 100 m:

$${\mathrm{d}}_{\mathrm{i},\mathrm{j}}=\mathrm{l}+\mathrm{r}+\mathrm{D}$$

(17)

In Equation (17), D represents the straight portion of the Dubins curve between adjacent task points in the task planning trajectory, and l and r are the arc lengths between adjacent task points and the straight line, respectively.

The optimization objective of the UAV swarm Dubins multiple traveling salesman model is to minimize the total flight distance of the UAV swarm. According to Equation (17), the objective function can be represented as follows:

$$\min \mathrm{L}={\sum }_{\mathrm{i}=1,\mathrm{j}=1,\mathrm{q}=1}^{\mathrm{n},\mathrm{n},\mathrm{m}}{\mathrm{d}}_{\mathrm{i},\mathrm{j}}\ast {\mathrm{k}}_{\mathrm{i},\mathrm{j},\mathrm{q}}$$

(18)

In Equation (18), ${\mathrm{k}}_{\mathrm{i},\mathrm{j},\mathrm{q}}$ represents whether UAV q has a planned task between target points i and j. The value of ${\mathrm{k}}_{\mathrm{i},\mathrm{j},\mathrm{q}}$ is either 0 or 1. The constraint functions are as follows:

$${\sum }_{\mathrm{i}=1,\mathrm{j}=1,\forall \mathrm{q}\in \mathrm{m}}^{\mathrm{n},\mathrm{n}}{\mathrm{k}}_{\mathrm{i},\mathrm{j},\mathrm{q}}\ge \mathrm{k}$$

(19)

$${\sum }_{\mathrm{i}=1,\mathrm{j}=1,\mathrm{q}=1}^{\mathrm{n},\mathrm{n},\mathrm{m}}{\mathrm{k}}_{\mathrm{i},\mathrm{j},\mathrm{q}}=\mathrm{n}$$

(20)

$${\sum }_{\mathrm{i}\in \mathrm{n},\mathrm{j}\in \mathrm{n},\mathrm{q}=1}^{\mathrm{m}}{\mathrm{k}}_{\mathrm{i},\mathrm{j},\mathrm{q}}=1$$

(21)

$${\mathrm{d}}_{1,\mathrm{w}}=\left|\left({\mathrm{x}}_1,{\mathrm{y}}_1\right),\left({\mathrm{x}}_{\mathrm{w}},{\mathrm{y}}_{\mathrm{w}}\right)\right|=0,\forall (\mathrm{x},\mathrm{y})\in \mathrm{T}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}$$

(22)

Equation (19) is the number of targets assigned to each drone, which ensures that each drone is assigned no less than k targets; Equation (20) represents the total number of targets allocated to the drone cluster, which ensures that all target points are allocated; Equation (21) is the number of times the target point is assigned to the drone, and this constraint ensures that each point is only assigned once; and Equation (22) represents the decision criterion for the starting and ending points. $\left({\mathrm{x}}_{\mathrm{w}},{\mathrm{y}}_{\mathrm{w}}\right)$ denotes the coordinates of the last target point planned for a single UAV, indicating that each aircraft must start and end at the same point.

Dijkstra algorithm is a classic method for solving the MTSP (multiple traveling salesman problem), but it is only suitable for scenarios with a small number of samples. The computational complexity increases exponentially with the number of samples, making it unsuitable for large-scale UAV formations in the current context. The discrete DTSP method (DDTSP) is another classical approach for handling the MTSP problem. In DDTSP, for each target, h possible angles are considered, and the interval [0, 2π) is uniformly divided into h sub-intervals. The resulting discrete DTSP (DDTSP) can be transformed into an integer linear programming (ILP) problem. DDTSP method is equivalent to the generalized asymmetric traveling salesman problem with n clusters and 1 + h (n − 1) cities (Noon Bean, 1991), which requires a small value of h. However, the minimum time function of the Dubins model is non-continuous. Therefore, to approximate the optimal solution of DTSP, a larger value of h is required [33]. When the number of task points is small, DDTSP can provide the optimal path, but there are problems of high computational load and slow computational efficiency in large-scale task points and multiple traveling merchants [34,35]. In order to balance the number of task points and computational efficiency, the fourth section proposes a DTSP processing method based on ETSP classification, which simplifies the computational workload while ensuring the shortest path. The specific data comparisons are shown in the fourth section.

3.2. UAV Formation EDTSP Model

In the scenario of multiple unmanned aerial vehicles (UAVs) tasked with multiple objectives, direct application of the DDTSP algorithm for solving the shortest path is deemed imprudent. Hence, a preprocessing step for the sequential arrangement of mission objectives becomes imperative. This paper introduces a preprocessing DTSP model based on the ETSP algorithm called the EDTSP model. The specific approach of this model is outlined as follows.

Let the set of target points be denoted as $\mathrm{T}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}=\left[({\mathrm{x}}_1,{\mathrm{y}}_1),({\mathrm{x}}_2,{\mathrm{y}}_2),\cdots ,({\mathrm{x}}_{\mathrm{n}},{\mathrm{y}}_{\mathrm{n}})\right]\in {\mathrm{R}}^2$, and the ensemble of UAV formations as U. We conceptualize the DTSP problem as an ETSP problem for solving and computing the mission target points, obtaining subsets ${\mathrm{S}\mathrm{u}\mathrm{b}}_{\mathrm{i}}=\left\{({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}})\right\}\in \mathrm{T}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}$. These obtained subsets of targets serve as the target point set to be solved. We use the compensatory look-ahead algorithm proposed in Section 2 of this paper to determine the heading angle size for each target point, completing the path planning. The planning results of several subsets constitute the overall path planning. Pre-grou** the original task target points with the ETSP algorithm reduces the dimensionality in the DTSP solving process, thereby reducing computational complexity.

Subsequently, within each subset ${\mathrm{S}\mathrm{u}\mathrm{b}}_{\mathrm{i}}$, we use the compensatory look-ahead algorithm to solve the DTSP problem for the target sequence, thereby deriving UAV formation planning paths that meet the specified requirements. The specific process of the model is shown in Figure 5. To illustrate our model method, let us illustrate our model approach with the following example: Suppose we have a task point set consisting of 3n target points, and we have m UAVs for surveillance, with each UAV assigned to n targets. We use the ETSP algorithm to plan paths for the target set, resulting in three subsets of targets. After allocating the subsets, under the constraints of Dubins curves and UAV performance, where $\mathsf{\theta }$ belongs to $\left[-1/\mathsf{\gamma },1/\mathsf{\gamma }\right]$, we solve the TSP problem to obtain the optimal trajectories. We compute the results for the three subsets and combine them to form the overall planning results. If we were to directly solve the DTSP model using bio-inspired algorithms, we would need to consider $\left(3\mathrm{n}\right)!$ cases. However, by preprocessing the target point set with the ETSP model, we only need to consider $3\ast \mathrm{n}!$ cases, greatly reducing the computational complexity.

4. The Analysis of the Simulations

In the third and fourth sections, we introduced a compensatory look-ahead algorithm based on four-point optimization and a path-planning model for multi-UAV formations in multi-tasking scenarios. To prove the accuracy of our algorithm and model, we randomly distributed a certain number of target points in both dense (1000 m × 1000 m) and sparse (5000 m × 5000 m) regions. Using various algorithms, we generated the shortest planned paths.

Simulation 1: Using a single UAV, taking off from any location, we executed trajectory planning for ten task points in the dense region. Upon completing the tasks, the UAV returned to the departure point. We compared the itineraries of planning paths utilizing the compensation look-ahead algorithm, 2-opt look-ahead algorithm, alternating algorithm, and the ETSP model. Comparing the search efficiency of the four models.

Figure 6 describes the planned trajectories utilizing the compensatory look-ahead algorithm, 2-opt look-ahead algorithm, alternating algorithm, and the ETSP model. The straight-line trajectory based on the ETSP algorithm serves as the reference standard under ideal conditions. It can be seen that the model trajectory with the compensatory look-ahead algorithm is notably smoother. Through the comparison of total mileage in

Table 1. Total mileage of path planning in dense areas using different algorithms.

	Compensation Look-Ahead Algorithm	2-Opt Look-Ahead Algorithm	Alternating Algorithm	ETSP
Total mileage (102 m)	36.7695	42.2248	48.0164	27.6005

, it can be affirmed that in dense regions, the compensatory look-ahead algorithm exhibits a minimum improvement of 12.9% compared to other algorithms. The increase in Euclidean distance is minimal, with a mere 33.2% increment. This is for the constraints imposed by the maximum curvature of the UAV in dense regions.

Table 2. Total mileage of path planning with different target points for different algorithms in dense regions.

Target Number (Total Mileage 102 m)	Compensation Look-Ahead Algorithm	2-Opt Look-Ahead Algorithm	Alternating Algorithm	ETSP
5	23.9851	24.8415	26.8562	22.3707
10	36.7695	42.2248	48.0164	27.6005
15	47.1584	55.5124	60.8612	35.2885
20	56.5147	64.3615	76.1574	37.7948
25	65.1845	74.6124	87.1451	41.1114

displays the path-planning total mileage for different algorithms with 5, 10, 15, 20, and 25 target points. Through this table, it is discernible that the compensatory look-ahead algorithm is consistently better than other algorithmic models. This further substantiates the efficacy of the compensatory look-ahead algorithm in optimizing planning paths within dense regions, aligning seamlessly with our requirements.

Simulation 2: Using a single unmanned aerial vehicle (UAV) from a designated location, we conducted trajectory planning for 10 target points within a sparse region. Upon completing the tasks, the UAV returned to the starting point. We compared the travel distances of the planned paths using the compensatory look-ahead algorithm, 2-opt look-ahead algorithm, alternating algorithm, and ETSP model. Comparing the search efficiency of the four models.

Figure 7 describes the planned paths using the compensatory look-ahead algorithm, 2-opt look-ahead algorithm, alternating algorithm, and ETSP model within the sparse region. The straight-line trajectory based on the ETSP algorithm serves as the reference standard under ideal conditions. In comparison to simulation 1, the disparities among various algorithms in simulation 2 are not substantial. The reason is the vast distances between adjacent points in the sparse region, where subtle variations in heading angles have less impact compared to dense areas. The comparison of total mileage across different algorithms in

Table 3. Total mileage of path planning in sparse areas using different algorithms.

	Compensation Look-Ahead Algorithm	2-Opt Look-Ahead Algorithm	Alternating Algorithm	ETSP
Total mileage (102 m)	36.7695	42.2248	48.0164	27.6005

reveals that the compensatory look-ahead algorithm still yields the shortest planned path. In contrast to other algorithms, it boasts an improvement of up to 3.1%, with the increase in the shortest Euclidean distance merely reaching 1.4%.

Simulation 3 involves three unmanned aerial vehicles (UAVs) taking off from a designated location and planning trajectories for 30 target points within a sparse region. After completing the tasks, the UAVs return to the starting point, a comparison of the travel distances of planned paths using the compensatory look-ahead algorithm, 2-opt look-ahead algorithm, alternating algorithm, and ETSP model. Comparing the search efficiency of the four models.

Figure 8 describes the trajectory planning for multiple unmanned aerial vehicles (UAVs) and various target points employing distinct algorithms for takeoff from the same points. The straight-line trajectory based on the ETSP algorithm serves as the reference standard under ideal conditions. From the illustration, it is evident that the compensatory look-ahead algorithm adeptly accomplishes the planning of paths for multiple UAVs and multiple targets. Notably, the total mileage of path planning for multiple UAVs and multiple targets using different algorithms is shown in

Table 4. The total mileage of path planning for multiple UAVs and multiple targets using different algorithms starting from the same point.

	Compensation Look-Ahead Algorithm	2-Opt Look-Ahead Algorithm	Alternating Algorithm	ETSP
Total mileage (102 m)	141.17	154.115	159.702	109.995

. The compensatory look-ahead algorithm demonstrates an improvement of no less than 8.4% in the context of the multi-tasking, multi-UAV path-planning model.

Simulation 4 involves three unmanned aerial vehicles (UAVs) taking off from different locations and planning trajectories for 30 target points within a sparse region. After completing the tasks, the UAVs return to the starting point, a comparison of the travel distances of planned paths using the compensatory look-ahead algorithm, 2-opt look-ahead algorithm, alternating algorithm, and ETSP model. Comparing the search efficiency of the four models.

Figure 9 shows the planned trajectory of multiple unmanned aerial vehicles (UAVs) taking off from different points using different algorithms. The straight-line trajectory based on the ETSP algorithm serves as the reference standard under ideal conditions. The three drones are takeoff at [0, 7.5], [12.5, 0], and [25, 7.5]. From the graph, it can be seen that the compensation look-ahead algorithm can complete the multi-objective path-planning task for multiple drones. Compared with the 2-opt look-ahead algorithm and the alternating algorithm, the planned trajectory is more in line with practical requirements, and the total distance of the planned path is the shortest. According to

Table 5. Total mileage of path planning for multiple UAVs and multiple targets using different algorithms starting from different points.

	Compensation Look-Ahead Algorithm	2-Opt Look-Ahead Algorithm	Alternating Algorithm	ETSP
Total mileage (102 m)	114.1742	124.1158	138.4113	94.7549

, it can be calculated that the compensation forward algorithm has a maximum improvement of 17.5% in the multi-task and multi-drone path-planning model compared to the alternating algorithm.

5. Conclusions

This article mainly studies the shortest-path planning problem of unmanned aerial vehicle formation under uncertain target sequences. Firstly, in order to improve the efficiency of collaborative search in drone clusters, a four-point optimized heading angle calculation method is proposed based on the horizon algorithm. On the basis of Ma et al.’s research, a general formula for calculating heading angle was proposed. Through the simulation data in the fourth section, it can be proven that our proposed compensation look-ahead algorithm reduces the situation where drones need to perform large-scale curve movements to reach the target point coordinates. Compared with existing algorithms, it has a maximum improvement of 12.9%. Secondly, for the large-scale formation problem, the ETSP algorithm is used to pre-allocate the large-scale target points to be allocated and a single drone, forming multiple subsets of target points. Then, the compensation foresight algorithm proposed in the second section is used to calculate each subset, ensuring that the computational scale of single-machine path planning is reduced while the scale of the target point set remains unchanged, achieving a balance between computational complexity and optimized path. The simulation data show that the EDTSP model proposed in this paper has at least 8.4% and 17.5% improvement in the sparse interval for takeoff from different points.