Nomenclature

$\boldsymbol{{A}}$

steady-state state matrix

$b$

Wingspan in m

$\overline{c}$

mean aerodynamic chord in m

$\boldsymbol{{C}}$

steady-state observation matrix

${C_{{D_0}}}$

drag force coefficient at zero lift

${C_{{L_{{\delta _e}}}}}$

derivative of lift force coefficient w.r.t. elevator deflection

${C_{{L_0}}}$

lift force coefficient at zero angle-of-attack

${C_{{L_\alpha }}}$

derivative of lift force coefficient w.r.t. angle-of-attack

${C_{{Y_{{\delta _r}}}}}$

derivative of side force coefficient w.r.t. rudder deflection

${C_{{Y_0}}}$

side force coefficient at zero sideslip angle

${C_{{Y_p}}}$

derivative of side force coefficient w.r.t. roll rate

${C_{{Y_r}}}$

derivative of side force coefficient w.r.t. yaw rate

${C_{{Y_\beta }}}$

derivative of side force coefficient w.r.t. sideslip angle

${C_{{l_{{\delta _a}}}}}$

derivative of rolling moment coefficient w.r.t. aileron deflection

${C_{{l_{{\delta _r}}}}}$

derivative of rolling moment coefficient w.r.t. rudder deflection

${C_{{l_0}}}$

rolling moment coefficient at zero sideslip angle

${C_{{l_p}}}$

derivative of rolling moment coefficient w.r.t. roll rate

${C_{{l_r}}}$

derivative of rolling moment coefficient w.r.t. yaw rate

${C_{{l_\beta }}}$

derivative of rolling moment coefficient w.r.t. sideslip angle

${C_{{m_{{\delta _e}}}}}$

derivative of pitching moment coefficient w.r.t. elevator deflection

${C_{{m_0}}}$

pitching moment coefficient at zero angle-of-attack

${C_{{m_q}}}$

derivative of pitching moment coefficient w.r.t. pitch rate

${C_{{m_\alpha }}}$

derivative of pitching moment coefficient w.r.t. angle-of-attack

${C_{{n_{{\delta _r}}}}}$

derivative of yawing moment coefficient w.r.t. rudder deflection

${C_{{n_0}}}$

Yawing moment coefficient at zero sideslip angle

${C_{{n_p}}}$

derivative of yawing moment coefficient w.r.t. roll rate

${C_{{n_r}}}$

derivative of yawing moment coefficient w.r.t. yaw rate

${C_{{n_\beta }}}$

derivative of yawing moment coefficient w.r.t. sideslip angle

${C_D}$

nondimensional drag force coefficient

${C_L}$

nondimensional lift force coefficient

${C_Y}$

nondimensional side force coefficient

${C_{\boldsymbol{{l}}}}$

nondimensional rolling moment coefficient

${C_{\boldsymbol{{m}}}}$

nondimensional pitching moment coefficient

${C_{\boldsymbol{{n}}}}$

nondimensional yawing moment coefficient

$\boldsymbol{\mathcal{F}}$

Fisher information matrix

$\boldsymbol{{F}}$

process noise distribution matrix

${F_t}$

thrust force in N

g

acceleration due to gravity in ${\textrm{m}}/{{\textrm{s}}^2}$

$\boldsymbol{\mathcal{G}}$

gradient vector

$G$

measurement noise distribution matrix

${I_{xx}}$

mass moment of inertia about body x-axis in ${\textrm{kg }}{{\textrm{m}}^2}$

${I_{xz}}$

product moment of inertia in body xz-plane in ${\textrm{kg }}{{\textrm{m}}^2}$

${I_{yy}}$

mass moment of inertia about body y-axis in ${\textrm{kg }}{{\textrm{m}}^2}$

${I_{zz}}$

mass moment of inertia about body z-axis in ${\textrm{kg }}{{\textrm{m}}^2}$

$\boldsymbol{{J}}$

cost function

$k$

induced drag force correction factor

$\boldsymbol{{K}}$

steady-state Kalman gain matrix

$\boldsymbol{{L}}$

likelihood function

m

mass of UAV in kg

${n_p}$

number of parameters

${n_x}$

number of state variables

${n_y}$

number of observation variables

$p$

roll rate in rad/s

$\boldsymbol{{P}}$

steady-state prediction error covariance matrix

q

pitch rate in rad/s

$r$

Yaw rate in ${\textrm{rad}}/{\textrm{s}}$

$\boldsymbol{{R}}$

measurement noise covariance matrix

S

Wing planform area in ${{\textrm{m}}^2}$

${t_k}$

discrete-time

$\boldsymbol{{u}}$

control input vector

$\overline{\boldsymbol{{u}}}$

discrete-time control input vector

$\boldsymbol{{v}}$

measurement noise vector

V

freestream velocity in m/s

$\boldsymbol{{V}}$

measurement noise matrix

$\boldsymbol{{w}}$

process noise vector

$\boldsymbol{{x}}$

state vector

$\widehat{\boldsymbol{{x}}}$

corrected state vector

$\widetilde{\boldsymbol{{x}}}$

predicted state vector

$\boldsymbol{{y}}$

output vector

