2.1 UAV model

In this paper, N UAVs in a two-dimensional space are taken into consideration. The dynamic model of each UAV is described as follows [Reference Qiu and Duan14]:

(1) \begin{equation}\begin{cases}\dot{P}_x^i=V_{xy}^i\cos\phi^i,\\[5pt] \dot{P}_y^i=V_{xy}^i\sin\phi^i,\\[5pt] \dot{V}_{xy}^i=u_x^i\cos\phi^i+u_y^i\sin\phi^i,\\[5pt] \dot{\phi}^i=u_y^i\cos\phi^i-u_x^i\sin\phi^i,\end{cases}\end{equation}

where, for UAV i, $P^i=(P_x^i, P_y^i)$ is the position vector, $V^i=(V_x^i, V_y^i)$ is the velocity vector, $V_{xy}^i$ is the horizontal airspeed, $\phi^i$ represents yaw angle, and $u^i = (u_x^i, u_y^i)$ is the control input vector. In addition, the following constraints must be met:

(2) \begin{equation}\begin{cases}V_{xy}^{\min} \leq V_{xy}^i \leq V_{xy}^{\max},\\[5pt] \left| \dot{\phi}^i \right| \leq \frac{n^{\max}g}{V_{xy}^i},\\[5pt] V_{xy}^i=V^e_{xy}, & \text{if $\left|V_{xy}^i-V^e_{xy}\right|<V_{xy}^{l},$}\\[5pt] \phi^i=\phi^e, & \text{if $\left|\phi^i-\phi^e\right|<\phi^{l},$}\end{cases}\end{equation}

where, $n^{\max}=10$ g is the maximum lateral overload, $g=10\mathrm{m/s}^2$ is the gravitational acceleration, $V_{xy}^{\min}=5\mathrm{m/s}$ , $V_{xy}^{\max}=15\mathrm{m/s}$ are the lower and upper limits of horizontal airspeed, respectively, $V_{xy}^{l}=0.25\mathrm{m/s}$ , $\phi^{l}=0.1$ rad are allowable control error of horizontal airspeed, yaw angle, respectively, and $\phi^e = a\tan \left( \frac{V_y^e}{V_x^e} \right)$ , $V^e = (V_x^e,V_y^e)=(10,0)\mathrm{m/s}$ , $V^e_{xy}$ are the desired yaw angle, desired velocity vector, desired horizontal airspeed, respectively.

2.2 Flocking model

To achieve velocity controllable UAV flocking, the following movement rules need to be guaranteed: (1) Each UAV stay close to its neighbours; (2) All UAVs fly at a common velocity; (3) UAVs avoid conflict with each other; (4) The flocking velocity of UAVs is controllable.

Therefore, control input of UAV i consists of four components: $u_f^i$ is flocking geometry control component, which regulates the distance between UAV i and its neighbour. $u_{av}^i$ is horizontal airspeed alignment control component, which regulates the velocity of UAV i to be consistent with its neighbour. $u_c^i$ is collision avoidance control component, which regulates the safe distance among UAVs. $u_{vf}^i$ is flocking velocity control component, which regulates the velocity of UAV i to be consistent with the desired flocking velocity. Details are as follows [Reference Qiu and Duan14]:

(3) \begin{equation}u_k^i=u_f^i+u_{av}^i+u_c^i+u_{vf}^i-V_k^i, \quad k=x,y,\end{equation}

where, $k=x, y$ denotes the corresponding components in the x, y directions.

The flocking geometry control component $u_f^i$ is as the following equation:

(4) \begin{equation}u_f^i=C_f \sum_{j \in \left\{ d^{ij} \leq D_c \right\}} W_j^i(P_k^j-P_k^i) ln\! \left( \frac{d^{ij}}{D_d} \right)\end{equation}

