Flight Control Design Based on Nonlinear Dynamic Inversion for Post-Stall Maneuvers

1. Introduction

The ability to perform transient maneuvers at low speeds and high angle of attack (AOA) is an important advantage in air combat [1,2,3]. However, for an aircraft that does not have the maneuverability in the post-stall region, once it approaches the stall AOA, the flow over the airfoil separates and the performance of the aerodynamic control surfaces drops sharply. On the other hand, the nonlinear effect induced by the inertial coupling moment [4,5,6] is particularly significant at high AOA. The first problem has been solved by using thrust vectoring control (TVC) to compensate the loss of aerodynamic control at high AOA [7,8]. The use of TVC technology can enhance the maneuverability and stability of the aircraft in stall conditions [2,9,10,11], and help increase the maneuverability of close air combat [12]. In addition, TVC can reduce the use of traditional control surfaces, such as horizontal tails and vertical tails, and realize the design of tailless aircraft [13,14].

The second problem raised by high AOA flight is the nonlinear effects, which bring great difficulties to the design of linear flight control. Linear controllers designed for a certain trim condition are often insufficient to handle a wider range of flight operations. Although with the help of robust controllers and gain scheduling controllers, the uncertain dynamics caused by nonlinear effects inevitably degrade the performance of linear flight control laws [15,16]. Besides linear control approaches, nonlinear control technologies have been developed in recent decades to design flight control systems, such as feedback linearization [17], backstepping control [18], sliding mode control [19], and nonlinear adaptive control [20]. Among them, nonlinear dynamic inversion (NDI) is one of the most widely used nonlinear control methods [21], especially in military applications. For example, the development of X-35, the predecessor of F-35, used NDI control law [22]. The main advantage of the NDI control method is that it can handle a wide range of flight operations, such as flying at high AOA [23,24], and directly incorporate nonlinear dynamics into the design of the controller [25,26]. In addition, NDI control can be applied to different aircraft, which can reduce the development time and cost of the control system.

The NDI control is a model-based control design, which employs an input-affine nonlinear model to describe the flight motion involving large attitude maneuvers and small control surface deflections by retaining the nonlinear characteristics of the flight dynamics, while allowing the control surface deflection to appear linearly in the state equations. A single NDI control law derived from the input-affine nonlinear model can cover the entire flight envelope without the need for gain scheduling. However, the existence of an exact NDI controller has its prerequisite that the number of control surfaces must be greater than the total number of state variables. To solve the problem of insufficient number of control surfaces, this paper divides the nine state variables into three loops based on their different time scales so that the angular-rate states constitute the inner loop, the attitude states constitute the outer loop, and the flight-path states constitute the outermost loop. Since a loop only involves three state variables, the corresponding NDI controller for each loop can be realized exactly by using three control inputs to achieve the goal of perfect tracking.

The proposed NDI flight control is validated numerically for an experimental F-16 installed with TVC in the post-stall flight. Since the existing F-16 only has aerodynamic control surfaces, when the AOA exceeds 35 degrees, it enters the post-stall region and loses its flight maneuverability. Therefore, although there have been many discussions in the literature about applying NDI control to F-16, they are all limited to flights below the stall AOA [23,27,28,29]. Here, we retain all the original design of the F-16 (including the engine), except for replacing its fixed nozzle by a thrust vectoring nozzle. What we want to test is whether the experimental F-16 can fly with high maneuverability in the post-stall region when the total thrust remains unchanged and only the direction of the thrust can be changed. The numerical simulations verify that the experimental F-16 under the NDI control can effectively track high AOA command for post-stall flight and successfully demonstrate the Herbst maneuver by performing the heading reversal with $70°$ AOA.

To facilitate high AOA turns, we first establish the equations of flight motion based on the wind-axis coordinates in Section 2. The obtained equations are further rewritten into an input-affine nonlinear form in Section 3 and separated into a three-layer loop structure to facilitate the independent design of the NDI controller for each layer. Section 4 designs the NDI controller for the angular-rate inner loop. The purpose is to determine the deflection angle of the control surfaces so that the aircraft's angular rate can track the command given by the outer loop. Section 5 designs the NDI controller for the attitude outer loop, whose purpose is to determine the angular-rate command so that the aircraft’s attitude tracks the command given by the outermost loop. In Section 6, the NDI controller developed previously is applied to an experimental F-16 to verify its high-AOA maneuverability in the post-stall region. The conclusion section discusses the modeling errors that may be induced by the input-affine nonlinear model with separable loop structure and proposes quantitative index for measuring the error and feasible methods to reduce the error.

2. Equations of Motion in Wind-Axis Coordinates

Many post-stall maneuvers require the aircraft to roll around the wind axis rather than the body axis. This is because when the aircraft rolls around the body axis, it will cause the interchange between the AOA $\alpha$ and the angle of sideslip $\beta$, as shown in Figure 1, where we can see that when the aircraft rotates 90 degrees around its body axis, the aircraft's AOA changes to the sideslip angle, and vice versa. This interchange between $\alpha$ and $\beta$ does not cause problem at small AOA, but at high AOA it would be unacceptable because the induced large sideslip angle would greatly degrade the rolling performance of the aircraft. To eliminate the kinematic coupling between $\alpha$ and $\beta$, an aircraft operating at high AOA is usually designed to roll around the wind axis and its flight motion is better described by the wind-axis coordinates than the body-axis coordinates.

To describe the attitude of the body-axis frame relative to the local horizontal frame, it is common to employ three Euler angles $(\varphi ,\theta ,\psi )$ with $\varphi$ being the rolling angle around the body axis. Instead of three Euler angles, here we use five angles to describe the aircraft’s attitude, namely AOA $\alpha$, sideslip angle $\beta$, bank angle $\mu$ (rotation angle around the wind axis), flight path angle $\gamma$ and heading angle $\chi$, of which the first three are used to represent the attitude of the aircraft’s body axis relative to the wind axis and the last two represent the attitude of the wind axis relative to the local horizontal coordinates.

The relative orientation between the three coordinate systems is illustrated in Figure 2. The origin of the local horizontal coordinates $\left({\widehat{x}}_E,{\widehat{y}}_E,{\widehat{z}}_E\right)$ is located on the earth's surface (ignoring the earth's rotation), where the ${\widehat{x}}_E$ axis points to the local north, the ${\widehat{y}}_E$ axis points to the local east, and the ${\widehat{y}}_E$ axis points to the center of the earth. The origin of the body-axis coordinates $\left({\widehat{x}}_b,{\widehat{y}}_b,{\widehat{z}}_b\right)$ is located at the center of mass of the aircraft, where the ${\widehat{x}}_b$ axis points to the direction of the nose along the symmetry line of the aircraft, the ${\widehat{y}}_b$ axis is perpendicular to ${\widehat{x}}_b$ and points to the right direction of the aircraft, and the ${\widehat{z}}_b$ axis points in the direction of ${\widehat{x}}_b\times {\widehat{y}}_b$.

The wind axis ${\widehat{x}}_w$ points to the velocity direction of the aircraft, so the aircraft’s velocity can be expressed as ${\widehat{V}}_T=V_T{\widehat{x}}_w$, where $V_T$ is the magnitude of the velocity and the velocity direction ${\widehat{x}}_w$ relative to the local horizontal coordinates is defined by $(\gamma ,\chi )$, where $\gamma$ is the flight path angle between the wind axis ${\widehat{x}}_w$ and the local horizontal plane (see Figure 2a), and $\chi$ is the heading angle between the projection of the wind axis on the local horizontal plane and the local north ${\widehat{x}}_E$ (see Figure 2c). By rotating an angle χ on the local horizontal plane, and then rotating an angle $\gamma$ on the vertical plane, we can transfer the local horizontal coordinates $\left({\widehat{x}}_E,{\widehat{y}}_E,{\widehat{z}}_E\right)$ to the wind axis coordinates $\left({\widehat{x}}_w,{\widehat{y}}_w,{\widehat{z}}_w\right)$, as expressed by Eq. (A.1). The change rate of the flight-path states $(V_T,\chi ,\gamma )$ satisfies the following equations (referring to Appendix A for the derivation):