$\widetilde{\boldsymbol{{{y}}}}$

predicted output vector

${\textbf{z}}$

measurement vector

Greek symbol

$\alpha $

angle-of-attack in rad

$\beta $

sideslip angle in rad

${\delta _a}$

aileron deflection angle in rad

${\delta _e}$

elevator deflection angle in rad

${\delta _r}$

rudder deflection angle in rad

$\theta $

pitch angle in rad

$\mu $

mean of a scalar Gaussian random variable

ρ

density of air in ${\textrm{kg}}/{{\textrm{m}}^3}$

$\sigma $

standard deviation of a scalar Gaussian random variable

$\phi $

roll angle in rad

M

mean of a Gaussian random vector

${\boldsymbol{\Sigma }}$

covariance of a Gaussian random vector

${\boldsymbol{\Theta }}$

parameter vector

2.0 Mathematical modelling

The flight data for this research was acquired using two indigenously designed propeller-driven wing-alone UAVs, namely CDFP and CDRW [Reference Saderla, Kim and Ghosh3]. CDFP, with a gross take-off mass of 3.5 kg, features a cropped delta wing with a flat-plate aerofoil, while CDRW, with a gross take-off mass of 3.6 kg, features a reflex aerofoil. The UAVs house their separate fuselages, and they have an all-moving vertical tail, which acts as both a vertical stabiliser and rudder. Both UAVs do not have a dedicated horizontal tail, but they feature elevons that act as the control surface to induce pitching and rolling moments independently [Reference Saderla, Dhayalan and Ghosh35]. Both CDFP and CDRW, whose prototypes are given in Fig. 1, were developed, instrumented and tested at the Flight Laboratory of IIT Kanpur. The standard geometrical parameters of CDFP and CDRW have been listed in Table 1.

Table 1. Standard geometric parameters of CDFP and CDRW

Figure 1. Cropped Delta wing UAVs [Reference Saderla, Dhayalan and Ghosh35].

Figure 2. A schematic representation of the cropped Delta UAV [Reference Kumar, Saderla and Kim36].

The UAVs can be treated as rigid bodies and their six DOF dynamics can be described using differential equations, which are nonlinear and coupled in nature. These governing equations of motion are decoupled, retaining their original non-linearity, under appropriate assumptions based on flight manoeuvres to represent longitudinal and lateral-directional motion independently. The simplifying assumptions include the thrust from the propulsion system being perfectly aligned with the body’s x-axis and the angle-of-attack $\alpha $ being small. A schematic representation of the UAV is given in Fig. 2. The description of the symbols used in the forthcoming sections is given in nomenclature. The governing equations for longitudinal motion used in the current research are given by Equation (1) as follows [Reference Schetz, Klein and Morelli25]:

(1) \begin{align}\dot V & = - \frac{{\rho S{V^2}}}{{2m}} C_D + g\sin\! ( {\alpha - \theta } ) + \frac{{{F_t}}}{m}\cos\! ( \alpha ) \nonumber\\ \dot \alpha & = - \frac{{\rho SV}}{{2m}}{C_L} + \frac{g}{V}\cos\! ( {\alpha - \theta } ) - \frac{{{F_t}}}{{mV}}\sin\! ( \alpha ) + q \nonumber\\ \dot q & = \frac{{\rho S\overline{c}{V^2}}}{{2{I_{yy}}}}{C_m} \nonumber\\ \dot \theta & = q \end{align}

The lateral-directional dynamic equations in current research are given by Equation (2) as follows [Reference Schetz, Klein and Morelli25]:

(2) \begin{align}\dot \beta & = - \frac{{\rho SV}}{{2m}}{C_Y} - \frac{{{F_t}}}{{mV}}\sin\! ( \beta ) + \frac{g}{V}\sin\! ( \phi ) - r \nonumber\\ \dot p & = \frac{1}{2}\left( {\frac{{\rho S{V^2}b\left( {{I_{zz}}{C_l} + {I_{xz}}{C_n}} \right)}}{{{I_{xx}}{I_{zz}} - I_{xz}^2}}} \right) \nonumber\\ \dot r & = \frac{1}{2}\left( {\frac{{\rho S{V^2}b\left( {{I_{xz}}{C_l} + {I_{xx}}{C_n}} \right)}}{{{I_{xx}}{I_{zz}} - I_{xz}^2}}} \right) \nonumber\\ \dot \phi & = p \end{align}

The aerodynamic coefficients involved in the dynamic equations, if the UAV is operated at a low angle-of-attack, can be modelled using the linear equations as follows [Reference Saderla, Kim and Ghosh3, Reference Bryan37]:

(3) \begin{align}{C_L} & = {C_{{L_0}}} + {C_{{L_\alpha }}}\alpha + {C_{{L_q}}}\frac{{q\overline{c}}}{{2V}} + {C_{{L_{{\delta _e}}}}}{\delta _e} \nonumber\\ {C_D} & = {C_{{D_0}}} + kC_L^2 \nonumber\\ {C_m} & = {C_{{m_0}}} + {C_{{m_\alpha }}}\alpha + {C_{{m_q}}}\frac{{q\overline{c}}}{{2V}} + {C_{{m_{{\delta _e}}}}}{\delta _e} \nonumber \\ {C_Y} & = {C_{{Y_0}}} + {C_{{Y_\beta }}}\beta + {C_{{Y_p}}}\frac{{pb}}{{2V}} + {C_{{Y_r}}}\frac{{rb}}{{2V}} + {C_{{Y_{{\delta _r}}}}}{\delta _r} \nonumber\\ {C_l} & = {C_{{l_0}}} + {C_{{l_\beta }}}\beta + {C_{{l_p}}}\frac{{pb}}{{2V}} + {C_{{l_r}}}\frac{{rb}}{{2V}} + {C_{{l_{{\delta _a}}}}}{\delta _a} + {C_{{l_{{\delta _r}}}}}{\delta _r} \nonumber\\ {C_n} & = {C_{{n_0}}} + {C_{{n_\beta }}}\beta + {C_{{n_p}}}\frac{{pb}}{{2V}} + {C_{{n_r}}}\frac{{rb}}{{2V}} + {C_{{n_{{\delta _r}}}}}{\delta _r} \end{align}

Equation (1) can be rewritten in the following vector form,

(4) \begin{align}\dot{\boldsymbol{{x}}}_{\boldsymbol{{L}}}(t) = {\boldsymbol{{f}}_{\boldsymbol{{L}}}}( {{\boldsymbol{{x}}_{\boldsymbol{{L}}}}(t),{\boldsymbol{{u}}_{\boldsymbol{{L}}}}(t),{{\boldsymbol{\Theta }}_{\boldsymbol{{L}}}}} ) \end{align}

In Equation (4),

\begin{align*}{\boldsymbol{{x}}_{\boldsymbol{{L}}}} = \left[ {V,{\textrm{ }}\alpha ,{\textrm{ }}q,{\textrm{ }}\theta } \right] \end{align*}

\begin{align*}{\boldsymbol{{u}}_{\boldsymbol{{L}}}} = \left[ {{\delta _e},{\textrm{ }}{F_t}} \right] \end{align*}

\begin{align*}{{\boldsymbol{\Theta }}_{\boldsymbol{{L}}}} = \left[ {{C_{{D_0}}},k,{C_{{L_0}}},{C_{{L_\alpha }}},{C_{{L_q}}},{C_{{L_{{\delta _e}}}}},{C_{{m_0}}},{C_{{m_\alpha }}},{C_{{m_q}}},{C_{{m_{{\delta _e}}}}}} \right] \end{align*}

are the state variables, control variables, and parameters for the longitudinal dynamics, respectively.

Equation (2) can be rewritten in the following vector form,

(5) \begin{align}{\dot{\boldsymbol{{x}}}_{\boldsymbol{{LD}}}}(t) = {\boldsymbol{{f}}_{\boldsymbol{{LD}}}}\left( {{\boldsymbol{{x}}_{\boldsymbol{{LD}}}}(t),{\boldsymbol{{u}}_{\boldsymbol{{LD}}}}(t),{{\boldsymbol{\Theta }}_{\boldsymbol{{LD}}}}} \right) \end{align}

In Equation (5),

\begin{align*}{\boldsymbol{{x}}_{\boldsymbol{{LD}}}}(t) = [ {\beta ,p,r,\phi } ] \end{align*}

\begin{align*}{\boldsymbol{{u}}_{\boldsymbol{{LD}}}}(t) = [ {{\delta _a},{\delta _r},{F_t}} ] \end{align*}

\begin{align*}{{\boldsymbol{\Theta }}_{\boldsymbol{{LD}}}}(t) = \left[ {{C_{{Y_0}}},{C_{{Y_\beta }}},{C_{{Y_p}}},{C_{{Y_r}}},{C_{{Y_{{\delta _r}}}}},{C_{{l_0}}},{C_{{l_\beta }}},{C_{{l_p}}},{C_{{l_r}}},{C_{{l_{{\delta _a}}}}},{C_{{l_{{\delta _r}}}}},{C_{{n_0}}},{C_{{n_\beta }}},{C_{{n_p}}},{C_{{n_r}}},{C_{{n_{{\delta _r}}}}}} \right] \end{align*}

are the state variables, control variables and parameters for the lateral-directional dynamics, respectively.

In general, the nonlinear equations of motion describing the longitudinal and lateral-directional dynamics of a UAV can be represented in state space [Reference Jategaonkar12] as follows in Equation (6),

(6) \begin{align}\dot{\boldsymbol{{x}}}(t) & = \boldsymbol{{f}}( {\boldsymbol{{x}}(t),\boldsymbol{{u}}(t),{\boldsymbol{\Theta }}} ) \nonumber\\ \boldsymbol{{y}}(t) & = \boldsymbol{{g}}( {\boldsymbol{{x}}(t),\boldsymbol{{u}}(t),{\boldsymbol{\Theta }}} ) \end{align}

