2.1 Design of formation control and agent tracking

The motion dynamics of each agent in a multi-agent system, designed to establish and maintain a V-shaped formation, are described using the simplified single-integrator model [Reference Kahagh, Pazooki and Haghighi7]. This V-formation is flexible, allowing the apex angle (α) to vary; in the limiting case where the angle increases to 180°, the formation degenerates into a linear structure. In this model, the rate of change of each agent’s position is directly determined by its control signal, as shown in Equation (1):

(1) \begin{equation}\dot{p}_{i}=u_{i}\end{equation}

where $p_{i}$ is the position vector of agenti in 2D or 3D space, and $u_{i}$ is the control input signal applied to agent $i$ .

This mathematical expression assumes homogeneous agents and disregards their physical dimensions, providing a suitable framework for analysing and developing formation control algorithms.

The control signal $u_{i}$ , which is central to this framework, is designed to guide each agent toward its desired position within the final formation. To achieve an efficient and applicable control law for UAVs, the interaction between agents is defined by the two-state function shown in Equation (2) [Reference Kahagh, Pazooki and Haghighi7]:

(2) \begin{equation} u_{i}= \begin{cases} -\displaystyle\sum\limits_{j=1}^{N}w_{i}({p_{ij}}-d_{ij})^{\beta } e_{ij} \qquad \qquad \qquad \qquad \ p_{ij}\geq d_{ij}\\[12pt] \displaystyle\sum\limits_{j=1}^{N}w_{i}(d_{ij}-p_{ij})^{\beta }\left({1 \over (p_{ij}-d_{c})^{2}}+1\right) e_{ij} \qquad \ p_{ij}\lt d_{ij} \end {cases} \end{equation}

where $p_{ij}=||p_{i}-p_{j}||$ and $e_{ij}={p_{i}-p_{j} \over p_{ij}}$ .

Figure 1.

The desired positions of agents in the target formation.

In Equation (2), $d_{ij}$ is the desired distance between agents $i$ and $j$ in the defined formation; $w_{i}$ is a weighting coefficient used to regulate the interaction intensity; and $d_{c}$ is a critical distance representing the non-violable clearance around each agent. The parameter $\beta$ is incorporated to achieve smooth transitions in the inter-agent force intensities, a feature essential for implementing advanced controls, such as sliding mode control on pre-determined paths.

To facilitate the coordinated group movement and guidance of the agents, a virtual leader (agent N + 1) is introduced, which dictates the overall path trajectory [Reference Kahagh, Pazooki and Haghighi7]. The interaction between the leader and the followers is established using potential functions, similar to the inter-agent communication (Equation 2), but with distinct weighting coefficients. Consequently, the complete system, comprising N follower agents and one leader, forms an N + 1 agent system structured in a modified Y-shaped formation. This structure allows for the definition of integrated geometric and control relations between the leader and all followers.

The control mechanism operates by having the leader move along a pre-determined path, generating an effective force field in the state space. This field attracts the followers, compelling them to track the leader while maintaining the defined geometric relationships. As the leader moves, the entire group transitions cohesively along the intended path without compromising the formation structure.

Therefore, the position of each agent i at any moment, based on Equation (1), is given by Equation (3):

(3) \begin{align}\dot{p}_{i}=\sum\limits_{j=1,j\neq i}^{N+1}u_{j} \end{align}

To effectively utilise the control signal defined in Equation (2), the desired distance between every pair of agents in the pre-defined formation (illustrated in Figure 1) must be calculated. A systematic approach to determine the desired geometric positions of each agent in a V-shaped formation involves defining the desired distance between any pair of agents i and j, which forms the basis of the potential functions in the formation control [Reference Kahagh, Pazooki and Haghighi7]. For the inter-agent distance $d_{ij}$ , Equation (4) is used:

(4) \begin{align}& if\left(i,j\leq {N+1 \over 2} \ or \ i,j\geq {N+1 \over 2}\right)d_{ij}=\left| i-j\right| \times D \nonumber\\&\quad else \,d_{ij}=\sqrt{\left(\left(| i-j| \right)\times D\times \sin \left({\alpha \over 2}\right)\right)^{2}+\left(\left(N+1-(i+j)\right)\times D\times \mathit{\cos } \left({\alpha \over 2}\right)\right)^{2}}\end{align}

Furthermore, the desired distance between the leader (agent N + 1) and any follower agent k is determined as:

(5) \begin{align}d_{(N+1)k}&=\sqrt{\left(\left| {(N+1) \over 2}-k\right| \times D\times cos\left({\alpha \over 2}\right)+D_{L}\right)^{2}+\left(\left| {(N+1) \over 2}-k\right| \times D\times \sin \left({\alpha \over 2}\right)\right)^{2}}\nonumber\\ k&=1,2,\ldots ,N\end{align}