$${\dot{V}}_T={\displaystyle \frac{1}{m}}\left[\left(-D+\overline{Y}\sin \beta -mg\sin \gamma \right)+\left(T_x\cos \beta \cos \alpha +T_y\sin \beta +T_z\cos \beta \sin \alpha \right)\right],$$

(1)

$$\begin{array}{cc}\dot{\chi }=& {\displaystyle \frac{1}{m{V}_T\cos \gamma }}[\left(L\sin \mu +\overline{Y}\cos \mu \cos \beta \right)+T_x\left(\sin \mu \sin \alpha -\cos \mu \sin \beta \cos \alpha \right)\hfill \\ & +T_y\left(\cos \mu \cos \beta \right)-T_z(\cos \mu \sin \beta \sin \alpha +\sin \mu \cos \alpha )],\hfill \end{array}$$

(2)

$$\begin{array}{cc}\dot{\gamma }=& {\displaystyle \frac{1}{m{V}_T}}[\left(L\cos \mu -\overline{Y}\sin \mu \cos \beta -mg\cos \gamma \right)+T_x\left(\sin \mu \sin \beta \cos \alpha +\cos \mu \sin \alpha \right)\hfill \\ & -T_y\left(\sin \mu \cos \beta \right)+T_z(\sin \mu \sin \beta \sin \alpha -\cos \mu \cos \alpha )].\hfill \end{array}$$

(3)

The lift $L$, the drag $D$, and the lateral force $\overline{Y}$ are defined on the wind-axis coordinates and are related to the force components $(X,Y,Z)$ based on the body axis via the following relation,

$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\overline{q}S\left[\begin{array}{c}{C}_{XT}\left(\alpha ,\beta ,p,q,r,{\delta }_e,{\delta }_a,{\delta }_r\right)\\ C_{YT}\left(\alpha ,\beta ,p,q,r,{\delta }_e,{\delta }_a,{\delta }_r\right)\\ C_{ZT}\left(\alpha ,\beta ,p,q,r,{\delta }_e,{\delta }_a,{\delta }_r\right)\end{array}\right]=\left[\begin{array}{ccc}-\cos \alpha \cos \beta & -\cos \alpha \sin \beta & \sin \alpha \\ -\sin \beta & \cos \beta & 0\\ -\sin \alpha \cos \beta & -\sin \alpha \sin \beta & -\cos \alpha \end{array}\right]\left[\begin{array}{c}D\\ \overline{Y}\\ L\end{array}\right],$$

(4)

where $\overline{q}=\rho {V_T}^2/2$ is the dynamic pressure, $S$ is the area of the main wing, $\rho$ is the atmospheric density, and $C_{XT}$, $C_{YT}$, $C_{ZT}$ are the aerodynamic force coefficients of the three directions of the aircraft body axis, which can be read from the aero data tables of F-16 according to the deflection angles of the elevator ${\delta }_e$, the aileron ${\delta }_a$, and the rudder ${\delta }_r$ (see Appendix A for details). The aero data tables are based on the wind tunnel test results of F-16 at NASA Langley and Ames Research center. The database is valid for $-20°\le \alpha \le 90°$, $-30°\le \beta \le 30°$ and $V_T\le 06\mathrm{M}\mathrm{a}\mathrm{c}\mathrm{h}$, which covers the post-stall flight region considered here.

The thrust components ${\left[T_x,T_y,T_z\right]}^T$ in Eq. (1) – Eq. (3) are generated by the thrust vectoring nozzle:

$$\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]=\left[\begin{array}{c}T\cos {\delta }_z\cos {\delta }_y\\ T\sin {\delta }_y\\ -T\sin {\delta }_z\cos {\delta }_y\end{array}\right],$$

(5)

where $T$ is the total thrust of the aircraft, ${\delta }_y$ is the left-right deflection angle of the nozzle (defining it as positive to the left), and ${\delta }_z$ is the up-down deflection angle of the nozzle (defining downward deflection as positive), as shown by Figure 2b. When the thrust vector deflects to the left, a rightward thrust $T_y$ will be generated along the positive ${\widehat{y}}_b$ direction of the body axis. When the nozzle deflects downward, an upward thrust $T_z$ will be generated along the negative ${\widehat{z}}_b$, yielding a negative sign of $T_z$.

The wind-axis frame is an intermediate frame between the local horizontal frame and the aircraft body-axis frame. By a sequence of rotations, $\mu \to \beta \to \alpha$ as shown by Eq. (A2), the wind-axis coordinates $\left({\widehat{x}}_w,{\widehat{y}}_w,{\widehat{z}}_w\right)$ can be transferred into the body-axis coordinates $\left({\widehat{x}}_b,{\widehat{y}}_b,{\widehat{z}}_b\right)$. The governing equations for $(\alpha ,\beta ,\mu )$ can be established as follows (see Appendix A for the derivation):

$$\dot{\alpha }=q-\left(p\cos \alpha +r\sin \alpha \right)\tan \beta +{\displaystyle \frac{1}{m{V}_T\cos \beta }}\left[\left(-L+mg\cos \gamma \cos \mu \right)+\left(-T_x\sin \alpha +T_z\cos \alpha \right)\right],$$

(6)

$$\dot{\beta }=-r\cos \alpha +p\sin \alpha +{\displaystyle \frac{1}{m{V}_T}}\left[\left(\overline{Y}\cos \beta +mg\cos \gamma \sin \mu \right)+\left(-T_x\sin \beta \cos \alpha +T_y\cos \beta +T_z\sin \beta \sin \alpha \right)\right],$$

(7)

$$\begin{gathered} \dot{\mu }=\left(p\cos \alpha +r\sin \alpha \right)\sec \beta +{\displaystyle \frac{L}{m{V}_T}}\left(\tan \gamma \sin \mu +\tan \beta \right)+{\displaystyle \frac{\left(\overline{Y}+T_y\right)}{m{V}_T}}\tan \gamma \cos \mu \cos \beta -{\displaystyle \frac{g}{V_T}}\cos \gamma \cos \mu \tan \beta \\ +{\displaystyle \frac{\left(T_x\sin \alpha -T_z\cos \alpha \right)}{m{V}_T}}\left(\tan \gamma \sin \mu +\tan \beta \right)-{\displaystyle \frac{\left(T_x\cos \alpha +T_z\sin \alpha \right)}{m{V}_T}}\tan \gamma \cos \mu \sin \beta . \end{gathered}$$

(8)

It can be seen above that the attitude change rates $(\dot{\alpha },\dot{\beta },\dot{\mu })$ depend on the body-axis angular rates $(p,q,r)$, which in turn depend on the external moments $\mathbf{M}={\left[M_x,M_y,M_z\right]}^T$ as

$$\dot{p}=\left(C_{1}r+C_{2}p\right)q+C_3{M}_x+C_4{M}_z,$$

(9)

$$\dot{q}=C_{5}pr-C_6\left(p^2-r^2\right)+C_7{M}_y,$$

(10)

$$\dot{r}=\left(C_{8}p-C_{2}r\right)q+C_4{M}_x+C_9{M}_z,$$

(11)