Equation (6) does not necessarily capture the dynamics perfectly because of the inaccuracies in the model arising due to unmodelled dynamics and uncertainties in the parameters. These uncertainties in the model can be quantified by modelling the process noise. This research assumes that the process noise follows the Gaussian distribution with zero mean and identical covariance. During flight tests, the variables of interest are measured using sensors. However, all sensors are not perfect. The data captured by the sensors contain undesirable frequencies, which is termed measurement noise. In this research, it is assumed that the measurement noise is Gaussian with zero mean and identical covariance.

When measurement and process noise are added to Equation (6), it becomes [Reference Morelli and DeLoach5]

(7) \begin{align}\dot{\boldsymbol{{x}}}(t)& = \boldsymbol{{f}}( {\boldsymbol{{x}}(t),\boldsymbol{{u}}(t),{\boldsymbol{\Theta }}} ) + \boldsymbol{{Fw}}(t) \nonumber\\ \boldsymbol{{y}}(t)& = \boldsymbol{{g}}( {\boldsymbol{{x}}(t),\boldsymbol{{u}}(t),{\boldsymbol{\Theta }}} ) \nonumber\\ \boldsymbol{{z}}( {{t_k}} ) & = \boldsymbol{{y}}( {{t_k}} ) + \boldsymbol{{Gv}}( {{t_k}} ) \end{align}

F and G are process and measurement noise distribution matrices, respectively.

Since the noises are assumed to be random, w and v can be treated as random vectors, which means that the system states x is also a random vector, complicating the parameter estimation process because the state cannot be obtained by direct integration of state equations. It must be estimated using a suitable minimum variance estimator, as described in Section 3.

3.0 Parameter estimation methodology

3.1 Multivariate Gaussian distribution

From [Reference Bertsekas and Tsitsiklis38], if $X$ is a Gaussian random variable that takes on values from a discrete set $s\{ {{x_1},{x_2}, \ldots ,{x_m}} \}$ with mean $\mu $ and variance ${\sigma ^2}$, then its probability distribution function ${p_X}( x )$ is given by,

(8) \begin{align}{p_X}( {X = x} ) = \frac{1}{{\sqrt {2\pi {\sigma ^2}} }}\exp \left( { - \frac{{{{( {x - \mu } )}^2}}}{{2{\sigma ^2}}}} \right) \end{align}

If $\boldsymbol{{X}}$ is a Gaussian random vector with mean $\boldsymbol{{M}}$ and covariance ${\boldsymbol{\Sigma }}$, of dimension $n \times 1$ where each of its elements ${X_i}$ take on values from discrete sets ${s_i}\{ {{x_{i1}},{x_{i2}}, \ldots ,{x_{im}}} \}$ and if all ${X_i}$ are independent of each other, the joint probability distribution function ${p_{\boldsymbol{{X}}}}(\boldsymbol{{x}})$ is given by,

(9) \begin{align}{p_{\boldsymbol{{X}}}}( {\boldsymbol{{X}} = \boldsymbol{{x}}} ) = \frac{1}{{{{( {\sqrt {2\pi } } )}^n}}} \times \frac{1}{{\sqrt {{\textrm{det}}( {\boldsymbol{\Sigma }} )} }} \times \exp \left( { - \frac{1}{2}{{\left( {\boldsymbol{{x}} - \boldsymbol{{M}}} \right)}^T}{{\boldsymbol{\Sigma }}^{ - 1}}( {\boldsymbol{{x}} - \boldsymbol{{M}}} )} \right) \end{align}

If $\boldsymbol{{X}}$ is a Gaussian random matrix of dimension $n \times N$ where each of its $N$ columns ${X_i}$ are independent random vectors with identical variance $\boldsymbol{{M}}$, the joint probability distribution ${p_{\boldsymbol{{X}}}}( \boldsymbol{\mathcal{X}} )$ is given by,

(10) \begin{align}{p_{\boldsymbol{{X}}}}( {\boldsymbol{{X}} = \boldsymbol{\mathcal{X}}} ) = {\left( {{{( {2\pi } )}^n}{\textrm{det}}( {\boldsymbol{\Sigma }} )} \right)^{ - \frac{N}{2}}} \times \exp \left( { - \frac{1}{2}\mathop \sum \nolimits_{i = 1}^N {{\left( {{\boldsymbol{{x}}_i} - {\boldsymbol{\mu} _i}} \right)}^T}{{\boldsymbol{\Sigma }}^{ - 1}} ( {{\boldsymbol{{x}}_i} - {\boldsymbol{\mu} _i}} )} \right) \end{align}

3.2 Conventional filter error method (FEM)

The need for a state estimator is described briefly in this subsection. Since the dynamics of the UAV is nonlinear in both longitudinal and lateral cases, EKF has been preferred as the state estimator. In this research, the steady-state EKF algorithm has been used because it dramatically reduces the computation costs if the system is linear time-invariant (LTI) and the deviations from the nominal trajectory are minor. The steady-state EKF algorithm consists of prediction and correction steps given by [Reference Jategaonkar12].

Prediction step