where, $C_f=0.1$ denotes strength coefficient of the flocking geometry control component, $W_j^i$ denotes influence weight of UAV j to UAV i for flocking geometry control within [0,1], $d^{ij}=\sqrt{(P_x^i-P_x^j)^2 + (P_y^i-P_y^j)^2}$ denotes the horizontal distance between UAV i and j, $D_c=20$ m denotes the maximum horizontal communication distance, $D_d=10$ m denotes the desired horizontal distance between UAVs.

The horizontal airspeed alignment control component $u_{av}^i$ is as the following equation:

(5) \begin{equation}u_{av}^i=C_{av} \sum_{j \in \left\{ d_{ij} \leq D_c \right\}} W_j^i(V_k^j-V_k^i)\end{equation}

where, $C_{av}=0.1$ denote strength coefficient of the horizontal airspeed alignment control component.

The collision avoidance control component $u_c^i$ is as the following equation:

(6) \begin{equation}u_c^i=C_c \sum_{j \in \left\{ d^{ij} \leq D_{l1} \right\}} \left( \frac{1}{d^{ij}}-\frac{1}{D_{l1}} \right)^2 \frac{P_k^i-P_k^j}{d^{ij}}\end{equation}

where, $C_c=10^5$ denote strength coefficient of the collision avoidance control component, $D_{l1}=2$ m denotes the minimum distance between UAVs to avoid collision.

The flocking velocity control component $u_{vf}^i$ is as the following equation:

(7) \begin{equation}u_{vf}^i=C_{vf}W_i^iefv^i_k\end{equation}

where, $C_{vf}=1$ denotes strength coefficient of the flocking velocity control component, $W_i^i$ denotes influence weight of UAV i to UAV i for flocking velocity control within $[1,1.1]$ , $efv^i=(efv^i_x,efv^i_y)=(V^{e}_{xy} \cos\delta^i, V^{e}_{xy} \sin\delta^i)$ denotes the desired flocking velocity vector of UAV i, and $\delta^i$ denotes the desired flocking yaw angle of UAV i.

2.3 Obstacle avoidance and velocity tracking model

Because the magnitude of the desired velocity remains constant throughout the flight, obstacle avoidance and velocity tracking can be achieved by setting an appropriate desired flocking yaw angle $\delta^i$ which is calculated as follows:

(8) \begin{equation}\delta^{i}=a\tan \left( \frac{P^{core}_y + C_{\delta}\sin(\delta^{core}) - P^i_y}{P^{core}_x + C_{\delta}\cos(\delta^{core}) - P^i_x} \right)\end{equation}

(9) \begin{equation}P^{core}_k=\frac{1}{N}\sum_{i=1}^{N} P_k^i, \quad k=x,y,\end{equation}

where $P^{core}=(P^{core}_x,P^{core}_y)$ is position vector of flocking core calculated by Equation (9), $C_{\delta}$ is the parameter that controls the geometry of flocking, $\delta^{core}$ is the desired flocking yaw angle of flocking core, which is directly set as the reference desired flocking yaw angle in the velocity tracking process, but its equation is as follows in the obstacle avoidance process:

(10) \begin{equation}\delta^{core}=\begin{cases}a\tan \left( \dfrac{P^A_y - P^{core}_y}{P^A_x - P^{core}_x} \right), & \text{if $S_o \ne \emptyset,$}\\[9pt] a\tan \left( \dfrac{V_y^e}{V_x^e} \right), & \text{else,}\end{cases}\end{equation}

where $S_o$ is the set of obstacles within maximum perceived distance ahead of the UAV flocking in the desired velocity direction, as shown in Fig. 1. A is the point closest to the flocking core among the two edge points projected in the vertical direction of $V^e$ by the envelope of the obstacle closest to the flocking core in $S_o$ , $D_{l2}=10$ m is the minimum distance allowed between the UAV and obstacle, $R_{formation}$ is radius of the UAV flocking in the desired flocking yaw angle direction. Note that although the UAV flocking bypass the current obstacle, it cannot bypass the next obstacle until it passes through the current obstacle. The maximum perceived distance of UAV is 105.

Figure 1. Obstacle avoidance model.