In Equations (4) and (5), D represents the desired lateral distance between adjacent agents, $\alpha$ is the apex angle of the V-shaped formation, and $D_{L}$ is the desired distance between the central agent ((N + 1)/2) and the leader (N + 1).

This approach allows for precise definition of the geometric formation with high tunability; by changing the input parameters $(D,D_{L},\alpha )$ , different forms of V-shaped formations with variable angles can be created. Equations (2) through (5) can be used for any number of agents to design the UAV formation flight trajectory and determine each agent’s position. The only required adjustment is to modify the weight coefficient of the inter-agent forces due to the change in the problem’s scale as the number of agents increases.

2.2 Obstacle crossing mechanism

Figure 2 illustrates the obstacle crossing mechanism [Reference Kahagh, Pazooki, Haghighi and Asadi19]. $p_{i}$ represents the virtual path of a UAV that would collide with an obstacle having radius $R_{O}$ and centre $O$ at points $B_{1}$ and $B_{2}$ . Point H is located at the midpoint of these two collision points, and the centre $O$ defines a circle with radius $R_{a}$ , known as the activation circle. When the UAV’s distance from the obstacle’s centre is less than this radius (i.e. the UAV enters the activation circle), the obstacle crossing mechanism is activated, causing the UAV to deviate from the virtual path to safely bypass the obstacle.

Figure 2.

Obstacle crossing mechanism [Reference Kahagh, Pazooki, Haghighi and Asadi19].

Considering a case where a larger number of UAVs travelling along the virtual path $p(i)$ collide with an obstacle (as shown in Figure 3), each UAV theoretically collides with the obstacle at two points. The midpoint $H(i)$ is calculated for each virtual path using Equation (6).

Figure 3.

Virtual path of multiple UAVs colliding with an obstacle.

(6) \begin{equation}H{(i)}={B_{1}(i)+B_{2}(i) \over 2}\end{equation}

To pass an obstacle, each UAV can deviate from its virtual path in one of two directions, following a curve before returning to the virtual path, as shown in Figure 2. However, only one of these two directions provides the optimal path and a shorter curve for obstacle clearance, unless the virtual path passes precisely through the obstacle’s centre, in which case both directions are optimal.

Considering $H(i)=\Big[\begin{array}{c} x_{H}(i)\\ y_{H}(i) \end{array}\Big]$ and $O=\Big[\begin{array}{c} x_{O}\\ y_{O} \end{array}\Big]$ , paths where $y_{H}(i)\geq y_{O}$ are included in Direction Group 1 and paths where $y_{H}(i)\lt y_{O}$ are included in Direction Group 2 (Figure 3). To optimise the obstacle crossing curve, the distance of point $H(i)$ from the obstacle centre $O$ is calculated for the paths in each direction. The path corresponding to the smaller distance, based on $[\! \begin{array}{cc} min(||H_{1}(i)-O||) {} min(||H_{2}(i)-O||) \end{array} \!]$ is selected. Thus, only one path is selected from each direction, resulting in a total of two path curve optimisations.

The optimisation of the selected paths is performed only for the points that are within a certain range of the obstacle. This is done to improve simulation performance by reducing the search space and variables. For this purpose, the activation circle (Figure 2) is used, and only the path points within this circle are optimised. The centre $H(i)$ and radius $R_{a}(i)$ of this circle are calculated according to Reference [Reference Kahagh, Pazooki, Haghighi and Asadi19] using Equation (7).

(7) \begin{equation}R_{a}(i)=\sqrt{4{R}_{m}^{2}-(2R_{m}-R_{O}+||\overrightarrow {OH}||)^{2}+||\overrightarrow{OH}||^{2}}\end{equation}

where $\overrightarrow{OH}(i)=H(i)-O$ is the distance between $H(i)$ and $O$ , and $\vec {e}_{OH}(i)={\overrightarrow{OH}(i) \over ||{\overrightarrow{OH}}||},$ and $R_{m}$ is the minimum turning radius of the UAV, which depends on the minimum speed ( $v_{min}$ ) and the maximum bank angle ( $\phi _{max}$ ): $R_{m}={{v}_{min}^{2} \over g\mathit{\tan } \left(\phi _{max}\right)}$ .