where the constants $C_i$ are defined by the given moments of inertia of F-16 (see Appendix A). The external moment $\mathbf{M}$ experienced by the aircraft is composed of the aerodynamic moment ${\mathbf{M}}^{\left(A\right)}$ and the thrust moment ${\mathbf{M}}^{\left(T\right)}$. The aerodynamic moment $M^{\left(A\right)}$ is determined by the aerodynamic coefficients:

$${\mathbf{M}}^{\left(A\right)}=\left[\begin{array}{c}{M}_x^{\left(A\right)}\\ M_y^{\left(A\right)}\\ M_z^{\left(A\right)}\end{array}\right]=\left[\begin{array}{c}L\\ M\\ N\end{array}\right]=\overline{q}S\left[\begin{array}{c}b{C}_{lT}\left(\alpha ,\beta ,p,q,r,{\delta }_e,{\delta }_a,{\delta }_r\right)\\ \overline{C}{C}_{mT}\left(\alpha ,\beta ,p,q,r,{\delta }_e,{\delta }_a,{\delta }_r\right)\\ {bC}_{nT}\left(\alpha ,\beta ,p,q,r,{\delta }_e,{\delta }_a,{\delta }_r\right)\end{array}\right],$$

(12)

where $b$ is the span of the main wing, $\overline{C}$ is the average chord length, and $C_{lT}$, $C_{mT}$, and $C_{nT}$ are the aerodynamic moment coefficients in the three directions of the aircraft body axis. The thrust moment ${\mathbf{M}}^{\left(T\right)}$ is produced by the thrust about the center of gravity of the aircraft and can be expressed as

$${\mathbf{M}}^{\left(T\right)}=\left[\begin{array}{c}{M}_x^{\left(T\right)}\\ M_y^{\left(T\right)}\\ M_z^{\left(T\right)}\end{array}\right]=\left[\begin{array}{c}0\\ l_T{T}_z\\ -l_T{T}_y\end{array}\right]=\left[\begin{array}{c}0\\ {-l}_{T}T\sin {\delta }_z\cos {\delta }_y\\ -l_{T}T\sin {\delta }_y\end{array}\right],$$

(13)

where $l_T$ is the distance from the center of gravity to the nozzle, and $T_y$ and $T_z$ are given by Eq. (5).

