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

Preprint

Article

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

Altmetrics

Downloads

102

Views

100

Comments

This version is not peer-reviewed

Submitted:

06 August 2024

Posted:

07 August 2024

You are already at the latest version

Alerts

Abstract

The nonlinear dynamics exhibited by aircraft flying at high angle of attack (AOA) pose great challenges to the design of flight control. The control design based on nonlinear dynamic inversion (NDI) is to reversely find the required control input based on the current nonlinear dynamics of the aircraft and the control objectives to be achieved. Since it is derived from the nonlinear dynamics of the aircraft, a single NDI controller can cover the entire flight envelope without the need for controller switching or gain scheduling. This paper employs the layered structure of the flight control loop to implement the NDI controller through the combination of aerodynamic control and thrust vectoring control (TVC). The resulting NDI controller is applied to an experimental F-16 installed with TVC to verify its high AOA maneuverability in the post-stall region. The verification simulation shows that the experimental F-16 under the NDI control can effectively track the high AOA command and successfully demonstrate the Herbst maneuver by performing the heading reversal with 70° AOA.

Keywords:

Subject: Engineering - Aerospace Engineering

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

α

and the angle of sideslip

β

, 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

α

and

β

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

α

and

β

, 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

(ϕ, θ, ψ)

with

ϕ

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

α

, sideslip angle

β

, bank angle

μ

(rotation angle around the wind axis), flight path angle

γ

and heading angle

χ

, 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

({\hat{x}}_{E}, {\hat{y}}_{E}, {\hat{z}}_{E})

is located on the earth's surface (ignoring the earth's rotation), where the

{\hat{x}}_{E}

axis points to the local north, the

{\hat{y}}_{E}

axis points to the local east, and the

{\hat{y}}_{E}

axis points to the center of the earth. The origin of the body-axis coordinates

({\hat{x}}_{b}, {\hat{y}}_{b}, {\hat{z}}_{b})

is located at the center of mass of the aircraft, where the

{\hat{x}}_{b}

axis points to the direction of the nose along the symmetry line of the aircraft, the

{\hat{y}}_{b}

axis is perpendicular to

{\hat{x}}_{b}

and points to the right direction of the aircraft, and the

{\hat{z}}_{b}

axis points in the direction of

{\hat{x}}_{b} \times {\hat{y}}_{b}

The wind axis

{\hat{x}}_{w}

points to the velocity direction of the aircraft, so the aircraft’s velocity can be expressed as

{\hat{V}}_{T} = V_{T} {\hat{x}}_{w}

, where

V_{T}

is the magnitude of the velocity and the velocity direction

{\hat{x}}_{w}

relative to the local horizontal coordinates is defined by

(γ, χ)

, where

γ

is the flight path angle between the wind axis

{\hat{x}}_{w}

and the local horizontal plane (see Figure 2a), and

χ

is the heading angle between the projection of the wind axis on the local horizontal plane and the local north

{\hat{x}}_{E}

(see Figure 2c). By rotating an angle χ on the local horizontal plane, and then rotating an angle

γ

on the vertical plane, we can transfer the local horizontal coordinates

({\hat{x}}_{E}, {\hat{y}}_{E}, {\hat{z}}_{E})

to the wind axis coordinates

({\hat{x}}_{w}, {\hat{y}}_{w}, {\hat{z}}_{w})

, as expressed by Eq. (A.1). The change rate of the flight-path states

(V_{T}, χ, γ)

satisfies the following equations (referring to Appendix A for the derivation):

{\dot{V}}_{T} = \frac{1}{m} [(- D + \bar{Y} \sin β - m g \sin γ) + (T_{x} \cos β \cos α + T_{y} \sin β + T_{z} \cos β \sin α)],

(1)

\begin{matrix} \dot{χ} = & \frac{1}{m V_{T} \cos γ} [(L \sin μ + \bar{Y} \cos μ \cos β) + T_{x} (\sin μ \sin α - \cos μ \sin β \cos α) \\ + T_{y} (\cos μ \cos β) - T_{z} (\cos μ \sin β \sin α + \sin μ \cos α)], \end{matrix}

(2)

\begin{matrix} \dot{γ} = & \frac{1}{m V_{T}} [(L \cos μ - \bar{Y} \sin μ \cos β - m g \cos γ) + T_{x} (\sin μ \sin β \cos α + \cos μ \sin α) \\ - T_{y} (\sin μ \cos β) + T_{z} (\sin μ \sin β \sin α - \cos μ \cos α)] . \end{matrix}

(3)

The lift

L

, the drag

D

, and the lateral force

