2.1 Controller architecture

Figure 2 presents the controller architecture for the pitch rate control. The control elements include a baseline PID controller, an RLS-based DI controller and a NN adaptive controller. The stability analysis is detailed in the appendix.

Figure 2. Pitch rate model reference adaptive controller using PID, RLS dynamic inversion and adaptive NN.

Controlling the pitch rate $q$ for longitudinal control aims to track the desired short period performances [Reference Stevens, Lewis and Johnson3]. Therefore, the proposed control algorithm is designed to track a reference signal  ${q_{ref}}$ .

2.2 Controller elements

Nonlinear Dynamic Inversion Controller

For linear systems, DI is straightforward to execute [Reference Stevens, Lewis and Johnson3]. For nonlinear systems, the inversion is more difficult, as nonlinear systems may include unmodeled dynamics. According to Ref. [Reference Stevens, Lewis and Johnson3], most aircraft dynamics can be represented in affine form, with respect to the control input. Therefore, Equation (1) can be rewritten as:

(2) \begin{align}\begin{array}{l}{\dot {\mathbf{x}}} = {\mathbf{f}}( {\mathbf{x}} ) + {\mathbf{g}}( {\mathbf{x}} )\eta \\y = q = {\mathbf{Cx}}\end{array} \end{align}

where ${\rm{f}}$ and ${\mathbf{g}}$ are nonlinear functions, and ${\rm{y}}$ is the output. Hence, $\rm{C} = [ 0\quad {}0\quad {}1\quad {}0]$ for pitch rate output.

Inverting Equation (1), the control input may be written as:

(3) \begin{align}{\rm{\eta }} = \aleph ( {\dot {\mathbf{x}},\mathbf{x}} ) \end{align}

where $\aleph = {{\unicode{x1D6D7}}^{ - 1}}$ is the inverse function of ${\unicode{x1D6D7}}$ . Next, let’s express ${\rm{\eta }} = {{\rm{\delta }}_{\rm{e}}}$ according to the control algorithm in Fig. 2. Considering Equations (1) and (2):

(4) \begin{align}{\dot {\mathbf{x}}} = {\unicode{x1D6D7}}( {{\mathbf{x}},{\rm{\eta }}} ) = {\rm{f}}( {\mathbf{x}} ) + {\mathbf{g}}( {\mathbf{x}} ){\rm{\eta }} \end{align}

Developing Equation (4) using the first-order linear approximation around ( ${{\mathbf{x}}_0},{\eta _0})$ , we obtain:

(5) \begin{align}{\unicode{x1D6D7}}( {{\mathbf{x}},{\rm{\eta }}} ) = {\unicode{x1D6D7}}( {{{\mathbf{x}}_0},{{\rm{\eta }}_0}} ) + {{{\unicode{x1D6D7}'}}_{\rm{x }}}{( {{\mathbf{x}},{\rm{\eta }}} )_{{{\mathbf{x}}_0},{{\rm{\eta }}_0}}}( {{\mathbf{x}} - {{\mathbf{x}}_0}} ) + {{\unicode{x1D6D7}^\prime}_{\rm{\eta }}}{( {{\mathbf{x}},{\rm{\eta }}} )_{{{\mathbf{x}}_0},{{\rm{\eta }}_0}}}( {{\rm{\eta }} - {{\rm{\eta }}_0}} ) \end{align}

where ${{\unicode{x1D6D7}}}( {{{\mathbf{x}}_0},{\eta _0}} )$ is the value at ( ${{\mathbf{x}}_{\mathbf{0}}},{\eta _0})$ , ${{{\unicode{x1D6D7}'}}_x}{( {{\mathbf{x}},\eta } )_{{{\mathbf{x}}_0},{\eta _0}}}$ and ${{{\unicode{x1D6D7}'}}_\eta }{( {{\mathbf{x}},\eta } )_{{{\mathbf{x}}_0},{\eta _0}}}$ are the partial derivatives of ${{\unicode{x1D6D7}}}$ with respect to ${\mathbf{x}}$ and $\eta $ , respectively. We may formulate the terms in Equation (5), such as:

(6.1) \begin{align}{{\unicode{x1D6D7}}}( {{{\mathbf{x}}_0},{\eta _0}} ) = {\mathbf{f}}( {{{\mathbf{x}}_0}} ) + {\mathbf{g}}( {{{\mathbf{x}}_0}} )\eta \end{align}

(6.2) \begin{align} {{{\unicode{x1D6D7}'}}_x}( {{\mathbf{x}},\eta } ) = \frac{{\partial {\mathbf{f}}( {\mathbf{x}} )}}{{\partial {\mathbf{x}}}} + \frac{{\partial {\mathbf{g}}( {\mathbf{x}} )}}{{\partial {\mathbf{x}}}}\eta \end{align}

(6.3) \begin{align}{{{\unicode{x1D6D7}'}}_\eta }( {{\mathbf{x}},\eta} ) = \frac{{\partial {\rm{f}}( {\mathbf{x}} )}}{{\partial {\rm{\eta }}}} + \frac{{\partial \left[ {{\mathbf{g}}( {\mathbf{x}} ){\rm{\eta }}} \right]}}{{\partial {\rm{\eta }}}} = {\mathbf{g}}( {\mathbf{x}} ) \end{align}

Considering the following changes, so that ${\boldsymbol{\Delta }}{\mathbf{x}} = \left( {{\mathbf{x}} - {{\mathbf{x}}_0}} \right){\rm{\;}}$ and ${\boldsymbol{\Delta }}\eta = \left( {\eta - {\eta _0}} \right)$ , and substituting Equation (6) into Equation (5), the following expression is obtained:

(7) \begin{align}{{\unicode{x1D6D7}}}( {{\mathbf{x}},\eta } ) = {\mathbf{f}}( {{{\mathbf{x}}_0}} ) + {\mathbf{g}}( {{{\mathbf{x}}_0}} )\eta + {\left( {\frac{{\partial {\mathbf{f}}( {\mathbf{x}} )}}{{\partial {\mathbf{x}}}} + \frac{{\partial {\mathbf{g}}( {\mathbf{x}} )}}{{\partial {\mathbf{x}}}}\eta } \right)_{{{\mathbf{x}}_0},{\eta _0}}}{\mathbf{\Delta x}} + {\mathbf{g}}{( {\mathbf{x}} )_{{{\mathbf{x}}_0},{\eta _0}}}{\boldsymbol{\Delta }}\eta \end{align}

Let us express $\dot x$ in Equation (4) for small perturbations around $\left( {{{\mathbf{x}}_0},{\rm{\;}}{\eta _0}} \right)$ :

(8) \begin{align}{\dot {\mathbf{x}}} = {{\dot {\mathbf{x}}}_0} + {{\Delta \dot {\mathbf{x}}}} = {{\unicode{x1D6D7}}}\left( {{\mathbf{x}},\eta } \right) \end{align}

From Equations (7) and (8), and considering ${{\dot {\mathbf{x}}}_0} = {\mathbf{f}}( {{{\mathbf{x}}_0}} ) + {\mathbf{g}}( {{{\mathbf{x}}_0}} )\eta $ , ${\boldsymbol{\Delta }}{\dot {\mathbf{x}}}$ can be written such as:

(9) \begin{align}{\boldsymbol{\Delta }}{\dot {\mathbf{x}}} = {\left( {\frac{{\partial {\mathbf{f}}( {\mathbf{x}} )}}{{\partial {\mathbf{x}}}} + \frac{{\partial {\mathbf{g}}( {\mathbf{x}} )}}{{\partial {\mathbf{x}}}}\eta } \right)_{{{\mathbf{x}}_0},{\eta _0}}}{\boldsymbol{\Delta }}{\mathbf{x}} + {\mathbf{g}}{( {\mathbf{x}} )_{{{\mathbf{x}}_0},{\eta _0}}}{\boldsymbol{\Delta }}\eta \end{align}

Making the following changes:

(10) \begin{align}{{\mathbf{F}}_{{{\mathbf{x}}_0},{\eta _0}}} = {\left( {\frac{{\partial {\mathbf{f}}\left( {\mathbf{x}} \right)}}{{\partial {\mathbf{x}}}} + \frac{{\partial {\mathbf{g}}\left( {\mathbf{x}} \right)}}{{\partial {\mathbf{x}}}}\eta } \right)_{{{\mathbf{x}}_0},{\eta _0}}};{\rm{\;}}{{\mathbf{G}}_{{{\mathbf{x}}_0},{\eta _0}}} = {\mathbf{g}}{\left( {\mathbf{x}} \right)_{{{\mathbf{x}}_0},{\eta _0}}} \end{align}

Equation (9) can be expressed in a simplified form:

(11) \begin{align}{\boldsymbol{\Delta }}{\dot {\mathbf{x}}} = {{\mathbf{F}}_{{{\mathbf{x}}_0},{\eta _0}}}{\mathbf{\Delta x}} + {{\mathbf{G}}_{{{\mathbf{x}}_0},{\eta _0}}}{\boldsymbol{\Delta }}\eta \end{align}

Recalling from Equation (2), ${\boldsymbol{\Delta }}\dot y$ may be computed as follows:

(12) \begin{align}{{\Delta \dot {\mathbf{y}}}} = {{C\Delta \dot {\mathbf{x}}}} \end{align}

Substituting ${{\Delta \dot {\mathbf{x}}}}$ from Equation (11) into Equation (12), we obtain:

(13) \begin{align}{{\Delta \dot {\mathbf{y}}}} = {\mathbf{C}}{{\mathbf{F}}_{{{\mathbf{x}}_0},{{\rm{\eta }}_0}}}{{\boldsymbol{\Delta}} \mathbf{x}} + {\mathbf{C}}{{\mathbf{G}}_{{{\mathbf{x}}_0},{{\rm{\eta }}_0}}}{\rm{\Delta \eta }} \end{align}

From Equation (13), ${\boldsymbol{\Delta }}\eta $ may be obtained by inversion as follows:

(14) \begin{align}{\boldsymbol{\Delta }}\eta = {\left( {{\mathbf{C}}{{\mathbf{G}}_{{{\mathbf{x}}_0},{\eta _0}}}} \right)^{ - 1}}\left( {{\rm{\;\Delta }}\dot y - {\mathbf{C}}{{\mathbf{F}}_{{{\mathbf{x}}_0},{\eta _0}}}{\boldsymbol{\Delta} \mathbf{x}}} \right) \end{align}

The control task consists of tracking a reference signal $r$ . Consequently, the error dynamics $e$ may be determined as follows:

(15.1) \begin{align}\quad\, e = r - y \end{align}

(15.2) \begin{align} \Rightarrow {\rm{\;}}\dot e = \dot r - \dot y \end{align}

where $y{\rm{\;}}$ is the output variable.

Taking $\dot y \approx {\boldsymbol{\Delta }}\dot y$ for small discrete time steps and recalling from Equation (15.2), for which $\dot y = \dot r - \dot e$ , then ${\boldsymbol{\Delta }}\eta $ may be written as such:

(16) \begin{align}{\boldsymbol{\Delta }}\eta = {\left( {{\mathbf{C}}{{\mathbf{G}}_{{{\mathbf{x}}_0},{\eta _0}}}} \right)^{ - 1}}\left( {{\rm{\;}}\dot r + \nu - {\mathbf{C}}{{\mathbf{F}}_{{{\mathbf{x}}_{\mathbf{0}}},{\eta _0}}}{\mathbf{\Delta x}}} \right) \end{align}

where $\nu $ is chosen so that the following relationship is fulfilled [Reference Stevens, Lewis and Johnson3]:

(17) \begin{align}\dot e = - \nu \end{align}

Thus, the error can be stabilised by choosing a suitable form of $\nu $ .

Online parameter estimation

Finding analytical solutions of ${{\mathbf{F}}_{{{\mathbf{x}}_0},{\eta _0}}}$ and ${{\mathbf{G}}_{{{\mathbf{x}}_0},{\eta _0}}}{\rm{\;}}$ given in Equation (10) can be tasking. Using real-time estimations, such as the RLS method, is an alternative. The goal is to approximate ${{\mathbf{F}}_{{{\mathbf{x}}_0},{\eta _0}}}$ and ${{\mathbf{G}}_{{{\mathbf{x}}_0},{\eta _0}}}$ based on local state observations. The RLS formulation can be found in Ref. [Reference Haykin15]. This technique allows us to estimate the states and control variations based on their increments [Reference Mahadi, Ballal, Moinuddin and Al-Saggaf16, Reference Mohseni and Bernstein17].

Transforming Equation (13) into discrete times yields the following representation:

(18) \begin{align} {\boldsymbol{\Delta }}{y_{k + 1}} = {\mathbf{C}}{{\tilde {\mathbf{F}}}_k}{\boldsymbol{\Delta }}{{\mathbf{x}}_k} + {\mathbf{C}}{{\tilde {\mathbf{G}}}_k}{\boldsymbol{\Delta }}{\eta _k} \end{align}

where ${{\tilde {\mathbf{F}}}_k} $ and ${{\tilde {\mathbf{G}}}_k} $ are discrete form of ${{\mathbf{F}}_{{{\mathbf{x}}_0},{\eta _0}}}$ and ${{\mathbf{G}}_{{{\mathbf{x}}_0},{\eta _0}}}$ .

${\boldsymbol{\Delta }}{y_{k + 1}}{\rm{\;}}$ corresponds to the output variation, and depends on the current states increment ${\boldsymbol{\Delta }}{{\mathbf{x}}_k}$ and input increment ${\boldsymbol{\Delta }}{\eta _k}$ . ${{\tilde {\mathbf{F}}}_k}$ and ${{\tilde {\mathbf{G}}}_k}$ may then be replaced by their estimates ${{\hat {\mathbf{F}}}_k}$ and ${{\hat {\mathbf{G}}}_k}$ so that:

(19) \begin{align} {\boldsymbol{\Delta }}{\hat y_{k + 1}} = {\mathbf{C}}[ {\begin{array}{*{20}{c}}{{\boldsymbol{\Delta }}{{\mathbf{x}}_k}}\, {}{{\boldsymbol{\Delta }}{\eta _k}}\end{array}} ]\left[ {\begin{array}{*{20}{c}}{{\hat {\mathbf{F}}}_k^T}\\{{\hat {\mathbf{G}}}_k^T}\end{array}} \right] \end{align}

Let us consider the vector of the estimated parameter ${\rm{\;}}{{\hat {\boldsymbol{\Theta}}}_k} = \left[ {\begin{array}{*{20}{c}}{{\hat {\mathbf{F}}}_k^T}\\{{\hat {\mathbf{G}}}_k^T}\end{array}} \right] $ and the vector of past increments ${{\mathbf{X}}_k} = {\left[ {{\boldsymbol{\Delta }}{{\mathbf{x}}_k}{\rm{\;\Delta }}{\eta _k}} \right]^T}$ . Then, from Equation (19), ${\boldsymbol{\Delta }}{\hat y_{k + 1}}$ can be compacted as follows:

(20) \begin{align}{\boldsymbol{\Delta }}{\hat y_{k + 1}} = {\mathbf{CX}}_k^T{{\hat {\boldsymbol{\Theta}}}_k} \end{align}

The following defines the error matrix ${\epsilon _k}$ as:

(21) \begin{align}{\epsilon _k} = {\boldsymbol{\Delta }}{{\mathbf{x}}_k} - {\boldsymbol{\Delta }}{{\hat {\mathbf{x}}}_k} \end{align}

From Ref. [Reference Haykin15], ${{\hat {\boldsymbol{\Theta}}}_k}$ can be numerically computed such as:

(22) \begin{align} {{\hat {\boldsymbol{\Theta}}}_k} = {{\hat {\boldsymbol{\Theta}}}_{k - 1}} + \frac{{{{\mathbf{L}}_{k - 1{\rm{\;}}}}{{\mathbf{X}}_k}}}{{\kappa + {\mathbf{X}}_k^T{{\mathbf{L}}_{k - 1}}{{\mathbf{X}}_k}}}{\epsilon _k} \end{align}

where the covariance matrix ${{\mathbf{L}}_k}$ can be expressed as follows:

(23) \begin{align}{{\mathbf{L}}_k} = \frac{1}{\kappa }\left( {{{\mathbf{L}}_{k - 1}} - \frac{{{{\mathbf{L}}_{k - 1}}{{\mathbf{X}}_k}{\mathbf{X}}_k^T{{\mathbf{L}}_{k - 1}}}}{{\kappa + {\mathbf{X}}_k^T{{\mathbf{L}}_{k - 1}}{{\mathbf{X}}_k}}}} \right) \end{align}

${{\hat {\boldsymbol{\Theta}}}_k}$ and ${{\mathbf{L}}_k}$ are obtained using the matrix inversion lemma [Reference Haykin15]. ${{\hat {\mathbf{F}}}_k}$ and ${{\hat {\mathbf{G}}}_k}$ entirely depend on the state variations and control inputs. The initial parameters ${{\boldsymbol{\Theta }}_0}$ , ${{\mathbf{L}}_0}$ and forgetting factor $\kappa \;$ must be set beforehand.

Algorithm 1 presents the update of the covariance matrix ${{\mathbf{L}}_k}$ . The covariance matrix is updated when the estimation error ${\epsilon _k}$ reaches a threshold value ${\boldsymbol{\epsilon} _{treshold}}$ .

Algorithm 1 RLS algorithm and covariance matrix update

Proportional Integral Derivative (PID) controller

From Equation (14), the DI requires a linear feedback controller $\nu $ [Reference Stevens, Lewis and Johnson3]. Its value should be chosen so that stability is ensured according to Lyapunov criteria detailed in the appendix. For this reason, a basic PID controller is sufficient and may be expressed in Laplace form such as [Reference Stevens, Lewis and Johnson3]:

(24) \begin{align}\nu = {\rm{\;}}\left( {{K_P} + \frac{{{K_I}}}{s} + {K_D}s} \right)e{\rm{\;}} \end{align}

where ${K_P},{K_I}{\rm{\;}}$ and ${K_D}$ are the PID gains.

Neural network adaptive controller

Using the estimations in Equation (20) frequently results in prediction errors in real-world applications [Reference Ge, Hang, Lee and Zhang24–Reference Emami, Castaldi and Banazadeh26]. Including prediction errors ${\epsilon _1}$ into (14), the input $\Delta {\eta _k}$ can be written as:

(25) \begin{align} \Delta {\eta _k} = {\big( {\mathbf{C}{{\hat {\mathbf{G}}}_k}} \big)^{ - 1}}\big[ {\Delta {y_{k + 1}} - \mathbf{C}{{\hat {\mathbf{F}}}_k}{\boldsymbol{\Delta}} \mathbf{x} - {\epsilon _1}} \big]\; \end{align}

From Equations (16) and (17), $\Delta {\eta _k}$ in Equation (25) can be written:

(26) \begin{align}{\boldsymbol{\Delta }}{\eta _k} = {\big( {{\mathbf{C}}{{{\hat {\mathbf{G}}}}_k}} \big)^{ - 1}}\big[ {\dot r + v - {\mathbf{C}}{{{\hat {\mathbf{F}}}}_k}{\boldsymbol{\Delta }}{\mathbf{x}} - {\epsilon _1}} \big] \end{align}

The variables ${\mathbf{C}}{{\hat {\mathbf{G}}}_k}$ and ${\mathbf{C}}{{\hat {\mathbf{F}}}_k}{\boldsymbol{\Delta }}{\mathbf{x}}$ in Equation (26) represent the aircraft’s local dynamics. To rectify the uncertainties and estimation errors for ${{\hat {\mathbf{F}}}_k}$ and ${{\hat {\mathbf{G}}}_k}$ , the estimation error can be compensated with an adaptive term ${\nu _{ad}}$ so that:

(27) \begin{align}{\epsilon _1} = {\nu _{ad}} \end{align}

An FFNN was used to train ${\nu _{ad}}\;$ as seen in Fig. 3(b).

Figure 3. Representation of a feedforward neural network.

The effectiveness of the FFNN in Fig. 3 is determined by its hyperparameters, including the number of neurons and hidden layers, the activation function $\boldsymbol{\sigma} $ and the weights ${\mathbf{W}}$ and ${\mathbf{V}}$ . The output from the NN may be written such as [Reference Stevens, Lewis and Johnson3]:

(28) \begin{align}{\hat \nu _{ad}} = {{\mathbf{W}}^{\mathbf{T}}}\sigma ( {{{\mathbf{V}}^{\mathbf{T}}}{\bar {\mathbf{x}}}} ){\rm{\;}} \end{align}

where $\sigma \in {{\mathbb E}^{N \times N}}$ , with $N$ – the number of neurons in the hidden layer, ${\bar {\mathbf{x}}} = {[ {{\delta _e}{\rm{\;}}{q_{ref}}{\rm{\;}}q} ]^{\mathbf{T}}}$ , with $ {\bar {\mathbf{x}}} \in {{\mathbb E} ^{m \times 1}}$ denotingthe input vector, and ${\mathbf{W}} \in {{\mathbb E}^{N \times n}}{\rm{\;}}$ and ${\mathbf{V}} \in {{\mathbb E}^{m \times N}}$ , with $n$ representing the number of network outputs. The activation function in the middle layer corresponds to a tangent sigmoid function such as:

(29) \begin{align}{\sigma _i}( z ) = \frac{1}{{1 + {e^{ - z}}{\rm{\;}}}} \end{align}

where $i$ denotes the neuron’s rank. An affine function was used in the input and output layers.

The adaptive signal ${\hat \nu _{ad}}{\rm{\;}}$ is multiplied by a gain ${K_{NN}}$ , which allows us to control the amplitude of the adaptive term ${\hat \nu _{ad}}$ . Consequently:

(30) \begin{align}{\nu _{ad}} = {K_{NN}}{\hat \nu _{ad}} = {K_{NN}}{{\mathbf{W}}^T}\sigma ( {{{\mathbf{V}}^T}{\bar {\mathbf{x}}}} ) \end{align}

Adaptation laws

The following adaptation laws taken from the backpropagation principle [Reference Narendra and Annaswamy38–Reference Ghazi and Botez45, Reference Campbell, Nguyen, Kaneshige and Krishnakumar47–Reference Kim and Calise52] are used to update the NN hyperparameters:

(31) \begin{align}\begin{array}{ccccc}{\dot {\mathbf{W}}} {} = - \left[ {\left( {{\boldsymbol{\sigma }} - {\mathbf{\sigma '}}{{\mathbf{V}}^T}{\bar {\mathbf{x}}}} \right){e^T} + \lambda \left\| e \right\|{\mathbf{W}}} \right]{{\boldsymbol{\Gamma }}_{\rm{W}}}\\[5pt]{\dot {\mathbf{V}}} {} = - {{\boldsymbol{\Gamma }}_{\boldsymbol{V}}}\left[ {{\bar {\mathbf{x}}}{e^T}{{\mathbf{W}}^T}{\boldsymbol{\sigma '}} + \lambda \left\| e \right\|{\mathbf{V}}} \right]\end{array} \end{align}

where $e$ is the tracking error, ${{\boldsymbol{\Gamma }}_{\rm{W}}}$ and ${{\boldsymbol{\Gamma }}_{\mathbf{V}}}$ are learning rate matrices, $\sigma '$ is the derivative of $\sigma $ , and $\lambda $ is a scalar that is equal to 1. In Ref. [Reference Narendra and Annaswamy38], the additional terms $ - \lambda \| e \|{\mathbf{W}}{{\boldsymbol{\Gamma }}_{\rm{W}}}$ and $ - \lambda \left\| e \right\|{\mathbf{V}}{{\boldsymbol{\Gamma }}_{{\mathbf{V}}.}}$ The adaptation law shown in Equation (31) ensures robustness and stability. A fixed learning rate ${\boldsymbol{\Gamma }}$ was considered for diagonal matrices ${{\boldsymbol{\Gamma }}_{\rm{W}}}{\rm{\;}}$ and ${{\boldsymbol{\Gamma }}_{\mathbf{V}}}{\rm{\;}}$ for the input and output weights, respectively. The backpropagation algorithm is considered as the steepest descent approximation when ${\boldsymbol{\Gamma }}$ is set to constant.

2.3 Control law

The overall control law may be obtained by substituting Equations (24), (27) and (30) into Equation (26), such that:

(32) \begin{align}{\boldsymbol{\Delta }}{\eta _k} = {\big( {{\mathbf{C}}{{{\hat {\mathbf{G}}}}_k}} \big)^{ - 1}}\left[ {{\rm{\;}}\dot r + \left( {{K_P} + \frac{{{K_I}}}{s} + {K_D}s} \right)e - {\mathbf{C}}{{{\hat {\mathbf{F}}}}_k}{\boldsymbol{\Delta }}{{\mathbf{x}}_k} - {K_{NN}}{{\mathbf{W}}^T}\sigma \left( {{{\mathbf{V}}^T}{\bar {\mathbf{x}}}} \right)} \right] \end{align}

2.4 Reference model

The proposed algorithm consists of a model reference adaptive controller and controls the pitch rate $q$ , which is a primary state variable for short period motion. The short period and can be approximated by a second-order transfer function:

(33) \begin{align}\frac{{{q_{ref}}}}{{{q_{cmd}}}} = \frac{{\omega _{sp}^2{\rm{\;}}( {1 + {T_{{\theta _2}}}s} )}}{{{s^2} + 2{\xi _{sp}}{\omega _{sp}}s + \omega _{sp}^2}} \end{align}

where $s{\rm{\;}}$ is the Laplace variable. The ${q_{ref}}$ signal consists of the desired second order pitch rate signal, while ${q_{cmd}}$ is a step input signal.

The design requirements for the damping ratio ${\xi _{sp}}$ and natural frequency ${\omega _{sp}}$ are given by the flight qualities in MIL-STD-1797A [36]. Equation (33) is used to model the pitch rate reference model as depicted in Fig. 2. The zero of Equation (33) is defined as $ - 1/{T_{{\theta _2}}}$ , where ${T_{{\theta _2}}}$ is the short period time constant.

2.5 Flight quality requirements

According to the MIL-STD 1797A [36], the Cessna Citation X can be categorised as a Class II aircraft. The proposed longitudinal controller was designed for Category B cruise operations. Table 3 establishes the Level 1 flight quality (FQ) requirements for short period dynamics.

Table 3. Level 1 longitudinal flight quality requirements for category B cruise flights

2.6 Performance metrics

The overall sum squared error denoted as overall SSE was used for optimal tuning of the controller for a given number of flight conditions, so that:

(34) \begin{align}overall{\rm{\;}}SSE = \mathop \sum \limits_{i = 1}^{i = N} SS{E_i} \end{align}

where $N$ is the number of flight conditions.

The SSE for a single flight condition is computed as follows:

(35) \begin{align}SSE = \mathop \sum \limits_{t = 0}^{t = T} {( {{q_t} - {q_{re{f_t}}}} )^2} \end{align}

where ${q_t}$ and ${q_{re{f_t}}}$ correspond to the pitch rate and the pitch rate reference at a given time $t$ , respectively, and $T$ is the simulation time with a sample time ${\boldsymbol{\Delta }}t = 0.02{\rm{\;}}$ sec.

3.1 Controller tuning parameters

The controller in Fig. 2 was tuned through sensitivity-based tuning. First, the aircraft is trimmed for all cruise conditions over the operating envelope, for given altitude $h$ , calibrated airspeed CAS, mass $m$ and centre of gravity ${x_{cg}}$ . Second, the PID gains ${K_P}$ , ${K_I}$ and ${K_D}$ are tuned so that convergence is fulfilled. Third, the RLS inversion controller is tuned, for which the covariance matrix ${{\mathbf{L}}_0}$ , the parameter vector ${{\boldsymbol{\Theta }}_0}$ and the forgetting factor $\kappa $ are initialised. ${{\mathbf{L}}_{\mathbf{k}}}$ is then iteratively updated following Algorithm 1. Then, the neural network learning rate ${\boldsymbol{\Gamma }}$ and gain ${K_{NN}}$ are set. The NN weights are initialised from 0. Finally, robustness tests are performed. The control parameters are then adjusted throughout the process. The overall control parameter values are summed up in Table 4.

Table 4. Design parameters for the pitch rate controller

3.2 Simulation for cruise at 35,000 ft and 290 kn

Figure 4 shows the pitch rate control at 10,668 m [35,000 ft] and 537 km/h [290 kn]. The hyperparameter variations are shown in Fig. 5. Figure 6 shows the RLS parameters, including the covariance matrix ${{\mathbf{L}}_k}$ in Equation (23) and the prediction error ${\epsilon _k}$ in Equation (21). Figure 7 shows the RLS state estimations.

Figure 4. Pitch rate and elevator responses for cruises at 35,000 ft and 290 kt.

Figure 5. Hyperparameter variation during cruise flight.

Figure 6. (a) Covariance matrix ${{\mathbf{L}}_k}$ update and (b) prediction error during cruise.

Figure 7. Longitudinal states estimation using the RLS algorithm for cruise flight at 10,668 m AGL and 537 km/h.

Figure 8. Selected cruise conditions for $m = 15,000{\rm{\;}}$ kg and CoG ${x_{cg}} = 25{\rm{\% }}$ MAC.

3.3 Simulation and validation over the flight envelope

3.3.2 NN and RLS parameter sensitivity analysis

The effect of the NN and RLS parameters on the pitch rate reference tracking are shown in Figs. 9 and 10.

Figure 9. Pitch rate response for different learning rates ${\boldsymbol{\Gamma }}$ and forgetting factors $\kappa {\rm{\;}}$ at all cruise conditions.

Figure 10. Pitch rate response for three different neural networks gains ${\mathbf{K}_{\mathbf{NN}}}$ values for all cruise conditions.

3.3.3 Controller tuning under wind gusts and turbulences

Figure 11 shows the overall SSE in Equation (34) for $N = 10$ flight conditions. The simulations in Fig. 11(b) include wind gusts with ${u_{wind}}{\rm{\;}}$ = 40 m/s [78 kn] and ${w_{wind}} = 5$ m/s [9.7 kn]. The cruise flights in Fig. 11(c) were simulated using the Dryden turbulence model [Reference Stevens, Lewis and Johnson3] with moderate intensity. From Fig. 11(b) and Fig. 11(c), an optimal compromise of ${\boldsymbol{\Gamma }}$ and ${K_{NN}}\;$ were found and summed up in Table 4.

Figure 11. Sum squared error contour plot for 10 flight conditions for (a) ideal conditions, (b) wind gust and (c) moderate turbulence.

Figure 12. Pitch rate response for different configurations for all cruise conditions: (a) PID controller, (b) PID – RLS Dynamic Inversion controller and (c) PID – RLS Dynamic Inversion – NN controller.

3.4 Comparison of controller configurations

In Fig. 12, the tracking error is smaller when the adaptive controllers are successively added. This shows the contribution of each control element and the controller’s ability to adapt across the flight envelope.

3.5 Flight quality validation

The pitch-rate controller was designed to satisfy the FQ requirements specified in MIL-STD-1797A for Class II, Category B aircraft. In this study, the reference model parameters (in Table 5) were selected to introduce a small overshoot of approximately 20% and steady state relative error below 1.2% in the pitch-rate response. This design choice was made to ensure compliance with the dropback requirement specified in Table 3 and defined in MIL-STD-1797A [36], rather than representing a formal Level 1 requirement. The FQ validation was carried out by assessing the short-period damping ratio ${\xi _{sp}}$ , natural frequency ${\omega _{sp}}$ , and control anticipation parameter (CAP) across the flight envelope. As shown in Fig. 13, all identified short-period parameters and CAP values remained within Level 1 boundaries, confirming that the aircraft meets the required flying qualities. The controller’s transient performances for all flight conditions are summarised in Table 6.

Table 5. Pitch rate reference model dynamic characteristics for level 1 flight quality

Table 6. Reference signal and pitch rate dynamic performance results

Figure 13. Short period flight qualities for (a) damping and frequency requirements and (b) CAP requirements.

The aircraft flight qualities were assessed in Fig. 13(a) and Fig. 13(b). The dynamic parameters ${\xi _{sp}}$ , ${\omega _{sp}}$ and ${T_{{\theta _2}}}$ in Equation (33) have been identified for each flight condition using system identification. For each flight condition, optimisation errors less than 2% were obtained, small enough to validate the identification. The CAP in Fig. 13(b) was also determined [Reference Hodgkinson37] and computed as follows:

\begin{align}CAP \approx \frac{{\omega _{sp}^2}}{{\frac{{{V_0}}}{g}\frac{1}{{{T_{{\theta _2}}}}}}} \end{align}

where ${V_0}$ is the airspeed, $g$ is the gravitational constant, equal to 9.81 ${\rm{m}}/{{\mathbf{s}}^2}$ [32.18 ${\rm{ft}}/{{\mathbf{s}}^2}$ ], ${\omega _{sp}}$ is the short period frequency and ${T_{{\theta _2}}}$ is the short period time constant.

The Level 1 FQ boundaries given by MIL-STD-1797A [36] are depicted in Fig. 13(a) and 13(b). The FQ boundaries show the FAA requirements for optimal dynamic parameters ${\xi _{sp}}$ , ${\omega _{sp}}$ , ${T_{{\theta _2}}}$ and CAP parameter. The pitch rate reference model dynamic parameters according to Table 6 are also represented. As seen in Fig. 13(a), the controller’s dynamic parameters are close to the reference model and lying within the Level 1 region for all 64 flight conditions. This shows that the adaptive controller’s dynamic performance is close to the targeted reference model and respects the FAA Level 1 criteria for the overall flight envelope. Furthermore, the CAP parameters for all flight conditions are also lying within the Level 1 region as seen in Fig. 13(b). The CAP parameters are gathered far from the Level 1 boundary, showing dynamic stability margin.

4.1 Wind gust and turbulence tests

The controller was tested against wind gusts in Fig. 14 and turbulence in Fig. 15 [Reference Gao and Liu53].

Figure 14. Pitch rate and elevator responses under tailwind gust.

Figure 15. Turbulence profiles at different altitudes, pitch rate and elevator responses.

The tailwind gust was simulated for ${u_{wind}} = 40$ m/s and ${w_{wind}} = 5$ m/s at $t = 5$ sec. The lateral wind component ${v_{wind}}$ was set to $0$ m/s. Figure 15 illustrates the controller’s behaviour in the presence of turbulence. Dryden turbulence model was simulated for moderate turbulence in Simulink.

4.2 Actuator fault tests

Finally, the flight controller was validated for actuator faults [Reference Hussein, Zhang and Liu54]. Fifty per cent loss-of-effectiveness (LOE) on the elevators was simulated in Fig. 16(a), while white noises on the actuators are shown in Fig. 16(b). The actuator faults were added to the elevator’s effectiveness $C_l^{{\delta _e}}$ .

Figure 16. Pitch rate, elevator and adaptive element responses for (a) 50% LOE and (b) actuator noise.