4.1 UAV flocking obstacle avoidance and velocity tracking based on MOPSO

UAV flocking obstacle avoidance and velocity tracking based on MOPSO is a special case of UAV flocking control framework (see Algorithm 1). It uses MOPSO [Reference Coello Coello and Lechuga32] to calculate the optimal influence weight vector $\hat{W}^i$ , The specific steps of the MOPSO are as follows:

Input: Objective function $obj^i_1$ and $obj^i_2$ , maximum generation $M_G$ , number of particles $N_P$ , and the upper and lower bounds of $L_g^h$ and $S_g^h$ are the same as those of influence weight vector. Where, $L_g^h$ and $S_g^h$ are position vector and velocity vector of the generation g of particle h respectively, $h=1,2,\dots, N_P$ , $g=1,2,\dots, M_G$ , $L_g^h$ is a feasible influence weight vector.

Initialisation: Initialise the $L_g^h$ and $S_g^h$ of $N_P$ particles within their upper and lower bounds respectively, current iteration $g = 1$ .

Step 1: Calculate the objective function $obj^i_1$ and $obj^i_2$ of all the particles by Equations (1)–(8) and (12)–(13).

Step 2: Rank the particles positions in the set $L_g$ by Pareto sorting scheme, store the positions of particles in the Pareto frontier in the repository REP.

Step 3: $pbest^h = L_1^h$ .

Step 4: gbest is selected randomly from REP.

Step 5: Update the velocity of each particle by the following equation.

(14) \begin{equation}S_{g+1}^{h}=aS_{g}^{h}+r_{1}^T(pbest^h-L_{g}^h)+r_{2}^T(gbest-L_{g}^h)\end{equation}

Where, $a = 0.4$ is the inertia weight, $r_{1}$ and $r_{2}$ are random number vectors uniformly distributed in the interval [0,1].

Step 6: Update the position of each particle by the following equation.

(15) \begin{equation}L_{g+1}^{h}=L_{g}^h+V_{g+1}^{h}\end{equation}

Step 7: Initialises the state of the particle if its position exceeds the upper and lower bounds.

Step 8: Calculate the objective function $obj^i_1$ and $obj^i_2$ of all the particles by Equations (1)–(8) and (12)–(13).

Step 9: If $L_{g+1}^{h}$ is better than $pbest^h$ , then $pbest^h = L_{g+1}^{h}$ .

Step 10: Rank the particles positions in the set $L_{g+1}$ and REP by Pareto sorting scheme, store the positions of particles in the Pareto frontier in the repository REP.

Step 11: The current iteration $g=g+1$ . If $g \le M_G$ , then return to Step 4. Otherwise, output the optimal influence weight vector which is the position of particle in the REP whose objective function $obj^i_2$ has the smallest value.

4.2 UAV flocking obstacle avoidance and velocity tracking based on NN

UAV flocking obstacle avoidance and velocity tracking based on NN is a special case of UAV flocking control framework(see Algorithm 1). It uses NN to calculate the optimal influence weight vector $\hat{W}^i$ , and the NN is generated as follows.

Generate training samples: The UAV flocking velocity tracking algorithm based on MOPSO is repeated $N_{repeat}=100$ times to obtain the training sample set $\{input^i,label^i\}_{i=1}^{N_{repeat} \times \frac{T_{\max}}{ts}}$ . Where, $input^i$ is the input vector i and expressed by Equation (17), $label^i$ is the label vector i and expressed by Equation (18), $P_x^i$ , $P_y^i$ are random numbers within [0,30], [75,150] respectively, but horizontal distance between any two UAVs should be less than maximum horizontal communication distance of $D_c$ , number of UAVs $N=5$ , maximum simulation time $T_{\max}=50$ s, sampling time ts=0.5s, $C_{\delta} = \infty$ , $\delta^{core}$ is expressed by Equation (16).

