3.1 Maximum likelihood estimation

A cost function, $J$ , is defined to optimise the parameters using the maximum likelihood function for unknown parameter vector, $\Theta $ , and the measurement noise covariance matrix, $R$ , as follows^{(Reference Jategaonkar2,Reference Verma and Peyada45,Reference Verma and Peyada46)} :

(5) \begin{equation}J(\Theta ,R) = \frac{1}{2}\sum\limits_{i = 1}^N {{E_i}^T{R^{ - 1}}{E_i}} + \frac{N}{2}\ln \{\det \,(R)\} + \frac{{{n_y}N}}{2}\ln \left( {2\pi } \right)\end{equation}

where ${E_i}$ denotes the residual error, expressed as ${E_i} = \left[ {Z(i) - Y(i)} \right]$ .

In the optimisation procedure, it is assumed that the input data samples are well generated using the aerodynamic parameters and a priori measured motion and control variables for prediction of the observations. Hence, the measurement noise covariance matrix is known before the optimisation, which can be expressed as Equation (6) by applying the partial derivative of the cost function, $J\left( {\Theta ,R} \right)$ w.r.t $R$ :

(6) \begin{equation}R = \frac{1}{N}\sum\limits_{i = 1}^N {{E_i}{E_i}^T} \end{equation}

The cost function $J(\Theta )$ can thus be expressed as

(7) \begin{equation}J(\Theta ) = \frac{1}{2}{n_y}N + \frac{N}{2}\ln \{\det \,(R)\} + \frac{{{n_y}N}}{2}\ln \left( {2\pi } \right)\end{equation}

The equation (7) includes constant terms, so without affecting the minimum value, the cost function, $J(\Theta )$ can be written in the simplest form as follows:

(8) \begin{equation} J(\Theta ) = \det \,(R)\end{equation}

For optimisation of the unknown parameter vector $\Theta $ , the cost function, $J(\Theta )$ is minimised by applying the GN method and the parameter update vector, $\Delta \Theta$ in the k ^th iteration is expressed as

(9) \begin{equation} \Delta \Theta = - {\left[ {{{\left( {\frac{{{\partial ^2}J}}{{\partial {\Theta ^2}}}} \right)}_k}} \right]^{ - 1}}{\left( {\frac{{\partial J}}{{\partial \Theta }}} \right)_k}\end{equation}

where the gradient vector, $\left( {\frac{{\partial J}}{{\partial \Theta }}} \right)$ of the k ^th iteration is given as ${\left( {\frac{{\partial J}}{{\partial \Theta }}} \right)_k} = - \sum\limits_{i = 1}^N {{{\left[\! {\frac{{\partial Y(i)}}{{\partial \Theta }}} \!\right]}^T}{R^{ - 1}}{E_i}}$ ; the information matrix, $\left( {\frac{{{\partial ^2}J}}{{\partial {\Theta ^2}}}} \right)$ of the k ^th iteration is given as ${\left( {\frac{{{\partial ^2}J}}{{\partial {\Theta ^2}}}} \right)_k} = \sum\limits_{i = 1}^N {{{\left[ \!{\frac{{\partial Y(i)}}{{\partial \Theta }}} \!\right]}^T}{R^{ - 1}}\left[\! {\frac{{\partial Y(i)}}{{\partial \Theta }}} \!\right]} \,$ , which is also known as the Hessian matrix. The response gradient matrix, $\left( {{{\partial Y(i)} \mathord{\left/{\vphantom {{\partial Y(i)} {\partial \Theta }}} \right.} {\partial \Theta }}} \right)$ is computed by applying forward difference approximation theory.

The parameter vector of the (k + 1)^th iteration, ${\Theta _{k + 1}}$ is computed by using the parameter vector of the k ^th iteration, ${\Theta _k}$ as follows:

(10) \begin{equation} {\Theta _{k + 1}} = {\Theta _k} + \Delta \Theta \end{equation}

3.2 Recurrent extreme learning machine-based parameter optimisation