(11) \begin{align} \widetilde{\boldsymbol{{x}}} ( {{t_{k + 1}}} ) & = \widehat{\boldsymbol{{x}}} ( {{t_k}} ) + \mathop \int\limits_{{t_k}}^{{t_{k + 1}}} \boldsymbol{{f}}( {x(t),\overline{\boldsymbol{{u}}}({{t_k}}),{\boldsymbol{\Theta }}} ) {\textrm{d}}t,\quad \widehat{\boldsymbol{{x}}}( {{t_0}} ) = {\boldsymbol{{x}}_{\boldsymbol{0}}} \nonumber\\ \widetilde{\boldsymbol{{y}}} ({{t_k}} ) & = \boldsymbol{{g}}( {\widetilde{\boldsymbol{{x}}}( {{t_k}} ), \overline{\boldsymbol{{u}}}( {{t_k}} ), {\boldsymbol{\Theta }}} ) \end{align}

Correction step

(12) \begin{align}\widehat{\boldsymbol{{x}}} ( {{t_k}} ) = \widetilde{\boldsymbol{{x}}} ( {{t_k}} ) + \boldsymbol{{K}}[ {\boldsymbol{{z}}( {{t_k}} ) - \widetilde{\boldsymbol{{y}}}( {{t_k}} )} ] \end{align}

The steady-state Kalman gain $\boldsymbol{{K}}$ in Equation (12) is obtained by solving the steady-state Ricatti equation given by Equation (13) for steady-state prediction error covariance matrix $\boldsymbol{{P}}$ as follows,

(13) \begin{align} \boldsymbol{{AP}} + {\boldsymbol{{PA}}^{\boldsymbol{{T}}}} & \! - \frac{1}{{\Delta t}}\boldsymbol{{P}}{\boldsymbol{{C}}^{\boldsymbol{{T}}}}{\boldsymbol{{R}}^{ - \boldsymbol{1}}}\boldsymbol{{CP}} + \boldsymbol{{F}}{\boldsymbol{{F}}^{\boldsymbol{{T}}}} = \boldsymbol{0} \nonumber\\ & \qquad\boldsymbol{{K}} = \boldsymbol{{P}}{\boldsymbol{{C}}^{\boldsymbol{{T}}}}{R^{ - \boldsymbol{1}}} \end{align}

The State and Observation matrices in Equation (13) are obtained using Newton’s Central Difference scheme. Proper tuning of F and G matrices decides the success of EKF. This adaptive filtering problem is solved using the combined formulation technique proposed by Maine and Iliff in Ref. [Reference Maine and Iliff15]. The EKF by itself does not do parameter estimation. It gives the best estimate of the state based on current parameters. The goodness of this generated estimate is measured with a cost function. This cost function must be minimised for the parameters using an optimisation algorithm. In conventional EKF implementation, Gauss-Newton (GN) method is used for carrying out the optimisation.

Using Equations (7) and (10), the likelihood function $L( {\boldsymbol{{z}}{\textrm{|}}{\boldsymbol{\Theta }}} ),$ which captures the probability of observing the data $\boldsymbol{{z}}$ given the parameters ${\boldsymbol{\Theta }}$ for the case of $\boldsymbol{{z}}$ being drawn from Gaussian distribution is given by,

(14) \begin{align}L ( {\boldsymbol{{z}}|{\boldsymbol{\Theta }},\boldsymbol{{R}}} ) = {( {{{( {2\pi } )}^n}{\textrm{det}}(\boldsymbol{{R}})} )^{ - \frac{N}{2}}} \times \exp \left( { - \frac{1}{2}\mathop \sum \nolimits_{i = 1}^N {{( {{\boldsymbol{{z}}_i} - {{\widetilde{\boldsymbol{{y}}}}_i}} )}^T}{\boldsymbol{{R}}^{ - \boldsymbol{1}}}( {{\boldsymbol{{z}}_i} - {{\widetilde{\boldsymbol{{y}}}}_i}} )} \right) \end{align}

It is possible to convert the parameter estimation problem into a nonlinear optimisation problem constrained by dynamic equations. If the measurement z and covariance R are fixed, when the optimal value of parameters ${{\boldsymbol{\Theta }}^*}$ is used in Equation (14), $L$ will have its maximum value, which is the maximum likelihood estimate (MLE) for $\boldsymbol{{z}}$. In this research, the cost $J$ is taken as the negative logarithm of Equation (14), as given below,

(15) \begin{align}J( {\boldsymbol{{z}}{\textrm{|}}{\boldsymbol{\Theta }},\boldsymbol{{R}}} ) = \frac{1}{2}\mathop \sum \nolimits_{i = 1}^N {( {{\boldsymbol{{z}}_i} - {{\widetilde{\boldsymbol{{y}}}}_i}} )^T} {\boldsymbol{{R}}^{ -\boldsymbol{1}}}( {\boldsymbol{{z}}_i} - {{\widetilde{\boldsymbol{{y}}}_i}} ) + \frac{N}{2}\log\! ( {{\textrm{det}}(\boldsymbol{{R}})} ) + \frac{{N{n_y}}}{2}{\textrm{log}}( {2\pi } ) \end{align}

If the value of $\boldsymbol{{R}}$ is unknown, it is taken to be the value that minimises Equation (15) when ${\boldsymbol{\Theta }}$ is a known value as proposed in the two-step relaxation strategy [Reference Jategaonkar12]. It is found by partially differentiating Equation (15) with respect to $\boldsymbol{{R}}$ and equating it to $\boldsymbol{0}$ as follows,

