Robust nonlinear control of close formation flight

TL;DR

This paper presents a robust nonlinear control method for close formation flight of aircraft, ensuring precise position tracking despite aerodynamic disturbances and uncertainties, with validated numerical simulations.

Contribution

It introduces a novel control algorithm combining baseline control laws and disturbance observers for robust formation flight control.

Findings

Position error remains bounded under disturbances.

Enhanced drag reduction achieved in simulations.

Control algorithm stabilizes nonlinear formation dynamics.

Abstract

This paper investigates the robust nonlinear close formation control problem. It aims to achieve precise position control at dynamic flight operation for a follower aircraft under the aerodynamic impact due to the trailing vortices generated by a leader aircraft. One crucial concern is the control robustness that ensures the boundedness of position error subject to uncertainties and disturbances to be regulated with accuracy. This paper develops a robust nonlinear formation control algorithm to fulfill precise close formation tracking control. The proposed control algorithm consists of baseline control laws and disturbance observers. The baseline control laws are employed to stabilize the nonlinear dynamics of close formation flight, while the disturbance observers are introduced to compensate system uncertainties and formation-related aerodynamic disturbances. The position control performance can be guaranteed within the desired boundedness to harvest enough drag reduction for a follower aircraft in close formation using the proposed design. The efficacy of the proposed design is demonstrated via numerical simulations of close formation flight of two aircraft.

Tables

Table 1. Table 1: Aircraft geometry and mass parameters

Parameter	Symbol	Value	Unit	Parameter	Symbol	Value	Unit
Wing area	$S$	$27.87$	$m^2$	Vertical tail height	$h_t$	$3.05$	$m$
Vertical tail area	$S_v$	$5.09$	$m^2$	Quarter-chord sweep angle	${\Lambda }_s$	$0.57$	$rad$
Horizontal tail area	$S_h$	$10.034$	$m^2$	Dihedral angle	${\Lambda }_d$	$0$	$rad$
Wing span	$b$	$9.14$	$m$	Gross mass	$m$	$9295.44$	$kg$
Tail wing span	$b_t$	$5.49$	$m$	Roll moment of inertia	$I_x$	$12874.8$	$kg\cdot m^2$
Mean aerodynamic chord	$\overline{c}$	$3.45$	$m$	Pitch moment of inertia	$I_y$	$75673.6$	$kg\cdot m^2$
Root chord	$c_r$	$5.02$	$m$	Yaw moment of inertia	$I_z$	$85552.1$	$kg\cdot m^2$
Tip chord	$c_t$	$1.07$	$m$	Product moment of inertia	$I_{xz}$	$1331.4$	$kg\cdot m^2$

Table 2. Table 2: Aircraft aerodynamic coefficients

Coefficient	Symbol	Value	Coefficient	Symbol	Value
Zero-lift drag coefficient	$C_{D_0}$	$0.02$	Rolling moment coefficient	$C_{{ℒ}_{{\delta }_r}}$	$0.02636$
Oswald efficiency number	$e_o$	$0.663$	Pitching moment coefficient	$C_{{ℳ}_0}$	$-0.02029$
Section lift curve slope	$C_{l_{\alpha }}$	$5.3$	Pitch stiffness	$C_{{ℳ}_{\alpha }}$	$0.0466$
Lift coefficient	$C_{L_0}$	$0.05$	Pitch damping	$C_{{ℳ}_q}$	$-5.159$
Wing lift curve slope	$C_{L_{\alpha }}$	$5.3$	Pitching moment coefficient	$C_{{ℳ}_{{\delta }_e}}$	$-0.60123$
Vertical tail efficiency factor	$c_{\eta }$	$0.95$	Weathercock stability derivative	$C_{{𝒩}_{\beta }}$	$0.2993$
Dihedral effect	$C_{{ℒ}_{\beta }}$	$-0.1059$	Cross coupling derivative	$C_{{𝒩}_p}$	$0.02678$
Roll damping coefficient	$C_{{ℒ}_p}$	$-0.4127$	Yaw damping coefficient	$C_{{𝒩}_r}$	$-0.36988$
Vertical tail effect coefficient	$C_{{ℒ}_r}$	$0.0625$	Yawing moment coefficient	$C_{{𝒩}_{{\delta }_a}}$	$-0.03349$
Rolling moment coefficient	$C_{{ℒ}_{{\delta }_a}}$	$-0.1463$	Yawing moment coefficient	$C_{{𝒩}_{{\delta }_r}}$	$-0.081159$

Table 3. Table 3: Command filter parameters for 𝐋 c subscript 𝐋 𝑐 \mathbf{L}_{c} and χ ˙ r subscript ˙ 𝜒 𝑟 \dot{\chi}_{r}

Natural frequency	Value	Damping ratio	Value	Natural frequency	Value	Damping ratio	Value
${\omega }_{l_x}$	5	${\zeta }_{l_x}$	1	${\omega }_{l_y}$	5	${\zeta }_{l_y}$	1
${\omega }_{l_z}$	5	${\zeta }_{l_z}$	1	${\omega }_{{\chi }_r}$	5	${\zeta }_{{\chi }_r}$	1

Table 4. Table 4: Outerloop control parameters

Parameter	Symbol	Value	Parameter	Symbol	Value	Parameter	Symbol	Value
Time constant	${𝒯}_{Wx}$	$0.8$	Control gain	$K_x$	$0.3$	Control gain	$C_y$	${10}^{-4}$
Time constant	${𝒯}_{Wy}$	$0.5$	Control gain	$K_z$	$0.2$	Control gain	$C_z$	$5\times {10}^{-7}$
Time constant	${𝒯}_{Wz}$	$0.4$	Control gain	$K_V$	$1.75$	Natural frequency	${\omega }_V$	$8$
Time constant	${𝒯}_V$	$0.25$	Control gain	$K_{\gamma }$	$0.75$	Natural frequency	${\omega }_{\gamma }$	$8$
Time constant	${𝒯}_{\gamma }$	$0.2$	Control gain	$K_{\chi }$	$1.75$	Damping ratio	${\zeta }_V$	$1$
Time constant	${𝒯}_{\chi }$	$0.2$	Control gain	$C_x$	${10}^{-5}$	Damping ratio	${\zeta }_{\gamma }$	$1$

Table 5. Table 5: Innerloop control parameters

Parameter	Symbol	Value	Parameter	Symbol	Value	Parameter	Symbol	Value
Time constant	${𝒯}_{\mu }$	$0.8$	Control gain	$K_p$	$12$	Natural frequency	${\omega }_q$	$5$
Time constant	${𝒯}_{\alpha }$	$0.8$	Control gain	$K_q$	$7.5$	Natural frequency	${\omega }_r$	$5$
Time constant	${𝒯}_{\beta }$	$0.8$	Control gain	$K_r$	$7.5$	Damping ratio	${\zeta }_{\mu }$	$1$
Time constant	${𝒯}_p$	$0.02$	Control gain	$c_p$	${10}^{-5}$	Damping ratio	${\zeta }_{\alpha }$	$1$
Time constant	${𝒯}_q$	$0.02$	Natural frequency	${\omega }_{\alpha }$	$8$	Damping ratio	${\zeta }_p$	$1$
Time constant	${𝒯}_r$	$0.02$	Natural frequency	${\omega }_{\mu }$	$8$	Damping ratio	${\zeta }_q$	$1$
Control gain	$K_{\mu }$	$5$	Control gain	$c_r$	${10}^{-5}$	Damping ratio	${\zeta }_r$	$1$
Control gain	$K_{\alpha }$	$5$	Control gain	$c_q$	${10}^{-5}$
Control gain	$K_{\beta }$	$5$	Natural frequency	${\omega }_p$	$25$

Equations

$$\Vert \mathbf{x}\left(t_0\right)\Vert \le {\mathcal{A}}_a\Rightarrow \Vert \mathbf{x}\left(t\right)\Vert \le {\mathcal{A}}_b\text{, }\forall t\ge t_0+\mathcal{T}$$

$${\mathcal{T}}_d\dot{\hat{d}}=-\hat{d}+d$$

$$x_d=x_l+l_x\text{,}{y}_d=y_l+l_y\text{,}\text{and}{z}_d=z_l+l_z$$

$${\dot{x}}_d=V_l\cos {\gamma }_l\cos {\chi }_l+{\dot{l}}_x\text{, }{\dot{y}}_d=V_l\cos {\gamma }_l\sin {\chi }_l+{\dot{l}}_y\text{,}\text{and}{\dot{z}}_d=-V_l\sin {\gamma }_l+{\dot{l}}_z$$

$$\left[\begin{array}[]{c}\dot{\mathscr{S}}_{c}\\ \ddot{\mathscr{S}}_{c}\end{array}\right]=\underbrace{\left[\begin{array}[]{cc}0\text{,}&1\\ -\omega_{\mathscr{S}}^{2}\text{,}&-2\zeta_{\mathscr{S}}\omega_{\mathscr{S}}\end{array}\right]}_{\mathbf{A}_{\mathscr{S}}}\left[\begin{array}[]{c}{\mathscr{S}}_{c}\\ \dot{\mathscr{S}}_{c}\end{array}\right]+\left[\begin{array}[]{c}0\\ \omega_{\mathscr{S}}^{2}\end{array}\right]\mathscr{S}\left(t\right)$$

$$\Vert {\mathbf{e}}_{\mathcal{S}}{\Vert }_2\le \sqrt{\frac{{\lambda }_{max}\left({\mathbf{P}}_{\mathcal{S}}\right)}{{\lambda }_{min}\left({\mathbf{P}}_{\mathcal{S}}\right)}}{e}^{-\frac{t}{2{\lambda }_{max}\left({\mathbf{P}}_{\mathcal{S}}\right)}}|\dot{\mathcal{S}}\left(0\right)|+\left(1-e^{-\frac{t}{2{\lambda }_{max}\left({\mathbf{P}}_{\mathcal{S}}\right)}}\right)\frac{2{\lambda }_{max}^2\left({\mathbf{P}}_{\mathcal{S}}\right)\Vert \ddot{\mathcal{S}}+2{\zeta }_{\mathcal{S}}{\omega }_{\mathcal{S}}\dot{\mathcal{S}}{\Vert }_{\infty }}{{\lambda }_{min}\left({\mathbf{P}}_{\mathcal{S}}\right)}$$

$$\tilde{\mathcal{S}}=\mathcal{O}\left(1/{\omega }_{\mathcal{S}}\right)\text{ and }\dot{\tilde{\mathcal{S}}}/{\omega }_{\mathcal{S}}=\mathcal{O}\left(1/{\omega }_{\mathcal{S}}\right)\text{vspace}-3mm$$