The aircraft’s motion is governed by the external forces and moments due to the aerodynamic effects, engine thrust and gravity, which is represented appropriately by the motion variables such as the true velocity of the aircraft, angle-of-attack, side-slip angle, attitude angles, angular rates, linear accelerations etc., and the control variables using the elevator, aileron and rudder deflections and the engine throttle setting. Such a concept has been applied previously in time-series neural modelling and parameter estimation for longitudinal and lateral directional motions^{(Reference Peyada and Ghosh33,Reference Verma and Peyada45,Reference Verma and Peyada46)} . The current research paper investigates the use of a recurrent neural model in the optimisation of the aerodynamic and stall parameters of the quasi-steady stall manoeuvre of an aircraft. Such a stall phenomenon is categorised as a large-amplitude manoeuvre that exhibits hysteresis in the aerodynamic forces and moments with a considerable amount of turbulence near stall. In the present study, an RELM network is used for the nonlinear mapping between the input and output variables of dimension $X \in {R^{{n_x}}}$ and $Y \in {R^{{n_y}}}$ , respectively. The network consists of a FFNN with feedback connections from the output nodes, adopting the structure of a conventional Jordan artificial neural network (ANN) as shown in Fig. 1. Let us assume a single hidden layer with neurons of dimension $H \in {R^{{n_h}}}$ . In general, the output of such a network can be expressed mathematically using the time delays m and n in the input and output variables, respectively, as follows:

(11) \begin{equation} y(i + 1) = f(x(i),...,x(i - m),y(i),...,y(i - n))\end{equation}

Figure 1. Structure of a typical recurrent neural model.

The main objective of the present investigation is to use this network to optimise the aerodynamic parameters, hence the time delays in the input and output variables are restricted so as to avoid overfitting of the network. Most commonly, ANNs are trained using back-propagation algorithms such as steepest descent, Levenberg–Marquardt and scaled conjugate gradients, or using global optimisation algorithms such as the genetic algorithm, particle swarm optimisation, artificial bee colony optimisation etc., which may lead to overfitting and also may consume more time in some cases. Such issues are overcome by using the nonlinear mapping suggested by Huang et al. for a single hidden layer neural network, known as an ELM network^{(Reference Huang, Zhu and Siew38–Reference Huang, Zhou, Ding and Zhang40)}. The learning procedure for an ELM network is employed in the present study. The input variables of the network are considered to be the coefficients of the aerodynamic forces and moments $({C_D},{C_L},{C_m})$ of the i ^th instant, which are derived from the measured and geometric quantities of the aircraft as follows:

(12) \begin{equation} {C_D}(i) = - {C_X}(i)\cos (\alpha (i)) - {C_Z}(i)\sin (\alpha (i))\end{equation}

(13) \begin{equation} {C_L}(i) = {C_X}(i)\sin (\alpha (i)) - {C_Z}(i)\cos (\alpha (i))\,\,\,\end{equation}

(14) \begin{equation} {C_m}(i) = [{I_y}\dot q(i)\, + ({p^2}(i) - {r^2}(i)){I_{xz}} - p(i)r(i)({I_z} - {I_x}) - T(i){z_{ENCG}}]/(\bar{q}(i)S\bar{c})\end{equation}

where ${C_X},{C_z}$ are the body force coefficients, expressed as follows:

(15) \begin{equation} \begin{array}{l}{C_X}(i) = {{\,m{a_X}(i)} \mathord{\left/ {\vphantom {{\,m{a_X}(i)} {\bar{qS}}}} \right. } {\bar{qS}}}\\[3pt] {C_Z}(i)\, = {{\,m{a_Z}(i)} \mathord{\left/{\vphantom {{\,m{a_Z}(i)} {\bar{qS}}}} \right.} {\bar{qS}}}\end{array}\end{equation}

In the supervised learning of the RELM network, the motion variables $(\alpha ,\theta ,q,V,{a_x},{a_z})$ are considered as (i + 1)^th instant target variables, while for their feedback to the input nodes, they are considered at the i ^th instant. The input–output variables have different physical significance, hence their normalisation is essential to enhance the nonlinear mapping for a given flight dataset. Generally, the input and output variables are normalised within a specified range of $({\bar{x}_{\max}},{\bar{x}_{\min }})$ as follows^{(Reference Sola and Sevilla42)}:

(16) \begin{equation} {\bar{x}_i} = {\bar{x}_{i,\min }} + \frac{{({{\bar{x}}_{i,\max }} - {{\bar{x}}_{i,\min }})}}{{({x_{i,\max }} - {x_{i,\min }})}}({x_i} - {x_{i,\min }})\end{equation}

In the supervised learning of the network, the input and output samples of the i ^th instant are processed from the context nodes to the neurons of the hidden layer through the connecting elements of the input weight matrix, $\bar{V}$ with dimension $({n_x} + {n_y} + 1) \times {n_h}$ . Log sigmoid-based transfer functions are used in the neurons of the hidden layer, thus the output of the j ^th neuron is given by the expression