(16) \begin{equation}\delta^{core}=\begin{cases}\dfrac{\pi}{3} \sin\left( \dfrac{2 \pi t}{0.8T_{\max}} \right), & \text{if $t \le 0.8T_{\max},$}\\[12pt] a\tan \left( \dfrac{V^e_y}{V^e_x} \right), & \text{if $t > 0.8T_{\max},$}\end{cases}\end{equation}

(17) \begin{equation}\begin{aligned}input^i &=\bigg( \underbrace{P_x^1,P_y^1,V^1_x,V^1_y,efv^1_x,efv^1_y}_{State \, of \, UAV \, 1},\\[5pt]&\mathrel{\phantom{=}} \phantom{{}\bigg( } \underbrace{P_x^2,P_y^2,V^2_x,V^2_y,efv^2_x,efv^2_y}_{State \, of \, UAV \, 2},\\[5pt]&\mathrel{\phantom{=}} \phantom{{}\bigg( } \dots, \\[5pt]&\mathrel{\phantom{=}} \phantom{{}\bigg( } \underbrace{P_x^N,P_y^N,V^N_x,V^N_y,efv^N_x,efv^N_y}_{State \, of \, UAV \, N} \bigg)\end{aligned}\end{equation}

(18) \begin{equation}\begin{aligned}label^i &=\bigg( \underbrace{\hat{W}^1_1,\hat{W}^1_2,\dots,\hat{W}^1_N}_{Optimal \, influence \, weight \, of \, UAV \, 1},\\[5pt]&\mathrel{\phantom{=}} \phantom{{}\bigg( } \underbrace{\hat{W}^2_1,\hat{W}^2_2,\dots,\hat{W}^2_N}_{Optimal \, influence \, weight \, of \, UAV \, 2},\\[5pt]&\mathrel{\phantom{=}} \phantom{{}\bigg(Optima} \dots, \\[5pt]&\mathrel{\phantom{=}} \phantom{{}\bigg( } \underbrace{\hat{W}^N_1,\hat{W}^N_2,\dots,\hat{W}^N_N}_{Optimal \, influence \, weight \, of \, UAV \, N} \bigg)\end{aligned}\end{equation}

Neural network structure: NN has four hidden layers, each layer has 6 $\times$ N neurons, and the dimensions of the input and output layers are 6 $\,\times\,$ N and N $\,\times\,$ N, respectively. Its structure is shown in Fig. 2. Where, the activation function of hidden layer, output layer are shown in Equations (19), (20), respectively.

(19) \begin{equation}Tansig(o) = \frac{2}{1+exp(-2o)}-1\end{equation}

(20) \begin{equation}Purelin(o) = o\end{equation}

Training neural network: Each dimension of the input data is normalised to $[-1, 1]$ before training, the loss function is the mean square error function is expressed by Equation (21), the training function is trainscg (scaled conjugate gradient backpropagation), epochs=1,000, and The ratio between training set, validation set and test set is $7:1.5:1.5$ .

(21) \begin{equation}mse = \frac{1}{N_d}\sum_{k=1}^{N_d}(label_k - output_k)^2\end{equation}

Figure 2. Neural network structure.

4.3 Analysis of generalisation capability of UAV flocking control based on NN

We discuss the generalisation capability of UAV flocking control based on NN from the following two aspects:

Position of UAVs: Deviation occurs in the calculation result when UAV position exceeds min-max normalisation range of the NN input vector. This problem can be solved by moving all UAVs in the flocking to min-max normalisation range without changing their relative positions.

Number of UAVs: The input vector of NN only supports the state of k UAVs, and can not be directly applied to larger scale UAV flocking control. This problem can be solved by distributed computing, that is, each UAV only calculates its next state according to the nearest k-1 UAV and its own state.

5.1 Simulation of UAV flocking obstacle avoidance and velocity tracking

In the simulation, five UAVs form a flocking, they first spend 50s crossing obstacles, then track the velocity for 50s. Initial parameters of UAVs and obstacles are shown in Table 1, and other initial parameters are as follows: maximum generation MG = 58, number of particles NP = 20, maximum simulation time $T_{\max}=100$ s, sampling time ts=0.5s, $C_{\delta} = \infty$ , and calculate $\delta^{core}$ by Equation (22).

(22) \begin{equation}\delta^{core}=\begin{cases}-\dfrac{\pi}{3} \sin\left( \dfrac{2 \pi (t-50)}{30} \right), & \text{if $t > 50 \, {\rm and} \, t \le 80,$}\\[12pt] a\tan \left( \dfrac{V^e_y}{V^e_x} \right), & \text{if $t > 80,$}\end{cases}\end{equation}

Figure 3 depicts the simulation results of five UAVs flocking obstacle avoidance and velocity tracking based on NN where, $RD_{ij}$ denotes the relative distance between UAV i and UAV j, $D_{l1}$ is the minimum distance to avoid collision between UAVs. As shown in Fig. 3(a), five UAVs successfully completed obstacle avoidance and velocity tracking tasks and flew out of a smooth trajectory curve. As shown in Fig. 3(b), the difference of horizontal airspeed among UAVs was less than 1.5480m/s, and the horizontal airspeed of the UAV flocking quickly converges to the desired horizontal airspeed after fluctuating within the allowable range. As shown in Fig. 3(c), the UAV flocking has a good velocity tracking effect, and a slight yaw angle difference. As shown in Fig. 3(d), relative distance between UAVs is converges rapidly and fluctuates less and always greater than the minimum distance to avoid collision between UAVs. As shown in Fig. 3(e) and (f), the geometry and velocity of the UAV flocking converge rapidly and stably.

Table 1. Initial parameters of the UAVs and obstacles

Figure 3. Simulation results of 5 UAV flocking obstacle avoidance and velocity tracking based on NN: (a) trajectory curves, (b) horizontal airspeed curves, (c) yaw angle curves, (d) relative distance curves, (e) $\sum obj_1^i$ curves, and (f) $\sum obj_2^i$ curves.

To illustrate the superiority of the UAV flocking obstacle avoidance and velocity tracking based on NN, we compare it with the UAV flocking obstacle avoidance and velocity tracking based on the MOPSO. The statistical results of the two algorithms are shown in Table 2. From the data in the table, the NN-based algorithm achieves faster convergence of horizontal airspeed and yaw angle, especially the simulation time has been greatly reduced. The simulation results of the UAV flocking obstacle avoidance and velocity tracking based on the MOPSO is shown in Fig. 4. From the comparison between Figs 3 and 4, it can be seen more intuitively that the experimental results based on NN have smaller fluctuation, faster convergence, and better obstacle avoidance and velocity tracking effects.

Table 2. Statistical results of two algorithms for 5 UAVs

^a Convergence time of horizontal airspeed.

^b Convergence time of yaw angle.

^c Total simulation time.

Figure 4. Simulation results of 5 UAV flocking obstacle avoidance and velocity tracking based on MOPSO: (a) trajectory curves, (b) horizontal airspeed curves, (c) yaw angle curves, (d) relative distance curves, (e) $\sum obj_1^i$ curves, and (f) $\sum obj_2^i$ curves.

5.2 Simulation of generalisation capability of UAV flocking control based on NN

In the simulation, we adopt the neural network for 5 UAVs in the previous experiment to verify its generalization capability in a flocking of 49 UAVs, which first track the velocity for 60s, then spend 60s crossing obstacles, and finally track the velocity for 60s. The initial velocity of the UAV and radius of obstacle is the same as that shown in Table 1, and the positions of UAVs and obstacles have been moved away from where they were during training NN (due to the large amount of data, no specific data is given here). Other initial parameters are as follows: maximum simulation time $T_{\max}=180$ s, sampling time ts=0.5s, calculate $C_\delta$ by Equation (23), calculate $\delta^{core}_1$ for the first velocity tracking by Equation (24), calculate $\delta^{core}_2$ for the second velocity tracking by Equation (25).

(23) \begin{equation}C_\delta=\begin{cases}80, & \text{if $S_o \ne \emptyset,$}\\[5pt] \infty, & \text{else,}\end{cases}\end{equation}

(24) \begin{equation}\delta^{core}_1=\begin{cases}\dfrac{\pi}{3} \sin\left( \dfrac{2 \pi t}{40} \right), & \text{if $t \le 40,$}\\[9pt] a\tan \left( \dfrac{V^e_y}{V^e_x} \right), & {\text{if} \, t > 40 \,\,\, \text{and} \,\,\, t < 50,}\end{cases}\end{equation}

(25) \begin{equation}\delta^{core}_2=\begin{cases}-\dfrac{\pi}{3} \sin\left( \dfrac{2 \pi (t-120)}{40} \right), & {\text{if} \, t > 120 \,\,\, \text{and} \,\,\, t \le 160,}\\[12pt] a\tan \left( \dfrac{V^e_y}{V^e_x} \right), & \text{if $t > 160,$}\end{cases}\end{equation}

Figure 5. Simulation results of 49 UAV flocking obstacle avoidance and velocity tracking based on NN: (a) trajectory curves, (b) horizontal airspeed curves, (c) yaw angle curves, (d) relative distance curves, (e) $\sum obj_1^i$ curves, and (f) $\sum obj_2^i$ curves.

Figure 5 depicts the simulation results of 49 UAV flocking obstacle avoidance and velocity tracking based on NN. where, $RD_{\min}$ is the minimum relative distance between UAVs. As shown in Fig. 5(a), the UAV flocking can safely pass through the obstacle area through geometric transformation, and the UAV position exceeds min-max normalisation range of the NN input vector, but it has no adverse impact on flight. As shown in Fig. 5(b), the horizontal airspeed of the UAV flocking quickly converges to the desired horizontal airspeed after fluctuating within the allowable range. As shown in Fig. 5(c), the UAV flocking has a good velocity tracking effect, and the yaw angle quickly converges to the desired yaw angle after fluctuations. As shown in Fig. 5(d), during the whole flight, the minimum relative distance between UAVs is always greater than the minimum distance to avoid collision between UAVs. As shown in Fig. 5(e) and (f), the geometry and velocity of the UAV flocking converge rapidly and stably. To sum up, UAV flocking controller based on NN has good generalisation capability and can control the flocking geometry to a certain extent.

To illustrate the superiority of the UAV flocking obstacle avoidance and velocity tracking based on NN, we compare it with the UAV flocking obstacle avoidance and velocity tracking based on the MOPSO. The statistical results of the two algorithms are shown in Table 3. From the data in the table, the NN-based algorithm achieves faster convergence of horizontal airspeed, especially the simulation time has been greatly reduced. The simulation results of the 49 UAV flocking obstacle avoidance and velocity tracking based on the MOPSO is shown in Fig. 6. From the comparison between Figs 5(a) and 6(a), it can be seen that in the NN-based experimental results the UAV is further away from the obstacle and safer. The comparison of Figs 5(b) and 6(b) shows that the horizontal airspeed deviates from the desired value in a shorter time, smaller amplitude, more regular, and achieves convergence earlier in the NN-based experimental results. Contrast Fig. 5(c) with Fig. 6(c) indicates that the NN-based method has a fast response to changes in yaw angle. In conclusion, the experimental results based on NN have smaller fluctuation, faster convergence, and better obstacle avoidance and velocity tracking effects.

Table 3. Statistical results of two algorithms for 49 UAVs

^a Convergence time of horizontal airspeed.

^b Convergence time of yaw angle.

^c Total simulation time.

Figure 6. Simulation results of 49 UAV flocking obstacle avoidance and velocity tracking based on MOPSO: (a) trajectory curves, (b) horizontal airspeed curves, (c) yaw angle curves, (d) relative distance curves, (e) $\sum obj_1^i$ curves, and (f) $\sum obj_2^i$ curves.