$$\left\{\begin{array}[]{l}\dot{{x}}_{r}=\dot{x}_{l}+\dot{l}_{cx}\text{,}\qquad\dot{{y}}_{r}=\dot{y}_{l}+\dot{l}_{cy}\sin{\chi_{r}}\text{,}\qquad\dot{{z}}_{r}=\dot{z}_{l}+\dot{l}_{cz}\\ V_{r}=\sqrt{\dot{x}_{r}^{2}+\dot{y}_{r}^{2}+\dot{z}_{r}^{2}}\text{, }\quad\gamma_{r}=-\sin^{-1}\left(\frac{\dot{z}_{r}}{V_{r}}\right)\text{, }\quad\chi_{r}=\chi_{l}+\sin^{-1}\left(\frac{-\dot{l}_{cx}\sin{\chi_{l}}+\dot{l}_{cy}\cos{\chi_{l}}}{V_{r}\cos{\gamma_{r}}}\right)\end{array}\right.$$

$$\left\{\begin{array}[]{lcl}\dot{x}_{f}&=&V_{f}\cos{\gamma_{f}}\cos{\chi_{f}}+W_{x}\\ \dot{y}_{f}&=&V_{f}\cos{\gamma_{f}}\cos{\chi_{f}}+W_{y}\\ \dot{z}_{f}&=&-V_{f}\sin{\gamma_{f}}+W_{z}\end{array}\right.\qquad\left\{\begin{array}[]{lcl}\dot{V}_{f}&=&\frac{1}{m_{f}}\left(T\cos{\alpha_{f}}\cos{\beta_{f}}-\overline{D}\right)-g\sin{\gamma_{f}}+d_{V}\\ \dot{\gamma}_{f}&=&\frac{\left(\overline{L}+T\sin{\alpha_{f}}\right)\cos{\mu_{f}}}{m_{f}V_{f}}-\frac{g}{V_{f}}\cos{\gamma_{f}}+d_{\gamma}\\ \dot{\chi}_{f}&=&\frac{\left(\overline{L}+T\sin{\alpha_{f}}\right)\sin{\mu_{f}}}{m_{f}V_{f}\cos{\gamma_{f}}}+d_{\chi}\end{array}\right.$$

$$\left\{\begin{array}[]{lcl}d_{V}&=&-\frac{D-\overline{D}+\Delta D}{m_{f}}-\frac{\dot{W}_{W_{x}}}{m_{f}}\\ d_{\gamma}&=&\frac{\left(L-\overline{L}+\Delta L\right)\cos{\mu_{f}}-\left(Y+\Delta Y-T\cos{\alpha_{f}}\sin{\beta_{f}}\right)\sin{\mu_{f}}}{m_{f}V_{f}}-\frac{\dot{W}_{W_{y}}\sin{\mu_{f}}+\dot{W}_{W_{z}}\cos{\mu_{f}}}{m_{f}V_{f}}\\ d_{\chi}&=&\frac{\left(L-\overline{L}+\Delta L\right)\sin{\mu_{f}}+\left(Y+\Delta Y-T\cos{\alpha_{f}}\sin{\beta_{f}}\right)\cos{\mu_{f}}}{m_{f}V_{f}\cos{\gamma_{f}}}+\frac{\dot{W}_{W_{y}}\cos{\mu_{f}}-\dot{W}_{W_{z}}\sin{\mu_{f}}}{m_{f}V_{f}\cos{\gamma_{f}}}\\ \end{array}\right.$$

$$\hat{\mathbf{W}}={\mathbf{\lambda }}_W+{\mathcal{T}}_W^{-1}{\mathbf{X}}_{\mathbf{P}}\text{, }{\dot{\mathbf{\lambda }}}_W=-{{\mathcal{T}}^{-1}}_W{\mathbf{\lambda }}_W-{{\mathcal{T}}^{-1}}_W\left({\mathcal{T}}_W^{-1}{\mathbf{X}}_{\mathbf{P}}+{\mathbf{U}}_{\mathbf{P}}\right)$$

$${\dot{\tilde{W}}}_x=-{\tilde{W}}_x/{\mathcal{T}}_{Wx}\text{, }{\dot{\tilde{W}}}_y=-{\tilde{W}}_y/{\mathcal{T}}_{Wy}\text{,}\text{and }{\dot{\tilde{W}}}_z=-{\tilde{W}}_z/{\mathcal{T}}_{Wz}$$

$${\dot{x}}_f={\hat{V}}_f\cos {\hat{\gamma }}_f\cos {\hat{\chi }}_f\text{,}{\dot{y}}_f={\hat{V}}_f\cos {\hat{\gamma }}_f\sin {\hat{\chi }}_f\text{,}\text{and}{\dot{z}}_f=-{\hat{V}}_f\sin {\hat{\gamma }}_f$$

$$\small\left\{\begin{array}[]{ll}\widehat{V}_{f}&=\sqrt{\left({V}_{f}\cos{{\gamma}_{f}}\cos{{\chi}_{f}}+\widehat{W}_{x}\right)^{2}+\left({V}_{f}\cos{{\gamma}_{f}}\cos{{\chi}_{f}}+\widehat{W}_{y}\right)^{2}+\left(-{V}_{f}\sin{{\gamma}_{f}}+\widehat{W}_{z}\right)^{2}}\\ \widehat{\gamma}_{f}&=-\sin^{-1}\left(\frac{\widehat{W}_{z}-{V}_{f}\sin{\gamma_{f}}}{\widehat{V}_{f}}\right)\\ \widehat{\chi}_{f}&=\chi_{f}+\sin^{-1}\left(\frac{\widehat{W}_{y}\cos{\chi_{f}}-\widehat{W}_{x}\sin{\chi_{f}}}{\widehat{V}_{f}\cos{\widehat{\gamma}_{f}}}\right)\end{array}\right.$$

$$e_x=\cos {\hat{\chi }}_f{x}_e+\sin {\hat{\chi }}_f{y}_e\text{, }{e}_y=-\sin {\hat{\chi }}_f{x}_e+\cos {\hat{\chi }}_f{y}_e\text{, }{e}_z=z_e$$

$$\left\{\begin{array}[]{lcl}\dot{e}_{x}&=&\widehat{V}_{f}\cos{\widehat{\gamma}_{f}}-{V}_{r}\cos{{\gamma}_{r}}\cos{e_{\chi}}+\dot{\widehat{\chi}}_{f}e_{y}\\ \dot{e}_{y}&=&{V}_{r}\cos{{\gamma}_{r}}\sin{e_{\chi}}-\dot{\widehat{\chi}}_{f}e_{x}\\ \dot{e}_{z}&=&-\widehat{V}_{f}\sin{\widehat{\gamma}_{f}}+{V}_{r}\sin{{\gamma}_{r}}\end{array}\right.$$

$$V_d=\left(-K_x{e}_x+V_r\cos {\gamma }_r\cos e_{\chi }\right)/\cos {\hat{\gamma }}_f-\delta V\text{, }{\gamma }_d={\sin }^{-1}\left(\left(K_z{e}_z+V_r\sin {\gamma }_r+{\hat{W}}_z\right)/V_f\right)$$

$$\left[\begin{array}[]{c}\dot{\mathscr{S}}_{c}\\ \ddot{\mathscr{S}}_{c}\end{array}\right]=\left[\begin{array}[]{cc}0\text{,}&1\\ -\omega_{\mathscr{S}}^{2}\text{,}&-2\zeta_{\mathscr{S}}\omega_{\mathscr{S}}\end{array}\right]\left[\begin{array}[]{c}\mathscr{S}_{c}\\ \dot{\mathscr{S}}_{c}\end{array}\right]+\left[\begin{array}[]{c}0\\ \omega_{\mathscr{S}}^{2}\end{array}\right]\mathscr{S}_{d}\text{, }\quad\mathscr{S}\in\left\{V\text{, }\gamma\right\}$$

$$\left\{\begin{array}[]{lcl}\dot{e}_{V}&=&\underbrace{\frac{T\cos{\alpha_{f}}\cos{\beta_{f}}-\overline{D}}{m_{f}}-g\sin{\gamma_{f}}}_{u_{V}}-\dot{V}_{c}+d_{V}\\ \dot{e}_{\gamma}&=&\underbrace{\frac{\left(\overline{L}_{0}+\overline{L}_{\alpha}\alpha_{f}+T\sin{\alpha_{f}}\right)\cos{\mu_{f}}-m_{f}g\cos{\gamma_{f}}}{m_{f}V_{f}}}_{u_{\gamma}}-\dot{\gamma}_{c}+d_{\gamma}\\ \dot{e}_{\chi}&=&\underbrace{\frac{\left(\overline{L}_{0}+\overline{L}_{\alpha}\alpha_{f}+T\sin{\alpha_{f}}\right)\sin{\mu_{f}}}{m_{f}V_{f}\cos{\gamma_{f}}}}_{u_{\chi}}-\dot{\widehat{\chi}}_{r}+{d}_{\chi}\end{array}\right.$$

$$u_V^d=u_{V0}^d-{\hat{d}}_V+{\dot{V}}_c\text{,}{u}_{\gamma }^d=u_{\gamma 0}^d-{\hat{d}}_{\gamma }+{\dot{\gamma }}_c\text{,}\text{and}{u}_{\chi }^d=u_{\chi 0}^d-{\hat{d}}_{\chi }+{\dot{\hat{\chi }}}_r\text{vspace}-3mm$$

$$u_{V0}^d=-K_V{e}_V-\frac{c_V{\epsilon }_x\cos {\hat{\gamma }}_f}{H}\text{,}{u}_{\gamma 0}^d=-K_{\gamma }{e}_{\gamma }\text{,}\text{and}{u}_{\chi 0}^d=-K_{\chi }{e}_{\chi }-\frac{c_{\chi }{e}_y{V}_r\cos {\gamma }_r\cos \left(\frac{e_{\chi }}{2}\right)}{H}$$

$${\dot{\xi }}_x=-K_x{\xi }_x+\left(V_c-V_d\right)\cos {\hat{\gamma }}_f\text{,}{\dot{\xi }}_z=-K_z{\xi }_z+V_f\left(\sin {\gamma }_d-\sin {\gamma }_f\right)\text{vspace}-3mm$$

$$\left\{\begin{array}[]{lcl}\dot{V}_{f}&=&-K_{V}e_{V}-c_{V}\frac{\varepsilon_{x}\cos\widehat{\gamma}_{f}}{H}+\dot{V}_{c}-\widehat{d}_{V}+d_{V}\\ \dot{\gamma}_{f}&=&-K_{\gamma}e_{\gamma}+\dot{\gamma}_{c}-\widehat{d}_{\gamma}+d_{\gamma}\\ \dot{\chi}_{f}&=&-K_{\chi}\sin\left(\frac{e_{\chi}}{2}\right)-\frac{c_{\chi}e_{y}V_{r}\cos{\gamma_{r}}\cos\left(\frac{e_{\chi}}{2}\right)}{H}+\dot{\widehat{\chi}}_{r}-\widehat{{d}}_{\chi}+{d}_{\chi}\end{array}\right.$$

$${\hat{\mathbf{D}}}_{\mathbf{D}}={\mathbf{\lambda }}_{\mathbf{D}}+{\mathcal{T}}_{\mathbf{D}}^{-1}{\mathbf{X}}_{\mathbf{D}}\text{, }{\dot{\mathbf{\lambda }}}_{\mathbf{D}}=-{\mathcal{T}}_{\mathbf{D}}^{-1}{\mathbf{\lambda }}_{\mathbf{D}}-{\mathcal{T}}_{\mathbf{D}}^{-1}\left({\mathcal{T}}_{\mathbf{D}}^{-1}{\mathbf{X}}_{\mathbf{D}}+{\mathbf{U}}_{\mathbf{D}}-{\hat{\mathbf{D}}}_{\mathbf{D}}\right)$$

$$\left\{\begin{array}[]{lcl}T_{c}&=&\frac{m_{f}\left(u_{V}^{d}+g\sin{\gamma_{f}}\right)+\overline{D}}{\cos{\alpha_{f}}\cos{\beta_{f}}}\\ \alpha_{d}&=&\frac{m_{f}V_{f}\sqrt{\left(u_{\gamma}^{d}+\frac{g}{V}\cos{\gamma_{f}}\right)^{2}+\left(u_{\chi}^{d}\right)^{2}\cos^{2}{\gamma_{f}}}-T\sin{\alpha_{f}}-\overline{L}_{0}}{\overline{L}_{\alpha}}\\ \mu_{d}&=&\tan^{-1}\left(\frac{m_{f}V_{f}u_{\chi}^{d}\cos{\gamma_{f}}}{m_{f}V_{f}u_{\gamma}+m_{f}g\cos{\gamma_{f}}}\right)\end{array}\right.$$

$$\left\{\begin{array}[]{lcl}\dot{\varepsilon}_{x}&=&\widehat{V}_{f}\cos{\widehat{\gamma}_{f}}-{V}_{r}\cos{{\gamma}_{r}}\cos{e_{\chi}}-\dot{\xi}_{x}+\dot{\widehat{\chi}}_{f}e_{y}\\ \dot{e}_{y}&=&{V}_{r}\cos{{\gamma}_{r}}\sin{e_{\chi}}-\dot{\widehat{\chi}}_{f}\varepsilon_{x}\\ \dot{\varepsilon}_{z}&=&-\widehat{V}_{f}\sin{\widehat{\gamma}_{f}}+{V}_{r}\sin{{\gamma}_{r}}-\dot{\xi}_{z}\end{array}\right.$$

$$\underset{t\to \infty }{\lim }{e}_x\to {\xi }_x\text{, }\underset{t\to \infty }{\lim }{e}_z\to {\xi }_z\text{, }\text{ and }\underset{t\to \infty }{\lim }{e}_y\to 0\text{vspace}-3mm$$

$$\underset{t\to \infty }{\lim }{\xi }_x\le \frac{1}{K_x}\Vert V_c-V_d{\Vert }_{\infty }\text{and }\underset{t\to \infty }{\lim }{\xi }_x\le \frac{1}{K_z}\Vert V_f{\Vert }_{\infty }\Vert \sin {\gamma }_d-\sin {\gamma }_f{\Vert }_{\infty }\text{vspace}-3mm$$

$$\dot{{\tilde{d}}_V}=-\frac{{\tilde{d}}_V}{{\mathcal{T}}_V}\text{,}\dot{{\tilde{d}}_{\gamma }}=-\frac{{\tilde{d}}_{\gamma }}{{\mathcal{T}}_{\gamma }}\text{,}\text{and}\dot{{\tilde{d}}_V}=-\frac{{\tilde{d}}_{\chi }}{{\mathcal{T}}_{\chi }}\text{vspace}-3mm$$

$$\mathbb{V}=H+\frac{{\epsilon }_z^2}{2}+8c_{\chi }^{-1}{\sin }^2\left(\frac{e_{\chi }}{4}\right)+c_V^{-1}\frac{e_V^2}{2}+\frac{e_{\gamma }^2}{2}+c_{\chi }^{-1}\frac{{\tilde{d}}_{\chi }^2}{2}+c_V^{-1}\frac{{{\tilde{d}}_V}^2}{2}+\frac{{\tilde{d}}_{\gamma }^2}{2}-1$$

Full text

Nomenclature

[TABLE]

[TABLE]

I Introduction

Close formation flight is inspired by migratory birds who adopt a “V-shape” formation flight when migrating from one habitat to another Lissaman1970Science ; Weimerskirch2001Nature ; Portugal2014Nature . In close formation, a follower aircraft, holding a close relative position to a leader aircraft, flies in the upwash wake region of the trailing vortices induced by the leader aircraft as shown in Figure 1, by which the follower aircraft reduces its drag and thus saves fuels. Drag reduction of close formation flight has been demonstrated by simulations Halaas2014SciTech ; Kent2015JGCD ; Zhang2017JA , observed by wind tunnel experiments Bangash2006JA ; Cho2014JMST , and confirmed by flight tests Ray2002AFMCE ; Pahle2012AFMC ; Bieniawski2014ASM ; Hanson2018AFMC .

Despite its benefits, close formation flight is challenging in terms of the accuracy and robustness requirement for guidance and control Zhang2017JA . The position control accuracy must be guaranteed within the consideration of system uncertainties and formation-related aerodynamic disturbances. Yet, different control algorithms have been proposed for close formation flight. Most of them are focused on level and straight flight with constant speeds Schumacher2000ACC ; Pachter2001JGCD ; Dogan2005JGCD ; Almeida2015GNC . Two different linear strategies were applied, namely formation holding control and formation tracking control. Both of them are limited. Formation-holding control assumes a follower aircraft is initially well-trimmed at its optimal position in close formation, such as the proportional-integral (PI) formation control Pachter2001JGCD , the close formation control by the linear-quadratic regulator (LQR) Dogan2005JGCD , and the linear model predictive control (MPC)-based control Almeida2015GNC . Linear formation-tracking control doesn’t require the follower aircraft to be initially located at its desired position in close formation Binetti2003JGCD ; Lavretsky2003GNC ; Chichka2006JGCD ; Zhang2016GNC ; Zhang2017GNC ; Zhang2018AESCTE , but they are not guaranteed to address complex nonlinear aircraft dynamics at dynamic flight operation. Additionally, linear control methods will experience dramatic performance degradation or even fail to stabilize the system, when being applied to nonlinear systems. Robust nonlinear formation control is, therefore, more preferable to accommodate close formation flight at dynamic operation.

Nonlinear close formation control is challenging. Contrary to linear cases, model uncertainties and aerodynamic disturbances are less predictable and harder to be described at nonlinear scenarios, making close formation control more difficult. Early robust nonlinear close formation control was investigated using sliding model control Singh2000IJRNC or high order sliding mode (HOSM) control Galzi2007ACC . Both of the two methods only focus on outer-loop design, and they require the vortex-upper bounds of induced forces or their derivatives to satisfy certain limits to guarantee the stability. The nonlinear robust formation design including both inner-loop and outer-loop control was reported in Brodecki2015JGCD using an incremental nonlinear dynamic inversion (INDI) method, but this method is not robust to model uncertainties and external disturbances. Therefore, present nonlinear methods either rely on unknown model information to ensure both stability and robustness like the sliding mode control Galzi2007ACC , or not robust to model uncertainties and external disturbances at general dynamic operations, such as the INDI-based control Brodecki2015JGCD . Therefore, it is still an open issue for nonlinear robust close formation control with certain performance guaranty only using available model information.

In this paper, we investigate the robust nonlinear control problem for close formation flight at dynamic operation. The control design is presented under a leader-following architecture. The fundamental objective is to secure highly precise position control for close formation flight at dynamic flight operation design with the consideration of system uncertainties and aerodynamic impact caused by trailing vortices of leader aircraft a robust nonlinear controller. Control robustness will be one of the critical concerns which significantly affects the possible accuracy for close formation flight as it is subject to system uncertainties and aerodynamic disturbances. A robust nonlinear formation controller, which consists of baseline controllers and disturbance observers, is proposed in this paper. The baseline controllers are designed based on a command filtered backstepping technique to stabilize the nominal nonlinear dynamics of an aircraft in close formation Farrell2005JGCD ; Farrell2009TAC , whereas the disturbance observers could estimate and compensate for system uncertainties and formation-related aerodynamic disturbances by purely observing system inputs and available states. In the proposed design, the follower aircraft is required to track its dynamic optimal relative position to a leader aircraft in the inertial frame under different flight maneuvers. Both the inner-loop and outer-loop control will be studied in this paper, which makes the formation control design more reliably but also more difficult. The assumption on a well-designed inner-loop controller in Zhang2018TIE is, therefore, removed in this paper. The proposed design is capable of achieving highly accurate and efficiently robust control performance without using any gradient or boundary information of formation aerodynamic disturbances. Position tracking errors will be ultimately bounded. The final boundaries could be regulated by choosing different control parameters.

The rest of the paper is organized as follows. Section II presents some preliminaries, while Section III formulates reference trajectories. In Section IV, robust nonlinear control is reported and analyzed. Numerical simulations are given in Section V. Conclusions are in Section VI.

II Preliminaries

Some preliminaries are provided, which will be used for the design and analysis in the sequel.

Definition 1** (Definition 4.6 Khalil2002Book ).**

A system $\dot{\mathbf{x}}=\mathbf{f}(t\text{, }\mathbf{x})$ is uniformly ultimately bounded if there are positive constants $\mathcal{A}_{b}$ and $\mathcal{A}_{c}$, there exists $\mathcal{T}=\mathcal{T}\left(\mathcal{A}_{a}\text{, }\mathcal{A}_{b}\right)$ for any $\mathcal{A}_{a}\in\left(0\text{, }\mathcal{A}_{c}\right)$, such that

[TABLE]

Lemma 1**.**

Let $d\left(t\right)$ be a bounded signal whose derivative $\dot{d}\left(t\right)$, namely $d\left(t\right)\in\mathscr{L}_{\infty}$ and $\dot{d}\left(t\right)\in\mathscr{L}_{\infty}$. Assume $\widehat{d}\left(t\right)$ is the estimation of $d\left(t\right)$ through a first-order filter as shown in (2).

[TABLE]

where $\mathcal{T}_{d}>0$ is a time constant. Define $\widetilde{d}=\widehat{d}-d$ as the estimation error. If $\widehat{d}\left(0\right)=0$,

$\widetilde{d}$* are globally bounded by $|\widetilde{d}|\leq\max\left\{|{d}\left(0\right)|\text{, }\mathcal{T}_{d}\|\dot{d}\|_{\infty}\right\}$;* 2. 2.

if $|{d}\left(0\right)|\neq\mathcal{T}_{d}\|\dot{d}\|_{\infty}$, there exist any small positive constant $\epsilon$ and time $t_{\epsilon}$ such that $|\widetilde{d}|<\mathcal{T}_{d}\|\dot{d}\|_{\infty}+\epsilon$ for all $t>t_{\epsilon}$, where $t_{\epsilon}=\max\left\{0\text{, }\mathcal{T}_{d}\ln\left(\frac{||{d}\left(0\right)|-\mathcal{T}_{d}\|\dot{d}\|_{\infty}|}{\epsilon}\right)\right\}$; 3. 3.

if $\lim_{t\to\infty}\dot{d}=0$, $\lim_{t\to\infty}\widetilde{d}=0$.

III Formulation of reference trajectories at dynamic operation

In this section, a motion planner is designed for follower aircraft at close formation. According to Zhang2017JA , the optimal relative position in close formation is static in the wind frame of the leader aircraft. Assume $\left[r_{x}\text{, }r_{y}\text{, }0\right]^{T}$ is the static optimal relative position in the wind frame of the leader aircraft, where $r_{x}$ ranges from $-2b$ to $-10b$ and $r_{y}$ is around $\pm 0.95b$ with $b$ denoting the wing span. When flying at close formation, the reference position of a follower aircraft in the inertial frame is

[TABLE]

where $x_{l}$, $y_{l}$, and $z_{l}$ are position coordinates of the leader aircraft in the inertial frame, and $\left[l_{x}\text{, }l_{y}\text{, }l_{z}\right]^{T}=\mathbf{C}_{IW}\left(\mu_{l},\gamma_{l},\chi_{l}\right)\left[r_{x}\text{, }r_{y}\text{, }0\right]^{T}$ where $\mu_{l}$, $\gamma_{l}$, and $\chi_{l}$ are the bank, flight path, and heading angles of the leader aircraft. Differentiating (3) yields

[TABLE]

where $V_{l}$, $\gamma_{l}$, and $\chi_{l}$ are the airspeed, flight path angle, and heading angle of leader aircraft, respectively. At dynamic operation, $l_{x}$, $l_{y}$, and $l_{z}$ are time-varying, but their derivatives cannot be accurately computed. Hence, in the design, we introduce a command filter (5) to get the command signals $l_{ci}$ and $\dot{l}_{ci}$ ($i\in\left\{x\text{, }y\text{, }z\right\}$). Let $\mathscr{S}\left(t\right)$ be a smooth signal, so the command filter is

[TABLE]

where $\omega_{\mathscr{S}}>0$ is the natural frequency and $\zeta_{\mathscr{S}}>0$ is the damping ratio. Let $\tilde{\mathscr{S}}=\mathscr{S}_{c}-\mathscr{S}$ and $\mathbf{e}_{\mathscr{S}}=\left[\tilde{\mathscr{S}}\text{, }\dot{\tilde{\mathscr{S}}}\right]^{T}$. If $\mathscr{S}_{c}\left(0\right)=\mathscr{S}\left(0\right)$ and $\dot{\mathscr{S}}_{c}\left(0\right)=0$, Lemma 2 exists for any bounded signal $\mathscr{S}\left(t\right)$.

Lemma 2**.**

The estimator (5) is input-to-state stable with respect to $\mathscr{S}\left(t\right)$. If both ${\dot{\mathscr{S}}}\left(t\right)$ and ${\ddot{\mathscr{S}}}\left(t\right)$ are bounded, $\tilde{\mathscr{S}}$ is uniformly and ultimately bounded, and the following inequality exists for $\mathbf{e}_{\mathscr{S}}\left(t\right)$.

[TABLE]

where $\lambda_{max}\left(\cdot\right)$ and $\lambda_{min}\left(\cdot\right)$ are the maximum and minimum eigenvalues of a matrix, respectively, $\mathbf{P}_{\mathscr{S}}$ is positive definite such that $\mathbf{P}_{\mathscr{S}}\mathbf{A}_{\mathscr{S}}+\mathbf{A}^{T}_{\mathscr{S}}\mathbf{P}_{\mathscr{S}}=-\mathbf{I}$. Furthermore, there exist

[TABLE]

where $\mathcal{O}\left(\cdot\right)$ is an order of magnitude notation Khalil2002Book .

In real flight, $\dot{\gamma}_{l}$, $\dot{\chi}_{l}$, and $\dot{\mu}_{l}$ and their derivatives are all bounded, so (5) is a valid choice to obtain $l_{ci}$ and $\dot{l}_{ci}$ ($i\in\left\{x\text{, }y\text{, }z\right\}$). In addition, Lemma 2 indicates that the command signals $l_{ci}$ and $\dot{l}_{ci}$ ($i\in\left\{x\text{, }y\text{, }z\right\}$) could be ensured to be arbitrarily close to their corresponding desired values $l_{i}$ and $\dot{l}_{i}$ by choosing proper command filter parameters. Note that $\mathscr{S}_{c}\left(0\right)=\mathscr{S}\left(0\right)$ is needed to avoid the peaking phenomenon (Page 613, Khalil2002Book ). If $\mathscr{S}_{c}\left(0\right)\neq\mathscr{S}\left(0\right)$, $\mathscr{S}_{c}$, the signal $\dot{\mathscr{S}}_{c}$ will transiently peak to $\mathcal{O}\left({\omega_{\mathscr{S}}}\right)$ before it exponentially decays, resulting in the so-called peaking phenomenon due to $\mathscr{S}_{c}\left(0\right)\neq\mathscr{S}\left(0\right)$. Without loss of generality, the following assumption is introduced.

Assumption 1**.**

The attitude signals $\gamma_{l}$, $\chi_{l}$, and $\mu_{l}$ and their derivatives are all bounded.

In light of (5), the command position for a follower aircraft in close formation is ${x}_{r}=x_{l}+{l}_{cx}$, ${y}_{r}=y_{l}+{l}_{cy}$, and ${z}_{r}=z_{l}+{l}_{cz}$, and accordingly,

[TABLE]

IV Robust nonlinear formation control design

The proposed design in this section can be easily extended to the case with more than three aircraft, even though it is discussed under the leader-follower architecture with two aircraft. In the proposed design, command filtered backstepping technique is employed, which avoids the analytic calculation of time derivatives of intermediate virtual inputs Farrell2005JGCD ; Sonneveldt2007JGCD ; Farrell2009TAC ; Sonneveldt2009JGCD . As shown in Figure 2, the entire design consists of two major loops: an outer loop for formation position control and an inner loop for attitude control. The outer-loop control allows a follower aircraft to track the planned motion by (6), and generates command thrust $T_{c}$, desired bank angle $\mu_{d}$, and desired angle of attack $\alpha_{d}$. The inner-loop control stabilizes follower aircraft’s attitudes to their desired values $\mu_{d}$ and $\alpha_{d}$ from the outer-loop control, while holding zero sideslip angle $\beta_{f}$.

IV.1 Outer-loop formation position control

Let $\overline{D}$ and $\overline{L}$ be the nominal values of the drag $D$ and lift $L$, respectively. They are obtained by either available aerodynamic data or certain analytical models Morelli1998ACC . The sideslip angle $\beta_{f}$ is negligibly small, as it is always stabilized to be zero. Accordingly, the side force $Y$ is small and taken as a model uncertainty. The outer-loop dynamics used for control design are

[TABLE]

where $x_{f}$, ${y}_{f}$, and ${z}_{f}$ are follower position coordinates in the inertial frame, $V_{f}$ is the airspeed, $\gamma_{f}$ and $\chi_{f}$ are the flight path and heading angles, $T$ is the thrust, $W_{x}$, $W_{y}$, and $W_{z}$ are induced wake velocities, and $d_{V}$, $d_{\gamma}$, and $d_{\chi}$ are the augmentation of system uncertainties and disturbances.

[TABLE]

where $\dot{W}_{W_{x}}$, $\dot{W}_{W_{y}}$, and $\dot{W}_{W_{z}}$ are the wake velocity derivatives denoted in the wind frame of follower aircraft, $\Delta L$, $\Delta D$, and $\Delta Y$ are the vortex-induced forces. According to Zhang2017JA , $W_{x}$, $W_{y}$, and $W_{z}$ are bounded, and have much slower dynamics in comparison with aircraft speed and attitudes, so their derivatives are relatively small. Furthermore, the following assumption is introduced.

Assumption 2**.**

Induced wake velocities $W_{x}$, $W_{y}$, and $W_{z}$ are all bounded, and furthermore, they are piecewise constant, namely $\dot{W}_{{x}}=0$, $\dot{W}_{{y}}=0$, and $\dot{W}_{{z}}=0$.

The following nonlinear disturbance observer is employed.

[TABLE]

where $\boldsymbol{\mathcal{T}}_{W}=diag\left\{\mathcal{T}_{Wx}\text{, }\mathcal{T}_{Wy}\text{, }\mathcal{T}_{Wz}\right\}>0$, $\boldsymbol{\lambda}_{W}=\left[\lambda_{Wx}\text{, }\lambda_{Wy}\text{, }\lambda_{Wz}\right]^{T}$, $\mathbf{X_{P}}=\left[x_{f}\text{, }y_{f}\text{, }z_{f}\right]^{T}$, $\mathbf{{U}_{P}}=\left[V_{f}\cos{\gamma_{f}}\cos{\chi_{f}}\text{, }V_{f}\cos{\gamma_{f}}\sin{\chi_{f}}\text{, }-V_{f}\sin{\gamma_{f}}\right]^{T}$, $\widehat{\mathbf{W}}=\left[\widehat{{W}}_{x}\text{, }\widehat{{W}}_{y}\text{, }\widehat{{W}}_{z}\right]^{T}$ where $\widehat{{W}}_{x}$, $\widehat{{W}}_{y}$, and $\widehat{{W}}_{z}$ are estimates of $W_{x}$, $W_{y}$, and $W_{z}$, respectively. It is chosen that $\boldsymbol{\lambda}_{W}\left(0\right)=-\boldsymbol{\mathcal{T}}_{W}^{-1}\mathbf{X_{P}}\left(0\right)$. Let $\widetilde{W}_{x}=\widehat{W}_{x}-W_{x}$, $\widetilde{W}_{y}=\widehat{W}_{y}-W_{y}$, and $\widetilde{W}_{z}=\widehat{W}_{z}-W_{z}$. Under Assumption 2, one has

[TABLE]

According to Lemma 1, $\widetilde{W}_{x}$, $\widetilde{W}_{y}$, and $\widetilde{W}_{z}$ can be made arbitrarily small by choosing sufficiently small time constants, even if $\dot{W}_{{x}}\neq 0$, $\dot{W}_{{y}}\neq 0$, and $\dot{W}_{{z}}\neq 0$. To simplify the analysis, it is assumed that $\widetilde{W}_{x}=\widetilde{W}_{y}=\widetilde{W}_{z}=0$. In light of (9), one has

[TABLE]

where

[TABLE]

Let $x_{e}=x_{f}-{x}_{r}$, $y_{e}=y_{f}-{y}_{r}$, and $z_{e}=z_{f}-{z}_{r}$. Transform $x_{e}$, $y_{e}$, and $z_{e}$ into a new frame.

[TABLE]

We have

[TABLE]

where $e_{\chi}=\widehat{\chi}_{f}-\chi_{r}$. The desired velocity and flight path angle are shown in (13).

[TABLE]

where $K_{x}\text{, }K_{z}>0$ are control parameters, and $\delta V=\widehat{V}_{f}-V_{f}$. The desired signals $V_{d}$ and $\gamma_{d}$ are passed through a command filter to obtain $V_{c}$, $\gamma_{c}$, and their rates. Hence,

[TABLE]

Let $e_{V}=V_{f}-V_{c}$ and $e_{\gamma}=\gamma_{f}-\gamma_{c}$. Note that $\overline{L}=\overline{L}_{0}+\overline{L}_{\alpha}\alpha_{f}$ where $\overline{L}_{0}$ is the lift at $\alpha_{f}=0$ and $\overline{L}_{\alpha}$ is the lift derivative with respect to the angle of attack. According to (7), one has

[TABLE]

where $u_{V}$, $u_{\gamma}$, and $u_{\chi}$ are intermediate control inputs, $\dot{\widehat{\chi}}_{r}$ is the estimation of $\dot{\chi}_{r}$ by passing $\chi_{r}$ through a 2nd-order filter similar to (5). Two more uncertainty terms ${\dot{\widetilde{\chi}}}_{r}$ and $\dot{\widetilde{\chi}}_{f}$ are included in ${d}_{\chi}$ in (15), where ${\dot{\widetilde{\chi}}}_{r}={\dot{\widehat{\chi}}}_{r}-\dot{\chi}_{r}$ and $\dot{\widetilde{\chi}}_{f}=\frac{d}{dt}\sin^{-1}\left(\frac{{\widehat{W}}_{y}\cos{\chi_{f}}-\widehat{W}_{x}\sin{\chi_{f}}}{\widehat{V}_{f}\cos{\widehat{\gamma}_{f}}}\right)$. Hence, ${d}_{\chi}$ in (15) is re-defined to be ${d}_{\chi}=\frac{\left(L-\overline{L}+\Delta L\right)\sin{\mu_{f}}+\left(Y+\Delta Y-\cos{\alpha_{f}}\sin{\beta_{f}}\right)\cos{\mu_{f}}}{m_{f}V_{f}\cos{\gamma_{f}}}+\dot{\widetilde{\chi}}_{r}+\dot{\widetilde{\chi}}_{f}+\frac{\dot{W}_{W_{y}}\cos{\mu_{f}}-\dot{W}_{W_{z}}\sin{\mu_{f}}}{m_{f}V_{f}\cos{\gamma_{f}}}$. Note that $\chi_{r}$ is a smooth signal with bounded derivatives. According to Lemma 2, ${\dot{\widetilde{\chi}}}_{r}$ and its derivative are bounded, so ${\dot{\widetilde{\chi}}}_{f}$ and its derivative are also bounded.

Assumption 3**.**

The uncertainties and disturbances $d_{V}$, $d_{\gamma}$, $d_{\chi}$, and their derivatives are bounded.

The following virtual inputs $u_{V}^{d}$, $u_{\gamma}^{d}$, and $u_{\chi}^{d}$ are choosing for the error systems (12) and (15).

[TABLE]

where $\widehat{d}_{V}$, $\widehat{d}_{\gamma}$, and $\widehat{{d}}_{\chi}$ are estimates of ${d}_{V}$, ${d}_{\gamma}$, and ${{d}}_{\chi}$, respectively, and $u_{V0}^{d}$, $u_{\gamma 0}^{d}$, and $u_{\chi 0}^{d}$ are

[TABLE]

where $K_{V}$, $K_{\gamma}$, $K_{\chi}$, $c_{V}$, $c_{\chi}>0$ are control parameters, $H=\sqrt{\varepsilon_{x}^{2}+e_{y}^{2}+1}$, $\varepsilon_{x}=e_{x}-\xi_{x}$, and $\varepsilon_{z}=e_{z}-\xi_{z}$, where $\xi_{x}$ and $\xi_{z}$ in (18) are used to counteract the estimation errors in the filter (14).

[TABLE]

Shown in Figure 3 is the command filter and auxiliary system for speed control. If one chooses $u_{V}={u}_{V}^{d}$, $u_{\gamma}={u}_{\gamma}^{d}$, and $u_{\chi}={u}_{\chi}^{d}$, there exists

[TABLE]

Based on (19), the nonlinear disturbance observer is

[TABLE]

where $\boldsymbol{\mathcal{T}_{D}}=diag\left\{\mathcal{T}_{V}\text{, }\mathcal{T}_{\gamma}\text{, }\mathcal{T}_{\chi}\right\}>0$ is a positive definite constant matrix, $\boldsymbol{\lambda_{D}}=\left[\lambda_{V}\text{, }\lambda_{\gamma}\text{, }\lambda_{\chi}\right]^{T}$, $\mathbf{X_{D}}=\left[V_{f}\text{, }\gamma_{f}\text{, }\chi_{f}\right]^{T}$, $\mathbf{U_{D}}=\left[{u}_{V0}^{d}\text{, }u_{\gamma 0}^{d}\text{, }u_{\chi 0}^{d}\right]^{T}$, and $\mathbf{\widehat{D}_{D}}=\left[\widehat{d}_{V}\text{, }\widehat{d}_{\gamma}\text{, }\widehat{d}_{\chi}\right]^{T}$. The uncertainty and disturbance estimates $\widehat{d}_{V}$, $\widehat{d}_{\gamma}$, and $\widehat{{d}}_{\chi}$ need to be fed back to the estimator for the next estimation updates. Combining (17) and (20), one is able to get $u_{V}^{d}$, $u_{\gamma}^{d}$, and $u_{\chi}^{d}$. Hence,

[TABLE]

The entire outer-loop control structure is illustrated in Figure 4. The following assumption is introduced for ${d_{V}}$, ${d_{\gamma}}$, and ${d_{\chi}}$ for the stability analysis.

Assumption 4**.**

${d_{V}}$, ${d_{\gamma}}$, and ${d_{\chi}}$ have slow dynamics, and furthermore, $\dot{d}_{V}\simeq 0$, $\dot{d}_{\gamma}\simeq 0$, and $\dot{d}_{\chi}\simeq 0$.

The following error dynamics will be obtained.

[TABLE]

Assume $\alpha_{f}$ and $\mu_{f}$ are able to be rapidly stabilized by its inner-loop attitude controller to their desired values $\alpha_{d}$ and $\mu_{d}$, respectively. The following theorem holds.

Theorem 1**.**

If Assumptions 3 and 4 hold, and $K_{V}$, $K_{\gamma}$, $K_{\chi}$, $\mathcal{T}_{V}$, $\mathcal{T}_{\gamma}$, $\mathcal{T}_{\chi}$, $c_{V}$, $c_{\chi}>0$, $0<K_{x}<2\zeta_{V}\omega_{V}$, and $0<K_{z}<2\zeta_{\gamma}\omega_{\gamma}$, the proposed outer-loop formation controller given by (16), (17), and (20) will stabilize the outer-loop formation error system composed of (15) and (22), and

[TABLE]

where $\xi_{x}$ and $\xi_{z}$ will be ultimately arbitrarily bounded by control parameters $\omega_{V}$, $K_{x}$, and $\omega_{\gamma}$, $K_{z}$, respectively. If $\xi_{x}\left(0\right)=\xi_{z}\left(0\right)=0$, there exist

[TABLE]

where $V_{c}-V_{d}$ will exponentially converge to $\mathcal{O}\left(\frac{1}{\omega_{V}}\right)$ and $\lim_{t\to\infty}\gamma_{f}-\gamma_{d}\to\mathcal{O}\left(\frac{1}{\omega_{\gamma}}\right)$. Therefore, by tuning $\omega_{V}$, $K_{x}$, and $\omega_{\gamma}$, $K_{z}$, the ultimate boundaries of ${e}_{x}$ and ${e}_{z}$ could be regulated accordingly.

Proof.

Let $\widetilde{d}_{i}=\widehat{d}_{i}-d_{i}$ ($i\in\left\{V\text{, }\gamma\text{, }\chi\right\}$). If Assumption 4 holds, and $\alpha_{f}$ and $\mu_{f}$ are rapidly stabilized to be $\alpha_{d}$ and $\mu_{d}$, respectively, it is easily to obtain

[TABLE]

Based on (15), (22), and (24), Theorem 1 is proven in two steps. The first step shows $\lim_{t\to\infty}\varepsilon_{x}\to 0$, $\lim_{t\to\infty}e_{y}\to 0$, and $\lim_{t\to\infty}\varepsilon_{z}\to 0$, and thus $\lim_{t\to\infty}e_{x}\to\xi_{x}$ and $\lim_{t\to\infty}e_{z}\to\xi_{z}$. The second step demonstrates that $\xi_{x}$ and $\xi_{z}$ are uniformly ultimately bounded, which implies $e_{x}$ and $e_{z}$ are ultimately bounded. The ultimate boundaries of $e_{x}$ and $e_{z}$ are related to control gains.

To show $\lim_{t\to\infty}\varepsilon_{x}\to 0$, $\lim_{t\to\infty}e_{y}\to 0$, and $\lim_{t\to\infty}\varepsilon_{z}\to 0$, we choose

[TABLE]

where $H=\sqrt{\varepsilon_{x}^{2}+e_{y}^{2}+1}$. Differentiating $\mathbb{V}$ yields

[TABLE]

Substituting (16) and (17) for corresponding terms in (25) yields

[TABLE]

If $K_{V}\text{, }K_{\chi}\text{, }K_{\gamma}\text{, }\mathcal{T}_{V}\text{, }\mathcal{T}_{\gamma}\text{, }\mathcal{T}_{\chi}>0$, one has

[TABLE]

Hence, there exists $\dot{\mathbb{V}}_{1}\leq 0$ by choosing $K_{x}>0$, $K_{z}>0$, $c_{V}>0$, and $c_{\chi}>0$. Therefore, $\mathbb{V}$ is a non-increasing positive definite function, which implies $\mathbb{V}(\infty)\leq\mathbb{V}(0)$ is a finite constant and $\lim_{t\to\infty}\int^{t}_{0}\dot{\mathbb{V}}_{1}(\tau)d\tau=\mathbb{V}(\infty)-\mathbb{V}(0)$ exists and is finite. Thus, $\mathbb{V}$ is bounded and has a finite limit as $t\to\infty$, and $\varepsilon_{x}$, $e_{y}$, $\varepsilon_{z}$, $e_{\chi}$, $e_{V}$, $e_{\gamma}$, $\widetilde{d}_{V}$, $\widetilde{d}_{\chi}$, and $\widetilde{d}_{\gamma}$ are all bounded as well. Furthermore, $\ddot{\mathbb{V}}_{1}$ is also bounded, as it is a function of $\varepsilon_{x}$, $e_{y}$, $\varepsilon_{z}$, $e_{\chi}$, $e_{V}$, $e_{\gamma}$, $\widetilde{d}_{V}$, $\widetilde{d}_{\chi}$, and $\widetilde{d}_{\gamma}$. Hence, $\dot{\mathbb{V}}_{1}$ is uniformly continuous. According to Barbǎlat’s lemma Ioannou1995Book , $\lim_{t\to\infty}\dot{\mathbb{V}}_{1}\to 0$, so $\lim_{t\to\infty}\varepsilon_{x}\to 0\text{, }\varepsilon_{z}\to 0\text{, }\lim_{t\to\infty}e_{\chi}\to 0\text{, }\lim_{t\to\infty}e_{V}\to 0\text{, }\lim_{t\to\infty}e_{\gamma}\to 0$, $\lim_{t\to\infty}\widetilde{d}_{V}\to 0$, $\lim_{t\to\infty}\widetilde{d}_{\chi}\to 0$, $\lim_{t\to\infty}\widetilde{d}_{\gamma}\to 0$. It is obviously that $\lim_{t\to\infty}e_{x}\to\xi_{x}$ and $\lim_{t\to\infty}e_{z}\to\xi_{z}$. One could apply Barbǎlat’s lemma to $\dot{e}_{\chi}$ to show $\lim_{t\to\infty}e_{y}\to 0$. Differentiating $\dot{e}_{\chi}$ by time will yield

[TABLE]

Obviously, $\ddot{\chi}_{e}$ is bounded according to Lemma 1 and Assumptions 1 and 3. Considering $\int\dot{e}_{\chi}dt=e_{\chi}$ and $\lim_{t\to\infty}e_{\chi}\to 0$, we have $\dot{\psi}_{e}\to 0$ as $t\to\infty$ in accordance with Barbǎlat’s lemma. According to (15), (16), and (17), $e_{y}=-\frac{\dot{e}_{\chi}+K_{\chi}e_{\chi}+\widetilde{d}_{\chi}}{c_{\chi}V_{r}\cos{\gamma_{r}}\cos{\frac{e_{\chi}}{2}}}H$. Since $\lim_{t\to\infty}\dot{e}_{\chi}\to 0$, $\lim_{t\to\infty}e_{\chi}\to 0$, $\lim_{t\to\infty}\widetilde{d}_{\chi}\to 0$, $c_{\chi}V_{r}\cos{\gamma_{r}}\cos{\frac{e_{\chi}}{2}}>0$, and $H>0$, it is readily to conclude that $\lim_{t\to\infty}e_{y}\to 0$.

The second step shows that $\xi_{x}$ and $\xi_{z}$ are uniformly ultimately bounded. From (13), one has

[TABLE]

Substituting the equation above for $V_{d}$ in (14) and (18) will generate

[TABLE]

The characteristic equation of (27) is $s^{3}+2\zeta_{V}\omega_{V}s^{2}+\omega_{V}^{2}s+K_{x}\omega_{V}^{2}=0$. According to the Routh–Hurwitz stability criterion, $\mathbf{A}_{v}$ is ensured to be Hurwitz all the time, if $0<K_{x}<2\zeta_{V}\omega_{V}$. Therefore, the system (27) is input-to-state stable with respect to $\frac{V_{r}\cos{\gamma_{r}}\cos{e_{\chi}}-K_{x}\varepsilon_{x}}{\cos{\widehat{\gamma}_{f}}}+\delta V$. According to the previous analysis, $\varepsilon_{x}$ will asymptotically converge to [math]. By Assumptions 1 and 2, both $\delta V$ and ${V_{r}\cos{\gamma_{r}}\cos{\xi_{e}}}$ are uniformly bounded. Since $\widehat{\gamma}_{f}\in\left(-\frac{\pi}{2}\text{, }\frac{\pi}{2}\right)$, $\frac{V_{r}\cos{\gamma_{r}}\cos{e_{\chi}}-K_{x}\varepsilon_{x}}{\cos{\widehat{\gamma}_{f}}}+\delta V$ is uniformly bounded, which implies $\xi_{x}$ is uniformly bounded.

The final boundary of $\xi_{x}$ is related to $\omega_{V}$ and $K_{x}$. If $\xi_{x}\left(0\right)=0$, there exists ${\xi}_{x}\leq\frac{1}{K_{x}}\left(1-e^{-K_{x}t}\right)\|V_{c}-V_{d}\|_{\infty}$, so $\lim_{t\to\infty}{\xi}_{x}\leq\frac{1}{K_{x}}\|V_{c}-V_{d}\|_{\infty}$. Therefore, by tuning $K_{x}$, the ultimate boundary of ${\xi}_{x}$ could be changed. In addition, $\xi_{x}$ is bounded-input-bounded-output with respect to $V_{c}-V_{d}$, and Lemma 2 indicates that $V_{c}-V_{d}$ will exponentially converge to $\mathcal{O}\left(\frac{1}{\omega_{V}}\right)$. The ultimate boundary of $\xi_{x}$ can be, therefore, altered by changing $\omega_{V}$ as well. Since $\lim_{t\to\infty}e_{x}\to\xi_{x}$, the ultimate boundary of $e_{x}$ could thus be regulated by changing $K_{x}$ and $\omega_{V}$.

According to (18), the dynamics of $\xi_{z}$ are input-to-state stable with respect to $V_{f}\left(\sin\gamma_{d}-\sin\gamma_{f}\right)$. Since $V_{f}\left(\sin\gamma_{d}-\sin\gamma_{f}\right)$ is bounded, $\xi_{z}$ is ultimately bounded. According to our previous analysis, one knows that $\lim_{t\to\infty}\gamma_{f}\to\gamma_{c}$, while $\gamma_{c}$ is from a 2nd-order command filter (14). As $\gamma_{c}-\gamma_{d}=\mathcal{O}\left(\frac{1}{\omega_{\gamma}}\right)$, $\lim_{t\to\infty}\gamma_{f}-\gamma_{d}\to\mathcal{O}\left(\frac{1}{\omega_{\gamma}}\right)$, and thus $\sin\gamma_{d}-\sin\gamma_{f}$ will eventually be limited by certain small boundaries related to $\omega_{\gamma}$. Therefore, $\xi_{z}$ is uniformly ultimately bounded. Choosing $\xi_{z}\left(0\right)=0$ implies ${\xi}_{z}\leq\frac{1}{K_{z}}\left(1-e^{-K_{z}t}\right)\|V_{f}\|_{\infty}\|\sin{\gamma_{d}}-\sin{\gamma_{f}}\|_{\infty}$, so $\lim_{t\to\infty}{\xi}_{z}\leq\frac{1}{K_{z}}\|V_{f}\|_{\infty}\|\sin{\gamma_{d}}-\sin{\gamma_{f}}\|_{\infty}$. As $\lim_{t\to\infty}e_{x}\to\xi_{x}$, $\lim_{t\to\infty}e_{z}\to\xi_{z}$, and both $\xi_{x}$ and $\xi_{z}$ are uniformly ultimately bounded, $e_{x}$ and $e_{z}$ will be ultimately bounded, respectively. ∎

IV.2 Inner-loop attitude control

The inner-loop dynamics of a follower aircraft in close formation is

[TABLE]

where $\mu_{f}$ is the bank angle, $\alpha_{f}$ is the angle of attack, $\beta_{f}$ is the sideslip angle, $p$, $q$, and $r$ are angular rates in the body frame, $\mathcal{L}$, $\mathcal{M}$, and $\mathcal{N}$ are moments, while $\Delta\mathcal{L}$, $\Delta\mathcal{M}$, and $\Delta\mathcal{N}$ are the moment disturbances induced by trailing vortices. The presented attitude controller will rapidly stabilize $\alpha_{f}$ and $\mu_{f}$ to their desired values $\alpha_{d}$ and $\mu_{d}$, respectively. Meanwhile, the sideslip angle $\beta_{f}$ will be kept to be [math]. Regarding the first objective, a command filter is employed again to estimate the derivatives of $\alpha_{d}$ and $\mu_{d}$, respectively.

[TABLE]

Define $e_{\alpha}=\alpha_{f}-\alpha_{d}$ and $e_{\mu}=\mu_{f}-\mu_{d}$. Let $\mathbf{e}_{\Theta}=\left[e_{\mu}\text{, }e_{\alpha}\text{, }\beta_{f}\right]^{T}$, $\boldsymbol{\Theta}_{d}=\left[\mu_{d}\text{, }\alpha_{d}\text{, }0\right]^{T}$, $\boldsymbol{\Theta}_{c}=\left[\mu_{c}\text{, }\alpha_{c}\text{, }0\right]^{T}$, $\boldsymbol{\Psi}=\left[\gamma_{f}\text{, }\chi_{f}\right]^{T}$, and $\boldsymbol{\Omega}=\left[p\text{, }q\text{, }r\right]^{T}$. According to (28), we have

[TABLE]

where $\boldsymbol{\dot{\Psi}}$ will be estimated by $\boldsymbol{\widehat{\dot{\Psi}}}=\left[u_{\gamma}+\widehat{d}_{\gamma}\text{, }u_{\chi}+\widehat{d}_{\chi}\right]^{T}$, and

[TABLE]

When Assumptions 3 and 4 hold, $\lim_{t\to\infty}\left(u_{\gamma}+\widehat{d}_{\gamma}\right)\to\dot{\gamma}$, $\lim_{t\to\infty}\left(u_{\chi}+\widehat{d}_{\chi}\right)\to\dot{\chi}$ as $\lim_{t\to\infty}\widehat{d}_{\gamma}\to{d}_{\gamma}$ and $\lim_{t\to\infty}\widehat{d}_{\chi}\to{d}_{\chi}$.

Since Assumption 4 is barely met, $\mathbf{H}\left(\boldsymbol{{\dot{\Psi}}}-\boldsymbol{\widehat{\dot{\Psi}}}\right)$ is treated as model uncertainties. Additionally, $\boldsymbol{\dot{\Theta}}_{d}$ is replaced by $\boldsymbol{\dot{\Theta}}_{c}$ in terms of (29), as it is unavailable. Eventually, the model uncertainties of (30) are $\mathbf{d}_{\Theta}=\mathbf{H}\left(\boldsymbol{{\dot{\Psi}}}-\boldsymbol{\widehat{\dot{\Psi}}}\right)+\left(\boldsymbol{\dot{\Theta}}_{c}-\boldsymbol{\dot{\Theta}}_{d}\right)$. Define $\mathbf{u}_{\Theta}=\mathbf{G}\boldsymbol{\Omega}+\mathbf{H}-\boldsymbol{\widehat{\dot{\Psi}}}$, so

[TABLE]

The desired intermediate virtual inputs to stabilize (31) are proposed in (32).

[TABLE]

where $\mathbf{K}_{\Theta}=diag\left\{K_{\mu}\text{, }K_{\alpha}\text{, }K_{\beta}\right\}>0$ is a gain matrix and $\mathbf{\widehat{d}}_{\Theta}$ is the estimation of $\mathbf{d}_{\Theta}$, which is from

[TABLE]

where$\boldsymbol{\mathcal{T}}_{\Theta}=diag\left\{\mathcal{T}_{\mu}\text{, }\mathcal{T}_{\alpha}\text{, }\mathcal{T}_{\beta}\right\}>0$ is a diagonal time constant matrix, and $\boldsymbol{\lambda}_{\Theta}=\left[\lambda_{\mu}\text{, }\lambda_{\alpha}\text{, }\lambda_{\beta}\right]^{T}$. Since $\mathbf{u}_{\Theta}=\mathbf{G}\boldsymbol{\Omega}+\mathbf{H}-\boldsymbol{\widehat{\dot{\Psi}}}$, the desired angular rates are given by

[TABLE]

where $\boldsymbol{\Omega}_{d}=\left[p_{d}\text{, }q_{d}\text{, }r_{d}\right]^{T}$. The commanded angular rates $p_{c}$, $q_{c}$, and $r_{c}$ are obtained by (35).

[TABLE]

where $\mathscr{S}_{c}\left(0\right)=\mathscr{S}_{d}\left(0\right)$ for $\mathscr{S}\in\left\{p\text{, }q\text{, }r\right\}$. Define $\mathbf{e}_{\Omega}=\boldsymbol{\Omega}-\boldsymbol{\Omega}_{c}$. According to (28), one has

[TABLE]

where $\boldsymbol{{\tau}}=\left[\mathcal{L}\text{, }\mathcal{M}\text{, }\mathcal{N}\right]^{T}$, $\boldsymbol{\Delta\tau}=\left[\Delta\mathcal{L}\text{, }\Delta\mathcal{M}\text{, }\Delta\mathcal{N}\right]^{T}$, and $\mathbf{I}$ is the inertia matrix of the aircraft. The control inputs are control surface deflections, including aileron deflection $\delta_{a}$, elevator deflection $\delta_{e}$, and rudder deflection $\delta_{r}$. Let $\boldsymbol{\delta_{u}}=\left[\delta_{a}\text{, }\delta_{e}\text{, }\delta_{r}\right]^{T}$, and $\boldsymbol{\tau}=\boldsymbol{\tau}_{0}+\mathbf{M_{\tau}}\boldsymbol{\delta_{u}}$, where $\mathbf{M_{\tau}}$ is the control derivative matrix and $\boldsymbol{\tau}_{0}$ is the torque vector at $\boldsymbol{\delta_{u}}=0$. Both $\boldsymbol{\tau}_{0}$ and $\mathbf{M_{\tau}}$ cannot be accurately obtained, so they are approximated using available aerodynamic data from wind tunnel tests. Let $\boldsymbol{\overline{\tau}}_{0}$ and $\mathbf{\overline{M}_{\tau}}$ be the approximate results of $\boldsymbol{\tau}_{0}$ and $\mathbf{M_{\tau}}$, respectively, so

[TABLE]

Let $\mathbf{u}_{\tau}=-\mathbf{I}^{-1}\boldsymbol{\Omega}\times\mathbf{I}\boldsymbol{\Omega}+\mathbf{I}^{-1}\left(\boldsymbol{\overline{\tau}}_{0}+\mathbf{\overline{M}_{\tau}}\boldsymbol{\delta_{u}}\right)=\left[u_{p}\text{, }u_{q}\text{, }u_{r}\right]^{T}$. The error model (36) is rewritten as

[TABLE]

where $\mathbf{d}_{\tau}=\mathbf{I}^{-1}\left(\boldsymbol{\tau}-\boldsymbol{\overline{\tau}}+\Delta\boldsymbol{\tau}\right)$ is the sum of model uncertainties and formation aerodynamic disturbances. Eventually, the control law for $\mathbf{u}_{\tau}$ is proposed to be

[TABLE]

where $\mathbf{K}_{\Omega}=diag\left\{K_{p}\text{, }K_{q}\text{, }K_{r}\right\}>0$ is a diagonal constant gain matrix, $\mathbf{C}_{\Omega}=diag\left\{c_{p}\text{, }c_{q}\text{, }c_{r}\right\}>0$ is a constant matrix, $\mathbf{\widehat{d}_{\tau}}$ is the estimation of $\mathbf{d_{\tau}}$, and $\boldsymbol{\varepsilon}_{\Theta}=\left[\varepsilon_{\mu}\text{, }\varepsilon_{\alpha}\text{, }\varepsilon_{\beta}\right]^{T}=\mathbf{e}_{\Theta}-\boldsymbol{\xi}_{\Theta}$ where $\boldsymbol{\xi}_{\Theta}$ is from

[TABLE]

The uncertainty and disturbance estimator for $\mathbf{d_{\tau}}$ is given by

[TABLE]

where $\boldsymbol{\mathcal{T}}_{\Omega}=diag\left\{\mathcal{T}_{p}\text{, }\mathcal{T}_{q}\text{, }\mathcal{T}_{r}\right\}>0$ is the time constant matrix. Let $\boldsymbol{\delta_{c}}=\left[\delta_{ac}\text{, }\delta_{ec}\text{, }\delta_{rc}\right]^{T}$ be the commanded control surface deflection vector. Eventually, we have

[TABLE]

The inner-loop controller is shown in Figure 5. The following assumption is introduced for the sake of stability analysis.

Assumption 5**.**

Both $\mathbf{d}_{\Theta}$ and $\boldsymbol{d}_{\tau}$ are bounded with slow dynamics, namely $\mathbf{\dot{d}}_{\Theta}\simeq 0$ and $\mathbf{\dot{d}}_{\tau}\simeq 0$.

Define $\mathbf{\widetilde{d}}_{\Theta}=\mathbf{\widehat{d}}_{\Theta}-\mathbf{d}_{\Theta}$ and $\mathbf{\widetilde{d}}_{\tau}=\mathbf{\widehat{d}}_{\tau}-\mathbf{d}_{\tau}$. Under Assumption 5, we have

[TABLE]

Lemma 3**.**

If Assumption 5 holds, and $K_{\mu}$, $K_{\alpha}$, $K_{\beta}$, $K_{p}$, $K_{q}$, $K_{r}$, $\mathcal{T}_{\mu}$, $\mathcal{T}_{\alpha}$, $\mathcal{T}_{\beta}$, $\mathcal{T}_{p}$, $\mathcal{T}_{q}$, $\mathcal{T}_{r}$, $c_{p}$, $c_{q}$, and $c_{r}>0$, $\boldsymbol{\varepsilon}_{\Theta}$ and ${\mathbf{e}}_{\Omega}$ will exponentially converge to zero, namely $\lim_{t\to\infty}{\mathbf{e}}_{\Theta}\to\boldsymbol{\xi}_{\Theta}$, and $\sigma_{c}-\sigma_{d}=\mathcal{O}\left(\frac{1}{\omega_{\sigma}}\right)$ and ${\xi}_{\sigma}=\mathcal{O}\left(\frac{1}{\omega_{\sigma}}\right)$ with $\sigma\in\left\{p\text{, }q\text{, }r\right\}$, so ${\mathbf{e}}_{\Theta}$ is ultimately bounded.

In real implementations, it is only required that $\mathbf{d}_{\Theta}$, $\boldsymbol{d}_{\tau}$, and their derivatives are bounded. If $\mathbf{\dot{d}}_{\Theta}\simeq 0$ and $\boldsymbol{\dot{d}}_{\tau}\simeq 0$ fail to exist, ${\xi}_{\sigma}=\mathcal{O}\left(\frac{1}{\omega_{\sigma}}\right)$ still holds uniformly. However, instead of achieving $\lim_{t\to\infty}\boldsymbol{\varepsilon}_{\Omega}\to 0$, we can only ensure $\|{\boldsymbol{\varepsilon}}_{\Theta}-\boldsymbol{\bar{e}}_{\Theta}\|=\mathcal{O}\left(\epsilon_{1}\right)$ and $\|\boldsymbol{e}_{\Omega}-\boldsymbol{\bar{e}}_{\Omega}\|=\mathcal{O}\left(\epsilon_{1}\right)$ where $\epsilon_{1}$ is a certain positive small value related to the time constants of the disturbance observers (33) and (41), and $\boldsymbol{\bar{e}}_{\Theta}$ and $\boldsymbol{\bar{e}}_{\Omega}$ are the attitude tracking errors from the standard backstepping design. The dynamics of $\boldsymbol{\bar{e}}_{\Theta}$ and $\boldsymbol{\bar{e}}_{\Omega}$ are shown in (44). Obviously, there exist $\lim_{t\to\infty}\boldsymbol{\bar{e}}_{\Theta}\to 0$ and $\lim_{t\to\infty}\boldsymbol{\bar{e}}_{\Omega}\to 0$. Since $\boldsymbol{\varepsilon}_{\Theta}=\boldsymbol{e}_{\Theta}-\boldsymbol{\xi}_{\Theta}$, $\boldsymbol{{e}}_{\Theta}$ will be ultimately bounded. This conclusion is summarized in Proposition 1.

[TABLE]

Proposition 1**.**

Assume $\mathbf{d}_{\Theta}$, $\boldsymbol{d}_{\tau}$, and their derivatives are bounded. The proposed inner-loop attitude control law composed of (32), (33), (34), (35), (39), and (41) will make $\|{\boldsymbol{\varepsilon}}_{\Theta}-\boldsymbol{\bar{e}}_{\Theta}\|=\mathcal{O}\left(\epsilon_{1}\right)$ and $\|\boldsymbol{e}_{\Omega}-\boldsymbol{\bar{e}}_{\Omega}\|=\mathcal{O}\left(\epsilon_{1}\right)$ uniformly hold, where $\epsilon_{1}=\max\left\{\mathcal{T}_{\mu}\text{, }\mathcal{T}_{\alpha}\text{, }\mathcal{T}_{\beta}\text{, }\mathcal{T}_{p}\text{, }\mathcal{T}_{q}\text{, }\mathcal{T}_{r}\right\}$. Furthermore, $\boldsymbol{e}_{\Theta}$ will be uniformly ultimately bounded, and $\|{\boldsymbol{e}}_{\Theta}-\boldsymbol{\bar{e}}_{\Theta}\|=\mathcal{O}\left(\epsilon_{2}\right)$ where $\epsilon_{2}=\max\left\{\epsilon_{1}\text{, }\frac{1}{\omega_{p}}\text{, }\frac{1}{\omega_{q}}\text{, }\frac{1}{\omega_{r}}\right\}$.

V Simulation verification

The proposed robust nonlinear close formation controller is validated based on two F-16 aircraft. A high-fidelity model presented in Russel2003Tech will be employed to simulate the nonlinear dynamics of a F-16 aircraft. The aerodynamic effects by close formation flight are characterized using the aerodynamic model developed in Zhang2017JA . Note that the aerodynamic data used to build up the nonlinear aircraft model are assumed to be unavailable for the control design. Instead, a global nonlinear parameter modeling technique in Morelli1998ACC is employed to calculate the necessary aerodynamic coefficients and $\overline{D}$, $\overline{L}_{0}$, $\overline{L}_{\alpha}$, $\boldsymbol{\overline{\tau}}_{0}$, and $\mathbf{\overline{M}_{\tau}}$. The aircraft geometry and mass parameters are listed in Table 1, while the necessary aerodynamic parameters are given in Table 2. The numerical simulations are carried out at two scenarios. In the first scenario, the robustness of the disturbance observer-based controller is verified by being compared with the control without disturbance observers (DO). In the second scenario, close formation flight is conducted at different velocities with the same group of control parameters to further confirm the efficacy of the proposed design.

V.1 Scenario 1: With/without disturbance observers

The initial conditions of the leader aircraft are $x_{l}\left(0\right)=45$ $m$, $y_{l}=-15$ $m$, $z_{l}=-5015$ $m$, $V_{l}=200$ $m/s$, $\gamma_{l}=\chi_{l}=0$ $deg$. According to the aerodynamic analysis in Zhang2017JA , the optimal relative position vector is selected to be $\left[-36\text{, }9\text{, }0\right]^{T}$ $m$. The parameters of the command filters for generating $\mathbf{L}_{c}$ and estimating are given in Table 3. The initial conditions for the follower aircraft are $x_{f}\left(0\right)=45$ $m$, $y_{f}\left(0\right)=-15$ $m$, $z_{f}\left(0\right)=-5015$ $m$, $V_{f}\left(0\right)=200$ $m/s$, $\gamma_{f}\left(0\right)=\chi_{f}\left(0\right)=\mu_{f}\left(0\right)=\beta_{f}\left(0\right)=0$ $deg$, $\alpha_{f}\left(0\right)=2.774$ $deg$, $p=q=r=0$ $rad/s$. The outerloop and innerloop control parameters are presented in Table 4 and 5, respectively. The formation trajectory in the inertial frame is illustrated in Figures 6. From [math] to $35$ seconds, the formation trajectory is at level and straight flight. From $35$ to $145$ seconds, aircraft in close formation are required to make a turn, and meanwhile reduce their altitudes. After $145$ seconds, level and straight trajectory is recovered. Shown in Figure 7 are the postures and relative positions of the leader and follower aircraft at four different time instants under the proposed robust nonlinear control. The follower aircraft is initially far away from its optimal position relative to the leader aircraft. Under the proposed robust nonlinear controller, the follower aircraft is able to quickly catch up the optimal relative position. Shown in Figure 8 are highlights of the top and front views of the relative positions between the leader and follower aircraft at $180$ seconds.

Two different nonlinear controllers are implemented in the first scenario. One is the proposed robust nonlinear formation controller, the other one is purely the baseline nonlinear formation controller without including disturbance observers (DO). The position tracking errors under the two controllers are shown in Figures 9, 10, and 11. According to Zhang2017JA , close formation flight will fail if the optimal lateral and vertical relative positions cannot be tracked with in $10\%$ wing span, while $90\%$ of the maximum drag reduction could be retained if the position tracking error is kept under $5\%$ wing span. For efficient close formation flight, it is required that both the lateral and vertical tracking errors are smaller than at least $5\%$ wing span. To show the validness of the proposed design, the regions covered by $5\%$ and $10\%$ wing span are highlighted in Figures 10 and 11. Obviously, the nonlinear control without including disturbance observers failed to achieve reasonable close formation flight, as the vertical position tracking errors are much far away from the region of interest. After incorporating disturbance observers (DO), both the lateral and vertical position tracking errors are confined to be smaller than $5\%$ wing span, so accurate close formation flight is fulfilled. Furthermore, position tracking errors under the proposed robust nonlinear controller would converge to zero when close formation is at level and straight flight. When the leader aircraft is taking maneuvers ($35$ to $145$ $s$), steady tracking errors are observed for lateral position tracking. The steady tracking errors are under $5\%$ wing span, which implies efficient close formation flight is still guaranteed. Shown in Figures 12, 13, and 14 are tracking errors of speed, flight path angle, and heading angle, respectively. The non-zero lateral steady tracking errors from $35$ to $145$ $s$ result from the non-zero tracking errors in heading angle control as shown in Figure 14. When leader aircraft is taking maneuvers, the heading angle will have non-zero tracking errors due to the difference between $\dot{\chi}_{r}$ and $\dot{\hat{\chi}}_{r}$. In real implementations, the derivative of the reference heading angle $\chi_{r}$ is always unavailable, so a second order filter is introduced to approximate $\dot{\chi}_{r}$. When the leader aircraft is under level and straight flight, $\chi_{r}$ is constant, which implies $\dot{\chi}_{r}\equiv 0$. In this case, $\dot{\hat{\chi}}_{r}$ could be ensured to be equal to $\dot{\chi}_{r}$, so asymptotical stability is able to be obtained as shown in Figures 9, 10, 11, and 14. However, when the leader aircraft is taking maneuvers, $\chi_{r}$ is not constant, and $\dot{\hat{\chi}}_{r}$ can only be guaranteed to converge to a certain value close to $\dot{\chi}_{r}$. The difference between $\dot{\hat{\chi}}_{r}$ and $\dot{\chi}_{r}$ might lead to the steady heading tracking errors from $35$ to $145$ $s$ as shown in Figure 14, which is reflected in lateral position tracking control as given in Figure 10.

The inner-loop state responses are summarized in Figure 15. The sideslip angle of the follower aircraft is always kept to be zero by the proposed robust nonlinear formation controller. Shown in Figure 16 are responses of control inputs. As mentioned before, the baseline formation controller without disturbance observers cannot achieve successful close formation flight. The follower aircraft under the baseline controller is eventually one span away from its optimal relative position to the leader aircraft, in which case the influence of the trailing vortices is quite small. Therefore, the steady performance of the follower aircraft by the baseline formation controller is similar to that of an aircraft at solo flight. Compared with the baseline controller, the robust nonlinear controller will eventually have $13.876\%$ decrease in throttle inputs, which implies that around $13.876\%$ energy saving could be obtained by close formation flight.

V.2 Scenario 2: Different flight speeds

The efficacy of the proposed robust nonlinear controller is further verified by running the close formation flight under different velocities. All the control parameters, the reference trajectories, and the initial conditions will be same as those at the first scenario. Without loss of generality, only position tracking errors and control input responses are given. As compared with a linear controller, nonlinear formation controller can be applied to much wider flight scenarios with stability and performance guaranty. This advantage of the nonlinear controller is demonstrated in Figure 17, 18, and 19. It is observed that the proposed robust nonlinear controller is able to ensure almost the same control performance under different velocities. The control input responses under different velocities are illustrated in Figure 20. Another interesting observation is that increasing speed will result in better performance in lateral position tracking as shown in Figure 18. Explaining this observation is difficult from the nonlinear perspective, but we could analyze it from a linear method at a special case. Assume $V_{f}=V_{r}=Const$, $\gamma_{f}=\gamma_{r}=0$, and $\chi_{r}=0$. The nonlinear closed outer-loop lateral dynamics are linearized about its equilibrium, which results in

[TABLE]

where $d_{\chi}$ denotes any uncertainties, disturbances, or inputs. The transfer function from $d_{\chi}$ to ${y}_{e}$ is

[TABLE]

According to the final value theorem, the increase of $V_{f}$ leads to smaller steady values in ${y}_{e}$ as illustrated in Figure 18.

VI Conclusions

The paper presented a robust nonlinear controller for autonomous close formation flight under different flight maneuvers. The proposed controller was developed by combining the command filtered backstepping method and the disturbance observation technique. Both inner-loop and outer-loop controllers were designed in this paper. Based on the proposed design, a follower aircraft is able to track its optimal relative position to a leader aircraft under different flight maneuvers. The proposed design was able to be extended to close formation flight of more than three aircraft, though it was described in the scenario of two-aircraft close formation. Enough robustness and high accuracy could be achieved by the presented design. Different numerical simulations were conducted to demonstrate the efficacy of the presented robust nonlinear close formation controller.

Appendix A: Proof of Lemma 1

For (2), the error dynamic equation is

[TABLE]

By solving the first-order differential equation (47), one has $\widetilde{d}\left(t\right)=e^{-\frac{t}{\mathcal{T}_{d}}}\widetilde{d}\left(0\right)+\int_{0}^{t}e^{-\frac{t-\tau}{\mathcal{T}_{d}}}\dot{d}\left(\tau\right)d\tau$. If choosing $\widehat{d}\left(0\right)=0$, we have $\widetilde{d}\left(0\right)=-{d}\left(0\right)$. If $\dot{d}\left(t\right)\in\mathscr{L}_{\infty}$,

[TABLE]

Obviously, $|\widetilde{d}|\leq\max\left\{|{d}\left(0\right)|\text{, }\mathcal{T}_{d}\|\dot{d}\|_{\infty}\right\}$ globally exists.

For any positive constant $\epsilon$, if $\mathcal{T}_{d}\|\dot{d}\|_{\infty}+\epsilon>|{d}\left(0\right)|$, there exists $|\widetilde{d}|<\mathcal{T}_{d}\|\dot{d}\|_{\infty}+\epsilon$ $\forall t_{\epsilon}>0$ according to the global boundedness of $|\widetilde{d}|$. If $\mathcal{T}_{d}\|\dot{d}\|_{\infty}+\epsilon<|{d}\left(0\right)|$, the first term on the right-hand side of (48) will decrease, while the second term will increase. At certain time $t_{\epsilon}$, there exists $e^{-\frac{t}{\mathcal{T}_{d}}}|{d}\left(0\right)|+\mathcal{T}_{d}\left(1-e^{-\frac{t}{\mathcal{T}_{d}}}\right)\|\dot{d}\|_{\infty}=\mathcal{T}_{d}\|\dot{d}\|_{\infty}+\epsilon$. After $t_{\epsilon}$, the right hand side of (48) will be smaller than $\mathcal{T}_{d}\|\dot{d}\|_{\infty}+\epsilon$. By simple mathematical calculation, one has $t_{\epsilon}=\mathcal{T}_{d}\ln\left(\frac{||{d}\left(0\right)|-\mathcal{T}_{d}\|\dot{d}\|_{\infty}|}{\epsilon}\right)$, so the second conclusion of Lemma 1 is obtained. In addition, the estimator error is input-to-state stable with respect to $\dot{d}$ according to (47). Hence, it is obvious that $\lim_{t\to\infty}\widetilde{d}=0$, if $\lim_{t\to\infty}\dot{d}=0$.

Appendix B: Proof of Lemma 2

The following error dynamics are easy to obtain.

[TABLE]

where $\mathbf{B}_{\mathscr{S}}=\left[0\text{, }-1\right]^{T}$. By choosing $\zeta_{\mathscr{S}}\text{, }\omega_{\mathscr{S}}>0$, $\mathbf{A}_{\mathscr{S}}$ is Hurwitz, which implies $\mathbf{P}_{\mathscr{S}}\mathbf{A}_{\mathscr{S}}+\mathbf{A}^{T}_{\mathscr{S}}\mathbf{P}_{\mathscr{S}}=-\mathbf{I}$ with $\mathbf{P}_{\mathscr{S}}>0$. Choose $\mathbb{V}=\mathbf{e}_{\mathscr{S}}^{T}\mathbf{P}_{\mathscr{S}}\mathbf{e}_{\mathscr{S}}$ as the Lyapunov function for (49), so

[TABLE]

If both $\ddot{\mathscr{S}}$ and $\dot{\mathscr{S}}$ are bounded, differentiating $\mathbb{V}$ with respect to time will yield

[TABLE]

Let $\mathbb{W}\left(t\right)=\sqrt{\mathbb{V}}$, so $\dot{\mathbb{W}}=\frac{\dot{\mathbb{V}}}{2\sqrt{\mathbb{V}}}$, when $\mathbb{V}\neq 0$. Hence, $\dot{\mathbb{W}}\leq-{\mathbb{W}}\left/{\left(2\lambda_{max}\left(\mathbf{P}_{\mathscr{S}}\right)\right)}\right.+{\lambda_{max}\left(\mathbf{P}_{\mathscr{S}}\right)}\|\ddot{\mathscr{S}}+2\zeta_{\mathscr{S}}\omega_{\mathscr{S}}\dot{\mathscr{S}}\|_{\infty}\left/{\sqrt{\lambda_{min}\left(\mathbf{P}_{\mathscr{S}}\right)}}\right.$. According to the comparison principle (Page 102, Khalil2002Book ), one can get

[TABLE]

According to (50), $\mathbb{W}\leq\sqrt{\lambda_{max}\left(\mathbf{P}_{\mathscr{S}}\right)}\|\mathbf{e}_{\mathscr{S}}\|_{2}$, so $\mathbb{W}\left(0\right)\leq\lambda_{max}\left(\mathbf{P}_{\mathscr{S}}\right)\|\mathbf{e}_{\mathscr{S}}\left(0\right)\|_{2}$. By setting $\mathscr{S}_{c}\left(0\right)=\mathscr{S}\left(0\right)$, and $\dot{\mathscr{S}}_{c}\left(0\right)=0$, one has $\|\mathbf{e}_{\mathscr{S}}\left(0\right)\|=|\dot{\mathscr{S}}\left(0\right)|$. Eventually,

[TABLE]

Therefore, $\mathbf{e}_{\mathscr{S}}$ is uniformly and ultimately bounded. With the consideration of $\|\tilde{\mathscr{S}}\|_{2}\leq\|\mathbf{e}_{\mathscr{S}}\|_{2}$, we are able to conclude that $\tilde{\mathscr{S}}$ is uniformly and ultimately bounded.

The second conclusion of Lemma 2 could be demonstrated by virtue of the perturbation theory. According to (49), the following singularly perturbed system is readily obtained.

[TABLE]

where $\tilde{\mathscr{S}}_{1}=\tilde{\mathscr{S}}$ and $\tilde{\mathscr{S}}_{2}=\frac{1}{\omega_{\mathscr{S}}}\dot{\tilde{\mathscr{S}}}$. Obviously, the effect of $\frac{1}{\omega_{\mathscr{S}}}\ddot{\mathscr{S}}+2\zeta_{\mathscr{S}}\dot{\mathscr{S}}$ will be diminished by increasing the natural frequency $\omega_{\mathscr{S}}$. According to the properties of singularly perturbed system, we readily have $\tilde{\mathscr{S}}_{1}=\mathscr{S}_{c}-\mathscr{S}=\mathcal{O}\left(\frac{1}{\omega_{\mathscr{S}}}\right)$ and $\dot{\tilde{\mathscr{S}}}_{2}=\frac{\dot{\tilde{\mathscr{S}}}}{\omega_{\mathscr{S}}}=\mathcal{O}\left(\frac{1}{\omega_{\mathscr{S}}}\right)$.

Appendix C: Proof of Lemma 3

In light of (31), (36), (32), (38), (39), (40), and (43), one has

[TABLE]

The stability of (51) is shown by picking the Lyapunov function $\mathbb{V}$ as below.

[TABLE]

The derivative of $\mathbb{V}$ is

[TABLE]

where $\mathbf{\widetilde{d}}_{\Theta}$ and $\mathbf{\widetilde{d}}_{\tau}$ are both vectors. Let $\mathbf{\widetilde{d}}_{\Theta}=\left[\widetilde{d}_{\mu}\text{, }\widetilde{d}_{\alpha}\text{, }\widetilde{d}_{\beta}\right]^{T}$ and $\mathbf{\widetilde{d}}_{\tau}=\left[\widetilde{d}_{p}\text{, }\widetilde{d}_{q}\text{, }\widetilde{d}_{r}\right]^{T}$, so

[TABLE]

It is easy to show that $\dot{\mathbb{V}}_{2}$ is negative definite, if all control parameters are chosen to be positive. Hence, $\boldsymbol{\varepsilon}_{\Theta}$ and ${\mathbf{e}}_{\Omega}$ will exponentially converge to zero. To show the second conclusion in Lemma 3, a singular perturbation model is established based on (35) and (40). To simplify the analysis complexity, assume that $\omega_{p}=\omega_{q}=\omega_{r}=\omega_{\Omega}$ and $\zeta_{p}=\zeta_{q}=\zeta_{r}=\zeta_{\Omega}$. Define a new variable $\boldsymbol{\bar{\Omega}}_{cdot}=\frac{\boldsymbol{\dot{\Omega}}_{c}}{\omega_{\Omega}}$. Therefore, the singular perturbation model is

[TABLE]

Obviously, the command filter system composed of (54b) and (54c) has much faster dynamics than the auxiliary system (54a), if $\omega_{\Omega}$ is chosen to be sufficiently large. The reduced system of (54) is given by $\boldsymbol{\dot{\bar{\xi}}}_{\Theta}=-\mathbf{K}_{\Theta}\boldsymbol{\bar{\xi}}_{\Theta}$ by setting $\omega_{\Omega}$ to be infinity, where $\boldsymbol{\bar{\xi}}_{\Theta}$ is the reduced system state vector. If $\boldsymbol{\xi}_{\Theta}\left(0\right)=\boldsymbol{\bar{\xi}}_{\Theta}\left(0\right)=0$, $\|\boldsymbol{\xi}_{\Theta}\|=\mathcal{O}\left(\frac{1}{\omega_{\Omega}}\right)$ will hold uniformly according to the Tikhonov’s theorem, which implies ${\xi}_{\sigma}=\mathcal{O}\left(\frac{1}{\omega_{\sigma}}\right)$. It is easy to get $\sigma_{c}-\sigma_{d}=\mathcal{O}\left(\frac{1}{\omega_{\sigma}}\right)$ according to Lemma 2.

Appendix D: Proof of Proposition 1

When $\mathbf{\dot{d}}_{\Theta}\neq 0$ and $\boldsymbol{\dot{D}}_{\tau}\neq 0$, the closed-loop error dynamics (51) will be rewritten as

[TABLE]

If the time constants $\boldsymbol{\mathcal{T}}_{\Theta}$ and $\boldsymbol{\mathcal{T}}_{\Omega}$ are chosen to be sufficiently small, (55) will be a standard perturbation model whose reduced system is given in (44). The reduced system (44) is apparently exponentially stable. Since both $\mathbf{\dot{d}}_{\Theta}$ and $\boldsymbol{\dot{D}}_{\tau}$ are bounded, their impact will be diminished by reducing $\boldsymbol{\mathcal{T}}_{\Theta}$ and $\boldsymbol{\mathcal{T}}_{\Omega}$, respectively. Since $\boldsymbol{\Omega}_{c}\left(0\right)=\boldsymbol{\Omega}_{c}\left(0\right)$ and $\boldsymbol{\xi}_{\Theta}\left(0\right)=0$, one has $\boldsymbol{{\varepsilon}}_{\Theta}\left(0\right)=\boldsymbol{\bar{e}}_{\Theta}\left(0\right)$ and $\boldsymbol{{e}}_{\Omega}\left(0\right)=\boldsymbol{\bar{e}}_{\Omega}\left(0\right)$. According to the Tikhonov’s theorem for a standard perturbation model, one is able to conclude that $\|{\boldsymbol{\varepsilon}}_{\Theta}-\boldsymbol{\bar{e}}_{\Theta}\|=\mathcal{O}\left(\epsilon_{1}\right)$ and $\|\boldsymbol{e}_{\Omega}-\boldsymbol{\bar{e}}_{\Omega}\|=\mathcal{O}\left(\epsilon_{1}\right)$ will uniformly hold. According to the definition of the oder of magnitude, it is easy to find that $\epsilon_{1}=\max\left\{\mathcal{T}_{\mu}\text{, }\mathcal{T}_{\alpha}\text{, }\mathcal{T}_{\beta}\text{, }\mathcal{T}_{p}\text{, }\mathcal{T}_{q}\text{, }\mathcal{T}_{r}\right\}$.

Furthermore, with the consideration of (54) and (55), we have

[TABLE]

[TABLE]

Notice that Eq. (57) will perform as fast dynamics, if $\omega_{\Omega}$ is chosen sufficiently large and $\boldsymbol{\mathcal{T}}_{\Theta}$ and $\boldsymbol{\mathcal{T}}_{\Omega}$ are chosen sufficiently small. The reduced system composed of (56) and (57) is still (44). By picking $\boldsymbol{\xi}_{\Theta}\left(0\right)=0$, one has $\boldsymbol{e}_{\Theta}\left(0\right)=\bar{e}_{\Theta}$, so $\|{\boldsymbol{e}}_{\Theta}-\boldsymbol{\bar{e}}_{\Theta}\|=\mathcal{O}\left(\epsilon_{2}\right)$ will uniformly hold, where $\epsilon_{2}=\max\left\{\epsilon_{1}\text{, }\frac{1}{\omega_{p}}\text{, }\frac{1}{\omega_{q}}\text{, }\frac{1}{\omega_{r}}\right\}$. In addition, $\lim_{t\to\infty}\boldsymbol{\bar{e}}_{\Theta}\to 0$, so $\boldsymbol{e}_{\Theta}$ will be ultimately bounded.

Acknowledgments

This research work presented in this paper is supported by Natural Science and Engineering Research Council of Canada (NSERC) Discovery Grant 227674.