(17) \begin{equation} {\bar{h}_j} = g\left( {\sum\limits_{i = 1}^{{n_x}} {{{\bar{x}}_i}{{\bar{v}}_{ij}}} + \sum\limits_{i = {n_x} + 1}^{{n_y}} {{{\bar{y}}_i}{{\bar{v}}_{ij}}} + {{\bar{v}}_{({n_x} + {n_y} + 1)j}}} \right)\end{equation}

where $i = 1,2,3,...,({n_x} + {n_y} + 1)$ ; $j = 1,2,3,...,{n_y}$ ; ${\bar{v}_{ij}}$ denotes the connecting element of $\bar{V}$ between the i ^th and j ^th neuron of the input and hidden layer, respectively; ${\bar{v}_{({n_x} + {n_y} + 1)j}}$ denotes the bias of the j ^th neuron of the hidden layer; $g( \cdot )$ denotes the transfer function relationship of the log sigmoid used by the hidden layer neuron.

Let us assume that the weight between the hidden and output layer is $\bar{W}$ . Therefore, the k ^th node of the network’s output layer can be approximated using the elements of the output weight, $\bar{W}$ as follows:

(18) \begin{equation} {\bar{y}_k} = \sum\limits_{j = 1}^{{n_h}} {{{\bar{h}}_j}{{\bar{w}}_{jk}}} \end{equation}

where $k = 1,2,3,...,{n_y}$ ; ${\bar{h}_j}$ denotes the output of the j ^th neuron of the hidden layer; ${\bar{w}_{jk}}$ denotes the element of the output weight $\bar{W}$ between the j ^th neuron and the k ^th node of the hidden and output layer, respectively.

Equations (17) and (18) represent a typical way of processing the data through a single hidden layer recurrent neural model. The output of the hidden layer neurons can be expressed in vector form as $\bar{H} = [ {{{\bar{h}}_1}}\ {{{\bar{h}}_2}}\ \cdots \ {{{\bar{h}}_{{n_h}}}} ]$ , thus the predicted output can be given as $\bar{Y} = \bar{H}\bar{W}$ . In the supervised learning of the network, the objective is to minimise the residual error between the predicted and measured output for output weight $\bar{W}$ as follows:

(19) \begin{equation} \mathop {min }\limits_{\bar{W}} ||\,\bar{H}\bar{W} - \bar{Z}||\end{equation}

Equation (19) becomes a set of linear equations $\bar{H}\bar{W} = \bar{Z}$ that provides a number of solutions. Huang et al.^{(Reference Huang, Zhu and Siew38–Reference Huang, Zhou, Ding and Zhang40)} suggested to use the minimum norm of $\bar{W}$ among all the solutions. The output weight can be estimated as follows:

(20) $$\bar H\bar W = \bar Z \to \hat \bar W = {H^\dag }\bar Z$$

where $${\bar H^\dag } = {({\bar H^T}\bar H)^{ - 1}}{\bar H^T}$$ is a Moore–Penrose generalised inverse matrix of $\bar{H}$ . The equation (20) can be determined analytically for given input weights, biases and number of hidden layer neurons. Bartlett^{(Reference Bartlett41)} pointed out that the magnitude of the network weights is more important than the value of the residual error of the training example in terms of performance generalisation. A smaller norm of the network weight would tend to indicate better generalisation. As analysed above, the learning procedure not only reaches a smaller value of the residual error but also a smaller norm of the output weight for given input weights and hidden layer neuron biases.

In this study, we use two statistical measures, namely the mean square error (MSE) and the coefficient of determination (R ²), to describe the training and generalisation abilities of the RELM network^{(Reference Sahin43,Reference Verma and Peyada46)} . The MSE is defined as the mean of the sum of the squared error between the predicted and measured quantities and is expressed as follows:

(21) \begin{equation} MSE = \,\frac{1}{{{n_y}N}}\sum\limits_{i = 1}^N {{{({Y_i} - {Z_i})}^2}} \end{equation}

The coefficient of determination (R ²) is a quantity used to assess the quality of the prediction by the neural model by illustrating how well the nonlinear map operates between the chosen input–output dataset. The value of R ² lies between 0.0 and 1.0. A value closer to 1.0 indicates good predictions for similar input data samples. It is expressed for the chosen network as follows:

(22) \begin{equation} {R^2} = \frac{1}{{{n_y}}}\left[ {1 - \frac{{\sum\limits_{i = 1}^N {{{({Z_i} - {Y_i})}^2}} }}{{\sum\limits_{i = 1}^N {{{({Z_i} - {Z_{avg.}})}^2}} }}} \right]\end{equation}

where ${Z_{avg.}}$ denotes the mean of the measured output data samples, $Z$ .