(16) \begin{align}\frac{{\partial J(\boldsymbol{{R}})}}{{\partial \boldsymbol{{R}}}} = \boldsymbol{0} \end{align}

From, [Reference Petersen and Pedersen39]

(17) \begin{align}\boldsymbol{{R}} = \frac{1}{N}\mathop \sum \nolimits_{i = 1}^N ( {{\boldsymbol{{z}}_i} - {{\widetilde{\boldsymbol{{y}}}}_i}} ){( {{\boldsymbol{{z}}_i} - {{\widetilde{\boldsymbol{{y}}}}_i}} )^T} \end{align}

Substituting Equation (17) in Equation (15),

(18) \begin{align}J( {\boldsymbol{{z}}{\textrm{|}}{\boldsymbol{\Theta }}} ) = \frac{{N{n_y}}}{2} + \frac{N}{2}\log\! ( {{\textrm{det}}(\boldsymbol{{R}})} ) + \frac{{N{n_y}}}{2}\log\! ( {2\pi } ) \end{align}

Since $J$ is to be optimised for ${\boldsymbol{\Theta }}$, for a given measurement z, Equation (18) can further be reduced into [Reference Jategaonkar12]

(19) \begin{align} J ( {\boldsymbol{\Theta }} ) = {\textrm{det}}(\boldsymbol{{R}}) \end{align}

GN optimisation algorithm is used to optimise Equation (19) w.r.t. ${\boldsymbol{\Theta }}$. The parameter ${\boldsymbol{\Theta }}$ iteratively converges towards the optimal solution ${{\boldsymbol{\Theta }}^*}$ using the gradient information of the cost function. The cost function $J ( {\boldsymbol{\Theta }} )$ can be linearly approximated [Reference Peyada and Ghosh28], and parameter correction ${\boldsymbol{\Delta }}{\boldsymbol{\Theta }}$ [Reference Jategaonkar12] is given as follows,

(20) \begin{align} {\boldsymbol{\Delta }}{\boldsymbol{\Theta }} = - {\left( {\frac{{{\partial ^2}J}}{{\partial {{\boldsymbol{\Theta }}^2}}}} \right)_i}^{ - 1}{\left( {\frac{{\partial J}}{{\partial {\boldsymbol{\Theta }}}}} \right)_i} \cong - {\boldsymbol{\mathcal{F}}^{ - 1}}\boldsymbol{\mathcal{G}} \end{align}

In Equation (20), the first derivative of $J$ w.r.t. ${\boldsymbol{\Theta }}$ represents the gradient $\boldsymbol{\mathcal{G}}$ and an approximation of the second derivate represents the Fisher information matrix $\boldsymbol{\mathcal{F}}$. They are evaluated as follows [Reference Jategaonkar12]:

(21) \begin{align} \boldsymbol{\mathcal{G}} & = - {\boldsymbol{{S}}}^{\boldsymbol{{T}}}{\boldsymbol{{R}}^{ - \boldsymbol{1}}}\mathop \sum \nolimits_{i = 1}^N [ {{\boldsymbol{{z}}_i} - {{\widetilde{\boldsymbol{{y}}}}_i}} ] \nonumber\\ \boldsymbol{\mathcal{F}} & = {\boldsymbol{{S}}^{\boldsymbol{{T}}}}{\boldsymbol{{R}}^{ - \boldsymbol{1}}}\boldsymbol{{S}} \end{align}

In Equation (21), $\boldsymbol{{S}}$ is a matrix of dimension (${n_y} \times {n_p}$) and it is known as the sensitivity matrix, which quantifies the effect of change in parameters on the estimated output vector $\widetilde{\boldsymbol{{y}}}$. Mathematically, it is given as,

(22) \begin{align} \boldsymbol{{S}} = \mathop \sum \limits_{i = 1}^N \frac{{\partial {{\widetilde{\boldsymbol{{y}}}}_{\boldsymbol{{i}}}}}}{{\partial {\boldsymbol{\Theta }}}} \end{align}

The response gradience $\dfrac{{\partial \widetilde{\boldsymbol{{y}}}}}{{\partial {\boldsymbol{\Theta }}}}$ in Equation (22) is obtained numerically by giving a small perturbation $\delta {\boldsymbol{\Theta }}$ to each element of the parameter vector and computing the corresponding perturbed response ${\widetilde{\boldsymbol{{y}}}}_{\boldsymbol{{p}}}$. The ${l^{th}}$ column of the response gradience is computed using the following linear approximation:

(23) \begin{align} {\left( {\frac{{\partial \widetilde{\boldsymbol{{y}}}}}{{\partial {\boldsymbol{\Theta }}}}} \right)_l} = \frac{{{{\widetilde{\boldsymbol{{y}}}}_{{\boldsymbol{{p}}_{\boldsymbol{{l}}}}}} - {{\widetilde{\boldsymbol{{y}}}}_l}}}{{\delta {\boldsymbol{\Theta }}}} \end{align}

Finally, the parameter update is given by,

