## 3 Formation Methods and Simulation
### 11 \(\times\) 11 Matlab Simulations of UAVs Magic Square Formation  
#### 3.3.1 UAV Model  
In real UAVs with different model parameters, there are multiple aerodynamic configurations, causing the variance of mathematical modeling. To simplify the algorithm of upper control, we suppose that the UAV internal-loop is controlled by autopilot. Thus, the model could be built with the UAV position and velocity external-loop model as the upper control algorithm, as shown in Figure 17.  
Figure 16: Reference for the UAV at square 1. Figure 17: UAV model. Figure 15: A visual reference topological structure diagram of a 7 \(\times\) 7 nested magic square UAV formation. The mathematical model of the \(UAV_{i}\) is expressed as Equation (8)  
\[ \begin{cases}\dot{x_{i}}=V_{i}\cos\gamma_{i}\cos\chi_{i}\\ \dot{y_{i}}=V_{i}\cos\gamma_{i}\sin\chi_{i}\\ \dot{h_{i}}=V\sin\gamma_{i}\\ \dot{V_{i}}=\dfrac{T_{i}-D_{i}}{m_{i}}-g\sin\gamma_{i}\\ \dot{\gamma_{i}}=\dfrac{L\cos\Phi_{i}-m_{i}g\cos\gamma_{i}}{m_{i}V_{i}}\\ \dot{x_{i}}=\dfrac{L_{i}\sin\Phi_{i}}{m_{i}V_{i}\cos\gamma_{i}}\end{cases}  \tag{8}\]  
where \(i=1,\cdots,N\). \(x_{i},y_{i}\), and \(h_{i}\) correspond to the down-range, cross-range, and altitude displacement. \(V_{i}\) refers to the airspeed of \(UAV_{i}\), \(\gamma_{i}\) is the plane path angle, and \(\chi_{i}\) represents the heading angle. \(T_{i}\) is the engine thrust, \(D_{i}\) refers to drag, \(m_{i}\) is the quality of \(UAV_{i}\), and \(g\) represents the gravity acceleration. Furthermore, \(L_{i}\) refers to lift, and \(\Phi_{i}\) is the bank angle. Equation (9) can be achieved with the transformation of the mathematical model. \[ \begin{cases}\dot{x_{i}}=u_{xi}\\ \dot{y_{i}}=u_{yi}\\ \ddot{z_{i}}=u_{hi}\end{cases}  \tag{9}\]  
\(u_{xi}\), \(u_{yi}\), and \(u_{hi}\) are the subjunctive control input, and the transformation relationship between the executive order and subjunctive control input can be expressed as Equation (10),  
\[ \begin{cases}\Phi_{i}=\arctan\left(\dfrac{u_{yi}\cos\chi_{i}-u_{xi}\sin\chi_{i }}{(u_{hi}+g)\cos\gamma_{i}-(u_{xi}\cos\chi_{i}+u_{yi}\sin\chi_{i})\sin\gamma_ {i}}\right)\\ L_{i}=m_{i}\dfrac{(u_{hi}+g)\cos\gamma_{i}-(u_{xi}\cos\chi_{i}+u_{yi}\sin\chi_ {i})\sin\gamma_{i}}{\cos\Phi_{i}}\\ T_{i}=m_{i}[(u_{hi}+g)\sin\gamma_{i}+(u_{xi}\cos\chi_{i}+u_{yi}\sin\chi_{i}) \cos\gamma_{i}]+D_{i}\end{cases}  \tag{10}\]  
where \(\tan(\chi_{i})=\dot{y_{i}}/\dot{x_{i}}\), and \(\sin(\gamma_{i})=\dot{h_{i}}/V_{i}\). Therefore, the subjunctive control input is designed as Equation (9), and the real input of the UVA could be calculated through Equation (10), which can be expressed as the state place:  
\[ \begin{cases}\dot{z_{i}}=Az_{i}+Bu_{i}\\ p_{i}=C_{p}z_{i}\\ v_{i}=C_{p}z_{i}\end{cases}  \tag{11}\]  
where \(z_{i}=[p_{i}T,v_{i}]^{T}\), \(p_{i}\) refers to the position vector, \(v_{i}\) is the speed vector, and \(u_{i}=[u_{xi}^{T},u_{yi}^{T},u_{hi}^{T}]^{T}\) shows the subjunctive control input. \[A_{i}= \begin{bmatrix}0&1\\ 0&0\end{bmatrix} \otimes I_{3},B_{i}= \begin{bmatrix}0\\ 1\end{bmatrix} \otimes I_{3},C_{p}= \begin{bmatrix}1&0\end{bmatrix} \otimes I_{3},C_{v}= \begin{bmatrix}0&1\end{bmatrix} \otimes I_{3} \tag{12}\]  
\(I_{3}\in R^{3\times 3}\) refers to the identity matrix, and \(\otimes\) is the Kronecker product. In Equation (10), the air resistance \(D_{i}\) can be expressed as Equation (13). \[D_{i}=0.5\rho(V_{i}-V_{wi})^{2}SC_{D0}+\dfrac{2k_{d}k_{n}^{2}L^{2}}{\rho(V_{i}- V_{wi})^{2}Sg^{2}} \tag{13}\]  
where \(\rho\) refers to the air density, \(C_{D0}\) represents the zero-lift drag coefficient, \(V_{wi}\) refers to gust, \(S\) is the wing area, \(k_{d}\) is the induced drag, and \(k_{n}\) refers to the load-factor effectiveness. The mathematical modeling of gust can be expressed as Equation (14). \[ \begin{cases}V_{wi}=\overline{V}_{wi}+\delta V_{wi}\\ \overline{V}_{wi}=0.215V_{m}\log_{10}(h_{i})+0.285V_{m}\end{cases}  \tag{14}\]  
where \(\overline{V}_{wi}\) is normal wind shear, \(V_{m}\) refers to the mean wind speed and \(\delta V_{wi}\) is the wind gust turbulence.