This optimisation methodology is applied herein to extract the aerodynamic parameters of a postulated aerodynamic model using the RELM network as shown in Fig. 2. The structure of the recurrent neural model is used for predictions of similar input samples. The input variables of the neural model are the coefficients of the aerodynamic forces and moment, which are computed in two ways. The first approach employs the computation based on the measurements of the linear and angular accelerations as expressed in Equations (12)–(14), as utilised in the training of the neural model. The other approach is to postulate them by applying Taylor’s series expansion theory on the governing functions of the coefficients with respect to the motion and control variables measured during flight. The latter approach is commonly used in the postulation of aerodynamic models and can be computed analytically for given values of the parameters. In solving the parameter identification problem, one must optimise the parameters of the aerodynamic model using some optimisation algorithm based on the residual error minimisation principle.

Figure 2. Flow chart of the RELM-GN optimisation methodology.

In the optimisation of the stability and control derivatives, denoted as $\Theta$ , the observations are predicted by propagating the analytically generated input samples through the RELM network as follows:

(23) \begin{equation} Y(i + 1) = f(X(i),Y(i),\bar{V},\bar{W},\Theta )\end{equation}

where $f( \bullet )$ denotes the nonlinear function relationship of the RELM network; $Y(i + 1)\,$ is a vector of the predicted observations of the (i + 1)^th instant; $Y(i)\,$ is a vector of the training observations of the i ^th instant; $X(i)$ is a vector of the input samples of the i ^th instant; $\bar{V},\bar{W}$ are the input and output weight matrix, respectively.

The coefficients of the aerodynamic forces and moment are computed analytically using the initial parameters $\left( {{\Theta _0}} \right)$ , and the generated input samples are propagated through the trained RELM network to obtain their corresponding predictions using Equation (23). Similarly, the response gradient matrix ${\left( {{{\partial Y(i)} \mathord{\left/{\vphantom {{\partial Y(i)} {\partial \Theta }}} \right.} {\partial \Theta }}} \right)_{ij}}$ is computed for each of the perturbed parameters ${(\Theta + \delta {\Theta _j})}$ by applying the forward difference approximation theory as follows:

(24) \begin{align} {\left( {\frac{{\partial Y(i)}}{{\partial \Theta }}} \right)_{ij}} \approx \frac{{f\left( {X\left( {i - 1} \right),Y\left( {i - 1} \right),\bar{V},\bar{W},\left( {\Theta + \delta {\Theta _j}} \right)} \right) - f\left( {X\left( {i - 1} \right),Y\left( {i - 1} \right),\bar{V},\bar{W},\Theta } \right)}}{{\delta {\Theta _j}}}\end{align}

where $i = 1,2,3,...,N;\ j = 1,2,3,...,{n_\Theta };\ \delta {\Theta _j}$ denotes a small perturbation in the j ^th parameter of $\Theta$ .

Furthermore, the cost function is computed by using Equation (8), and the convergence (CNVG) of the algorithm is checked. If the optimisation process does not converge, then the parameters are updated by following Equations (9) and (10) of the GN method. The optimisation process is continued until convergence occurs.

Finally, the statistical value of the optimised parameter in terms of the standard deviation (STD) is computed, given by the square root of the diagonal element of the estimation error covariance matrix, P, as follows:

(25) \begin{equation} {\sigma _{{\Theta _j}}} = \sqrt {{p_{jj}}} \end{equation}

where $j = 1,2,3,...,{n_\Theta };$ P is approximated as the inverse of the Hessian matrix as follows:

(26) \begin{equation} P \approx {\left\{ {\sum\limits_{i = 1}^N {{{\left[ {\frac{{\partial Y(i)}}{{\partial \Theta }}} \right]}^T}\left[ {{R^{ - 1}}} \right]} \left[ {\frac{{\partial Y(i)}}{{\partial \Theta }}} \right]} \right\}^{ - 1}}\end{equation}

3.3 Model validation

Generally, the model is validated by the quality of fit using one of the statistical approaches. Theil’s inequality coefficient (TIC) is one of these, being used to represent the goodness of fit in terms of the residuals^{(Reference Jategaonkar2)} and expressed as

(27) \begin{equation} U = \frac{{\sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{{\left[ {{z_i} - {y_i}} \right]}^2}} } }}{{\sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{z_i}^2} } + \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{y_i}^2} } }}\end{equation}

where $z$ and $y$ are the measured and model predicted responses, respectively. $U$ is a normalised index whose value lies in the range $\left[ {0,1} \right]$ . $U = 0$ corresponds to the case of equality, whereas $U = 1$ corresponds to the case of maximum inequality.