The optimal obstacle crossing curve is created by moving each virtual path point within the activation circle in the direction of the unit vector $\vec {e}_{OH}(i)$ (Equation (7)). This process generates the actual path points of the UAV. The number of virtual path points within the activation circle may be large, depending on the time step used for solving the dynamic equations. For this reason, a sampling rate is applied, selecting one point for every n virtual path points to define the real path for obstacle crossing. However, since this deviation from the virtual path must also be applied to the path of neighbouring UAVs to maintain formation, the number of reduced points must be returned to the initial quantity, which is achieved by interpolating the reduced points.

2.3 Other equations

In some cases, it is necessary to determine the projection of a point along a vector onto a circular circumference. Let the equation of the line passing through point P in the direction of the unit vector $\hat{u}={u \over \left|\left|u\right|\right|}$ be given by Equation (8), where t is the line parameter, and the equation of the circle with cente C and radius R be given by Equation (9):

(8) \begin{equation}L(t)=P+t.\hat{u}\end{equation}

(9) \begin{equation}\left|\left|X-C\right|\right|^{2}=R^{2}\end{equation}

By substituting Equation (8) into Equation (9) and considering $P-C=\sigma$ , Equation (10) is obtained:

(10) \begin{equation}||\sigma +t.\hat{u}||^{2}=R^{2}\end{equation}

Expanding and simplifying this equation yields the standard quadratic equation in Equation (11):

(11) \begin{equation}t^{2}+2(\sigma .\hat{u})t+(||\sigma ||^{2}-R^{2})=0\end{equation}

If the required solution to this equation is $t_{\textit{required}}$ , the point projected onto the circle is given by Equation (12):

(12) \begin{equation}P_{\textit{projected}}=P+t_{\textit{required}}.\hat{u}\end{equation}

In some cases, it is also necessary to transfer the coordinates of the path points from the original coordinate system to a rotated coordinate system that has an angle θ with respect to the original coordinate system.

The coordinate transformation matrix from the original frame to the rotated frame ( ${C}_{m}^{r}$ ) is expressed by Equation (13):