(24) \begin{align} {{\boldsymbol{\Theta }}_{{\textbf{i}} + \boldsymbol{1}}} = {{\boldsymbol{\Theta }}_{\textbf{i}}} + {\boldsymbol{\Delta} \boldsymbol{\Theta }} \end{align}

3.3. FEM augmented with particle swarm optimisation (FEM-PSO)

The cost function in Equation (19) is optimised using the PSO algorithm. The PSO belongs to a class of computational techniques inspired by biological systems. This evolutionary computational technique performs stochastic optimisation based on the social behaviour of particles. The main advantage of using the PSO algorithm over the conventional Gauss-Newton (GN) algorithm is that the PSO does not require gradience computation at each iteration which significantly reduces the computation cost of the onboard computer. The working of the PSO algorithm for parameter estimation is briefly explained as follows:

Initialise  The search space $S$ where the optimal solution ${{\boldsymbol{\Theta }}^*}$ is likely to be observed is chosen as $\left[ {{\boldsymbol{\Theta }}_{\boldsymbol{{min}}}^{\boldsymbol{{i}}}, {\boldsymbol{\Theta }}_{\boldsymbol{{max}}}^{\boldsymbol{{i}}}} \right]$ where ${\boldsymbol{\Theta }}_{\boldsymbol{{min}}}^{\boldsymbol{{i}}}$ and ${\boldsymbol{\Theta }}_{\boldsymbol{{max}}}^{\boldsymbol{{i}}}$ represent the maximum and minimum limits of the solution. The initial position of all the particles ${{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}( 0 )$ are chosen randomly within $S.$ The initial velocity of all the particles ${\boldsymbol{\Delta }}{{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(0)$ are chosen to be $0$. Initial value of personal best position for each particle ${\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}^{\boldsymbol{{i}}}(0)$ is chosen to be their respective initial positions since the particles have not started moving yet. The initial global best position ${{\boldsymbol{\Theta }}_{\boldsymbol{{gb}}}}(0)$ is chosen to be the personal best position that minimises the cost function given in Equation (19).

(25) \begin{align} {{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(0) &\, :\!=\, {\textrm{rand}}\left[ {{\boldsymbol{\Theta }}_{\boldsymbol{{min}}}^{\boldsymbol{{i}}},{\textrm{ }}{\boldsymbol{\Theta }}_{\boldsymbol{{max}}}^{\boldsymbol{{i}}}} \right],\forall {\textrm{ }}i \in [ {1,2, \ldots ,{N_p}} ] \nonumber\\ {\boldsymbol{\Delta }}{{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(0)&\, :\!=\, \boldsymbol{0} \nonumber\\ {\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}^i(0)&\, :\!=\, {{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(0) \nonumber\\ {{\boldsymbol{\Theta }}_{\boldsymbol{{gb}}}}(0)&\, :\!=\, \mathop {{\textrm{argmin}}}\limits_{{{\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}}} J\left( {{\boldsymbol{\Theta }}_{pb}^i} \right) \end{align}

Update  The position of the particles is updated iteratively, followed by its velocity update to give a new generation of particles. The velocity is updated based on the relative position of the particle from its personal best position and global best position in the current iteration. The following equations list the velocity and position updates of the particles, respectively, for a given generation $t$,

(26) \begin{align} {\boldsymbol{\Delta }}{{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}( {t + 1} ) = w( {t + 1} ){\boldsymbol{\Delta }}{{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(t) + {c_1}\boldsymbol{{a}}_{\boldsymbol{1}}^{\boldsymbol{{T}}}\left[ {{\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}^{\boldsymbol{{i}}}(t) - {{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(t)} \right] + {c_2}\boldsymbol{{a}}_{\boldsymbol{2}}^{\boldsymbol{{T}}}\left[ {{{\boldsymbol{\Theta }}_{\boldsymbol{{gb}}}}(t) - {{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(t)} \right] \end{align}

(27) \begin{align} {{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}( {t + 1} ) = {{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}(t) + {\boldsymbol{\Delta }}{{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}( {t + 1} ) \end{align}

In Equation (26), $w, c_1$, and ${c_2}$ are the hyperparameters of the PSO algorithm and are called the inertial coefficient, personal cognitive coefficient, and social cognitive coefficient, respectively. In this research, the values of ${c_1}$ and ${c_2}$ are chosen to be 2 as given in [Reference Kennedy and Eberhart32]. The value of $w$ is chosen to decrease with each iteration by the rule $w( {t + 1} ) = 0.99 w(t)$ to encourage the particles to follow the optimal direction as given in [Reference Kumar, Saderla and Kim36]. ${\boldsymbol{{a}}_{\boldsymbol{1}}}$ and ${\boldsymbol{{a}}_{\boldsymbol{2}}}$ in Equation (26) are random vectors of appropriate dimensions.

It is to be asserted that ${{\boldsymbol{\Theta }}^{\textbf{i}}}( {t + 1} ) \in S, \forall \in [ {1, 2, \ldots ,{N_p}} ]$. If the ${i^{th}}$ particle is not in $S$, it is to be re-initialised. The personal best position of the particles and the global best position are updated iteratively as follows:

(28) \begin{align} {\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}^{\boldsymbol{{i}}}( {t + 1} )& = \begin{cases}{{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}( {t + 1} ), & \ \ \text{if}\ J\left( {{{\boldsymbol{\Theta }}^{\boldsymbol{{i}}}}( {t + 1} )} \right) \le J\left( {{\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}^{\boldsymbol{{i}}}(t)} \right)\\ {\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}^{\boldsymbol{{i}}}(t), & \ \ \text{Otherwise}\end{cases} \nonumber\\ {{\boldsymbol{\Theta }}_{\boldsymbol{{gb}}}}( {t + 1} ) &= \begin{cases}{\boldsymbol{\Theta }}_{\boldsymbol{{pb}}}^{\boldsymbol{{i}}}( {t + 1} ), & \text{if}\ \ J\left( {{\boldsymbol{\Theta }}_{{{\boldsymbol{{pb}}}}}^{\boldsymbol{{i}}}( {t + 1} )} \right) \le J\left( {{{\boldsymbol{\Theta }}_{\boldsymbol{{gb}}}}(t)} \right)\\ {{\boldsymbol{\Theta }}_{\boldsymbol{{gb}}}}(t),& \text{Otherwise}\end{cases} \end{align}

Convergence  If $t \geqslant {t_{max}}$, the iteration is stopped and ${{\boldsymbol{\Theta }}_{\boldsymbol{{gb}}}}$ is taken as ${{\boldsymbol{\Theta }}^*}$.

Figure 3. PSO algorithm for parameter estimation.

Figure 4. Comparison between actual flight data and simulated data for CDFP in the longitudinal flight regime.

Figure 5. Comparison between flight data and simulated data for CDRW in the longitudinal flight regime.

The above description has been summarised in Fig. 3 as a flowchart.

4.0 Flight data gathering

Flight data gathering (FDG) can be defined as acquiring the time history of motion variables and the corresponding applied control inputs. FDG is such an essential step in the parameter estimation process that the accuracy and reliability of the estimated parameters directly depend on the quality of the gathered data. The amount of information content is a direct measure of the quality of flight data, and it can be held high only when the control input sufficiently excites the dynamics of the UAV independently. For this purpose, the UAV should house an industrial-grade data acquisition system with an inertial measurement unit (IMU) to capture motion variables, standard quality pressure sensors to indirectly measure velocity, and flow angle sensors to measure the angle-of-attack and sideslip angles. The technical details of the data acquisition system and the nature of the flight tests performed with CDFP and CDRW are presented briefly in Ref. [Reference Saderla, Dhayalan and Ghosh35]. The generated flight data are labelled with the UAV name followed by the nature of excitation – L for longitudinal and LD for lateral directional, followed by the serial number of the flight test. Therefore, CDFP_Li would mean that it is the ${i^{\textrm{th}}}$ CDFP dataset where its longitudinal mode is excited.

A total of 16 compatible flight data sets, 4 for longitudinal parameter estimation and 4 for lateral-directional parameter estimation for CDFP and CDRW each, are used in the current research. From Figs 4 and 5, it can be noticed that sinusoidal elevator inputs are applied to the UAVs for the excitation of longitudinal dynamics. Since the maximum angle-of-attack observed in the manoeuvres performed by the UAVs is close to ${10^ \circ }$ and from the wind tunnel experiments carried out using the above UAVs [Reference Saderla40], it was observed that the variation of the aerodynamic coefficients was linear when the angle-of-attack varied between $ - 5^\circ $ and $10^\circ $. Hence the use of the linear aerodynamic model presented in Equation (3) is suitable for the aerodynamic characterisation of the above UAVs. From Figs 7 and 8, it can be referred that the aileron is subjected to sinusoidal inputs varied between $ \pm 10^\circ $ to excite the lateral-directional dynamics of the UAVs. It can be observed that the rudder inputs are negligible in most cases because the UAVs have a high aspect ratio all-movable vertical tail rudder, which indicates high control power, and the slightest deflection of the rudder can make the task of controlling the UAV difficult for a ground stationed pilot. Therefore, lateral-directional dynamics are excited through aileron deflection alone [Reference Saderla, Kim and Ghosh3].

Table 2. Estimated longitudinal parameters for CDFP_L1, CDFP_L2

Note: Values in the () denote Cramer-Rao lower bounds.

Table 3. Estimated longitudinal parameters for CDFP_L3, CDFP_L4

Table 4. Estimated longitudinal parameters for CDRW_L1, CDRW_L2

Table 5. Estimated longitudinal parameters for CDRW_L3, CDRW_L4

Figure 6. Scatter plot comparing the estimated longitudinal parameters with wind tunnel estimates.

Figure 7. Comparison between flight data and simulated data for CDFP in the lateral-directional flight regime.

Figure 8. Comparison between flight data and simulated data for CDRW in the lateral-directional flight regime.

Table 6. Estimated lateral-directional parameters for CDFP_LD1, CDFP_LD2

Table 7. Estimated lateral-directional parameters for CDFP_LD3, CDFP_LD4

Table 8. Estimated lateral-directional parameters for CDRW_LD1, CDRW_LD2

Table 9. Estimated lateral-directional parameters for CDRW_LD3, CDRW_LD4