The flight motion under the wind-axis coordinates is described by nine states, for which the flight-path states $(V_T,\chi ,\gamma )$ are given by Eq. (1) – Eq. (3), the attitude states $(\alpha ,\beta ,\mu )$ are given by Eq. (6) – Eq. (8), and the angular-rate state $(p,q,r)$ are given by Eq. (9) – Eq. (11). As shown in Figure 3, The output of the nine states is governed by the five control inputs ${\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$, which generate the aerodynamic force $F^{\left(A\right)}$ by Eq. (4), the aerodynamic moment $M^{\left(A\right)}$ by Eq. (12), the thrust force $F^{\left(T\right)}$ by Eq. (5), and the thrust moment $M^{\left(T\right)}$ by Eq. (13). By solving the nine equations of motion subjected to the applied forces and moments, we obtain the position, velocity and attitude of the aircraft relative to the local horizontal coordinates at any instant $t$. Compared with the body-axis based state variables $(\psi ,\theta ,\varphi ,u,v,w,p,q,r)$, the wind-axis based state variables $(V_T,\chi ,\gamma ,\alpha ,\beta ,\mu ,p,q,r)$ are more suitable for designing flight control law of aircraft maneuvering at high AOA.

3. NDI Flight Control Law Design

The equations of motion of the nine states $\mathbf{x}={\left[V_T,\chi ,\gamma ,\mu ,\alpha ,\beta ,p,q,r\right]}^T$ derived previously can be expressed as an input-affine nonlinear form:

$$\begin{array}{ccc}\dot{\mathbf{x}}=\mathbf{f}\left(\mathbf{x}\right)+\mathbf{g}\left(\mathbf{x}\right)\mathbf{u},& \mathbf{x}\left(t_0\right)={\mathbf{x}}_0,& \mathbf{y}=\mathbf{h}\left(\mathbf{x}\right)\end{array},$$

(14)

where $\mathbf{f}\left(\mathbf{x}\right)$ is the vector function with 9 elements and $\mathbf{g}\left(\mathbf{x}\right)$ is the matrix function with dimension $9\times 5$, and $\mathbf{y}=\mathbf{h}\left(\mathbf{x}\right)$ is the system’s output vector. The control input $\mathbf{u}={\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$ consists of the aerodynamic deflection angles $({\delta }_a,{\delta }_e,{\delta }_r)$ and the nozzle deflection angles $({\delta }_y,{\delta }_z)$. Due to the limitations of the practical deflection mechanism, the deflection angle of each control surface is small so that the aerodynamic/thrust force and moment can all be expanded as a linear function of $\mathbf{u}$; however, the small angle assumption is not valid for the state $\mathbf{x}$. This situation is reflected in the input-affine nonlinear form (14), where the nonlinear nature of the state $\mathbf{x}$ is retained, while the control input $\mathbf{u}$ appears linearly.

The NDI control architecture is shown in Figure 4. The essence of the NDI control design is to determine the control signal $\mathbf{u}$ through the dynamic inverse process to transform the nonlinear system into an equivalent linear system. After the process of feedback linearization, an auxiliary linear controller is then designed to reduce the tracking error between the output $\mathbf{y}$ and the command ${\mathbf{y}}_{\mathbf{c}}$.

The control input $\mathbf{u}$ that linearizes the nonlinear system (14) can be found as,

$$\mathbf{u}={\mathbf{g}}^{-1}\left(\mathbf{x}\right)\left(\mathbf{v}-\mathbf{f}\left(\mathbf{x}\right)\right).$$

(15)

If this control is used in Eq. (14), we obtain the desired result:

$$\dot{\mathbf{x}}=\mathbf{f}\left(\mathbf{x}\right)+\mathbf{g}\left(\mathbf{x}\right){\mathbf{g}}^{-1}\left(\mathbf{x}\right)\left(\mathbf{v}-\mathbf{f}\left(\mathbf{x}\right)\right)=\mathbf{v}.$$

(16)

Accordingly, the effect of the NDI control amounts to making the inner closed-loop system behave like a linear system as shown by Figure 4. After the feedback linearization is achieved, the new control input $\mathbf{v}$ can then be designed by any linear control method to drive the system output $\mathbf{y}$ to follow the target output ${\mathbf{y}}_{\mathbf{c}}$.

The existence of the ideal NDI controller (15) depends on the invertibility of the matrix function $\mathbf{g}\left(\mathbf{x}\right)$. Unfortunately, the inverse of the matrix $\mathbf{g}\left(\mathbf{x}\right)$ with dimension $9\times 5$ considered here does not exist. The mathematical reason is that the number of its rows (number of equations) is greater than the number of its columns (number of unknowns), and the physical reason is that there are 9 desired states in Eq. (14) to be tracked, but only 5 degrees of freedom of control are available. Apparently, the desired state rate $\mathbf{v}={\dot{\mathbf{x}}}_{\mathrm{d}\mathrm{e}\mathrm{s}}$ can only be tracked partially and the best control we can design is to minimize the tracking error by replacing the inverse ${\mathbf{g}}^{-1}$ in Eq. (15) with its pseudo inverse,

$${\mathbf{g}}^{\#}={\left({\mathbf{g}}^T\mathbf{g}\right)}^{-1}{\mathbf{g}}^T,$$

(17)

where ${\mathbf{g}}^{\#}$ is the solution minimizing the distance ${‖\mathbf{g}{\mathbf{g}}^{\#}-{\mathbf{I}}_9‖}_2$, where ${\mathbf{I}}_9$ is the $9\times 9$ identity matrix.

For flight control applications, the pseudo-inverse approach is not the only way to solve the NDI control problem. Although the mathematical relationship in Eq. (14) indicates that the control $\mathbf{u}$ affects the 9 states simultaneously, the degree of influence on each state is very different. Due to the inherent dynamic characteristics of the aircraft, the first response to the change in the control surface deflections is the aircraft's body-axis angular rates ${[p,q,r]}^T$. When the aircraft’s angular rate changes, its attitude ${[\alpha ,\beta ,\mu ]}^T$ will change accordingly. The heading angle $\chi$ and flight path angle $\gamma$ of the aircraft will not respond until the attitude of the aircraft changes. The sequence of response is that the control $\mathbf{u}$ affects ${[p,q,r]}^T$ first, then ${[p,q,r]}^T$ affects ${[\alpha ,\beta ,\mu ]}^T$, and finally ${[\alpha ,\beta ,\mu ]}^T$ affects $(V_T,\chi ,\gamma )$. Therefore, the response speed of the 9 states to the control $\mathbf{u}$ is that the angular-rate state ${[p,q,r]}^T$ has the fastest response, the attitude state ${[\alpha ,\beta ,\mu ]}^T$ is slower, and the flight-path state $(V_T,\chi ,\gamma )$ is the slowest.

The concept of NDI design reverses the above sequence of dynamic response. It starts from the desired flight-path state $(V_T,\chi ,\gamma )$ to find reversely the required attitude state ${[\alpha ,\beta ,\mu ]}^T$, from the desired attitude state to find reversely the required angular-rate state ${[p,q,r]}^T$, and finally from the desired angular-rate state to find reversely the required control surface deflection $\mathbf{u}$.

Based on the fact that the three sets of states have three different time scales, the entire flight control loop can be separated into three layers of loop, as shown in Figure 5, where the inner loop is the angular-rate control loop regarding the control of the fast states ${[p,q,r]}^T$, the outer loop is the attitude control loop regarding the control of the slow states ${[\alpha ,\beta ,\mu ]}^T$, and the outermost loop is the flight-path control loop regarding the control of the slowest state $(V_T,\chi ,\gamma )$. The current article focuses on the NDI control design for the inner loop and the outer loop.

Since the nine states are in three different loops according to the difference in their time scales, independent NDI controller can be designed for each loop, which handles only part of the nine states and makes the pseudo-inverse operation unnecessary. In the inner loop of Figure 5, only three angular-rate states ${[p,q,r]}^T$ are concerned and the inner NDI control problem is to find reversely the control surface deflection $\mathbf{u}$ to achieve the desired angular-rate response ${[p_c,q_c,r_c]}^T$. Similarly, in the outer loop of Figure 5, only three attitude states ${[\alpha ,\beta ,\mu ]}^T$ are concerned and the outer NDI control problem is to find reversely the required angular-rate input ${[p_c,q_c,r_c]}^T$ to achieve the desired attitude response ${[{\alpha }_c,{\beta }_c,{\mu }_c]}^T$.

4. Angular-Rate NDI Control

The NDI control design in the inner loop is to determine the control surface deflection $\mathbf{u}$ so that the response of aircraft’s angular rates ${[p,q,r]}^T$ follows the command ${[p_c,q_c,r_c]}^T$, as shown in Figure 6.

To design the inner-loop NDI controller $\mathbf{u}$, the fully nonlinear equations (9) – (11) must be transformed into an affine nonlinear form, wherein the control input $\mathbf{u}={\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$ appears linearly. The main task involved in this transformation is to expand the aerodynamic/thrust moments in Eq. (9) – Eq. (11) as a linear function of $\mathbf{u}$, which yields the following affine nonlinear form:

$$\left[\begin{array}{c}\dot{p}\\ \dot{q}\\ \dot{r}\end{array}\right]={\mathbf{f}}_1\left(\mathbf{x}\right)+{\mathbf{g}}_1\left(\mathbf{x}\right)\mathbf{u}=\left[\begin{array}{c}{f}_p\left(\mathbf{x}\right)\\ f_q\left(\mathbf{x}\right)\\ f_r\left(\mathbf{x}\right)\end{array}\right]+\left[\begin{array}{ccccc}{g}_{p_{{\delta }_a}}& 0& g_{p_{{\delta }_r}}& g_{p_{{\delta }_y}}& 0\\ 0& g_{q_{{\delta }_e}}& 0& 0& g_{q_{{\delta }_z}}\\ g_{r_{{\delta }_a}}& 0& g_{r_{{\delta }_r}}& g_{r_{{\delta }_y}}& 0\end{array}\right]\left[\begin{array}{c}{\delta }_a\\ {\delta }_e\\ {\delta }_r\\ {\delta }_y\\ {\delta }_z\end{array}\right],$$

(18)

where explicit expressions for the elements in the vector ${\mathbf{f}}_1\left(\mathbf{x}\right)$ and the matrix ${\mathbf{g}}_1\left(\mathbf{x}\right)$ are given by Eq. (A17) in Appendix A. Although the three angular rates ${[p,q,r]}^T$ are the dominant states in Eq. (18), all the nine states in $\mathbf{x}$ are involved in the expressions of ${\mathbf{f}}_1\left(\mathbf{x}\right)$ and ${\mathbf{g}}_1\left(\mathbf{x}\right)$; however, as measured under the time scale of ${[p,q,r]}^T$, the remaining six states in $\mathbf{x}$, i.e., ${\left[V_T,\chi ,\gamma ,\alpha ,\beta ,\mu \right]}^T$, can be treated as constants in the inner loop.

The control targets in the inner loop are the three desired angular rates ${[p_c,q_c,r_c]}^T$, but there are 5 independent control inputs $\mathbf{u}={\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$ can be manipulated. Consequently, many possible combinations of the five control inputs can be employed to achieve the desired angular rates, and an appropriate control allocation is required to find the best combination. Mathematically, it means that the matrix ${\mathbf{g}}_1\left(\mathbf{x}\right)$ in Eq. (18) is (right) invertible, i.e., ${\mathbf{g}}_1{\mathbf{g}}_1^{-1}={\mathbf{I}}_3$, but its inverse ${\mathbf{g}}_1^{-1}$ is not unique. Due to the non-uniqueness of ${\mathbf{g}}_1^{-1}$, all the possible solutions of $\mathbf{u}={\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$ can be expressed as

$$\mathbf{u}=\mathbf{N}\overline{\mathbf{u}}⟶\left[\begin{array}{c}{\delta }_a\\ {\delta }_e\\ {\delta }_r\\ {\delta }_y\\ {\delta }_z\end{array}\right]=\left[\begin{array}{ccc}{N}_{aP}& N_{aQ}& N_{aR}\\ N_{eP}& N_{eQ}& N_{eR}\\ N_{rP}& N_{rQ}& N_{rR}\\ N_{yP}& N_{yQ}& N_{yR}\\ N_{zP}& N_{zQ}& N_{zR}\end{array}\right]\left[\begin{array}{c}{\delta }_P\\ {\delta }_Q\\ {\delta }_R\end{array}\right],$$

(19)

where $\mathbf{N}$ is the control allocation matrix, whose elements can be arbitrarily assigned as long as ${\mathbf{g}}_1\mathbf{N}$ is nonsingular. The effective control $\overline{\mathbf{u}}$ consists of three effective control surface deflections, ${\delta }_P$, ${\delta }_Q$, and ${\delta }_R$, required for controlling the rolling rate $p$, the pitching rate $q$, and the yawing $r$, respectively.

The elements of matrix $\mathbf{N}$ indicate how each control surface deflection is allocated to the control of the three angular rates. For example, the first equation in Eq. (19) reads,

$${\delta }_a=N_{aP}{\delta }_P+N_{aQ}{\delta }_Q+N_{aR}{\delta }_R,$$

(20)

which means that $N_{aP}$ percentage of the total aileron deflection ${\delta }_a$ is allocated to the rolling rate control, $N_{aQ}$ percentage to the pitching rate control, and $N_{aR}$ percentage to the rolling rate control. It is worth noting that ${\delta }_y$ and ${\delta }_z$ are the deflection angles of the thrust vectoring nozzle and the last two rows of $\mathbf{N}$ tell us how the total thrust is allocated to the control of the three angular rates. Regarding the setting of the elements of matrix $\mathbf{N}$, we must also consider the angle and rate limitations of the deflection mechanism imposed on different control surfaces. There are many control allocation methods, which can be used to determine the matrix $\mathbf{N}$, such as daisy chain, explicit ganging, pseudo control or control selector, and pseudo inverse [30,31].

The effective control $\overline{\mathbf{u}}$ is to be determined so that the aircraft’s angular acceleration ${[\dot{p},\dot{q},\dot{r}]}^T$ converges to the command ${[{\dot{p}}_c,{\dot{q}}_c,{\dot{r}}_c]}^T$. The ideal effective control $\overline{\mathbf{u}}$ can be derived by letting ${[\dot{p},\dot{q},\dot{r}]}^T={[{\dot{p}}_c,{\dot{q}}_c,{\dot{r}}_c]}^T$ and $\mathbf{u}=\mathbf{N}\overline{\mathbf{u}}$ in Eq. (18) to obtain

$$\overline{\mathbf{u}}={\left({\mathbf{g}}_1\mathbf{N}\right)}^{-1}\left(\left[\begin{array}{c}{\dot{p}}_c\\ {\dot{q}}_c\\ {\dot{r}}_c\end{array}\right]-\left[\begin{array}{c}{f}_p\left(\mathbf{x}\right)\\ f_q\left(\mathbf{x}\right)\\ f_r\left(\mathbf{x}\right)\end{array}\right]\right)⟶\mathbf{u}=\mathbf{N}\overline{\mathbf{u}}=\mathbf{N}{\left({\mathbf{g}}_1\mathbf{N}\right)}^{-1}\left(\left[\begin{array}{c}{\dot{p}}_c\\ {\dot{q}}_c\\ {\dot{r}}_c\end{array}\right]-\left[\begin{array}{c}{f}_p\left(\mathbf{x}\right)\\ f_q\left(\mathbf{x}\right)\\ f_r\left(\mathbf{x}\right)\end{array}\right]\right),$$

(21)

where ${\mathbf{g}}_1\mathbf{N}$ is a nonsingular $3\times 3$ matrix and its inverse is defined uniquely. With the effective control $\overline{\mathbf{u}}$ given by Eq. (21), the actual control surface deflections $\mathbf{u}={\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$ are then determined by passing $\overline{\mathbf{u}}$ via the allocation matrix $\mathbf{N}$ as shown in Figure 6.

With NDI controller $\mathbf{u}$ given by Eq. (21), the angular-rate dynamics (18) behaves like a linear system, for which a linear controller is to be designed to drive the output of the inner loop to the desired output ${[p_c,q_c,r_c]}^T$. Referring to Figure 6, the angular acceleration command ${[{\dot{p}}_c,{\dot{q}}_c,{\dot{r}}_c]}^T$ is adjusted by the linear controller to drive the tracking error $\mathbf{e}={[p_c,q_c,r_c]}^T-{[p,q,r]}^T$ to zero. Here a simple PI controller is used to generate the command ${[{\dot{p}}_c,{\dot{q}}_c,{\dot{r}}_c]}^T$ according to the tracking error $\mathbf{e}$:

$$\begin{array}{c}{\dot{p}}_c={\omega }_p\left(p_c-p\right)+\int {\mathsf{\Omega }}_p\left(p_c-p\right)\mathrm{d}t,\\ \begin{array}{c}{\dot{q}}_c={\omega }_q\left(q_c-q\right)+\int {\mathsf{\Omega }}_q\left(q_c-q\right)\mathrm{d}t,\\ {\dot{r}}_c={\omega }_r\left(r_c-r\right)+\int {\mathsf{\Omega }}_r\left(r_c-r\right)\mathrm{d}t,\end{array}\end{array}$$

(22)

where the desired angular rates ${[p_c,q_c,r_c]}^T$ come from the output of the outer loop, ${[p,q,r]}^T$ are the actual angular rates measured by the angular-rate sensors, ${\left[{\omega }_p{\omega }_q{\omega }_r\right]}^T$ are the adjustable proportional gains, and ${\left[{\mathsf{\Omega }}_p{\mathsf{\Omega }}_q{\mathsf{\Omega }}_r\right]}^T$ are the adjustable integral gains.

5. Attitude NDI Control

The dominant states in the outer NDI control loop are the slow attitude states ${[\mu ,\alpha ,\beta ]}^T$. Under the time scale of the attitude states, the fast angular-rate states ${[p,q,r]}^T$ have converged to their steady states ${[p_c,q_c,r_c]}^T$ by the inner NDI control, while the change of the flight-path states $(V_T,\chi ,\gamma )$ in the outermost loop are so slow that they can be considered as constants in the outer loop. As mentioned earlier, the dominant response of the angular-rate dynamics ${[p,q,r]}^T$ is caused by the control surface deflections $\mathbf{u}$, while the dominant response of the attitude dynamics ${[\mu ,\alpha ,\beta ]}^T$ is caused by the angular-rate input ${[p,q,r]}^T$, rather than by the control surface deflections $\mathbf{u}$. With the above understanding, the response of the attitude dynamics to the body-axis angular rates can be derived from Eq. (6) – Eq. (8) by letting ${[p,q,r]}^T={[p_c,q_c,r_c]}^T$ as the steady angular-rate states and ignoring the secondary influences of the control surface deflections $\mathbf{u}$ on the attitude states. The result can be expressed as the following affine nonlinear form:

$$\left[\begin{array}{c}\dot{\mu }\\ \dot{\alpha }\\ \dot{\beta }\end{array}\right]={\mathbf{f}}_2\left({\mathbf{x}}_2\right)+{\mathbf{g}}_2\left({\mathbf{x}}_2\right){\mathbf{u}}_2=\left[\begin{array}{c}{f}_{\mu }\left({\mathbf{x}}_2\right)\\ f_{\alpha }\left({\mathbf{x}}_2\right)\\ f_{\beta }\left({\mathbf{x}}_2\right)\end{array}\right]+\left[\begin{array}{ccc}\cos \alpha \sec \beta & 0& \sin \alpha \sec \beta \\ -\cos \alpha \tan \beta & 1& -\sin \alpha \tan \beta \\ \sin \alpha & 0& -\cos \alpha \end{array}\right]\left[\begin{array}{c}{p}_c\\ q_c\\ r_c\end{array}\right],$$

(23)

where the reduced state ${\mathbf{x}}_2$ are defined as ${\mathbf{x}}_2={\left[V_T,\chi ,\gamma ,\alpha ,\beta ,\mu \right]}^T$ and ${\mathbf{u}}_2={[p_c,q_c,r_c]}^T$ is treated as the control input. The three nonlinear functions in Eq. (23) are given by,

$$f_{\alpha }\left({\mathbf{x}}_2\right)={\displaystyle \frac{1}{m{V}_T\cos \beta }}\left[\left(-L+mg\cos \gamma \cos \mu \right)-T\sin \alpha \right],$$

$$\begin{array}{c}{f}_{\beta }\left({\mathbf{x}}_2\right)={\displaystyle \frac{1}{m{V}_T}}\left[\overline{Y}\cos \beta +mg\cos \gamma \sin \mu -T\sin \beta \cos \alpha \right],\end{array}$$

$$\begin{array}{cc}{f}_{\mu }\left({\mathbf{x}}_2\right)=& {\displaystyle \frac{L}{m{V}_T}}\left(\tan \gamma \sin \mu +\tan \beta \right)+{\displaystyle \frac{1}{m{V}_T}}\left(\overline{Y}\tan \gamma \cos \mu \cos \beta -mg\cos \gamma \cos \mu \tan \beta \right)\hfill \\ & +{\displaystyle \frac{T}{m{V}_T}}\left(\sin \alpha \tan \gamma \sin \mu +\sin \alpha \tan \beta -\cos \alpha \tan \gamma \cos \mu \sin \beta \right).\hfill \end{array}$$

In Eq. (23), the angular-rate state ${\mathbf{u}}_2={[p_c,q_c,r_c]}^T$ is the control input to the system and the attitude state ${[\mu ,\alpha ,\beta ]}^T$ is the response to the input. The concept of NDI reverses the sequence of dynamic response by requiring the desired response and reversely finding the necessary input. As shown in Figure 7, the control input ${\mathbf{u}}_2={[p_c,q_c,r_c]}^T$ is to be determined reversely by the NDI design so that the change rate of the attitude state ${[\dot{\mu },\dot{\alpha },\dot{\beta }]}^T$ tracks the desired response ${[{\dot{\mu }}_c,{\dot{\alpha }}_c,{\dot{\beta }}_c]}^T$.

By letting ${[\dot{\mu },\dot{\alpha },\dot{\beta }]}^T=[{\dot{\mu }}_c,{\dot{\alpha }}_c,{\dot{\beta }}_c]$ in Eq. (23), the ideal control input ${\mathbf{u}}_2$ can be determined as

$${\mathbf{u}}_2\triangleq \left[\begin{array}{c}{p}_c\\ q_c\\ r_c\end{array}\right]={\mathbf{g}}_2^{-1}\left({\mathbf{x}}_2\right)\left(\left[\begin{array}{c}{\dot{\mu }}_c\\ {\dot{\alpha }}_c\\ {\dot{\beta }}_c\end{array}\right]-\left[\begin{array}{c}{f}_{\mu }\left({\mathbf{x}}_2\right)\\ f_{\alpha }\left({\mathbf{x}}_2\right)\\ f_{\beta }\left({\mathbf{x}}_2\right)\end{array}\right]\right),$$

(24)

where ${\mathbf{g}}_2\left({\mathbf{x}}_2\right)$ is the $3\times 3$ nonsingular matrix defined in Eq. (23) and its inverse can be expressed analytically as

$${\mathbf{g}}_2^{-1}={\left[\begin{array}{ccc}\cos \alpha \sec \beta & 0& \sin \alpha \sec \beta \\ -\cos \alpha \tan \beta & 1& -\sin \alpha \tan \beta \\ \sin \alpha & 0& -\cos \alpha \end{array}\right]}^{-1}=\left[\begin{array}{ccc}\cos \alpha \cos \beta & 0& \sin \alpha \\ \sin \beta & 1& 0\\ \sin \alpha \cos \beta & 0& -\cos \alpha \end{array}\right]$$

(25)

Since ${\mathbf{u}}_2$ is determined uniquely, there is no control allocation problem in the outer NDI loop. After the outer loop is linearized by the NDI controller ${\mathbf{u}}_2$, a PI controller is incorporated into the outer loop to attain the desired aircraft’s attitude ${[{\mu }_c,{\alpha }_c,{\beta }_c]}^T$ according to the attitude tracking error, as shown in Figure 7. Like the inner-loop PI structure in Eq. (22), proportional gains ${[{\omega }_{\mu },{\omega }_{\alpha },{\omega }_{\beta }]}^T$ and integral gains ${[{\mathsf{\Omega }}_{\mu },{\mathsf{\Omega }}_{\alpha },{\mathsf{\Omega }}_{\beta }]}^T$ are used in the outer-loop attitude controller as follows:

$$\begin{array}{c}{\dot{\mu }}_c={\omega }_{\mu }\left({\mu }_c-\mu \right)+\int {\mathsf{\Omega }}_{\mu }\left({\mu }_c-\mu \right)dt,\\ \begin{array}{c}{\dot{\alpha }}_c={\omega }_{\alpha }\left({\alpha }_c-\alpha \right)+\int {\mathsf{\Omega }}_{\alpha }\left({\alpha }_c-\mu \right)dt,\\ {\dot{\beta }}_c={\omega }_{\beta }\left({\beta }_c-\beta \right)+\int {\mathsf{\Omega }}_{\beta }\left({\beta }_c-\mu \right)dt,\end{array}\end{array}$$

(26)

where ${[{\mu }_c,{\alpha }_c,{\beta }_c]}^T$ are the desired attitudes and ${[\mu ,\alpha ,\beta ]}^T$ are the aircraft’s actual attitudes provided by the attitude sensors. The magnitude of the proportional gains ${[{\omega }_{\mu },{\omega }_{\alpha },{\omega }_{\beta }]}^T$ reflects the bandwidth of the attitude control in the outer loop and is much smaller than proportional gains ${[{\omega }_p,{\omega }_q,{\omega }_r]}^T$, which reflects the bandwidth of the angular-rate control in the inner loop.

6. Verification of NDI Control for an Experimental F-16 in Post-Stall Flight

The experimental F-16 retains all the original design of the existing F-16 (including the engine), except for replacing its fixed nozzle by a thrust vectoring nozzle. The affine nonlinear model (14) can be used to describe the nonlinear dynamics of the experimental F-16 at high angles of attack by noting that the small angular changes in the thrust direction of the experimental F-16 (within a range of $\pm 15$ degrees) meet the basic assumption of the affine nonlinear model within which the control input $\mathbf{u}$ must appear linearly.

The flight performance of the experimental F-16 with and without TVC can be compared easily in the framework of NDI control by using two different control allocation matrices:

$${\mathbf{N}}_{\mathrm{o}\mathrm{f}\mathrm{f}}=\left[\begin{array}{ccc}0.75& 0& 0.25\\ 0& 1& 0\\ 0.25& 0& 0.75\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],\hspace{1em}\hspace{1em}{\mathbf{N}}_{\mathrm{o}\mathrm{n}}=\left[\begin{array}{ccc}0.75& 0& 0.25\\ 0& 1& 0\\ 0.25& 0& 0.75\\ 0.25& 0& 0.5\\ 0& 0.5& 0\end{array}\right].$$

(27)

Clearly, the main difference in the above two matrices is their last two rows. According to the definition of the allocation matrix in Eq. (19), the last two rows involve the allocation of the ${\delta }_y$ and ${\delta }_z$ deflections of the thrust to the three body-axis angular rates ${[p,q,r]}^T$. The elements in the last two rows of the ${\mathbf{N}}_{\mathrm{o}\mathrm{f}\mathrm{f}}$ matrix are all zero, which means that the TVC is not activated, so that the thrust does not contribute to the rotation about the three axes. Matrix ${\mathbf{N}}_{\mathrm{o}\mathrm{n}}$ corresponds to the situation when TVC is activated, in which the left-right deflection (${\delta }_y$) of the thrust contributes 25% to the rolling rate $p$ and 50% to the yawing rate $r$, and the up-down deflection (${\delta }_z$) of the thrust contributes 50% to the pitching rate $q$. The setting of the control allocation matrix $\mathbf{N}$ is not limited to the form expressed by Eq. (27). We can set different matrix $\mathbf{N}$ according to different flight mission.

The 6D dynamics of the experimental F-16 is modelled by Eq. (1) – Eq. (3) and Eq. (6) – Eq. (11) with five control inputs $\mathbf{u}={\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$ and nine state outputs $\mathbf{x}={\left[V_T,\chi ,\gamma ,\mu ,\alpha ,\beta ,p,q,r\right]}^T$. With the measured state outputs $\mathbf{x}$, the control inputs $\mathbf{u}$ are determined by a sequence of dynamic inverse processes. Here we assume that the attitude command ${[{\mu }_c,{\alpha }_c,{\beta }_c]}^T$ has been designed by an outermost flight-path loop not mentioned here, and we will start with ${[{\mu }_c,{\alpha }_c,{\beta }_c]}^T$ via two dynamic inverse processes to determine the five control surface deflections. At every time step, the 6D simulation contains the following iterative procedures:

• Input the current measurement of the nine states $\mathbf{x}={\left[V_T,\chi ,\gamma ,\mu ,\alpha ,\beta ,p,q,r\right]}^T$.
• Compute the current values of ${\mathbf{f}}_1\left(\mathbf{x}\right)$ and ${\mathbf{g}}_1\left(\mathbf{x}\right)$ by Eq. (18), and the current values of ${\mathbf{f}}_2\left({\mathbf{x}}_2\right)$ and ${\mathbf{g}}_2\left({\mathbf{x}}_2\right)$ by Eq. (23) with the new measurement $\mathbf{x}$.
• Compute the desired change rate of the attitude ${[{\dot{\mu }}_c,{\dot{\alpha }}_c,{\dot{\beta }}_c]}^T$ by Eq. (26) with given command ${[{\mu }_c,{\alpha }_c,{\beta }_c]}^T$.
• Compute the desired body-axis angular rates ${[p_c,q_c,r_c]}^T$ by Eq. (24) with given ${[{\dot{\mu }}_c,{\dot{\alpha }}_c,{\dot{\beta }}_c]}^T$.
• Compute the desired angular acceleration ${[{\dot{p}}_c,{\dot{q}}_c,{\dot{r}}_c]}^T$ by Eq. (22) with given ${[p_c,q_c,r_c]}^T$.
• Compute the required control surface deflections $\mathbf{u}={\left[{\delta }_a,{\delta }_e,{\delta }_r,{\delta }_y,{\delta }_z\right]}^T$ by Eq. (21) with given ${[{\dot{p}}_c,{\dot{q}}_c,{\dot{r}}_c]}^T$ and control allocation matrix $\mathbf{N}$.
• Solve the 6D model with the new control surface deflections $\mathbf{u}$ for the new value of the states $\mathbf{x}$ at next time step.

The above iterative procedures are applied to the following three cases of NDI flight simulation to verify the maneuverability of the experimental F-16 in the high AOA region.

6.1. NDI Flight Control of the F-16 without TVC

The first flight simulation is to verify that the traditional F-16 without TVC can fly near its stall AOA under the NDI flight control. The advantage of NDI flight control is that a single NDI controller can cover the entire flight envelope, without the need to divide the flight envelope and linearize the flight dynamics for gain scheduling. The parameters used in the design of NDI control are listed below.

●

The thrust setting is $5000\mathrm{l}\mathrm{b}\mathrm{f}$ with initial speed $V_T=500\mathrm{f}\mathrm{t}/\mathrm{s}$, initial height $h_0=15000\mathrm{f}\mathrm{t}$, and initial angular rates $p_0=q_0=r_0=0\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$.

●

In the attitude loop, the proportional gains and the integral gains are set to be ${\omega }_{\mu }={\omega }_{\alpha }={\omega }_{\beta }=2\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$ and ${\mathsf{\Omega }}_{\mu }={\mathsf{\Omega }}_{\alpha }={\mathsf{\Omega }}_{\beta }=1\mathrm{r}\mathrm{a}\mathrm{d}/{\mathrm{s}}^2$.

●

In the angular-rate loop, the proportional gains and the integral gains are set to be ${\omega }_{\mu }={\omega }_{\alpha }={\omega }_{\beta }=10\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$ and ${\mathsf{\Omega }}_{\mu }={\mathsf{\Omega }}_{\alpha }={\mathsf{\Omega }}_{\beta }=4\mathrm{r}\mathrm{a}\mathrm{d}/{\mathrm{s}}^2$.

●

The TVC function is deactivated in this case by choosing the control allocation matrix ${\mathbf{N}}_{\mathrm{o}\mathrm{f}\mathrm{f}}$ in Eq. (20).

●

The sideslip angle command ${\beta }_c$ is set to zero, and the AOA command ${\alpha }_c$ and the bank angle command ${\mu }_c$ are given, respectively, as

$${\alpha }_c=\left\{\begin{array}{cc}0°& ,t\le 2\hfill \\ 0°+\left(35°-0°\right)\times (t-2)/6& ,2<t\le 8\hfill \\ 35°-\left(35°-0°\right)\times (t-8)/7& ,8<t\le 15\hfill \\ 0°& ,15<t\le 30\hfill \end{array}\right.$$

(28)

$${\mu }_c=\left\{\begin{array}{cc}0°& ,t\le 2\hfill \\ 0°+\left(100°-0°\right)\times (t-2)/7& ,2<t\le 9\hfill \\ 100°-\left(100°-0°\right)\times (t-9)/6& ,9<t\le 15\hfill \\ 0°& ,15<t\le 30\hfill \end{array}\right.$$

(29)

The attitude command requires the F-16 to maintain level flight for the first 2 seconds, and then in the next 13 seconds, the AOA command ${\alpha }_c$ increases linearly from $0°$ to the stall angle $35°$ and then decreases to $0°$. Meanwhile, the bank angle command ${\mu }_c$ increases linearly from $0°$ to $100°$ and then decreases to $0°$. The attitude commands are represented by the black curves in Figure 8a–c, while the actual attitudes of the F-16 under the control of NDI are presented by the red curve. From Figure 8 we can see that the F-16 flies closely following the attitude commands. The mismatch between the attitude response and the attitude command mainly occurs between $t=9\mathrm{s}$ and $t=11\mathrm{s}$, within which the sideslip angle $\beta \left(t\right)$ deviates from zero about 3 degrees, and slight overshoots occur in the response of $\alpha \left(t\right)$ and $\mu \left(t\right)$.

Comparing Figure 8e, we find that the deviation of the attitudes between $t=9\mathrm{s}$ and $t=11\mathrm{s}$ is primarily caused by the large aileron deflection angle ${\delta }_a$ required by the NDI controller (black curve), which has far exceeded its saturation angle ($\pm 20°$). The flight maneuvers performed by this numerical test are very close to the performance ceiling of the experimental F-16 without TVC. Although it can barely follow the attitude commands to fly, the control surfaces have been deflected to their maximum angles.

6.2. NDI Flight Control beyond Stall AOA with and without TVC

We proceed to test the experimental F-16 to maneuver beyond the stall AOA with and without TVC. To conduct this flight maneuvers, we increase the maximum value of the AOA command in Eq. (28) from $35°$ to $40°$, while other settings remain unchanged. The simulation results are shown in Figure 9, where the black dotted curve is the input command, the blue chain curve is the aircraft’s response without TVC, and the red curve is the aircraft’s response with TVC. The attitude responses are shown in Figure 9a–c and the corresponding control surface deflections are shown in Figure 9d–f.

From the trend of the blue chain curves, we can see that tracking a $40°$ AOA command exceeds the performance limit of the F-16 without TVC. It appears that the response of the three attitudes began to diverge at $t=8\mathrm{s}$, when the F-16 reached its maximum AOA and entered the stall region, causing the aerodynamic surfaces to lose their control force and moment. Meantime, the deflection of the control surfaces was saturated at the maximum angle, indicating that the F-16 had no extra control force to use.

On the other hand, the red solid curves in Figure 9 show the performance of the experimental F-16 installed with TVC under the NDI control with the control allocation matrix ${\mathbf{N}}_{\mathrm{o}\mathrm{n}}$. From the time responses of the attitude, we can see that the high AOA attitude command is tracked very well with neglectable deviations, and the control surface deflections are safely within their allowable ranges. Through the integration of the aerodynamic control and the thrust control under the NDI framework, thrust control successfully replaces the role of aerodynamic control at high AOA, allowing the experimental F-16 to fly with high maneuverability in the post-stall region.

6.3. Herbst Maneuver with NDI Flight Control

The Herbst maneuver, also known as a J-turn, is an air combat maneuver that uses post-stall technology such as TVC and advanced flight controls to achieve high AOA. The Herbst maneuver allows an aircraft to quickly reverse direction by rolling around the wind axis at high AOA. Apparently, the Herbst maneuver is not demonstrable by the existing F-16. The third numerical verification is to illustrate the Herbst maneuver performed by the experimental F-16 equipped with TVC and NDI control law. The Herbst maneuver is mainly composed of three flight phases:

●

Phase 1: The aircraft’s AOA increases to $40°$ from the level flight to start the heading reversal maneuver.

●

Phase 2: The aircraft rotates around the wind-axis at high AOA near $70°$ to reverse the heading.

●

Phase 3: The aircraft stops rotating and lowers its AOA, returning to the level flight.

The typical Herbst maneuver is prescribed by the following attitude commands ${[{\mu }_c,{\alpha }_c,{\beta }_c]}^T$ [32]:

$${\alpha }_c=\left\{\begin{array}{cc}7°+\left(33°-0°\right)\times (t-0)/5& ,t\le 5\hfill \\ 40°+\left(70°-40°\right)\times (t-5)/2.8& ,5<t\le 7.8\hfill \\ 70°& ,7.8<t\le 8.2\hfill \\ 70°-\left(70°-50°\right)\times \left(t-8.2\right)/1.8& ,8.2<t\le 10\hfill \\ 50°-\left(50°-40.68°\right)\times \left(t-10\right)/2.5& ,10<t\le 12.5\hfill \\ 40.68°-\left(40.68°-40°\right)\times \left(t-12.5\right)/1.5& ,12.5<t\le 14\hfill \\ 40°-\left(40°-10°\right)\times \left(t-14\right)/6& ,14<t\le 20\hfill \end{array}\right.$$

(30)

$${\mu }_c=\left\{\begin{array}{cc}0°+\left(20.6°-0°\right)\times \left(t-0\right)/2& ,t\le 2\hfill \\ 20.6°+\left(-19.5°-20.6°\right)\times \left(t-2\right)/3& ,2<t\le 5\hfill \\ -19.5°+\left(-100°+19.5°\right)\times \left(t-5\right)/2& ,5<t\le 7\hfill \\ -100°+\left(-120°+100°\right)\times \left(t-7\right)/1.7& ,7<t\le 8.7\hfill \\ -120°-\left(100°-120°\right)\times \left(t-8.7\right)/1.7& ,8.7<t\le 10.4\hfill \\ -100°-\left(40.1°-100°\right)\times \left(t-10.4\right)/2& ,10.4<t\le 12.4\hfill \\ -40.1°-\left(20°-40.1°\right)\times \left(t-14\right)/6& ,12.4<t\le 15\hfill \\ -20°-\left(10°-20°\right)\times \left(t-15\right)/5& ,15<t\le 20\hfill \end{array}\right.$$

(31)

and ${\beta }_c=0$. Roughly speaking, the period $0\mathrm{s}\le t<5\mathrm{s}$ is the preparation period of the Herbst maneuver, which brings the aircraft to $40°$ AOA to start the heading reversal and belongs to phase 1. The period $5\mathrm{s}\le t<15\mathrm{s}$ is the working period of the Herbst maneuver, which performs the heading reversal and belongs to phase 2. The period $15\mathrm{s}\le t<20\mathrm{s}$ is the recovery period from the Herbst maneuver to the level flight, which belongs to phase 3.

The Herbst maneuver performed by the experimental F-16 is simulated by repeating the above seven iterative procedures with given attitude commands ${[{\mu }_c,{\alpha }_c,{\beta }_c]}^T$ and with thrust setting $15200\mathrm{l}\mathrm{b}\mathrm{f}$. The response of the experimental F-16 to the commands of the Herbst maneuver is illustrated in Figure 10 and Figure 11. It appears that the flight of the experimental F-16 closely follows the prescribed Herbst maneuver with reasonable tracking errors. As can be seen, in phase 1 and phase 3, the major change in aircraft attitude is AOA due to the elevator deflection ${\delta }_r$ (see Figure 10d) and the thrust deflection ${\delta }_z$ (see Figure 11a), while the deflection angles of aileron ${\delta }_a$ and rudder ${\delta }_r$ in phase 1 and phase 3 are relatively small (see Figure 10e,f). During phase 2, the aircraft rotates around the wind axis by $120°$ caused by the significant deflections of ${\delta }_a$, ${\delta }_r$, and ${\delta }_y$; meanwhile, the AOA increases from $40°$ to a maximum of $70°$, and then decreases to $40°$.

It is worth noting that during phase 2, the deflection angles of the aileron ${\delta }_a$ and the rudder ${\delta }_r$ has reached their saturation values, and the reason why the aircraft can still complete the heading reversal is mainly due to the deflection of the thrust in the y direction (see Figure 11b). In addition, in phase 2, we cannot ignore the high aerodynamic drag caused by high AOA. From Figure 11c, we can see that due to the high drag effect, the aircraft speed profile has an obvious sagging region in phase 2, causing the aircraft speed to drop from 250 ft/s to 130 ft/s, and then gradually rise again.

Figure 11d–f shows the trajectory of the experimental F-16 performing the Herbst maneuver. After the Herbst maneuver is completed, we can see that in the vertical plane, the height of the aircraft drops about 1000ft, while on the horizontal plane, the turning radius of the aircraft is about 2000ft. In addition, we observe from Figure 11 that the turning point of the Herbst maneuver occurs approximately at $t=10\mathrm{s}$, when the aircraft reaches the highest point, the speed drops to the lowest point, and the bank angle turns to the maximum angle. Overall, the experimental F-16 has accomplished the goal of improving the existing F-16 to a certain extent by proving that the F-16 installed with a thrust vectoring nozzle and controlled by an NDI controller has the maneuverability in the post-stall region.

7. Conclusions

The feasibility of the proposed NDI controller depends on how separable the three-layer loop is and how accurate the input-affine nonlinear model is. An exact input-affine nonlinear model with fully separable three-layer loops reflects that the attitude command can be perfectly tracked, leading to a null response of the sideslip angle, $\beta \left(t\right)\equiv 0$. Therefore, the deviation of the sideslip angle from zero can be regarded as a quantitative measure for the modeling error of the input-affine nonlinear model and the degree of separability of the three-layer loops. Verification simulations show that the transient response of the sideslip angle does deviate from zero. When the AOA command is $40°$, the maximum deviation of the sideslip angle is about 1 degree (see the red curve in Figure 9b). When the AOA command is increased to $70°$, the maximum deviation of the sideslip angle is about 3 degrees (see the red curve in Figure 10b). Obviously, the input-affine nonlinear model exhibits a greater modeling error at higher AOA. The induced sideslip angle can be suppressed by suitably choosing the PI gains in the auxiliary linear controller accompanying the NDI control, or instead by designing a robust linear controller to account directly for the modeling error.