(13) \begin{equation}{C}_{m}^{r}=\left[\begin{array}{l@{\quad}c} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\!,\, {P}_{i}^{r}={C}_{m}^{r}{P}_{i}^{m}\end{equation}

To calculate the radius of rotation of the paths created by the algorithm, the method of a circle passing through three non-collinear points is used. The general equation of a circle in the Cartesian coordinate system, its radius and centre are expressed as follows:

(14) \begin{align}x^{2}+y^{2} + Dx \, + \, & Ey+F=0 \nonumber\\[4pt] \text{center: }O\left(-{D \over 2},-{E \over 2}\right),\,& \text{Radius: } r={1 \over 2}\sqrt{D^{2}+E^{2}-4F}\end{align}

Given three points $A(x_{1},y_{1}$ ), $B(x_{2},y_{2})$ and $C(x_{3},y_{3})$ , the matrix equation system $M\vec {\kappa }=\vec {b}$ is obtained:

(15) \begin{equation}\left[\begin{array}{c@{\quad}c@{\quad}c} x_{1} & {} y_{1} & 1\\ x_{2} & y_{2} {}& 1\\ x_{3} & y_{3} {}& 1 \end{array}\right]\left[\begin{array}{c} D\\ E\\ F \end{array}\right]=\left[\begin{array}{c} -\left({x_{1}}^{2}+{y_{1}}^{2}\right)\\[3pt] -\left({x_{2}}^{2}+{y_{2}}^{2}\right)\\[3pt] -\left({x_{3}}^{2}+{y_{3}}^{2}\right) \end{array}\right]\end{equation}

Solving this equation yields $\vec {\kappa }=\Bigg[\begin{array}{c} D\\ E\\ F \end{array}\Bigg]$ . Finally, the radius of the circle passing through these three points is obtained using Equation (14).

To quantitatively assess the smoothness of the generated paths for UAV navigation, we employ a normalised bending energy criterion. Given a discrete path consisting of $N$ sequential waypoints $\{p_{1}, p_{2},\ldots ,p_{N}\}\subset \mathbb{R}^{d}$ (with $d=2 \ or \ 3$ ), we first define the step vectors as $v_{i}=p_{i+1}-p_{i}$ for $i=2,\ldots, N-1$ . The smoothness metric is then formulated as:

(16) \begin{equation}S_{\textit{smooth}}={1 \over L_{\textit{total}}}\sum \limits_{i=2}^{N-1}{\left|\left|v_{i}-v_{i-1}\right|\right|^{2} \over \left|\left|v_{i-1}\right|\right|+\left|\left|v_{i}\right|\right|}\end{equation}

where $L_{\textit{total}}={\sum }_{i=2}^{N-1}||v_{i-1}||$ denotes the total arc length of the path. The numerator $||v_{i}-v_{i-1}||^{2}$ represents the squared variation in the velocity vectors, which is directly proportional to the magnitude of the centripetal acceleration at waypoint $p_{i}$ . The denominator $||v_{i-1}||+||v_{i}||$ serves as an adaptive weighting factor that naturally accounts for non-uniform step lengths, ensuring that abrupt directional changes occurring over longer segments are penalised more heavily, in accordance with the flight dynamics and manoeuvrability constraints of UAVs. The normalisation by $L_{\textit{total}}$ renders the metric scale-invariant, enabling fair comparisons across different scenarios with varying path lengths. This criterion is particularly well-suited for densely sampled paths (typically comprising 1,000–1,500 points), as it approximates the continuous curvature energy while exhibiting high numerical stability and low sensitivity to local noise.

3.1 Overview of the standard GWO algorithm

The grey wolf optimisation (GWO) algorithm is a swarm intelligence algorithm inspired by the social behaviour and collective hunting of grey wolves. It simulates the leadership hierarchy of wolves (including alpha, beta, delta and omega wolves) and models the process of encircling, hunting and attacking prey to find the best possible solution in the problem space [Reference Mirjalili, Mirjalili and Lewis20, Reference Faris, Aljarah, Al-Betar and Mirjalili21]. In this method, the alpha wolf is considered the current best solution, the beta and delta wolves are secondary advisors, and the other wolves are considered followers. The position of each wolf is updated at each iteration using adaptive parameters ( $a$ , A and C) to establish a balance between exploration and exploitation [Reference Mirjalili, Mirjalili and Lewis20, Reference Faris, Aljarah, Al-Betar and Mirjalili21]. The mathematical model of GWO is based on the following equations:

(17) \begin{equation}\vec {D}=\big| \vec {C}.{\overrightarrow{X_{p}}}(t)-\vec {X}(t)\big|\end{equation}

(18) \begin{equation}\vec {X}(t+1)={\overrightarrow{X_{p}}}(t)-\vec {A}.\vec {D}\end{equation}

where ${\overrightarrow{X_{p}}}$ is the prey position, $\vec {X}$ is the wolf position, and $\vec {A}$ and $\vec {C}$ are coefficient vectors that are dynamically adjusted during the optimisation process:

$\vec {A}=2a.{\overrightarrow{r_{1}}}-a$ and $\vec {C}=2.{\overrightarrow{r_{2}}}$ .

${\overrightarrow{r_{1}}}$ and ${\overrightarrow{r_{2}}}$ are random vectors in the interval [0, 1], and $a$ is a control parameter that decreases linearly from 2 to 0 during the iterations, ensuring a gradual transition from the exploration to the exploitation phase. The parameter $\vec {A}$ plays a key role in balancing exploration and exploitation. Values of $| \vec {A}| \gt 1$ promote an expanded search (exploration), while values of $| \vec {A}| \lt 1$ lead to a localised search (exploitation) [Reference Mirjalili, Mirjalili and Lewis20, Reference Faris, Aljarah, Al-Betar and Mirjalili21].

Pseudocode of the GWO algorithm

Algorithm Long description

The pseudocode outlines the Grey Wolf Optimizer algorithm. It starts with initializing input parameters like the number of search agents, maximum iterations, objective function, search space bounds, and problem dimension. The output is the best solution found. The algorithm involves randomly initializing the population of wolves within the search space, evaluating their fitness values, and determining the top three wolves. It then iteratively updates the control parameter, generates random vectors, computes coefficient vectors, estimates candidate positions, and updates wolf positions. The fitness values of all wolves are evaluated, and the top three wolves are updated accordingly. This process repeats until the maximum number of iterations is reached, after which the best solution is returned.

3.2 Decision variable definition and problem encoding

In this study, each obstacle-crossing path is considered as a `wolf’. This wolf is defined by a set of points from the real path (or obstacle-crossing curve) that is generated by adjusting the positions of corresponding points on the virtual path that lie within the activation circle. It is important to note that these real path points are not all the virtual points within the activation circle but are selected by applying the sampling rate described in Section 2. However, the real path points in each wolf are not directly represented as coordinates $\Big[\begin{array}{c} x\\ y \end{array}\Big]$ . Instead, a wolf’s characteristic is defined by the distance between each virtual path point and its corresponding real path point. Therefore, each wolf contains a $1\times n$ vector of distances between two points of virtual and real path points. These n distances serve as the decision variables in the optimisation problem. If $d_{i}$ is the distance of the i-th point on the virtual path from its corresponding point on the real path, the corresponding real path point can be calculated. This is achieved by shifting the virtual path point along the vector $\vec {e}_{OH}(i)$ by the distance $d_{i}$ , as illustrated in Figure 4 and shown in Equation (19).

(19) \begin{align}\vec {p}_{ri}=\vec {p}_{i}+d_{i}\vec {e}_{OH}(i) \\[-22pt] \nonumber \end{align}

$\vec {e}_{OH}(i)$ is calculated according to Equation (7). $\vec {p}_{ri}$ represents the optimised points on the real path or obstacle-crossing curve, and $\vec {{\rm p}}_{{\rm i}}$ represents the corresponding points on the virtual path. It should be noted that points ${\rm A}_{1}$ and ${\rm A}_{2}$ , which are the first and last points of the virtual path in the activation circle (Figure 4), will not be optimised.

Figure 4.

Specification of the decision variable and derivation of the obstacle-crossing curve.

The lower limit of the decision variables is zero, as indicated in Figure 4. The upper limit is defined by a condition where the actual path is tangent to the circular obstacle, which is given by $\left(R_{O}+{d_{c} \over 2}+S.M\right)-||{\overrightarrow{OH}}||+d_{0}$ . A safety margin (S.M) is considered to prevent the UAV’s wings from colliding with the obstacle. $d_{0}$ is a small value (on the order of a few meters) to ensure a small distance between the upper limit of the decision variables and the safe boundary of the obstacle, allowing the algorithm to converge to the optimal value in this area during the iteration process. The paths created by the algorithm must possess several important properties to be both optimal and dynamically executable by the UAV. For optimality, the path must be shorter in length and must not intersect the obstacle boundary. Naturally, the closer the actual path is to the virtual path, the shorter its length will be. However, an essential feature for the UAV’s dynamic performance is a suitable turning radius for the trajectory. Failing to consider this feature could result in finding shorter, collision-free paths that are dynamically impossible to fly. Therefore, the turning radius of the trajectory must not be less than the minimum turning radius calculated according to Equation (7). Based on these requirements, the cost function for minimisation is defined by Equation (20):

(20) \begin{equation}J=\omega _{l}.\sum\limits_{n=1}^{N_{p}}\left(\Delta L_{n}\right)+\omega _{tr}.N_{tr}+\omega _{o}.N_{o}\end{equation}

In the above equation, $N_{p}$ is the number of consecutive points in the path $p_{ri}$ within the activation circle, and $\Delta L_{n}$ is the distance between two consecutive points. $N_{tr}$ is the number of points on the real path that do not have a suitable turning radius, and $N_{o}$ is the number of points on the real path within the activation circle whose distance to the obstacle centre is smaller than $\left(R_{O}+{d_{c} \over 2}+S.M\right)$ . $\omega _{l}$ , $\omega _{tr}$ and $\omega _{o}$ are penalty coefficients for the path length, the path’s turning radius constraint and the path’s collision with the obstacle, respectively. Although a sampling rate is determined to increase the convergence speed by selecting fewer points for optimisation, the search space remains large, and various states can occur for each path point. To manage this, a constraint is considered such that the real path points that collide with the obstacle are transferred to the obstacle environment, using Equations (8) through (12). This ensures that, in the initial iterations of the optimisation algorithm, a collision-free path is created, requiring only the radius of rotation and the path length to converge to the desired values in subsequent iterations. Appropriate initialisation in creating the initial population of wolves can also accelerate the convergence speed. If the initial initialisation of the decision variables is performed to create paths with a reasonable radius of rotation or appropriate approximate curvature, it is naturally possible to manage the interaction of exploration and exploitation. This is achieved by defining a suitable function for the parameter $a$ in the GWO algorithm so that more iteration operations are spent on local search (exploitation), accelerating convergence to a solution close to the optimal one. Determining an appropriate value for the penalty coefficients in the cost function can also be useful in properly directing the responses. All these actions are implemented to constrain the problem and reduce the search space, guiding the solutions toward the appropriate direction at each iteration.

3.3 Improved GWO with exponential decay strategy

In the standard GWO algorithm, the balance between exploration and exploitation is governed by the control parameter $a$ , which decreases linearly from 2 to 0 over the course of iterations, as given by $a_{1}$ . This linear decay ensures a smooth and gradual transition from global search (exploration) to local search (exploitation) [Reference Mirjalili, Mirjalili and Lewis20, Reference Faris, Aljarah, Al-Betar and Mirjalili21]. However, the linear strategy does not account for the specific characteristics of the optimisation problem, such as the dimension of the search space or the complexity of the objective function. In problems with a large number of decision variables – such as the obstacle avoidance trajectory optimisation considered in this study, where the number of optimised points ranges from 30 to 50 – the linear decay may allocate insufficient iterations to the exploitation phase, resulting in suboptimal convergence near the end of the optimissation process.

To address this limitation, we consider the functions given in Equation (21) for the parameter $a$ , defined as:

(21) \begin{align} a_{1}&=2\left(1-{t \over T}\right) \nonumber\\[3pt] a_{2}&=2e^{\left(-5\left({\frac{t}{\vphantom{x} T}}\right)^{2}\right)} \nonumber\\[3pt] a_{3}&=2e^{\left(-10\left({\frac{t}{\vphantom{x}T}}\right)^{2}\right)}\nonumber\\[3pt] a_{4}&=2\left(1-{\ln (t) \over \ln \left(T\right)}\right) \nonumber\\[3pt] a_{5}&={2 \over t}+e^{-3\left({\frac{t}{\vphantom{x} T}}\right)}\left({\ln \left(t/T\right) \over \ln \left(1/T\right)}-{1 \over t}\right)\end{align}

In this equation, $t$ represents the current iteration, and $T$ represents the maximum number of iterations.

Among the candidate functions presented in Equation (21), the exponential decay $a_{3}$ was selected after systematic evaluation. As illustrated in Figure 5, this function exhibits a rapid initial decline, followed by a prolonged period of very small values. This characteristic fundamentally alters the search behaviour of the algorithm: the exploration phase is significantly shortened, while the exploitation phase is substantially extended. In practical terms, the algorithm spends fewer iterations on a coarse global search and devotes a larger proportion of the total iterations to fine-tuning the solutions within promising regions of the search space.

Figure 5.

Changes of ${\rm a}$ versus iteration.

The effectiveness of this strategy is rooted in the specific structure of the optimisation problem. By constraining the initial population to generate paths with reasonable curvature – achieved through appropriate initialisation and penalty terms in the cost function – the search space is effectively narrowed from the outset. Consequently, the algorithm does not require an extended exploration phase to identify promising regions; instead, it benefits from a more aggressive transition to exploitation, where fine adjustments of path points are performed to minimise path length while satisfying turning radius and collision avoidance constraints.

The superiority of the exponential decay over the linear strategy is twofold:

1. I. Convergence Speed: The rapid decline of ${\rm a}_{3}$ forces the coefficient vector $\vec {{\rm A}}$ to assume values $| \vec {{\rm A}}| \lt 1$ earlier in the iteration process. Since $| \vec {{\rm A}}| \lt 1$ corresponds to localised searching (exploitation), the algorithm converges to a near-optimal solution in fewer iterations compared to the linear case, where a decreases slowly and maintains $| \vec {{\rm A}}| \gt 1$ for a longer duration. This is particularly beneficial when the number of iterations is limited, as the algorithm can allocate the majority of its computational budget to refining solutions rather than exploring distant regions of the search space.

2. II. Solution Quality: The extended exploitation phase allows the algorithm to thoroughly search the local neighbourhood of the current best solutions. Given that the decision variables are continuous and the objective function is multimodal due to the trade-offs among path length, turning radius and collision avoidance, the prolonged local search enables the algorithm to escape minor local minima and converge to a higher-quality solution.

It is also worth noting that the exponential decay function $a_{3}$ was selected after systematic comparison with other candidate functions, including $a_{2}$ , $a_{4}$ and $a_{5}$ . While $a_{4}$ and $a_{5}$ exhibit an even faster initial decay, reaching very small values within fewer than 100 iterations, this rapid transition to the exploitation phase introduces several critical drawbacks. First, the premature loss of population diversity restricts the algorithm’s ability to explore distant regions of the search space. Second, this early convergence often traps the search agents in suboptimal local minima, resulting in lower-quality final solutions. Third, the insufficient exploration phase prevents the algorithm from adequately sampling various combinations of path-point displacements, which are essential for generating paths with appropriate curvature and turning radius. In contrast, $a_{3}$ provides a more balanced decay profile: it accelerates the transition to exploitation sufficiently to enhance convergence speed, yet maintains a brief but adequate exploration window that preserves population diversity and prevents premature stagnation.

Finally, since the points are reduced by the sampling rate, after the GWO algorithm finishes, the number of points must be increased back to the initial value so that other UAVs can manoeuvre according to the actual path created to maintain the formation. For this purpose, a spline is fitted to the optimised points, which simultaneously performs interpolation and ensures that the created actual path is a smooth curve with an appropriate turning radius.

6.1 Limitations of idealised assumptions

The simulation framework presented in Sections 3 and 4 employs simplifying assumptions common in initial algorithmic validation [Reference Kahagh, Pazooki and Haghighi7, Reference Kahagh, Pazooki, Haghighi and Asadi19]: (I) perfect inter-agent communication with negligible latency; (II) ideal state estimation without sensor noise; and (III) wind-free environmental conditions. While these assumptions facilitate algorithmic development, practical deployment in operational missions such as cooperative surveillance and SAR (Section 5) requires careful consideration of their effects on the tracking of the pre-computed optimised trajectories.

As described in detail in Section 5, the proposed methodology follows a hierarchical two-stage design:

1. I. Stage 1 (Pre-Flight - Offline Design):

   • • The virtual leader’s kinematic trajectory is defined using Equations (22), (23).

   • • The formation control law based on APF in Equation (2) establishes and maintains the desired V-shaped formation geometry with prescribed inter-agent distances $d_{ij}$ and apex angle $\alpha$ . At this stage, the nominal paths for all agents are determined.

   • • Collision status with known obstacles is evaluated for each agent’s nominal path, and the agents whose paths intersect the obstacle boundary are identified.

   • • The modified GWO algorithm optimises the obstacle-crossing trajectories for the affected agents by adjusting the deviation distances $d_{i}$ in Equation (19). The cost function in Equation (20) simultaneously minimises path length, ensures appropriate turning radius and prevents obstacle collisions. Importantly, the number of agents that collide with the obstacle is unrestricted – the algorithm handles scenarios ranging from a single affected agent to multiple UAVs on either side of the formation, as demonstrated in Figure 8.

   • • Spline interpolation generates smooth, dynamically feasible trajectories for all agents.

2. II. Stage 2 (In-Flight - Online Execution):

   • • The pre-computed optimal trajectories are uploaded to the UAVs as reference paths.

   • • Each agent employs an APF-based tracking controller individually to follow its assigned reference trajectory. This controller utilises the potential function framework to generate appropriate control inputs that guide each UAV along the pre-computed path while compensating for minor deviations.

Key Observation: As clearly stated in Section 5, the APF-based formation control law in Equation (2) is employed only during the pre-flight design phase to establish the nominal formation geometry and determine the desired positions of all agents relative to the virtual leader. During actual flight, the UAVs execute the pre-computed trajectories using individual APF-based tracking controllers, not a formation-level APF controller that actively adjusts the formation in response to environmental disturbances. Consequently, the kinematic model in Equation (1) describes the intended motion during the design phase, but the actual motion during flight is subject to tracking errors arising from environmental disturbances. The formation geometry and inter-agent distances $d_{ij}$ that were established during the offline design phase are expected to be maintained only if each UAV accurately tracks its assigned reference trajectory using its individual APF-based tracking controller.

6.1.1 Wind disturbances

During flight execution, wind causes the UAV to deviate from its pre-computed reference trajectory. The actual kinematic relationship differs from the design model:

(24) \begin{equation}{\dot{p}}_{i}^{\textit{actual}}={\dot{p}}_{i}^{\textit{reference}}+W(p_{i},t)\end{equation}

where ${\dot{p}}_{i}^{\textit{reference}}$ is the desired velocity from the pre-computed trajectory, and $W(p_{i},t)$ denotes the wind velocity vector. This mismatch introduces:

1. I. Tracking Errors: Wind causes the UAV to deviate from its pre-computed reference trajectory, potentially violating the desired inter-agent distances $d_{ij}$ that were established during the offline formation design using the APF method. While the individual APF-based tracking controller attempts to keep each UAV on its reference path, its performance may be degraded under strong wind conditions.

2. II. Impact on Obstacle Avoidance: The optimised obstacle-crossing trajectory was computed assuming no wind. Wind-induced drift may cause the UAV to approach the obstacle boundary more closely than intended, potentially compromising the safety margin $S.M$ defined in the cost function. The minimum turning radius $R_{m}={v}_{min}^{2}/(g\mathit{\tan } (\phi _{max}))$ may also become infeasible if wind assists or opposes the turn.

6.1.2 Sensor noise

During trajectory tracking, the UAV relies on onboard sensors for state estimation. The individual APF-based tracking controller uses these measurements as feedback. Table 6 presents the noise characteristics that may influence the trajectory tracking performance.

Table 6.

Sensor noise sources and their effects on reference trajectory following

Table 6. Long description

A table with three columns and four rows, including a header row. The columns are labeled Noise source, Typical magnitude, and Effect on trajectory tracking. The rows list different noise sources and their characteristics. Row 1: GPS position error, 1-5 m (standard), 0.1-1 m (RTK), Errors in position feedback for the APF tracking controller. Row 2: IMU velocity drift, 0.1-0.5 m/s, Accumulated position errors. Row 3: Airspeed measurement error, 0.5-2 m/s, Errors in wind estimation.

6.2 Proposed adaptations for real-world implementation

Since the GWO optimisation is performed offline (as described in Section 5) and the optimised trajectories are provided to the UAVs as reference paths, the proposed adaptations focus on two complementary aspects: (I) enhancing the offline optimisation to account for expected uncertainties, and (II) improving the online tracking through individual APF-based feedback control to ensure accurate execution of the pre-computed trajectories.

6.2.2 In-flight phase enhancements (online tracking control)

Following the offline GWO optimisation and spline interpolation (Section 3), the optimised trajectories are uploaded to the UAVs as reference paths. During flight, each UAV employs an individual APF-based tracking controller to follow its reference trajectory. The following enhancements ensure accurate tracking:

1. I. Wind Estimation and Feedforward Compensation: Since wind is the primary source of tracking errors, a real-time wind estimation filter can be integrated into each UAV’s tracking control loop:

   (28) \begin{equation}\hat{W}(t)=\dot{\hat{p}}_{i}(t)-v_{\textit{airspeed}}(t)\end{equation}
   where $\dot{\hat{p}}_{i}(t)$ is the estimated velocity from GPS, and $v_{\textit{airspeed}}(t)$ is the measured airspeed. The estimated wind is used to generate feedforward compensation for the individual APF-based tracking controller:
   (29) \begin{equation}{u}_{i}^{comp}={u}_{i}^{APF}-\hat{W}(t)\end{equation}
   where ${u}_{i}^{APF}$ is the control signal generated by the individual APF-based tracking controller. This reduces steady-state tracking errors without requiring gain increases that could excite unmodeled dynamics.

1. II. Kalman Filter for State Estimation: Implement a discrete-time Kalman filter for each agent:

   (30) \begin{equation}\hat{p}_{i}{\left(k|k\right)}=\hat{p}_{i}{\left(k|k-1\right)}+K_{k}{[z_{i}(k)-H\hat{p}_{i}(k|k-1)]}\end{equation}
   where $z_{i}(k)$ is the noisy measurement vector (GPS/IMU), and $K_{k}$ is the optimal Kalman gain. This provides smoothed position estimates to the individual APF-based tracking controller, reducing the effects of sensor noise on tracking performance.

1. III. Adaptive Formation Expansion: If the estimated tracking uncertainty exceeds a threshold during flight, the formation can be dynamically expanded by increasing the inter-agent distance $d_{ij}$ :

   (31) \begin{equation} {d}_{ij}^{\textit{adapt}}={d}_{ij}^{\textit{nominal}}+\delta _{\textit{uncertainty}}\end{equation}

This maintains safety without requiring online GWO re-optimisation, which is computationally infeasible for onboard processors. The individual APF-based tracking controllers will then operate with larger safe margins. Consistent with the two-stage design philosophy described in Section 5, the proposed implementation maintains the following separation (Table 7):

Table 7.

Summary of implementation architecture

Table 7. Long description

A table summarizing the implementation architecture for UAV missions. The table has four columns: Phase, Activity, Timing, and Responsibility. It contains ten rows detailing various activities and their respective timings and responsibilities. The phases are divided into Pre-Flight and In-Flight. Pre-Flight activities include virtual leader trajectory definition, formation design using APF, collision detection, GWO optimization, spline interpolation, and uploading optimized trajectories. In-Flight activities include wind estimation, Kalman filtering, and trajectory tracking. Each activity is assigned a specific responsibility such as Mission planner, GWO algorithm, Post-processing, Ground station, Individual APF-based tracking controller, and Each UAV’s onboard controller.

The proposed adaptations preserve the core contribution of this work (the modified GWO algorithm for offline optimisation of obstacle avoidance trajectories) while adding the robustness necessary for practical implementation. During flight, each UAV employs an individual APF-based tracking controller to follow its pre-computed reference trajectory. The in-flight enhancements (wind estimation with feedforward compensation, Kalman filtering and adaptive formation expansion) serve only to ensure accurate tracking of the pre-computed paths, not to modify the optimised trajectories online. The APF-based tracking controller operates independently for each agent, generating appropriate control inputs based on the difference between the agent’s current position and its desired position on the reference trajectory.

This separation of offline optimisation and online tracking control offers several practical advantages:

• • The computationally intensive GWO optimisation does not burden onboard processors.

• • Mission planners can validate trajectories before deployment.

• • Contingency paths can be pre-computed for multiple obstacle configurations.

• • The formation can adapt to uncertainty through expansion rather than replanning.

• • The number of affected agents does not affect the computational complexity of the individual tracking controllers.

• • Each UAV’s APF-based tracking controller operates independently, providing inherent scalability.