2.1 Problem formulation

A novel formation control algorithm is developed based on the virtual structure, leader-follower, and artificial potential fields methods to benefit from the advantages found in these approaches. All guidance and navigation commands are given to the leader and all group members follow the leader.

Figure 1. V-shape formation structure.

The minimum communication information of the members is applied, i.e. merely the relative position information in the inertial coordinate frame. A mobile coordinate frame is defined with its origin placed on the leader. In the two-dimensional space, one of the axes of the introduced frame is tangent to the leader’s movement path, while the other one is perpendicular to the path. The position of each member is defined virtually with the desired formation in the mobile coordinate frame. Subsequently, the geometric positions of each virtual member are transferred to the inertial coordinate frame using transformation matrices. Finally, a desired virtual structure is obtained in the inertial coordinate frame, to which, the real members should correspondingly approach. It is necessary to obtain the angle between the tangent line to the leader’s movement path with one of the inertial coordinate frame axes, as follows, to find the orientation angle of the mobile coordinate frame relative to the inertial coordinate frame (Fig. 1).

(1) \begin{align}\theta (t) = {\rm{atan}}\left({\frac{{{y_l}(t) - {y_l}({t - \Delta t})}}{{{x_l}(t) - {x_l}({t - \Delta t})}}}\right)\end{align}

$\theta(t)$ : The orientation angle of the mobile coordinate frame relative to the inertial coordinate frame at time t.

$\left[ {\begin{array}{*{20}{c}}{{x_l}\left( t \right)}\\{{y_l}\left( t \right)}\end{array}} \right]$ : The relative position vector of the leader in the inertial coordinate frame.

The focus is on the formation and maintenance of V-shaped formation, while following the virtual leader in the simulation process. The formation algorithm is designed such that it would be adjustable in terms of the number of group members and the value of the V-shaped angle. The formation of the straight-line shape can also be achieved by increasing the value of the V-shaped angle to $\;{180^{\circ}}$ .

The positions of all virtual members should be calculated first in relation to the mobile coordinate frame according to Fig. 1, which will be then transferred to the inertial coordinate frame, as follows.

(2) $$\eqalign{ & \bar p_i^L = \left[ {\matrix{ {D\;\left| {{{\left( {n + 1} \right)} \over 2} - i} \right|\cos \left( {\pi - {\alpha \over 2}} \right) - {D_L}} \cr {D\;\left| {{{\left( {n + 1} \right)} \over 2} - i} \right|\sin \left( {\pi - {\alpha \over 2}} \right)\;\;\;\;\;\;\;\;} \cr } } \right]\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;If\;\;i{\rm{ < }}{{\left( {n + 1} \right)} \over 2} \cr & \bar p_i^L = \left[ {\matrix{ {D\;\left| {{{\left( {n + 1} \right)} \over 2} + i} \right|\cos \left( {\pi - {\alpha \over 2}} \right) - {D_L}} \cr {D\;\left| {{{\left( {n + 1} \right)} \over 2} + i} \right|\sin \left( {\pi - {\alpha \over 2}} \right)\;\;\;\;\;\;\;\;} \cr } } \right]\;\;\;\;\;\;\;\;\;\;\;\;If\;\;i < {{\left( {n + 1} \right)} \over 2} \cr & \matrix{ {\bar p_i^L = \left[ {\matrix{ { - {D_L}} \cr 0 \cr } } \right]\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;If\;\;\;i = {{\left( {n + 1} \right)} \over 2}} \hfill \cr } \cr} $$

n: Number of agents.

D: Desired inter-agent distance.

${D_L}$ : Desired distance between middle agents ( $({n + 1})/2$ ) and leader $\left( {n + 1} \right)$ .

$\bar p_i^L$ : The position vector of the virtual member in the mobile coordinate frame. The index (L) is the symbol of the mobile coordinate frame, which is attached to the leader.

The position of each virtual member relative to the inertial coordinate frame can be calculated as follows.

(3) \begin{align}\left[ {\begin{array}{*{20}{c}}{\bar x_i^E}\\[5pt]{\bar y_i^E}\\[5pt]1\end{array}} \right] = \left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c}}{{\rm{cos}}\left( \theta \right)} & { - sin\left( \theta \right)} & {{x_l}}\\[5pt] {sin\left( \theta \right)} & {cos\left( \theta \right)} & {{y_l}}\\[5pt] 0 & 0 & 1\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{\bar x_i^L}\\[5pt] {\bar y_i^L}\\[5pt] 1\end{array}} \right]\end{align}