\bar{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,

[\begin{matrix} X \\ Y \\ Z \end{matrix}] = \bar{q} S [\begin{matrix} C_{X T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ C_{Y T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ C_{Z T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \end{matrix}] = [\begin{matrix} - \cos α \cos β & - \cos α \sin β & \sin α \\ - \sin β & \cos β & 0 \\ - \sin α \cos β & - \sin α \sin β & - \cos α \end{matrix}] [\begin{matrix} D \\ \bar{Y} \\ L \end{matrix}],

(4)

where

\bar{q} = ρ {V_{T}}^{2} / 2

is the dynamic pressure,

S

is the area of the main wing,

ρ

is the atmospheric density, and

C_{X T}

C_{Y T}

C_{Z T}

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

δ_{e}

, the aileron

δ_{a}

, and the rudder

δ_{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 ° \leq α \leq 90 °

- 30 ° \leq β \leq 30 °

and

V_{T} \leq 06 M a c h

, which covers the post-stall flight region considered here.

The thrust components

{[T_{x}, T_{y}, T_{z}]}^{T}

in Eq. (1) – Eq. (3) are generated by the thrust vectoring nozzle:

[\begin{matrix} T_{x} \\ T_{y} \\ T_{z} \end{matrix}] = [\begin{matrix} T \cos δ_{z} \cos δ_{y} \\ T \sin δ_{y} \\ - T \sin δ_{z} \cos δ_{y} \end{matrix}],

(5)

where

T

is the total thrust of the aircraft,

δ_{y}

is the left-right deflection angle of the nozzle (defining it as positive to the left), and

δ_{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

{\hat{y}}_{b}

direction of the body axis. When the nozzle deflects downward, an upward thrust

T_{z}

will be generated along the negative

{\hat{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,

μ \to β \to α

as shown by Eq. (A2), the wind-axis coordinates

({\hat{x}}_{w}, {\hat{y}}_{w}, {\hat{z}}_{w})

can be transferred into the body-axis coordinates

({\hat{x}}_{b}, {\hat{y}}_{b}, {\hat{z}}_{b})

. The governing equations for

(α, β, μ)

can be established as follows (see Appendix A for the derivation):

\dot{α} = q - (p \cos α + r \sin α) \tan β + \frac{1}{m V_{T} \cos β} [(- L + m g \cos γ \cos μ) + (- T_{x} \sin α + T_{z} \cos α)],

(6)

\dot{β} = - r \cos α + p \sin α + \frac{1}{m V_{T}} [(\bar{Y} \cos β + m g \cos γ \sin μ) + (- T_{x} \sin β \cos α + T_{y} \cos β + T_{z} \sin β \sin α)],

(7)

\dot{μ} = (p \cos α + r \sin α) \sec β + \frac{L}{m V_{T}} (\tan γ \sin μ + \tan β) + \frac{(\bar{Y} + T_{y})}{m V_{T}} \tan γ \cos μ \cos β - \frac{g}{V_{T}} \cos γ \cos μ \tan β + \frac{(T_{x} \sin α - T_{z} \cos α)}{m V_{T}} (\tan γ \sin μ + \tan β) - \frac{(T_{x} \cos α + T_{z} \sin α)}{m V_{T}} \tan γ \cos μ \sin β .

(8)

It can be seen above that the attitude change rates

(\dot{α}, \dot{β}, \dot{μ})

depend on the body-axis angular rates

(p, q, r)

, which in turn depend on the external moments

M = {[M_{x}, M_{y}, M_{z}]}^{T}

\dot{p} = (C_{1} r + C_{2} p) q + C_{3} M_{x} + C_{4} M_{z},

(9)

\dot{q} = C_{5} p r - C_{6} (p^{2} - r^{2}) + C_{7} M_{y},

(10)

\dot{r} = (C_{8} p - C_{2} r) 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

M

experienced by the aircraft is composed of the aerodynamic moment

M^{(A)}

and the thrust moment

M^{(T)}

. The aerodynamic moment

M^{(A)}

is determined by the aerodynamic coefficients:

M^{(A)} = [\begin{matrix} M_{x}^{(A)} \\ M_{y}^{(A)} \\ M_{z}^{(A)} \end{matrix}] = [\begin{matrix} L \\ M \\ N \end{matrix}] = \bar{q} S [\begin{matrix} b C_{l T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ \bar{C} C_{m T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ {b C}_{n T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \end{matrix}],

(12)

where

b

is the span of the main wing,

\bar{C}

is the average chord length, and

C_{l T}

C_{m T}

, and

C_{n T}

are the aerodynamic moment coefficients in the three directions of the aircraft body axis. The thrust moment

M^{(T)}

is produced by the thrust about the center of gravity of the aircraft and can be expressed as

M^{(T)} = [\begin{matrix} M_{x}^{(T)} \\ M_{y}^{(T)} \\ M_{z}^{(T)} \end{matrix}] = [\begin{matrix} 0 \\ l_{T} T_{z} \\ - l_{T} T_{y} \end{matrix}] = [\begin{matrix} 0 \\ {- l}_{T} T \sin δ_{z} \cos δ_{y} \\ - l_{T} T \sin δ_{y} \end{matrix}],

(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}, χ, γ)

are given by Eq. (1) – Eq. (3), the attitude states

(α, β, μ)

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

{[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{T}

, which generate the aerodynamic force

F^{(A)}

by Eq. (4), the aerodynamic moment

M^{(A)}

by Eq. (12), the thrust force

F^{(T)}

by Eq. (5), and the thrust moment

M^{(T)}

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

(ψ, θ, ϕ, u, v, w, p, q, r)

, the wind-axis based state variables

(V_{T}, χ, γ, α, β, μ, 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

x = {[V_{T}, χ, γ, μ, α, β, p, q, r]}^{T}

derived previously can be expressed as an input-affine nonlinear form:

\begin{matrix} \dot{x} = f (x) + g (x) u, & x (t_{0}) = x_{0}, & y = h (x) \end{matrix},

(14)

where

f (x)

is the vector function with 9 elements and

g (x)

is the matrix function with dimension

9 \times 5

, and

y = h (x)

is the system’s output vector. The control input

u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{T}

consists of the aerodynamic deflection angles

(δ_{a}, δ_{e}, δ_{r})

and the nozzle deflection angles

(δ_{y}, δ_{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

u

; however, the small angle assumption is not valid for the state

x

. This situation is reflected in the input-affine nonlinear form (14), where the nonlinear nature of the state

x

is retained, while the control input

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

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

y

and the command

y_{c}

The control input

u

that linearizes the nonlinear system (14) can be found as,

u = g^{- 1} (x) (v - f (x)) .

(15)

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

\dot{x} = f (x) + g (x) g^{- 1} (x) (v - f (x)) = 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

v

can then be designed by any linear control method to drive the system output

y

to follow the target output

y_{c}

The existence of the ideal NDI controller (15) depends on the invertibility of the matrix function

g (x)

. Unfortunately, the inverse of the matrix

g (x)

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

v = {\dot{x}}_{d e s}

can only be tracked partially and the best control we can design is to minimize the tracking error by replacing the inverse

g^{- 1}

in Eq. (15) with its pseudo inverse,

g^{#} = {(g^{T} g)}^{- 1} g^{T},

(17)

where

g^{#}

is the solution minimizing the distance

{‖g g^{#} - I_{9}‖}_{2}

, where

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

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

{[α, β, μ]}^{T}

will change accordingly. The heading angle

χ

and flight path angle

γ

of the aircraft will not respond until the attitude of the aircraft changes. The sequence of response is that the control

u

affects

{[p, q, r]}^{T}

first, then

{[p, q, r]}^{T}

affects

{[α, β, μ]}^{T}

, and finally

{[α, β, μ]}^{T}

affects

(V_{T}, χ, γ)

. Therefore, the response speed of the 9 states to the control

u

is that the angular-rate state

{[p, q, r]}^{T}

has the fastest response, the attitude state

{[α, β, μ]}^{T}

is slower, and the flight-path state

(V_{T}, χ, γ)

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}, χ, γ)

to find reversely the required attitude state

{[α, β, μ]}^{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

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

{[α, β, μ]}^{T}

, and the outermost loop is the flight-path control loop regarding the control of the slowest state

(V_{T}, χ, γ)

. 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

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

{[α, β, μ]}^{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

{[α_{c}, β_{c}, μ_{c}]}^{T}

4. Angular-Rate NDI Control

The NDI control design in the inner loop is to determine the control surface deflection

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

u

, the fully nonlinear equations (9) – (11) must be transformed into an affine nonlinear form, wherein the control input

u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{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

u

, which yields the following affine nonlinear form:

[\begin{matrix} \dot{p} \\ \dot{q} \\ \dot{r} \end{matrix}] = f_{1} (x) + g_{1} (x) u = [\begin{matrix} f_{p} (x) \\ f_{q} (x) \\ f_{r} (x) \end{matrix}] + [\begin{matrix} g_{p_{δ_{a}}} & 0 & g_{p_{δ_{r}}} & g_{p_{δ_{y}}} & 0 \\ 0 & g_{q_{δ_{e}}} & 0 & 0 & g_{q_{δ_{z}}} \\ g_{r_{δ_{a}}} & 0 & g_{r_{δ_{r}}} & g_{r_{δ_{y}}} & 0 \end{matrix}] [\begin{matrix} δ_{a} \\ δ_{e} \\ δ_{r} \\ δ_{y} \\ δ_{z} \end{matrix}],

(18)

where explicit expressions for the elements in the vector

f_{1} (x)

and the matrix

g_{1} (x)

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

x

are involved in the expressions of

f_{1} (x)

and

g_{1} (x)

; however, as measured under the time scale of

{[p, q, r]}^{T}

, the remaining six states in

x

, i.e.,

{[V_{T}, χ, γ, α, β, μ]}^{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

u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{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

g_{1} (x)

in Eq. (18) is (right) invertible, i.e.,

g_{1} g_{1}^{- 1} = I_{3}

, but its inverse

g_{1}^{- 1}

is not unique. Due to the non-uniqueness of

g_{1}^{- 1}

, all the possible solutions of

u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{T}

can be expressed as

u = N \bar{u} ⟶ [\begin{matrix} δ_{a} \\ δ_{e} \\ δ_{r} \\ δ_{y} \\ δ_{z} \end{matrix}] = [\begin{matrix} N_{a P} & N_{a Q} & N_{a R} \\ N_{e P} & N_{e Q} & N_{e R} \\ N_{r P} & N_{r Q} & N_{r R} \\ N_{y P} & N_{y Q} & N_{y R} \\ N_{z P} & N_{z Q} & N_{z R} \end{matrix}] [\begin{matrix} δ_{P} \\ δ_{Q} \\ δ_{R} \end{matrix}],

(19)

where

N

is the control allocation matrix, whose elements can be arbitrarily assigned as long as

g_{1} N

is nonsingular. The effective control

\bar{u}

consists of three effective control surface deflections,

δ_{P}

δ_{Q}

, and

δ_{R}

, required for controlling the rolling rate

p

, the pitching rate

q

, and the yawing

r

, respectively.

The elements of matrix

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,

δ_{a} = N_{a P} δ_{P} + N_{a Q} δ_{Q} + N_{a R} δ_{R},

(20)

which means that

N_{a P}

percentage of the total aileron deflection

δ_{a}

is allocated to the rolling rate control,

N_{a Q}

percentage to the pitching rate control, and

N_{a R}

percentage to the rolling rate control. It is worth noting that

δ_{y}

and

δ_{z}

are the deflection angles of the thrust vectoring nozzle and the last two rows of

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

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

N

, such as daisy chain, explicit ganging, pseudo control or control selector, and pseudo inverse [30,31].

The effective control

\bar{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

\bar{u}

can be derived by letting

{[\dot{p}, \dot{q}, \dot{r}]}^{T} = {[{\dot{p}}_{c}, {\dot{q}}_{c}, {\dot{r}}_{c}]}^{T}

and

u = N \bar{u}

in Eq. (18) to obtain

\bar{u} = {(g_{1} N)}^{- 1} ([\begin{matrix} {\dot{p}}_{c} \\ {\dot{q}}_{c} \\ {\dot{r}}_{c} \end{matrix}] - [\begin{matrix} f_{p} (x) \\ f_{q} (x) \\ f_{r} (x) \end{matrix}]) ⟶ u = N \bar{u} = N {(g_{1} N)}^{- 1} ([\begin{matrix} {\dot{p}}_{c} \\ {\dot{q}}_{c} \\ {\dot{r}}_{c} \end{matrix}] - [\begin{matrix} f_{p} (x) \\ f_{q} (x) \\ f_{r} (x) \end{matrix}]),

(21)

where

g_{1} N

is a nonsingular

3 \times 3

matrix and its inverse is defined uniquely. With the effective control

\bar{u}

given by Eq. (21), the actual control surface deflections

u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{T}

are then determined by passing

\bar{u}

via the allocation matrix

N

as shown in Figure 6.

With NDI controller

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

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

e

\begin{matrix} {\dot{p}}_{c} = ω_{p} (p_{c} - p) + \int Ω_{p} (p_{c} - p) d t, \\ \begin{matrix} {\dot{q}}_{c} = ω_{q} (q_{c} - q) + \int Ω_{q} (q_{c} - q) d t, \\ {\dot{r}}_{c} = ω_{r} (r_{c} - r) + \int Ω_{r} (r_{c} - r) d t, \end{matrix} \end{matrix}

(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,

{[ω_{p} ω_{q} ω_{r}]}^{T}

are the adjustable proportional gains, and

{[Ω_{p} Ω_{q} Ω_{r}]}^{T}

are the adjustable integral gains.

5. Attitude NDI Control

The dominant states in the outer NDI control loop are the slow attitude states

{[μ, α, β]}^{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}, χ, γ)

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

u

, while the dominant response of the attitude dynamics

{[μ, α, β]}^{T}

is caused by the angular-rate input

{[p, q, r]}^{T}

, rather than by the control surface deflections

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

u

on the attitude states. The result can be expressed as the following affine nonlinear form:

[\begin{matrix} \dot{μ} \\ \dot{α} \\ \dot{β} \end{matrix}] = f_{2} (x_{2}) + g_{2} (x_{2}) u_{2} = [\begin{matrix} f_{μ} (x_{2}) \\ f_{α} (x_{2}) \\ f_{β} (x_{2}) \end{matrix}] + [\begin{matrix} \cos α \sec β & 0 & \sin α \sec β \\ - \cos α \tan β & 1 & - \sin α \tan β \\ \sin α & 0 & - \cos α \end{matrix}] [\begin{matrix} p_{c} \\ q_{c} \\ r_{c} \end{matrix}],

(23)

where the reduced state

x_{2}

are defined as

x_{2} = {[V_{T}, χ, γ, α, β, μ]}^{T}

and

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_{α} (x_{2}) = \frac{1}{m V_{T} \cos β} [(- L + m g \cos γ \cos μ) - T \sin α],

\begin{matrix} f_{β} (x_{2}) = \frac{1}{m V_{T}} [\bar{Y} \cos β + m g \cos γ \sin μ - T \sin β \cos α], \end{matrix}

\begin{matrix} f_{μ} (x_{2}) = & \frac{L}{m V_{T}} (\tan γ \sin μ + \tan β) + \frac{1}{m V_{T}} (\bar{Y} \tan γ \cos μ \cos β - m g \cos γ \cos μ \tan β) \\ + \frac{T}{m V_{T}} (\sin α \tan γ \sin μ + \sin α \tan β - \cos α \tan γ \cos μ \sin β) . \end{matrix}

In Eq. (23), the angular-rate state

u_{2} = {[p_{c}, q_{c}, r_{c}]}^{T}

is the control input to the system and the attitude state

{[μ, α, β]}^{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

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{μ}, \dot{α}, \dot{β}]}^{T}

tracks the desired response

{[{\dot{μ}}_{c}, {\dot{α}}_{c}, {\dot{β}}_{c}]}^{T}

By letting

{[\dot{μ}, \dot{α}, \dot{β}]}^{T} = [{\dot{μ}}_{c}, {\dot{α}}_{c}, {\dot{β}}_{c}]

in Eq. (23), the ideal control input

u_{2}

can be determined as

u_{2} ≜ [\begin{matrix} p_{c} \\ q_{c} \\ r_{c} \end{matrix}] = g_{2}^{- 1} (x_{2}) ([\begin{matrix} {\dot{μ}}_{c} \\ {\dot{α}}_{c} \\ {\dot{β}}_{c} \end{matrix}] - [\begin{matrix} f_{μ} (x_{2}) \\ f_{α} (x_{2}) \\ f_{β} (x_{2}) \end{matrix}]),

(24)

where

g_{2} (x_{2})

is the

3 \times 3

nonsingular matrix defined in Eq. (23) and its inverse can be expressed analytically as

g_{2}^{- 1} = {[\begin{matrix} \cos α \sec β & 0 & \sin α \sec β \\ - \cos α \tan β & 1 & - \sin α \tan β \\ \sin α & 0 & - \cos α \end{matrix}]}^{- 1} = [\begin{matrix} \cos α \cos β & 0 & \sin α \\ \sin β & 1 & 0 \\ \sin α \cos β & 0 & - \cos α \end{matrix}]

(25)

Since

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

u_{2}

, a PI controller is incorporated into the outer loop to attain the desired aircraft’s attitude

{[μ_{c}, α_{c}, β_{c}]}^{T}

according to the attitude tracking error, as shown in Figure 7. Like the inner-loop PI structure in Eq. (22), proportional gains

{[ω_{μ}, ω_{α}, ω_{β}]}^{T}

and integral gains

{[Ω_{μ}, Ω_{α}, Ω_{β}]}^{T}

are used in the outer-loop attitude controller as follows:

\begin{matrix} {\dot{μ}}_{c} = ω_{μ} (μ_{c} - μ) + \int Ω_{μ} (μ_{c} - μ) d t, \\ \begin{matrix} {\dot{α}}_{c} = ω_{α} (α_{c} - α) + \int Ω_{α} (α_{c} - μ) d t, \\ {\dot{β}}_{c} = ω_{β} (β_{c} - β) + \int Ω_{β} (β_{c} - μ) d t, \end{matrix} \end{matrix}

(26)

where

{[μ_{c}, α_{c}, β_{c}]}^{T}

are the desired attitudes and

{[μ, α, β]}^{T}

are the aircraft’s actual attitudes provided by the attitude sensors. The magnitude of the proportional gains

{[ω_{μ}, ω_{α}, ω_{β}]}^{T}

reflects the bandwidth of the attitude control in the outer loop and is much smaller than proportional gains

{[ω_{p}, ω_{q}, ω_{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

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:

N_{o f f} = [\begin{matrix} 0.75 & 0 & 0.25 \\ 0 & 1 & 0 \\ 0.25 & 0 & 0.75 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], N_{o n} = [\begin{matrix} 0.75 & 0 & 0.25 \\ 0 & 1 & 0 \\ 0.25 & 0 & 0.75 \\ 0.25 & 0 & 0.5 \\ 0 & 0.5 & 0 \end{matrix}] .

(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

δ_{y}

and

δ_{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

N_{o f 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

N_{o n}

corresponds to the situation when TVC is activated, in which the left-right deflection (

δ_{y}

) of the thrust contributes 25% to the rolling rate

p

and 50% to the yawing rate

r

, and the up-down deflection (

δ_{z}

) of the thrust contributes 50% to the pitching rate

q

. The setting of the control allocation matrix

N

is not limited to the form expressed by Eq. (27). We can set different matrix

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

u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{T}

and nine state outputs

x = {[V_{T}, χ, γ, μ, α, β, p, q, r]}^{T}

. With the measured state outputs

x

, the control inputs

u

are determined by a sequence of dynamic inverse processes. Here we assume that the attitude command

{[μ_{c}, α_{c}, β_{c}]}^{T}

has been designed by an outermost flight-path loop not mentioned here, and we will start with

{[μ_{c}, α_{c}, β_{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 $x = {[V_{T}, χ, γ, μ, α, β, p, q, r]}^{T}$ .
Compute the current values of $f_{1} (x)$ and $g_{1} (x)$ by Eq. (18), and the current values of $f_{2} (x_{2})$ and $g_{2} (x_{2})$ by Eq. (23) with the new measurement $x$ .
Compute the desired change rate of the attitude ${[{\dot{μ}}_{c}, {\dot{α}}_{c}, {\dot{β}}_{c}]}^{T}$ by Eq. (26) with given command ${[μ_{c}, α_{c}, β_{c}]}^{T}$ .
Compute the desired body-axis angular rates ${[p_{c}, q_{c}, r_{c}]}^{T}$ by Eq. (24) with given ${[{\dot{μ}}_{c}, {\dot{α}}_{c}, {\dot{β}}_{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 $u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{T}$ by Eq. (21) with given ${[{\dot{p}}_{c}, {\dot{q}}_{c}, {\dot{r}}_{c}]}^{T}$ and control allocation matrix $N$ .
Solve the 6D model with the new control surface deflections $u$ for the new value of the states $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 l b f$ with initial speed $V_{T} = 500 f t / s$ , initial height $h_{0} = 15000 f t$ , and initial angular rates $p_{0} = q_{0} = r_{0} = 0 r a d / s$ .
●: In the attitude loop, the proportional gains and the integral gains are set to be $ω_{μ} = ω_{α} = ω_{β} = 2 r a d / s$ and $Ω_{μ} = Ω_{α} = Ω_{β} = 1 r a d / s^{2}$ .
●: In the angular-rate loop, the proportional gains and the integral gains are set to be $ω_{μ} = ω_{α} = ω_{β} = 10 r a d / s$ and $Ω_{μ} = Ω_{α} = Ω_{β} = 4 r a d / s^{2}$ .
●: The TVC function is deactivated in this case by choosing the control allocation matrix $N_{o f f}$ in Eq. (20).
●: The sideslip angle command $β_{c}$ is set to zero, and the AOA command $α_{c}$ and the bank angle command $μ_{c}$ are given, respectively, as

$α_{c} = \{\begin{matrix} 0 ° & , t \leq 2 \\ 0 ° + (35 ° - 0 °) \times (t - 2) / 6 & , 2 < t \leq 8 \\ 35 ° - (35 ° - 0 °) \times (t - 8) / 7 & , 8 < t \leq 15 \\ 0 ° & , 15 < t \leq 30 \end{matrix}$

(28)

$μ_{c} = \{\begin{matrix} 0 ° & , t \leq 2 \\ 0 ° + (100 ° - 0 °) \times (t - 2) / 7 & , 2 < t \leq 9 \\ 100 ° - (100 ° - 0 °) \times (t - 9) / 6 & , 9 < t \leq 15 \\ 0 ° & , 15 < t \leq 30 \end{matrix}$

(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

α_{c}

increases linearly from

0 °

to the stall angle

35 °

and then decreases to

0 °

. Meanwhile, the bank angle command

μ_{c}

increases linearly from

0 °

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 s

and

t = 11 s

, within which the sideslip angle

β (t)

deviates from zero about 3 degrees, and slight overshoots occur in the response of

α (t)

and

μ (t)

Comparing Figure 8e, we find that the deviation of the attitudes between

t = 9 s

and

t = 11 s

is primarily caused by the large aileron deflection angle

δ_{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 °

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 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

N_{o 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

{[μ_{c}, α_{c}, β_{c}]}^{T}

[32]:

α_{c} = \{\begin{matrix} 7 ° + (33 ° - 0 °) \times (t - 0) / 5 & , t \leq 5 \\ 40 ° + (70 ° - 40 °) \times (t - 5) / 2.8 & , 5 < t \leq 7.8 \\ 70 ° & , 7.8 < t \leq 8.2 \\ 70 ° - (70 ° - 50 °) \times (t - 8.2) / 1.8 & , 8.2 < t \leq 10 \\ 50 ° - (50 ° - 40.68 °) \times (t - 10) / 2.5 & , 10 < t \leq 12.5 \\ 40.68 ° - (40.68 ° - 40 °) \times (t - 12.5) / 1.5 & , 12.5 < t \leq 14 \\ 40 ° - (40 ° - 10 °) \times (t - 14) / 6 & , 14 < t \leq 20 \end{matrix}

(30)

μ_{c} = \{\begin{matrix} 0 ° + (20.6 ° - 0 °) \times (t - 0) / 2 & , t \leq 2 \\ 20.6 ° + (- 19.5 ° - 20.6 °) \times (t - 2) / 3 & , 2 < t \leq 5 \\ - 19.5 ° + (- 100 ° + 19.5 °) \times (t - 5) / 2 & , 5 < t \leq 7 \\ - 100 ° + (- 120 ° + 100 °) \times (t - 7) / 1.7 & , 7 < t \leq 8.7 \\ - 120 ° - (100 ° - 120 °) \times (t - 8.7) / 1.7 & , 8.7 < t \leq 10.4 \\ - 100 ° - (40.1 ° - 100 °) \times (t - 10.4) / 2 & , 10.4 < t \leq 12.4 \\ - 40.1 ° - (20 ° - 40.1 °) \times (t - 14) / 6 & , 12.4 < t \leq 15 \\ - 20 ° - (10 ° - 20 °) \times (t - 15) / 5 & , 15 < t \leq 20 \end{matrix}

(31)

and

β_{c} = 0

. Roughly speaking, the period

0 s \leq t < 5 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 s \leq t < 15 s

is the working period of the Herbst maneuver, which performs the heading reversal and belongs to phase 2. The period

15 s \leq t < 20 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

{[μ_{c}, α_{c}, β_{c}]}^{T}

and with thrust setting

15200 l b 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

δ_{r}

(see Figure 10d) and the thrust deflection

δ_{z}

(see Figure 11a), while the deflection angles of aileron

δ_{a}

and rudder

δ_{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

δ_{a}

δ_{r}

, and

δ_{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

δ_{a}

and the rudder

δ_{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 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,

β (t) \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.

Author Contributions

Conceptualization, methodology and writing, C.D.; software, and validation, S.Y.

Funding

This research was funded by Taiwan National Science and Technology Council under Award No. MOST 111-2221-E-006-215-MY3.

Appendix A. The Input-Affine Nonlinear Form of F-16 with TVC

A.1. Rotation Matrix between Three Coordinate Systems

The flight motion of an aircraft is generally described by the body-axis coordinates. However, when the aircraft flies in a high AOA region, it is more suitable to be described by the wind-axis coordinate system. The purpose of this appendix is to derive the equations of motion in the wind-axis coordinate system and transform it to the input-affine nonlinear form for the use in the main text.

Through the rotation matrices generated by the flight path angle

γ

and heading angle

χ

, the local horizontal coordinates can be converted into the wind-axis coordinates as follows,

[\begin{matrix} {\hat{x}}_{w} \\ {\hat{y}}_{w} \\ {\hat{z}}_{w} \end{matrix}] = [\begin{matrix} \cos γ & 0 & - \sin γ \\ 0 & 1 & 0 \\ \sin γ & 0 & \cos γ \end{matrix}] [\begin{matrix} \cos χ & \sin χ & 0 \\ - \sin χ & \cos χ & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} {\hat{x}}_{E} \\ {\hat{y}}_{E} \\ {\hat{z}}_{E} \end{matrix}] = R_{E}^{w} [\begin{matrix} {\hat{x}}_{E} \\ {\hat{y}}_{E} \\ {\hat{z}}_{E} \end{matrix}] .

(A.1)

The inverse transformation of

R_{E}^{w}

is given by

R_{w}^{E} = {(R_{E}^{w})}^{- 1} = {(R_{E}^{w})}^{T}

. On the other hand, the wind-axis coordinates

({\hat{x}}_{w}, {\hat{y}}_{w}, {\hat{z}}_{w})

can be converted into the body-axis coordinates

({\hat{x}}_{b}, {\hat{y}}_{b}, {\hat{z}}_{b})

through the rotation matrices generated by the bank angle

μ

, the sideslip angle

β

and the AOA

α

, as follows

[\begin{matrix} {\hat{x}}_{b} \\ {\hat{y}}_{b} \\ {\hat{z}}_{b} \end{matrix}] = [\begin{matrix} \cos α & 0 & - \sin α \\ 0 & 1 & 0 \\ \sin α & 0 & \cos α \end{matrix}] [\begin{matrix} \cos (- β) & \sin (- β) & 0 \\ - \sin (- β) & \cos (- β) & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} 1 & 0 & 0 \\ 0 & \cos μ & \sin μ \\ 0 & - \sin μ & \cos μ \end{matrix}] [\begin{matrix} {\hat{x}}_{w} \\ {\hat{y}}_{w} \\ {\hat{z}}_{w} \end{matrix}] = R_{w}^{b} [\begin{matrix} {\hat{x}}_{w} \\ {\hat{y}}_{w} \\ {\hat{z}}_{w} \end{matrix}] .

(A.2)

The inverse transformation of

R_{w}^{b}

is given by

R_{b}^{w} = {(R_{w}^{b})}^{- 1} = {(R_{w}^{b})}^{T}

A.2. Equations of Motion in Body-axis Coordinates

The velocity of the center of mass of the aircraft under the body-axis coordinates is denoted by

{[\begin{matrix} u & v & w \end{matrix}]}^{T}

, and the angular velocity of the aircraft rotating around the body axis is denoted by

{[\begin{matrix} p & q & r \end{matrix}]}^{T}

. The motion of the center of mass of the aircraft is described by the following equation,

m ({[\begin{matrix} \dot{u} \\ \dot{v} \\ \dot{w} \end{matrix}]}_{b} + [\begin{matrix} 0 & - r & q \\ r & 0 & - p \\ - q & p & 0 \end{matrix}] {[\begin{matrix} u \\ v \\ w \end{matrix}]}_{b}) = {[\begin{matrix} F_{x} \\ F_{y} \\ F_{z} \end{matrix}]}_{b},

(A.3)

where the subscript “b” represents the component under the body-axis coordinate system. On the other hand, the rotational motion of the aircraft around its body axis is described by the following equations,

\dot{p} = (C_{1} r + C_{2} p) q + C_{3} M_{x} + C_{4} M_{z},

(A.4)

\dot{q} = C_{5} p r - C_{6} (p^{2} - r^{2}) + C_{7} M_{y},

(A.5)

\dot{r} = (C_{8} p - C_{2} r) q + C_{4} M_{x} + C_{9} M_{z},

(A.6)

where the constants

C_{i}

are defined by

C_{1} = [(I_{y} - I_{z}) I_{z} - I_{x z}^{2}] / Γ, C_{2} = [({I_{x} - I}_{y} + I_{z}) I_{x z}] / Γ, C_{3} = I_{z} / Γ, C_{4} = I_{x z} / Γ, C_{5} = ({I_{Z} - I}_{x}) / I_{y}, C_{6} = I_{x z} / I_{y}, C_{7} = 1 / I_{y}, C_{8} = [I_{x} (I_{x} - I_{y}) + I_{x z}^{2}] / Γ, C_{9} = I_{x} / Γ, Γ = I_{x} I_{z} - I_{x z}^{2}

In Eq. (A.3), the external force

{[\begin{matrix} F_{x} & F_{y} & F_{z} \end{matrix}]}^{T}

experienced by the aircraft is composed of the aerodynamic force

F^{(A)}

, thrust

F^{(T)}

and gravity

F^{(G)}

{[\begin{matrix} F_{x} \\ F_{y} \\ F_{z} \end{matrix}]}_{b} = {[\begin{matrix} F_{x}^{(A)} \\ F_{y}^{(A)} \\ F_{z}^{(A)} \end{matrix}]}_{b} + {[\begin{matrix} F_{x}^{(T)} \\ F_{y}^{(T)} \\ F_{z}^{(T)} \end{matrix}]}_{b} + {[\begin{matrix} F_{x}^{(G)} \\ F_{y}^{(G)} \\ F_{z}^{(G)} \end{matrix}]}_{b} .

(A.7)

The aerodynamic force

F^{(A)}

is computed by:

{[\begin{matrix} F_{x}^{(A)} \\ F_{y}^{(A)} \\ F_{z}^{(A)} \end{matrix}]}_{b} = [\begin{matrix} X \\ Y \\ Z \end{matrix}] = \bar{q} S [\begin{matrix} C_{X T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ C_{Y T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ C_{Z T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \end{matrix}],

(A.8)

where the body force coefficients

C_{X T}

C_{Y T}

, and

C_{Z T}

are computed from the following linear expansions:

C_{X T} = C_{X} (α, β, δ_{e}) + δ C_{X L E F} (1 - \frac{δ_{L E F}}{25}) + \frac{q \bar{C}}{2 V_{T}} [C_{X q} (α) + δ C_{X q L E F} (α) (1 - \frac{δ_{L E F}}{25})],

\begin{matrix} C_{Y T} = & C_{Y} (α, β) + δ C_{Y L E F} (1 - \frac{δ_{L E F}}{25}) + [δ C_{Y_{δ_{a}}} + δ C_{Y_{{δ_{a}}_{L E F}}} (1 - \frac{δ_{L E F}}{25})] (\frac{δ_{a}}{21.5}) + δ C_{Y_{δ_{r}}} (\frac{δ_{r}}{30}) \\ + (\frac{r b}{2 V_{T}}) [C_{Y r} (α) + δ C_{Y_{r_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + (\frac{p b}{2 V_{T}}) [C_{Y p} (α) + δ C_{Y_{p_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})], \end{matrix}

C_{Z T} = C_{Z} (α, β, δ_{e}) + δ C_{Z L E F} (1 - \frac{δ_{L E F}}{25}) + \frac{q \bar{C}}{2 V_{T}} [C_{Z q} (α) + δ C_{Z q L E F} (α) (1 - \frac{δ_{L E F}}{25})] .

The force coefficients in the above formula can be directly read from the aerodynamic coefficient table of F-16 with the following definitions.

●: $δ C_{X L E F} = C_{X L E F} (α, β) - {C_{X} (α, β, δ_{e})}_{δ_{e} = 0}$
●: $δ C_{Y L E F} = C_{Y L E F} (α, β) - C_{Y} (α, β)$
●: $δ C_{Y_{δ_{a}}} = C_{Y_{δ_{a}}} (α, β) - C_{Y} (α, β)$
●: $δ C_{Y_{{δ_{a}}_{L E F}}} = C_{Y_{{δ_{a}}_{L E F}}} (α, β) - C_{Y_{r_{L E F}}} (α, β) - δ C_{Y_{δ_{a}}}$
●: $δ C_{Y_{δ_{r}}} = C_{Y_{δ_{r}}} (α, β) - C_{Y} (α, β)$
●: $δ C_{Z L E F} = C_{Z L E F} (α, β) - {C_{Z} (α, β, δ_{e})}_{δ_{e} = 0}$

The three body-axis components

{[T_{x}, T_{y}, T_{z}]}^{T}

of the thrust are given by Eq. (5). The three body-axis components of gravity can be obtained by converting the geocentric component to the body-axis components:

{[\begin{matrix} F_{x}^{(G)} \\ F_{y}^{(G)} \\ F_{z}^{(G)} \end{matrix}]}_{b} = R_{E}^{b} {[\begin{matrix} F_{x}^{(G)} \\ F_{y}^{(G)} \\ F_{z}^{(G)} \end{matrix}]}_{E} = R_{E}^{b} [\begin{matrix} 0 \\ 0 \\ m g \end{matrix}] .

(A.9)

In Eq. (A4) – Eq. (A6), the external moment

M

is composed of the aerodynamic moment

M^{(A)}

and the thrust moment

M^{(T)}

, where the thrust moment

M^{(T)}

is given by Eq. (13) and

M^{(A)}

is given by Eq. (12) rewritten as follows:

M^{(A)} = [\begin{matrix} L \\ M \\ N \end{matrix}] = \bar{q} S [\begin{matrix} b C_{l T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ \bar{C} C_{m T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \\ {b C}_{n T} (α, β, p, q, r, δ_{e}, δ_{a}, δ_{r}) \end{matrix}],

(A.10)

where the moment coefficients

C_{l T}

C_{m T}

, and

C_{n T}

involved in

M^{(A)}

can be computed from the following linear expansions:

\begin{matrix} C_{l T} = & C_{l} (α, β, δ_{e}) + δ C_{l L E F} (1 - \frac{δ_{L E F}}{25}) + [δ C_{l_{δ_{a}}} + δ C_{l_{{δ_{a}}_{L E F}}} (1 - \frac{δ_{L E F}}{25})] (\frac{δ_{a}}{21.5}) + δ C_{l_{δ_{r}}} (\frac{δ_{r}}{30}) \\ + (\frac{r b}{2 V_{T}}) [C_{l r} (α) + δ C_{l_{r_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + (\frac{p b}{2 V_{T}}) [C_{l p} (α) + δ C_{l_{p_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + δ C_{l_{β}} (α) β, \end{matrix}

\begin{matrix} C_{m T} = & C_{m} (α, β, δ_{e}) + C_{Z T} [x_{c g r} - x_{c g}] + δ C_{m L E F} (1 - \frac{δ_{L E F}}{25}) \\ + \frac{q \bar{C}}{2 V_{T}} [C_{m q} (α) + δ C_{m_{q_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + δ C_{m} (α), \end{matrix}

\begin{matrix} C_{n T} = & C_{n} (α, β, δ_{e}) + δ C_{n L E F} (1 - \frac{δ_{L E F}}{25}) + [δ C_{n_{δ_{a}}} + δ C_{n_{{δ_{a}}_{L E F}}} (1 - \frac{δ_{L E F}}{25})] (\frac{δ_{a}}{21.5}) + δ C_{n_{δ_{r}}} (\frac{δ_{r}}{30}) \\ + (\frac{r b}{2 V_{T}}) [C_{n r} (α) + δ C_{n_{r_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + (\frac{p b}{2 V_{T}}) [C_{n p} (α) + δ C_{n_{p_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] \\ - C_{Y T} [x_{c g r} - x_{c g}] \frac{\bar{C}}{b} + δ C_{n_{β}} (α) β . \end{matrix}

The moment coefficients in the above formula can be directly read from the aerodynamic coefficient table of F-16 with the following definitions.

●: $δ C_{l L E F} = C_{l L E F} (α, β) - {C_{l} (α, β, δ_{e})}_{δ_{e} = 0}$
●: $δ C_{l_{δ_{a}}} = C_{l_{δ_{a}}} (α, β) - {C_{l} (α, β, δ_{e})}_{δ_{e} = 0}$
●: $δ C_{l_{{δ_{a}}_{L E F}}} = C_{l_{{δ_{a}}_{L E F}}} (α, β) - C_{l_{L E F}} (α, β) - δ C_{l_{δ_{a}}}$
●: $δ C_{l_{δ_{r}}} = C_{l_{δ_{r}}} (α, β) - {C_{l} (α, β, δ_{e})}_{δ_{e} = 0}$
●: $δ C_{m L E F} = C_{m L E F} (α, β) - {C_{m} (α, β, δ_{e})}_{δ_{e} = 0}$
●: $δ C_{n L E F} = C_{n L E F} (α, β) - {C_{n} (α, β, δ_{e})}_{δ_{e} = 0}$
●: $δ C_{n_{δ_{a}}} = C_{n_{δ_{a}}} (α, β) - {C_{n} (α, β, δ_{e})}_{δ_{e} = 0}$
●: $δ C_{n_{{δ_{a}}_{L E F}}} = C_{n_{{δ_{a}}_{L E F}}} (α, β) - C_{n_{L E F}} (α, β) - δ C_{n_{δ_{a}}}$
●: $δ C_{n_{δ_{r}}} = C_{n_{δ_{r}}} (α, β) - {C_{n} (α, β, δ_{e})}_{δ_{e} = 0}$

The actuators driving the deflections of the aero/thrust control surfaces are modeled by the first-order dynamics with angle and rate limits specified by Table A1. More detailed parameters about F-16 can be found in Ref. [33].

Table A1. The actuator dynamics of the aero surface deflections

(δ_{a}, δ_{e}, δ_{r})

and the nozzle deflections

(δ_{y}, δ_{z})

with angle and rate limits.

Table A1. The actuator dynamics of the aero surface deflections

(δ_{a}, δ_{e}, δ_{r})

and the nozzle deflections

(δ_{y}, δ_{z})

with angle and rate limits.

Actuators	Transfer function	Angle limits	Rate limits
Elevator $δ_{e}$	$\frac{20.2}{s + 20.2}$	$\pm 25 °$	$\pm 60 ° / s$
Aileron $δ_{a}$	$\frac{20.2}{s + 20.2}$	$\pm 21.5 °$	$\pm 80 ° / s$
Rudder $δ_{r}$	$\frac{20.2}{s + 20.2}$	$\pm 30 °$	$\pm 120 ° / s$
Thrust $δ_{z}$	$\frac{20.2}{s + 20.2}$	$\pm 15 °$	$\pm 60 ° / s$
Thrust $δ_{y}$	$\frac{20.2}{s + 20.2}$	$\pm 15 °$	$\pm 60 ° / s$

A.3. Equations of Motion in Wind-Axis Coordinates

Eq. (A.3) is the equation of motion of the aircraft's center of mass in body-axis coordinates. Now we convert Eq. (A.3) into the equation of motion in wind-axis coordinates by the following substitutions:

{\hat{V}}_{b} = {[\begin{matrix} u \\ v \\ w \end{matrix}]}_{b} \to {\hat{V}}_{w} = {[\begin{matrix} V_{T} \\ 0 \\ 0 \end{matrix}]}_{w}, {\hat{Ω}}_{E}^{b} = {[\begin{matrix} p \\ q \\ r \end{matrix}]}_{b} \to {\hat{Ω}}_{E}^{w} = {[\begin{matrix} - \dot{χ} \sin γ \\ \dot{γ} \\ \dot{χ} \cos γ \end{matrix}]}_{w},

(A.11)

where the velocity of the center of mass

{\hat{V}}_{b}

in the body-axis coordinates is replaced by the velocity

{\hat{V}}_{w}

in the wind-axis coordinates, and the angular velocity

{\hat{Ω}}_{E}^{b}

of the body-axis coordinates is replaced by the angular velocity

{\hat{Ω}}_{E}^{w}

of the wind-axis coordinates. After the above substitutions, Eq. (A.3) becomes,

m ({[\begin{matrix} {\dot{V}}_{T} \\ 0 \\ 0 \end{matrix}]}_{w} + [\begin{matrix} 0 & - \dot{χ} \cos γ & \dot{γ} \\ \dot{χ} \cos γ & 0 & \dot{χ} \sin γ \\ - \dot{γ} & - \dot{χ} \sin γ & 0 \end{matrix}] {[\begin{matrix} V_{T} \\ 0 \\ 0 \end{matrix}]}_{w}) = {[\begin{matrix} F_{x} \\ F_{y} \\ F_{z} \end{matrix}]}_{w} .

(A.12)

The above equation can be further simplified to

m [\begin{matrix} {\dot{V}}_{T} \\ V_{T} \dot{χ} \cos γ \\ {- V}_{T} \dot{γ} \end{matrix}] = {[\begin{matrix} F_{x} \\ F_{y} \\ F_{z} \end{matrix}]}_{w} = R_{b}^{w} {[\begin{matrix} F_{x} \\ F_{y} \\ F_{z} \end{matrix}]}_{b} = R_{b}^{w} [\begin{matrix} T_{x} \\ T_{y} \\ T_{z} \end{matrix}] + R_{b}^{w} [\begin{matrix} X \\ Y \\ Z \end{matrix}] + R_{E}^{w} [\begin{matrix} 0 \\ 0 \\ m g \end{matrix}] .

(A.13)

Substituting the rotation matrix

R_{E}^{w}

from Eq. (A.1) and the rotation matrix

R_{b}^{w}

from Eq. (A.2) into the above equation, the governing equations for the aircraft’s velocity components

(V_{T}, χ, γ)

in the wind-axis coordinates can be obtained as shown in Eq. (1) – Eq. (3) in the main text.

Next, we will derive the motion of the aircraft body axis relative to the wind axis. Mathematically, the relative rotation between the three frames can be described by

{\hat{Ω}}_{E}^{b} = {\hat{Ω}}_{w}^{b} + {\hat{Ω}}_{E}^{w},

(A.14)

which means that the angular velocity

{\hat{Ω}}_{E}^{b}

of the body-axis frame relative to the local horizontal frame is equal to the angular velocity

{\hat{Ω}}_{w}^{b}

of the body-axis frame relative to the wind-axis frame plus the angular velocity

{\hat{Ω}}_{E}^{w}

of the wind-axis frame relative to the local horizontal frame. The expressions for

{\hat{Ω}}_{E}^{b}

and

{\hat{Ω}}_{E}^{w}

have been given by Eq. (A.11), and the angular velocity

{\hat{Ω}}_{w}^{b}

can be expressed in terms of

(\dot{α}, \dot{β}, \dot{μ})

{\hat{Ω}}_{w}^{b} = {[\begin{matrix} \dot{μ} \\ \dot{α} + \dot{β} \sin μ \\ - \dot{β} \cos μ \end{matrix}]}_{w} .

(A.15)

Substituting the three angular velocities into Eq. (A.14), we get the following result,

{[\begin{matrix} p \\ q \\ r \end{matrix}]}_{b} = R_{w}^{b} {([\begin{matrix} \dot{μ} \\ \dot{α} + \dot{β} \sin μ \\ - \dot{β} \cos μ \end{matrix}] + [\begin{matrix} - \dot{χ} \sin γ \\ \dot{γ} \\ \dot{χ} \cos γ \end{matrix}])}_{w} .

(A.6)

The rotation matrix

R_{w}^{b}

is introduced because the angular velocity vectors

{\hat{Ω}}_{w}^{b}

and

{\hat{Ω}}_{E}^{w}

in (A.14) are expressed in the wind-axis coordinate system, while the vector

{\hat{Ω}}_{w}^{b}

is expressed in the body-axis coordinate system. Substituting the rotation matrix

R_{w}^{b}

from Eq. (A.2) into Eq. (A.16) and using Eq. (2) for

\dot{χ}

and Eq. (3) for

\dot{γ}

, we obtain the expression for

(\dot{α}, \dot{β}, \dot{μ})

as shown in Eq. (6) – Eq. (8).

A.4. Transformation into Input-Affine Nonlinear Form

The rationality of the affine nonlinear form is based on the fact that the deflection angle of the control surface is limited to a small range due to the mechanical constraint. By taking the linear expansion of the aerodynamic moments

{[L, M, N]}^{T}

given by Eq. (A.10) with respect to the control input

u = {[δ_{a}, δ_{e}, δ_{r}, δ_{y}, δ_{z}]}^{T}

and the linear expansion of the thrust moments given by Eq. (13), we can transform Eq. (A.4) – Eq. (A.6) into the following input-affine nonlinear form:

[\begin{matrix} \dot{p} \\ \dot{q} \\ \dot{r} \end{matrix}] = f_{1} (x) + g_{1} (x) u = [\begin{matrix} f_{p} (x) \\ f_{q} (x) \\ f_{r} (x) \end{matrix}] + [\begin{matrix} g_{p_{δ_{a}}} & 0 & g_{p_{δ_{r}}} & g_{p_{δ_{y}}} & 0 \\ 0 & g_{q_{δ_{e}}} & 0 & 0 & g_{q_{δ_{z}}} \\ g_{r_{δ_{a}}} & 0 & g_{r_{δ_{r}}} & g_{r_{δ_{y}}} & 0 \end{matrix}] [\begin{matrix} δ_{a} \\ δ_{e} \\ δ_{r} \\ δ_{y} \\ δ_{z} \end{matrix}] .

(A.17)

The nonlinear functions in

f_{1} (x)

, which are independent of the control surface deflections, are found as

f_{p} (x) = (C_{1} r + C_{2} p) q + C_{3} \hat{l} + C_{4} \hat{n},

f_{q} (x) = C_{5} p r - C_{6} (p^{2} - r^{2}) + C_{7} \hat{m},

f_{r} (x) = (C_{8} p - C_{2} r) q + C_{4} \hat{l} + C_{9} \hat{n},

where

\hat{l}

m

, and

\hat{n}

are control-independent moments defined by

\hat{l} = \bar{q} S b \{C_{l} (α, β) + C_{l_{L E F}} (1 - \frac{δ_{L E F}}{25}) + (\frac{r b}{2 V_{T}}) [C_{l r} (α) + δ C_{l_{r_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + (\frac{p b}{2 V_{T}}) [C_{l p} (α) + δ C_{l_{p_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + δ C_{l_{β}} (α) β\},

\hat{m} = \bar{q} S \bar{C} \{C_{m} (α, β) + C_{Z T} [x_{c g r} - x_{c g}] + δ C_{m_{L E F}} (1 - \frac{δ_{L E F}}{25}) + \frac{q \bar{C}}{2 V_{T}} [C_{m q} (α) + δ C_{m_{q_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + δ C_{m} (α)\},

\hat{n} = \bar{q} S b \{C_{n} (α, β) + C_{n L E F} (1 - \frac{δ_{L E F}}{25}) + (\frac{r b}{2 V_{T}}) [C_{n r} (α) + δ C_{n_{r_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + (\frac{p b}{2 V_{T}}) [C_{n p} (α) + δ C_{n_{P_{L E F}}} (α) (1 - \frac{δ_{L E F}}{25})] + δ C_{n_{β}} (α) β - C_{Y T} (x_{c g r} - x_{c g}) \frac{\bar{C}}{b}\} .

The matrix

g_{1} (x)

contains the linear expansion coefficients of the aerodynamic/thrust moments with respect to the five control surface deflections. The elements in the first three columns of

g_{1} (x)

are the expansion coefficients regarding the aero surface deflections and are found as

g_{p_{δ_{a}}} = \bar{q} S b \{C_{3} [δ C_{l_{δ_{a}}} + δ C_{l_{{δ_{a}}_{L E F}}} (1 - \frac{δ_{L E F}}{25})] (\frac{1}{21.5}) + C_{4} [δ C_{n_{δ_{a}}} + δ C_{n_{{δ_{a}}_{L E F}}} (1 - \frac{δ_{L E F}}{25})] (\frac{1}{21.5})\},

g_{p_{δ_{r}}} = \bar{q} S b [C_{3} δ C_{l_{δ_{r}}} (\frac{1}{30}) + C_{4} δ C_{n_{δ_{r}}} (\frac{1}{30})],

g_{q_{δ_{e}}} = \bar{q} S \bar{C} C_{7} {\hat{m}}_{δ_{e}},

g_{r_{δ_{a}}} = \bar{q} S b \{C_{4} [δ C_{l_{δ_{a}}} + δ C_{l_{{δ_{a}}_{L E F}}} (1 - \frac{δ_{L E F}}{25})] (\frac{1}{21.5}) + C_{9} [δ C_{n_{δ_{a}}} + δ C_{n_{{δ_{a}}_{L E F}}} (1 - \frac{δ_{L E F}}{25})] (\frac{1}{21.5})\},

g_{r_{δ_{r}}} = \bar{q} S b [C_{4} δ C_{l_{δ_{r}}} (\frac{1}{30}) + C_{9} δ C_{n_{δ_{r}}} (\frac{1}{30})] .

The elements in the last two columns of

g_{1} (x)

are the expansion coefficients regarding the nozzle deflections and are found as

g_{p_{δ_{y}}} = - C_{4} l_{T} T, g_{q_{δ_{z}}} = - C_{7} l_{T} T, g_{r_{δ_{y}}} = - C_{9} l_{T} T .

References

Herrick, P.W. Propulsion influences on air combat. AIAA-85-1457, July 1985.
Costes, P. Investigation of thrust vectoring and post-stall capability in air combat. AIAA-88-4160-CP, 1988.
Nguyen, L.; Gilert, W. Impact of emerging technologies on future combat aircraft agility. AIAA-90-1304, May 1990.
Day, R.E. Coupling dynamics in aircraft: a historical perspective. Technical Report NASA-SP-532, NASA Dryden Flight Research Center Edwards, CA United States, March 1997.
Pahle, J.W.; Wichman, K.D.; Foster, J.V.; Bundick, W.T. An overview of controls and flying qualities technology on the F/A-18 high alpha research vehicle. Technical Report NASA-H-2123, NASA Dryden Flight Research Center Edwards, January 1996.
Nguyen, L.T.; Ogburn, M.E.; Gilbert, W.P.; Kibler, K.S.; Brown, P.W.; Deal, P.L. Simulator study of stall/post-stall characteristics of a fighter airplane with relaxed longitudinal static stability. Technical Report NASA-TP-1538, NASA Langley Research Center Hampton, VA, United States, December 1979.
Ikaza, D. Thrust vectoring nozzle for military aircraft engines. 22nd International Congress of Aeronautical Sciences 2000, No. 5.3.R1, p. 534.
Capone, F.J.; Mason, M.L.; Leavitt, L.D. An experimental investigation of thrust vectoring two-dimensional convergent-divergent nozzles installed in a twin-engine fighter model at high angles of attack. Technical Report NASA-TM-4155, NASA Langley Research Center Hampton, VA, United States, February 1990.
Huber, P.; Seamount, P. X-31 High Angle of Attack Control System Performance. Fourth High Alpha Conference 1994, vol. 2, NASA CP-10143, NASA Dryden Flight Research Center.
Canter, D. X-31 Post-Stall Envelope Expansion and Tactical Utility Testing. Fourth High Alpha Conference 1994, vol. 2, NASA CP-10143, NASA Dryden Flight Research Center.
Canter D. E.; Groves, A.W. X-31 post-stall envelope expansion and tactical utility testing. Biennial Flight Test Conference 1994, AIAA Paper 1994-2171, Hilton Head, SC, U.S.A.
Asbury, S.C.; Capone, F.J. Multi-axis Thrust-Vectoring Characteristics of a Model Representative of the F-18 High-Alpha Research Vehicle at Angles of Attack from 0 deg to 70 deg. Technical Report NASA-TP-3531, NASA Langley Research Center Hampton, VA, United States, December 1995.
Chen, H. Effectiveness of Thrust Vectoring Control for Longitudinal Trim of a Blended Wing Body Aircraft. Master Thesis, Delft University of Technology, 2015.
Brinker, J.; Wise, K. Flight Testing of a Reconfigurable Flight Control Law on the X-36 Tailless Fighter Aircraft. Journal of Guidance, Control, and Dynamics 2001, 24, pp. 903-909.
Balas, G.J. Flight control law design: An industry perspective. European Journal of Control 2003, 9, pp. 207-226. [CrossRef]
Adams, R.J.; Buffington, J.M.; Sparks, A.G.; Banda, S.S. Robust multivariable flight control. Springer-Verlag, London, 1994.
Wei, J.; Chen, H. Designing backstepping control system for hypersonic vehicle based on feedback linearization, International Journal of Aerospace Engineering 2015, Article ID 916328, 11 pages, . [CrossRef]
Li, C.Y.; Jing, W.X.; Gao, C.S. Adaptive Backstepping-based Flight Control System Using Integral Filters. Aerospace Science and Technology 2009, 13, pp. 105-113. [CrossRef]
Hess, R.A.; Wells, S.R. Sliding mode control applied to reconfigurable flight control design. Journal of Guidance, Control, and Dynamics 2003, 26, pp. 452-462. [CrossRef]
Rysdyk, R.; Calise, A.J. Robust Nonlinear Adaptive Flight Control for Consistent Handling Qualities. IEEE Transactions on Control Systems Technology 2005, 13, pp. 896-910. [CrossRef]
Balas, G.G.; Hodgkinson, J. Control design methods for good flying qualities. AIAA Atmospheric Flight Mechanics Conference 2009, AIAA paper 2009-6319, Chicago.
Walker, G.P.; Allen, D.A. X-35B STOVL Flight control law design and flying qualities. Biennial International Powered Lift Conference and Exhibit 2002, AIAA Paper 2002-6018, Williamsburg, Virginia.
Albostan, O.; Gökaşan, M. High Angle of Attack Maneuvering Control of F-16 Aircraft Based on Nonlinear Dynamic Inversion and Eigen-structure Assignment, European Conference for Aeronautics and Space Sciences, 2017.
Snell, S.A.; Enns, D.F.; Garrard, W.L. Nonlinear inversion flight control for a supermaneuverable aircraft. Journal of guidance, control, and dynamics 1992, 15, pp. 976-984. [CrossRef]
Miller, C.J. Nonlinear Dynamic Inversion Baseline Control Law: Architecture and Performance Predictions, AIAA Guidance, Navigation, and Control Conference, AIAA Paper no. 2011-6467, Portland Oregon, 2011.
Miller, C.J. Nonlinear dynamic inversion baseline control law: Flight-test results for the full-scale advanced systems testbed F/A-18 airplane. AIAA Guidance, Navigation, and Control Conference, AIAA Paper no. 2011-6468, Portland Oregon, 2011.
Lubas, C.; Paglione, P. Nonlinear dynamic inversion applied to a F-16 aircraft model. 20th International Congress of Mechanical Engineering, Gramado, Brazil, 2009.
Jia, Q.; Zhang, W.; Shi, J.; Hu, J. Maneuverable aircraft flight control using nonlinear dynamic inversion. 18th International Conference on Control, Automation and Systems, Gangwon, Korea, 2018.
Adams, R.J.; Banda, S.S. Robust Flight Control Design Using Dynamic Inversion and Structured Singular Value Synthesis. IEEE Transactions on Control Systems 1993, 1, pp. 80-92. [CrossRef]
Durham, W.; Bordignon, K.A.; Beck, R. Aircraft control allocation; John Wiley & Sons, 2017.
Johansen, T.A.; Fossen, T.I. Control allocation: A survey. Automatica 2013, 49, pp. 1087-1103.
Well, K.H.; Faber, B.; Berger, E. Optimale taktische Flugmanoever fuer ein kamfflugzeug der 90er Jahre. Interner Bericht A-52-79/6, DFVLR, Germany, Oct. 1979.
Sonneveldt, L. Nonlinear F-16 model description, Delft University of Technology, Netherlands, 2006.

Figure 1. The interchange between the AOA

α

and the side slip angle

β

via rolling around the body axis

{\hat{x}}_{b}

by an angle

ϕ = 90 °

Figure 1. The interchange between the AOA

α

and the side slip angle

β

via rolling around the body axis

{\hat{x}}_{b}

by an angle

ϕ = 90 °

Figure 2. (a, c) The relative orientation between body-axis coordinates

({\hat{x}}_{b}, {\hat{y}}_{b}, {\hat{z}}_{b})

, wind-axis coordinates

({\hat{x}}_{w}, {\hat{y}}_{w}, {\hat{z}}_{w})

and local horizontal coordinates

({\hat{x}}_{E}, {\hat{y}}_{E}, {\hat{z}}_{E})

. (b) Multi-axis thrust vectoring nozzle with deflection

δ_{y}

along the

y_{b}

axis and deflection

δ_{z}

along the

z_{b}

axis.

Figure 2. (a, c) The relative orientation between body-axis coordinates

({\hat{x}}_{b}, {\hat{y}}_{b}, {\hat{z}}_{b})

, wind-axis coordinates

({\hat{x}}_{w}, {\hat{y}}_{w}, {\hat{z}}_{w})

and local horizontal coordinates

({\hat{x}}_{E}, {\hat{y}}_{E}, {\hat{z}}_{E})

. (b) Multi-axis thrust vectoring nozzle with deflection

δ_{y}

along the

y_{b}

axis and deflection

δ_{z}

along the

z_{b}

axis.

Figure 3. The computation of forces and moments exerting on the aircraft. The surface deflections

[δ_{a}, δ_{e}, δ_{r}]

generated by the actuators produce the aerodynamic force

F^{(A)}

given by Eq. (4) and moment

M^{(A)}

given by Eq. (12). The thrust deflections

[δ_{z}, δ_{y}]

produce the thrust force

F^{(T)}

given by Eq. (5) and moment

M^{(T)}

given by Eq. (13). The gravity force

F^{(G)}

is determined by the mass and the attitude of the aircraft.

Figure 3. The computation of forces and moments exerting on the aircraft. The surface deflections

[δ_{a}, δ_{e}, δ_{r}]

generated by the actuators produce the aerodynamic force

F^{(A)}

given by Eq. (4) and moment

M^{(A)}

given by Eq. (12). The thrust deflections

[δ_{z}, δ_{y}]

produce the thrust force

F^{(T)}

given by Eq. (5) and moment

M^{(T)}

given by Eq. (13). The gravity force

F^{(G)}

is determined by the mass and the attitude of the aircraft.

Figure 4. Block diagram of the NDI control law design. An ideal NDI controller

u

is expected to make the inner closed-loop system (surrounded by the red dashed box) behave like a linear system

\dot{x} = v

Figure 4. Block diagram of the NDI control law design. An ideal NDI controller

u

is expected to make the inner closed-loop system (surrounded by the red dashed box) behave like a linear system

\dot{x} = v

Figure 5. The NDI flight control structure can be separated into three layers of loop. The inner loop regards the angular rate control of the fast states

{[p, q, r]}^{T}

, the outer loop regards the attitude control of the slow states

{[α, β, μ]}^{T}

, and the outermost loop regards the flight-path control of the slowest states

{[V_{T}, γ, χ]}^{T}

Figure 5. The NDI flight control structure can be separated into three layers of loop. The inner loop regards the angular rate control of the fast states

{[p, q, r]}^{T}

, the outer loop regards the attitude control of the slow states

{[α, β, μ]}^{T}

, and the outermost loop regards the flight-path control of the slowest states

{[V_{T}, γ, χ]}^{T}

Figure 6. The inner angular-rate control loop is designed to determine the three aero deflections

[δ_{a}, δ_{e}, δ_{r}]

and two nozzle deflections

[δ_{z}, δ_{y}]

such that the angular-rate response

{[p, q, r]}^{T}

of the aircraft matches the desired response

{[p_{c}, q_{c}, r_{c}]}^{T}

Figure 6. The inner angular-rate control loop is designed to determine the three aero deflections

[δ_{a}, δ_{e}, δ_{r}]

and two nozzle deflections

[δ_{z}, δ_{y}]

such that the angular-rate response

{[p, q, r]}^{T}

of the aircraft matches the desired response

{[p_{c}, q_{c}, r_{c}]}^{T}

Figure 7. The outer attitude control loop is designed to generate the angular-rate command

{[p_{c}, q_{c}, r_{c}]}^{T}

for the inner loop according to the tracking error between the attitude command

{[μ_{c}, α_{c}, β_{c}]}^{T}

and the actual attitude

{[μ, α, β]}^{T}

. The yellow block represents the inner NDI control loop shown in Figure 6.

Figure 7. The outer attitude control loop is designed to generate the angular-rate command

{[p_{c}, q_{c}, r_{c}]}^{T}

for the inner loop according to the tracking error between the attitude command

{[μ_{c}, α_{c}, β_{c}]}^{T}

and the actual attitude

{[μ, α, β]}^{T}

. The yellow block represents the inner NDI control loop shown in Figure 6.

Figure 8. The time responses of the attitude

[μ (t), α (t), β (t)]

and the control surface deflections

[δ_{a} (t), δ_{e} (t), δ_{r} (t)]

of the F-16 flying below stall AOA without TVC.

Figure 8. The time responses of the attitude

[μ (t), α (t), β (t)]

and the control surface deflections

[δ_{a} (t), δ_{e} (t), δ_{r} (t)]

of the F-16 flying below stall AOA without TVC.

Figure 9. The time responses of attitude

[μ (t), α (t), β (t)]

and control surface deflections

[δ_{a 0} (t), δ_{e} (t), δ_{r} (t)]

for the experimental F-16 with and without TVC when the AOA command exceeds the stall limit. The black dotted curves are the attitude commands, the blue dash dotted curves denote the responses of the F-16 without TVC, and the red solid curves denote the responses of the F-16 installed with TVC.

Figure 9. The time responses of attitude

[μ (t), α (t), β (t)]

and control surface deflections

[δ_{a 0} (t), δ_{e} (t), δ_{r} (t)]

Figure 10. The time responses of the Herbst maneuver for the F-16 installed with TVC and controlled by the NDI controller. The attitude responses show that the F-16 completed the Herbst maneuver in 20 seconds via a rolling angle

μ = 120 °

around the wind axis at

α = 70 °

, while keeping the side slip angle

β (t)

as small as possible. With the help of TVC, all three aero surfaces operate below their saturation values except for some short periods of time.

μ = 120 °

around the wind axis at

α = 70 °

, while keeping the side slip angle

β (t)

as small as possible. With the help of TVC, all three aero surfaces operate below their saturation values except for some short periods of time.

Figure 11. The subgraphs (a) and (b) show the deflection angles of the thrust below their saturation values during the operation of the Herbst maneuver. Subgraph (c) is the total speed profile of F-16 starting from 250 ft/s and accelerating to 300 ft/s to enter the Herbst maneuver, during which the F-16 first decelerates with increasing AOA and reaches its minimum speed of 130 ft/s near

t = 10 s

to reverse the heading, and then accelerates to recover to the level flight state. Subgraphs (d), (e) and (f) shows the trajectory of the F-16 performing the Herbst maneuver, indicating that the height of the F-16 drops about 1000ft, and the turning radius on the horizontal plane is about 2000ft.

t = 10 s

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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

Abstract

1. Introduction

2. Equations of Motion in Wind-Axis Coordinates

3. NDI Flight Control Law Design

4. Angular-Rate NDI Control

5. Attitude NDI Control

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

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

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

6.3. Herbst Maneuver with NDI Flight Control

7. Conclusions

Author Contributions

Funding

Appendix A. The Input-Affine Nonlinear Form of F-16 with TVC

A.1. Rotation Matrix between Three Coordinate Systems

A.2. Equations of Motion in Body-axis Coordinates

A.3. Equations of Motion in Wind-Axis Coordinates

A.4. Transformation into Input-Affine Nonlinear Form

References

MDPI Initiatives

Important Links

Subscribe