(4) \begin{align}\bar p_i^E = C_L^E{\rm{\;}}\bar p_i^L\end{align}

$\bar p_i^E = \left[ {\begin{array}{*{20}{c}}{\bar x_i^E}\\{\bar y_i^E}\end{array}} \right]$ : The position vector of the virtual member in the inertial coordinate frame. The index (E) is the symbol of the inertial coordinate frame.

$C_L^E$ : Homogeneous transformation matrix.

After obtaining the position of the virtual structure in the inertial coordinate frame, the real members are forced to be pushed towards the virtual members to pursue the virtual members by defining the potential functions between the real and virtual members. To this end, the distance error between the corresponding virtual and real members is defined as follows.

(5) \begin{align}{e_i} = {\bar p_i} - {p_i}\end{align}

${\bar p_i} = \left[ {\begin{array}{*{20}{c}}{\bar x_i^E}\\{\bar y_i^E}\end{array}} \right]$ and ${p_i} = \left[ {\begin{array}{*{20}{c}}{x_i^E}\\{y_i^E}\end{array}} \right]$ : The position vectors of virtual and real members in the inertial coordinate frame, respectively. For simplicity, the superscript E is omitted in the rest of the paper.

The general kinematic model of the members for an n-member group with relative positions of ${p_i}$ in the inertial coordinate frame is defined as follows:

(6) \begin{align}\left\{ {\begin{array}{*{20}{c}}{\;{{\dot p}_i} = {v_i}\;}\\[5pt] {{{\dot v}_i} = {u_i}}\end{array}\;} \right.\!\!\!\!,\;i = 1,2, \ldots ,n\end{align}

Where, ${u_i}$ is the total control input vector, which is a combination of the control inputs as is described in the continuation of Section 2 and Section 3.

2.2 Formation control

The first input is a function directing the real members of the group to the corresponding virtual members. This function is obtained from the gradient of the potential function, which is defined on the error variable ( ${e_i}$ ) as follows.

(7) \begin{align}{u_{1i}} = - \Delta {U_{1i}} - {W_{1vi}}\left( {{v_i} - {{\bar v}_i}} \right)\end{align}

${U_{1i}}$ : The potential function between real and virtual members.

The second term in Equation (7) is the damper term, which can be adjusted by the weighting coefficient ( ${W_{1vi}})$ . The terms ${v_i}$ and ${\bar v_i}$ are the velocity of real and virtual members, respectively.

In this paper, the gradient was obtained based on the distance between the group members because agents move in the direction of the negative gradient of the potential function.

The following potential function is applied to converge the group real members toward the virtual structure.

(8) \begin{align}{U_{1i}} = \frac{{{W_{1i}}}}{2}{\left\| {{{\bar p}_i} - {p_i}} \right\|^2}\end{align}

Where, ${W_{1i}}$ is the weighting coefficient of attractive force between real and virtual agents.

2.3 Collision avoidance

A repulsive potential function is defined to guarantee collision avoidance of the group members. The gradient of this function drastically increases by reducing the distance between the members. The function is merely activated in the neighborhood of each member in a certain radius. Since the dynamics of the group members are as dimensionless point masses, the dimensions of members should be considered in calculating the radius of repulsive potential activation. Thus, the desired radius is obtained by summing the diameters of the real members, the safety margin, and the distance required to create the formation as follows:

(9) \begin{align}{r_{ij}} = {d_r} + {d_s} + {d_{ij}}\end{align}

Where, ${d_r}$ is the largest dimension of real dynamics of the group members, ${d_s}$ denotes the safety margin, ${d_{ij}}$ shows the distance required between two members of i, j to keep the desirable position in the formation and ${r_{ij}}$ is the radius of repulsive force activation (Fig. 2).

Figure 2. The position of two members relative to each other.

Thus, the second input of the formation control system is the gradient of the repulsive potential function around each member, defined as follows.

(10) \begin{align}_{2i} = \left\{ {\begin{array}{*{20}{c@{\quad}c}}{\Delta {U_{2i}} + {F_{damp}}} & \|{{p_i} - {p_j}\| \le {r_{ij}}}\\[5pt] 0 & \|{{p_i} - {p_j}\| \gt {r_{ij}}}\end{array}} \right.\end{align}

where, ${U_2}$ is the repulsive potential function and ${F_{damp}}$ shows the damper of the relative movement between two members. ${p_i}\;$ and $\;\;{p_j}$ are position vectors of i and j agents.

(11) \begin{align}{F_{damp}} = \mathop \sum\nolimits_{j = 1,j \ne i}^n {W_{2vi}}\left( {{v_j} - {v_i}} \right)\end{align}

${W_{2vi}}$ is weighting coefficient. ${v_i}$ and ${v_j}$ are the velocity of agents i and j, respectively.

In this paper, a novel potential function is introduced, which efficiently creates repulsive forces between the group members. The proposed function (Equation (12)) has a proper efficiency in creating a repulsive force in avoiding obstacle collision. Due to its soft behaviour when interacting with other forces, chattering at the collision area is less. It is also not very sensitive to changes in weighting coefficients, so adjusting the weighting coefficients of this function is easy.

(12) \begin{align}{U_{2i}} = \mathop \sum \nolimits_{j = 1,j \ne i}^n {w_{2i}}(\;\left( \|{{p_i} - {p_j}\|_\sigma - dr\;} \right)csc{h^{ - 1}}\left( {\frac{{\|{p_i} - {p_j}\|_\sigma - dr\;}}{{{\alpha _1}}}} \right) + {\alpha _1}\;sin{h^{ - 1}}\left( {\frac{{\|{p_i} - {p_j}\|_\sigma - dr\;}}{{{\alpha _1}}}} \right)\end{align}

(13) \begin{align}\Delta {U_{2i}} = \mathop \sum \nolimits_{j = 1,j \ne i}^n {w_{2i}}\;\;csc{h^{ - 1}}\left( {\frac{{\|{p_i} - {p_j}\|_\sigma - dr\;}}{{{\alpha _1}}}} \right) \times \;\frac{{\left( {{p_i} - {p_j}} \right)\;}}{{\left( {1 + \varepsilon \|{p_i} - {p_j}\|_\sigma } \right)}}\end{align}

Where, $\;{w_{2i}}\;\;$ and ${\alpha _1}$ are weighting coefficient and a constant coefficient, respectively. The values of ${w_{2i}}\;\;$ and ${\alpha _1}$ are obtained by trial and error in the simulation process. Their values were given in Table 2. The function ( ${\rm{csch}}$ ) is the hyperbolic cosecant function.

Table 2 The simulation parameters

Note 1: $\|z\|_\sigma $ (the σ-norm of a vector) is a ${R^m} \to R \ge 0$ mapping, which is defined as follows [Reference Olfati-Saber15].

(14) \begin{align}\|z\|_\sigma = \frac{1}{\varepsilon }\left( {\sqrt {1 + \varepsilon {\|z\|^2}} - 1} \right),\varepsilon \gt 0\end{align}

$\varepsilon \;$ is the small positive number shown in Table 2.

Figure 3 displays the gradient behaviour of the potential function.

Figure 3. Repulsive function for ${d_r} = 10m$ , ${w_{2i}} = 1$ , $\alpha = 1000$ .

Thus, the formation control input is obtained regardless of the environmental disturbances as follows.

(15) \begin{align}{u_{fi}} = {u_{1i}} + {u_{2i}}\end{align}

The algorithm introduced up to this stage can be utilised in all robots, whose dynamics can be described as the model presented in Equation (6). Nevertheless, the physical and functional constraints of robots should be considered in practical applications.

3.1 Static obstacles

As illustrated in Fig. 5, to apply the obstacle avoidance algorithm in online flights, the obstacle avoidance mechanism is activated, when the UAV reaches the distance ${R_{am}}$ (Equation (27)) from the obstacle. A tangent line to the flight path is drawn from the time the UAV reaches the ${R_{am}}$ radius. This line is named ‘collision line’. If the collision line ( ${p_{tt}}$ ) intersects the obstacle circle (the circle that surrounds the obstacle with radius of ${R_o}$ and the centre point of O) at two points, the collision will definitely happen. The intersection points are called ${c_1}$ and ${c_2}$ .

Next, a line is sketched from point H, in the middle of points ${c_1}$ and ${c_2}$ , which passes through the centre of the obstacle circle (point O). This line is perpendicular to the collision line. Subsequently, a circle with a radius ${R_m}$ and centre ${O_{m1}}$ (relation (23)) is generated so that the circle centre is located on the line perpendicular to the collision line, which is also tangent to the obstacle circle according to Fig. 5. Continuing the process, two circles ${O_{m2}}$ and ${O_{m3}}$ , which are tangent to the collision line as well as to the circle ${O_{m1}}$ , are derived. Finally, three circular sectors from the three introduced circles of ${O_{m1}}$ , ${O_{m2}}$ , and ${O_{m3}}$ are selected.

These circular sectors are indeed the ideal and desirable flight path for the UAVs collision avoidance with the obstacle. This is because it is the closest path to the obstacle, which is compatible with the functional constraints of the UAVs.

The information required to estimate the positions of the mentioned circles is given as follows.

(22) \begin{align}\vec O\left( i \right) = \vec O - \vec H\left( i \right) \to \;\;\;{\vec e_{HO}}\left( i \right) = \frac{{\overrightarrow {HO} }}{\|{\overrightarrow {HO} }\|}\end{align}

(23) \begin{align}{\vec O_{m1}}\left( i \right) = \left( {{R_m} - {R_o}} \right)\;{\vec e_{HO}}\left( i \right) + \vec O\end{align}

(24) \begin{align}{\vec O_{m2}}\left( i \right) = - {R_m}\;{\vec e_{HO}}\left( i \right) + {\vec a_1}\left( i \right)\end{align}

(25) \begin{align}{\vec O_{m3}}\left( i \right) = - {R_m}\;{\vec e_{HO}}\left( i \right) + {\vec a_2}\left( i \right)\end{align}

Where, ${\vec a_1}$ & ${\vec a_2}$ represent the relative position vectors of the intersections of the collision line and a circle with radius ${R_a}$ and the centre of O (centre of the obstacle). $\vec O$ is the position vector of obstacle circle.

(26) \begin{align}{R_a}\left( i \right) = \sqrt {4{R_m}^2 - {{\left( {2{R_m} - Ro + \|\overrightarrow {HO}\| } \right)}^2} + {\|{\overrightarrow {HO} }\|^2}} \end{align}

${R_a}\left( i \right)$ is the instantaneous radius of the activation circle. The activation circle acts as a switch, which activates the OA mechanism when any UAV enters inside it.

(27) \begin{align}{R_{am}}\left( i \right) = \sqrt {\left( {4{R_o}{R_m} - {R_o}^2} \right)} \end{align}

${R_{am}}$ is the maximum value of ${R_a}$ , i.e., $\overrightarrow {HO} = 0$ . Meanwhile, the minimum distance requires sensing the obstacle for successful operation of AO mechanism. As shown in Equation (27), it is only dependent on ${R_o}$ and ${R_m}$ .

Three potential functions in the centre point of each circle $\;{\vec O_{m1}}$ , ${\vec O_{m2}}$ , and ${\vec O_{m3}}$ are defined for guiding the UAVs to follow the desired path. These functions are active only inside the circles ${O_1}$ , ${O_2}$ , and ${O_3}$ , respectively. The circles ${O_1}$ , ${O_2}$ , and ${O_3}$ are shown in Fig. 6 and Fig. 7 (not illustrated in Fig. 5 due to the large number of images). Note that the direction of the repulsive force in the circles ${\vec O_{m2}}$ and ${\vec O_{m3}}$ is inwards, while it is in the opposite direction in the circle ${O_{m1}}$ .

The repulsive potential functions of obstacle avoidance are defined as follows.

(28) $$\eqalign{ & p{O_{m1}} = {O_{m1}} - {p_i}{_\sigma } \cr & \matrix{ {{U_{o1i}} = \left\{ {\matrix{ {{w_{o1\;}}\left( {\;p{O_{m1}}\;csc{h^{ - 1}}\left( {{{p{O_{m1}}\;} \over {{\alpha _2}}}} \right) + {\alpha _2}\;sin{h^{ - 1}}\left( {{{p{O_{m1}}\;} \over {{\alpha _2}}}} \right)} \right)} & {p{O_{m1}} \le {R_m}} \cr 0 & {p{O_{m1}}{\rm{ < }}{R_m}} \cr } } \right.} \hfill \cr } \cr} $$

(29) \begin{align}{U_{o2i}} = \left\{ {\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{\dfrac{{{w_{o2}}}}{2}\;{{\left( {{p_i} - {O_{m2}}} \right)}^2}\;\;\;\;\;\;\;\;\;\;} \|{{p_i} - {O_{m2}}\| \ge {R_m}}\end{array}}\\[5pt] {\begin{array}{*{20}{c}}{0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \|{{p_i} - {O_{m2}}\| \lt {R_m}}\end{array}}\end{array}} \right.\end{align}

(30) \begin{align}{U_{o3i}} = \left\{ {\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{\dfrac{{{w_{o3}}}}{2}\;{{\left( {{p_i} - {O_{m3}}} \right)}^2}} {\;\;\;\;\;\;\;\;\;\;\;\|{p_i} - {O_{m3}}\| \ge {R_m}}\end{array}}\\[5pt]{\begin{array}{*{20}{c}}{0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\|{p_i} - {O_{m3}}\| \lt {R_m}}\end{array}}\end{array}} \right.\end{align}

The results of the obstacle avoidance forces, obtained from the gradient of the defined potential functions, can be calculated as follows.

(31) \begin{align}{u_{oi}} = \Delta {U_{o1i}} + \Delta {U_{o2i}} + \Delta {U_{o3i}}\end{align}

Figure 6 illustrates an instance of applying the introduced algorithm on a UAV.

Note 2: Since the presented algorithm relies on the collision line at every moment, these collision lines must be drawn continuously tangent to the ideal path (the path that does not sense the obstacle). Hence, an ideal path is required to be defined simultaneously with the real path. This ideal path will be the basis for designing the obstacle avoidance mechanism. This virtual ideal path is shown in the simulation section in green color.

3.2 Dynamic obstacles

For optimal obstacle avoidance of moving obstacles, whose positions change with time, appropriate predictions of their future positions should be done. Otherwise (i.e., when information is online and there is no information on the obstacle’s behaviour in future moments), a safety range should be considered to avoid obstacle collision in case of unexpected behaviours. As such, it is deviated from the ideal and optimal behaviour in dealing with a static obstacle. Nonetheless, it is worth noting that this behaviour should be at least as close as possible to the ideal behaviour.

When designing the dynamic obstacle avoidance strategy, in addition to considering UAV functional constraints, the state of the collision must also be predicted. Thus, as depicted in Fig. 8, the obstacle circle was enclosed inside an MTR circle to have the longest distance between the obstacle circle and the turn circle in line with the obstacle’s trajectory. In this way, the collision line hits the MTR circle before hitting the obstacle circle to activate the avoidance mechanism. The design steps of this mechanism are explained as follows.

The first step of avoidance mechanism is to determine the minimum distance from the obstacle, where the mechanism should be activated in order to implement the avoidance process based on the desired criteria (Fig. 8), and is calculated as follows.

(32) \begin{align}{\vec O_{m4}}\left( i \right) = \vec O + \left( {{R_m} - {R_o}} \right)\;{\vec e_{{v_o}}}{\rm{\;}}\end{align}

(33) \begin{align}{R_{dam}}\left( i \right) = {\left[ {{{\left( {sin\beta \;\left( {{R_m} - {R_o}} \right) + {R_m}} \right)}^2} + {{\left( {2{R_m} - cos\beta \;\left( {{R_m} - {R_o}} \right)} \right)}^2}} \right]^{0.5}}\end{align}

$\beta $ : The angle between the obstacle velocity vector and the UAV velocity vector.

${R_{dam}}\left( i \right)$ : The minimum distance that requires sensing the obstacle for successful operation of AO mechanism.

${\vec e_{{v_o}}}$ : The unit vector of the obstacle velocity vector.

Continuing the process, dynamic AO mechanism is presented in Fig. 9. Figure 5 is a special case of Fig. 9. In Fig. 5, because the obstacle is fixed, a line perpendicular to the velocity vector of the UAVs (instead of the obstacle velocity vector) is used to design the OA mechanism. Two potential functions (Equations (37) and (38)) are defined in the centre point of each circle $\;{\vec O_{m4}}$ and ${\vec O_{m5}}$ for guiding the UAVs to follow the desired path. The first function is active only inside the circle ${O_4}$ . The other required data are presented as follows:

(34) $$\eqalign{ & \overline {{a_2}{O_4}} \left( i \right) = \sin \beta \;{\overline {PO} _{m4}}\left( i \right) \cr & \overline {P{a_2}} \left( i \right) = {({\overline {PO} _{m4}}^2 - {\overline {{a_2}{O_{m4}}} ^2})^{0.5}} \cr & \overline {{a_1}{a_2}} \left( i \right) = {\left[ {4R_m^2 - {{\left( {{R_m} - {k_1}\;\overline {{a_2}{O_{m4}}} } \right)}^2}} \right]^{0.5}} \cr & \overline {{a_1}P} \left( i \right) = \overline {{a_1}{a_2}} \left( i \right) - {k_1}\;\overline {P{a_2}} \left( i \right) \cr & {{\vec a}_1}\left( i \right) = \vec P\left( i \right) - \overline {{a_1}P} \left( i \right)\;{{\vec e}_{{v_i}}}\;\left( i \right) \cr & \matrix{ {{R_{da\left( i \right)}} = \vec O - {{\vec a}_1}\left( i \right)} \hfill \cr } \cr} $$

${\vec e_{{v_i}}}$ : The unit vector of the UAV velocity vector.

${R_{da}}\left( i \right)$ is the instantaneous radius of the activation circle.

Note 3: It is notable that when $\overline {PO} (i) \gt {\overline {PO} _{m4}}(i)$ , ${k_1} = 1$ ; otherwise, ${k_1} = - 1$ .

The position of the MTR circle that surrounds the obstacle is as follows:

(35) \begin{align}{\vec O_{m5}}(i) = {\vec a_1} + {k_2}\left[ {\begin{array}{*{20}{c}}{{{\vec e}_{{v_i}}}(2)}\\{ - {{\vec e}_{{v_i}}}\left( 1 \right)}\end{array}} \right]{R_m}\end{align}

For optimisation of the OA mechanism behaviour (In order to minimise deviation from the ideal path), the instantaneous speed of the obstacle must also be available, till proportion to the relative speed changes (between UAV and obstacle), the radius of circle O_m4 and the weight coefficients are adjusted. Of course, to do this, it is necessary to add an algorithm to the current mechanism to apply the previous settings according to the working conditions.

There are times when the collision is unexpected, and the angle of impact is such that the UAV does not have enough time to bypass the obstacle from behind. This occurs when the collision line experiences the first encounter with the obstacle from the front and at a very close distance to the obstacle.

In other words, the following condition is met according to Fig. 10.

(36) \begin{align}\|p_i' - P(i)\| \le \|{a_1} - (i)\|,\,\overline{{a_2} {O_{m4}}} \gt \mu {R_m}\end{align}

Figure 10. Unexpected AO mechanism.

Where $p_i'$ is the i ^th ideal member of the group (Note 2) and $\mu $ denotes the impact sharpness coefficient, which is a function of the relative speed between the UAV and obstacle where $0 \lt \mu \lt 1$ .

The mechanism of obstacle avoidance for this condition is illustrated in Fig. 10. The governing relation in this part is similar to the previous section, except it is necessary to select ${k_1} = - 1$ and ${k_2} = - 1$ .

In this section, the repulsive potential functions of obstacle avoidance are considered as follows.

(37) \begin{align}{U_{o4i}} = \left\{ \begin{array}{c@{\quad}c} {{w_{o4}}\left( {p{O_{m4}}csc{h^{ - 1}}\left( \dfrac{p{O_{m4}}}{{\alpha _2}} \right) + {\alpha _2}sin{h^{ - 1}}\left( \dfrac{{{p{O_{m4}}_\sigma }}}{{{\alpha _2}}} \right)} \right)\,} & {p{O_{m4}} \le {R_m}} \\[12pt] {0} & {p{O_{m4}} \gt {R_m}} \end{array} \right.\end{align}

(38) \begin{align}{U_{o5i}} = \left\{ {\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{\frac{{{w_{o5}}}}{2}\;{{\left( {{p_i} - {O_{m5}}} \right)}^2}} {\;\;\;\;\;\;\;\;\;\|{p_i} - {O_{m5}}\| \ge {R_m}}\end{array}}\\[5pt] {\;\;\;\;\;\;\;\;\;\;\begin{array}{*{20}{c}}{0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \|{{p_i} - {O_{m5}}\| \lt {R_m}}\end{array}}\end{array}} \right.\end{align}

where ${w_{o4\;}}$ and ${w_{o5}}$ are weighting coefficient. ${\alpha _2}$ is a constant coefficient.

The results of the obstacle avoidance forces, obtained from the gradient of the defined potential functions, can be calculated as follows.

(39) \begin{align}{u_{oi}} = {\Delta _{o4i}} + \Delta {U_{o5i}}\end{align}

Ultimately, the output of the total algorithm can be summarised, i.e. the aggregation of formation and the obstacle avoidance algorithms as follows.

(40) \begin{align}{u_i} = {u_{fi}} + {u_{oi}}\end{align}

where ${u_i}$ is the total input to the group members’ dynamics, ${u_{fi}}$ denotes the formation control input regardless of environmental disturbances (Equation (15)) and ${u_{oi}}$ is the sum of the obstacle avoidance forces.

Figure 11. Encountering an obstacle from centre.

Figure 12. Encountering an obstacle from sides.

Figure 13. Encountering an obstacle with velocity of ${\left[ {0{\rm{\;}} - 15} \right]^T}$ .

Figure 14. Encountering an obstacle with velocity of ${\rm{\;}}{\left[ {7{\rm{\;}} - 15} \right]^T}$ .

Figure 15. Encountering an obstacle with velocity of ${\rm{\;}}{\left[ { - 7{\rm{\;}} - 15} \right]^T}$ .

Figure 16. Functional characteristic of UAVs in OA.

4.1 Static AO simulation

The simulation results are presented in the Matlab software. The functional constraints of Aerosonde [Reference Niculescu24], which is a UAV manufactured by AAI corporation of USA, are considered in the simulation process. The dynamics of the virtual leader is assumed as follows.

(41) $$\eqalign{ & {{\dot x}_l} = {v_l}\cos \left( {{\psi _l}} \right) \cr & {{\dot y}_l} = {v_l}\;{\rm{sin}}\left( {{\psi _l}} \right) \cr & \matrix{ {{{\dot \psi }_l} = {w_l}} \hfill \cr } \cr} $$

where ${p_l}\left( t \right) = \left[ {{x_l}\left( t \right),\;\;{y_l}\left( {} \right)} \right]$ : Indicating the position vector of the leader at time t in the inertial coordinate system.

${v_l}$ & ${\psi _l}$ : The leader’s linear velocity and orientation angle, respectively. The leader control input is assumed as follows.

\begin{align*}{v_l}\left( t \right) = 20\;m/s\end{align*}

\begin{align*}{w_l} = 0.07\;{\rm{sin}}\left( {0.2t} \right)\end{align*}

Note that the following functional constraints are always applied to UAVs in all stages of the simulation.

\begin{align*}12\;m/s \lt {v_i} \lt 38\;m/s\end{align*}

\begin{align*} - {30^{\circ}} \lt {\phi _i} \lt {30^{\circ}}\end{align*}

\begin{align*}0 \lt {T_i} \lt 200\end{align*}

The attractive forces between virtual and real members should be calculated based on Equation (7) as well as the repulsive forces around a member for anti-collision between agents from Equation (10) and repulsive forces between members of the group with the obstacle based on Equations (29), (30) and the parameters indicated in Table 2. Finally, the total input of the system is obtained by substitution in Equation (40).

The initial conditions for simulation began with random locations in the [300m × 300m] space and with the random angles.

In all simulation processes, merely agent No. 2 is illustrated due to the large number of images. Agent No. 2 has the highest deviation from the normal pass, and thus has the worst motion condition.

Figure 11 displays the simulation of the movement of three UAVs with a straight-line formation or V-shaped formation with an angle of ${180^{\circ}}$ . The simulations are performed at a constant altitude and the zero-path angle (γ = 0). It has been conducted from different angles and in the most critical states possible in colliding with obstacles for fully demonstrating the capability of the designed algorithm.

Based on the initial assumption of the problem (considering the obstacle within a circle), only two cases may happen with the collision of the UAV from any angle and direction. The first case involves a collision with the centre of the obstacle, which would be the worst possible situation as the highest deviation in the UAV trajectory can occur (Fig. 11). In the second case, the collision may occur from lateral points outside the centre of the circle. In such a case, the path deviation rate depends on the distance of the collision point from the centre of the circle (Fig. 12).

The trajectories of the UAVs with the straight-line formation in facing an obstacle with a diameter of 60m are displayed in Fig. 11a. As shown, the UAVs have maintained the desired formation while following the leader. The members of the group return to the desired geometric shape with the minimum disorganisation after colliding with the obstacle. This feature of the algorithm suggests the robustness of the formation structure against environmental conditions. The magnification of the obstacle collision location is shown in Fig. 11b. As shown, the UAVs change direction very gently with a smooth curvature to pass the obstacle, indicating the compatibility of the designed path with the flight capabilities of UAVs. According to Fig. 11c, the UAVs reach from random positions to the desired distance of the adjacent UAV (radius of repulsive force activation in formation keeping algorithm: ${r_{ij}} = 50m$ ) after about 20 seconds and well maintain the distance over time, while following the leader. These distances are temporarily interrupted in the face of obstacles; however, they immediately return to the desired state. In the meantime, the distance between the UAVs never approaches the critical distance ( ${d_r} = 10m$ ) due to using the unbounded (from left) inter-agent repulsive function. This feature highly matters in practical applications to observe safety and prevent the collision of UAVs with each other.

Another scenario is displayed in Fig. 12, where the obstacle is in the middle of the paths of two adjacent UAVs. Accordingly, the obstacle is well avoided and crossed in this case as well with no collision and turbulence on the flight path. As indicated, state variables and control inputs have also shown proper behaviour.

Figure 12a and 12b indicate the trajectories of UAVs and the magnification of the trajectories, respectively. As shown, smooth passes are also generated in this case. As shown in Fig. 12b and 12c, in the time interval of 30 to 40 seconds, the formation composition is disturbed due to facing the obstacle. The distance between the UAVs has been reduced below 50 metres. That is, the repulsive force is activated.

Figure 12c demonstrates the distance between UAV2 and UAV3, which changes up to 120m around the obstacle but the distance between UAV1 and UAV2 does not reach below the critical distance ( ${d_r} = 10m$ ). Figure 12d and 12e illustrate the magnitude of state variables of velocities and heading angle. Figure 12f, g and h depict the time history of control inputs of thrust force, bank angle, and lift force, respectively.

The scenarios presented are indeed the most complex possible states that may occur in the case of obstacles, especially when the UAVs themselves are turning. Accordingly, other conditions of collision will be available.

4.2 Dynamic AO simulation

Similar to the previous section, the formation flight of UAVs was simulated with the same dynamic in this part, where the leader’s movement occurs with the fixed angle of ${\psi _l} = 0.5\;\left( {rad} \right)$ .

Figure 13 reveals the trajectories of three UAVs with the straight-line formation when encountering an obstacle at a ${\left[ {0{\rm{\;}} - 15} \right]^T}$ speed. As observed, the group of UAVs accomplished its mission, while following the leader, keeping the formation with the least disorganisation, when encountering a dynamic obstacle.

Figure 13c presents the distances between the UAVs and the centre of the obstacle. As illustrated, the distance of none of the UAVs with the obstacle reached below 60 m (radius of the obstacle circle), meaning that they did not collide with the obstacle.

Figure 13d shows the distance between UAVs. As observed, while the distance between UAVs 2 and 3 is about 150 m when encountering obstacles, the distance between UAVs 1 and 2 is by no means less than the critical distance (10 m).

Figures 14 and 15 show the obstacle movements at speeds of ${\left[ {7{\rm{\;}} - 15} \right]^T}$ and ${\left[ { - 7{\rm{\;}} - 15} \right]^T}$ , respectively. Formation kipping and lack of colliding of members with the obstacle and each other are evident in the trajectories of UAVs, similar to the previous scenario.

The behaviour of state variables and control inputs is also shown in Figure 16 to display the functional characteristics of UAVs in encountering a dynamic obstacle. For the sake of brevity, we only mentioned the mode, where the obstacle moved at a speed of ${\left[ { - 7{\rm{\;}} - 15} \right]^T}$ . Notably, the performance characteristics of UAVs in other scenarios are similar.

Figure 16f, g and h show the control inputs of thrust force, bank angle, and lift force, respectively. As shown, all of the inputs are saturated based on dictated limitations at the beginning of the simulation section. The key point is the moments of collision with an obstacle (i.e. the interval of 22s to 32s), where the control inputs act with maximum power. This feature is a declaration of accomplishing a coordinated turn by minimum radius, which is in compliance with the principle of flight dynamics reported in basic literature [Reference Roskam25].