Aerodynamic characterisation of delta wing unmanned aerial vehicle using non-gradient-based estimator

Abstract

Aerodynamic characterisation from flight testing is an integral subroutine for evaluating a new flight vehicle’s aerodynamic performance, stability and controllability. The estimation of aerodynamic parameters from flight test data has extensively been explored, in the past, using estimation methods such as the equation error method, output error method and filter error method. However, in the current era, non-gradient-based estimation techniques are gaining attention from researchers due to their inherent data-driven optimisation capability to find the global best solution. In this paper, a novel non-gradient-based estimation method is proposed for the aerodynamic characterisation of unmanned aerial vehicles from flight data, which relies on the maximum likelihood method augmented with particle swarm optimisation. Flight data sets of a wing-alone unmanned aerial vehicle are used to demonstrate the capabilities of the proposed method in estimating aerodynamic derivatives. Estimates from the proposed method are corroborated with the wind tunnel test and output error method results. It has been observed that simulated flight vehicle responses using estimated parameters are in good agreement with measured data in most of the manoeuvers considered. Confidence in the estimates of linear and nonlinear aerodynamic parameters is well established with the lower limit of Cramer-Rao bounds, which are minimal. The proposed method also demonstrates good predictability of the quasi-steady stall aerodynamic model by estimating stall characteristic parameters such as aerofoil static stall characteristics parameter, hysteresis time constant and breakpoint. The overall performance of the proposed estimation method is on par with the output error method and is validated with the proof-of-match exercise.

Nomenclature

$a$

aerofoil static stall characteristics parameter

$b$

wing span in m

$\bar c$

mean aerodynamic chord in m

${C_D}$

nondimensional drag force coefficient

${C_{{D_0}}}$

drag force coefficient at zero lift

${C_{{D_X}}}$

derivative of drag force coefficient w.r.t X

${C_l}$

nondimensional rolling moment coefficient

${C_{{l_0}}}$

rolling moment coefficient at zero deg sideslip angle

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

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

${C_{{l_p}}}$

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

${C_{{l_r}}}$

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

${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}$

nondimensional lift force coefficient

${C_{{L_0}}}$

lift force coefficient at zero deg angle-of-attack

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

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

${C_{{L_{{\alpha ^2}}}}}$

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

${C_{{L_q}}}$

damping coefficient of lift force w.r.t pitch rate

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

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

${C_m}$

nondimensional pitching moment coefficient

${C_{{m_0}}}$

pitching moment coefficient at zero deg angle-of-attack

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

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

${C_{{m_q}}}$

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

${C_{{m_X}}}$

derivative of pitching moment coefficient w.r.t X

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

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

${C_n}$

nondimensional yawing moment coefficient

${C_{{n_0}}}$

yawing moment coefficient at zero deg sideslip angle

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

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

${C_{{n_p}}}$

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

${C_{{n_r}}}$

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

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

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

${C_Y}$

nondimensional side force coefficient

${C_{{Y_0}}}$

side force coefficient at zero deg sideslip angle

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

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

${C_{{Y_p}}}$

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

${C_{{Y_r}}}$

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

${C_{{Y_{{\delta _a}}}}}$

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

${F_t}$

thrust produce by engine in N

$g$

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

${I_X}$

moment of inertia along body x-axis in ${\rm{kg}}\hbox{-}{{\rm{m}}^2}$

${I_Y}$

moment of inertia along body y-axis in ${\rm{kg}}\hbox{-}{{\rm{m}}^2}$

${I_Z}$

moment of inertia along body z-axis in ${\rm{kg}}\hbox{-}{{\rm{m}}^2}$

${I_{XZ}}$

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

$J$

cost function

$k$

induced drag force correction factor

$m$

mass of UAV in kg

$p$

roll rate in rad/s

$q$

pitch rate in rad/s

$r$

yaw rate in rad/s

$S$

wing planform area in ${{\rm{m}}^2}$

$V$

air speed in m/s

$X$

nondimensional distance of flow separation point

Greek symbol

$\alpha $

angle-of-attack in rad

${\alpha ^{\star}}$

breakpoint

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

$\phi $

roll angle in rad

$\psi $

yaw angle in rad

$\rho $

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

$\tau $

hysteresis time constant

$\theta $

pitch angle in rad

$\Theta $

vector of unknown parameters

1.0 Introduction

Aerodynamic characterisation of flight vehicles is one of the keen interests of many researchers and, indeed, a subroutine for aircraft design development, simulation and control. Analytical, semi-empirical, computational, and experimental methods have evolved and adapted to address the aforementioned objectives. Though wind tunnel testing suffers wall/sting interference, scale factor, and Reynolds number duplication, it is still a high-fidelity tool to estimate static aerodynamic parameters [1, 2]. With the advent of micro-electro-mechanical sensors and actuators, instrumentation of small-scale unmanned aerial vehicles (UAVs) is made feasible, which aids in acquiring flight data during manoeuvers for aerodynamic characterisation. Parameter estimates obtained from flight tests enable researchers to find out static and dynamic derivatives with great confidence that overcome the limitations of wind tunnel testing and computational fluid dynamics techniques. Equation error eethod [3–5], output error method (OEM) [6–8], filter error method [9–11], neural networks and fuzzy logic-based estimation methods [12–15] are extensively used by many researchers to estimate the aerodynamic parameters of manned and unmanned aircraft in various flight regimes. The equation error method based on the least squares cost function demands relatively less computational power, which can handle the estimation of linear aerodynamic derivatives from flight data in the absence of both measurement and process noise [16]. Due to the scale and size of small UAVs, sensors onboard are prone to the interference of various subsystems, like propulsion units, which eventually leads to measurement noise. Output error method based on maximum likelihood estimator is proved suitable, even in the presence of measurement noise, for estimating linear and nonlinear aerodynamic parameters of UAVs from flight tests data pertaining to various flight regimes [17, 18]. However, the OEM only estimates system parameters deterministically if system dynamics are appropriately modeled. Filter error method is a class of OEM formulated using the extended Kalman filter, which can estimate nondimensional aerodynamic derivatives of UAVs effectively even in the presence of both process and measurement noise [19]. Though filter error method is an efficient method for the aerodynamic characterisation of UAVs, the necessity of accurate priori information of states and heavy computational resources for gradient calculations and a solution to the Riccati equation might have limited its application for UAVs. It is well observed that the OEM and filter error method are sensitive to a priori information about initial conditions for better convergence.

On the contrary, the estimation methods with data-based intelligence make convergence faster, even if initial guess values have a large offset from actual values. The neural network estimation method is used to estimate linear, nonlinear and near-stall flight regime’s aerodynamic parameters of mini UAVs [16, 18]. The idea behind this method lies in the fact that the trained model from measured data is used to replace the system’s dynamics, which can not be a generalised flight dynamic model. Moreover, the confidence in estimates heavily depends upon the training criteria of neural networks, which may be different for different data sets. The fuzzy logic-based estimation method is also used by many researchers to characterise flight vehicle aerodynamics. Unlike the OEM and filter error method, the fuzzy logic-based estimation method does not require a mathematical model of system dynamics [20, 21]. It is generally observed that the model-based estimation techniques suffer sensitivity issues with data and artificial intelligence estimation methods have limitations with the data training, and all these techniques require heavy computational capabilities for accurate estimation. Without compromising the confidence and consistency in the estimates, the aforementioned challenges can be addressed with a hybrid method that involves data search optimisation with classical cost functions.

Particle swarm optimisation (PSO), a non-gradient-based optimisation method developed by Dr. Eberhart and Kennedy in 1995, is a technique inspired by the social behaviour of birds or fish around a food source [22]. Right from its invention, this optimisation tool has been widely used for numerous applications in engineering and sciences. In aerospace engineering, PSO has been employed to carry out aircraft design optimisation [23], flight controller gains optimisation [24] and flight routes optimisation [25]. Recently, PSO has started getting attention from the research community as a tool for the aerodynamic characterisation of flight vehicles. The maximum likelihood method augmented with PSO (ML-PSO) has been used as a technique in estimating linear longitudinal aerodynamic parameters of vertical takeoff and landing aircraft, symmetrical projectile, and UAV from flight test data [26–29]. Though the results from ML-PSO are auspicious, its applicability in estimating aerodynamic parameters in the nonlinear and near-stall flight regimes still needs to be explored.

In this paper, a comprehensive study is conducted to demonstrate the applicability and effectiveness of the ML-PSO in estimating the aerodynamic parameters in longitudinal linear, longitudinal nonlinear, longitudinal near stall and lateral-directional flight regimes of a UAV with limited control surfaces. The UAV considered for the current research is designed with cropped delta planform and reflex aerofoil wing cross-section; it is named Cropped Delta Reflex Wing (CDRW) configuration. Predefined control inputs are used to excite dynamic modes of the UAV to generate various flight data sets pertaining to different flight envelopes. The estimated aerodynamic parameters using the ML-PSO method are corroborated with the OEM and wind tunnel results. Furthermore, a proof-of-match exercise is also performed to validate the considered aerodynamic model and the estimates. The rest of the paper is organised as follows. Section 2 presents the mathematical and aerodynamic modeling of the UAV. The mathematical formulation of the ML-PSO is given in Section 3. Details of the UAV and flight data sets used for the study are presented in Section 4. The outcomes of the implementation of the ML-PSO are given in Section 5. Finally, the advantages and limitations of the ML-PSO in estimating aerodynamic derivatives are presented in Section 6.

2.0 Mathematical modeling of UAV dynamics

Newtonian mechanics was used to formulate the rigid body dynamics of a UAV. In general, equations representing UAV dynamics are coupled in nature, which can be decoupled in longitudinal and lateral-directional cases with suitable assumptions based on flight manoeuvers for aerodynamic characterisation. Equations (1)–(4) and Equations (6)–(9) are derived for longitudinal and lateral-directional motion, respectively, by considering that manoeuvers performed are independent of each other.

(1) \begin{align} \dot V = - \frac{{\rho S{V^2}}}{{2m}}{C_D} + g\sin\!(\alpha - \theta ) + \frac{{{F_t}}}{m}\cos \alpha \end{align}

(2) \begin{align} \dot \alpha = - \frac{{\rho SV}}{{2m}}{C_L} + \frac{g}{V}\cos\!(\alpha - \theta ) - \frac{{{F_t}}}{{mV}}\sin \alpha + q \end{align}

(3) \begin{align} \dot q = \frac{{\rho S\bar c{V^2}}}{{2{I_Y}}}{C_m} \end{align}

(4) \begin{align} \dot \theta = q \end{align}

(5) \begin{align}\dot{\textbf{X}}_{L} = f({{\bf{X}}_L},{\Theta _L},u) \end{align}

where ${{\bf{X}}_L} = {\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}V & \alpha & q & \theta \end{array}} \right]^T}$ , ${\Theta _L}$ is the vector of longitudinal aerodynamic parameters based on the longitudinal flight regime and $u$ is the vector of control inputs.

(6) \begin{align} \dot \beta = - \frac{{\rho SV}}{{2m}}{C_Y} - \frac{{{F_t}}}{{mV}}\sin \beta + \frac{g}{V}\sin \phi - r \end{align}

(7) \begin{align} \dot p = \frac{1}{2}\rho S{V^2}b[{I_Z}{C_l} + {I_{XZ}}{C_n}]\frac{1}{{{I_X}{I_Z} - I_{XZ}^2}} \end{align}

(8) \begin{align} \dot r = \frac{1}{2}\rho S{V^2}b[{I_{XZ}}{C_l} + {I_X}{C_n}]\frac{1}{{{I_X}{I_Z} - I_{XZ}^2}} \end{align}

(9) \begin{align} \dot \phi = p \end{align}

(10) \begin{align} {\dot{\textbf{X}}_{LD}} = f({{\bf{X}}_{LD}},{\Theta _{LD}},u) \end{align}

where ${{\bf{X}}_{LD}} = {\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\beta & p & r & \phi \end{array}} \right]^T}$ , ${\Theta _{LD}}$ is the vector of lateral-directional aerodynamic parameters based on the lateral-directional flight regime and $u$ is the vector of control inputs.

In general, the aerodynamics of a UAV is a nonlinear function of the angle-of-attack ( $\alpha $ ) and angle of sideslip ( $\beta $ ). However, the aerodynamic model can be classified into linear, nonlinear and stall regime models based on flight conditions. Rigorous wind tunnel tests of CDRW were performed at National Wind-tunnel Facility, Indian Institute of Technology Kanpur, to identify the aerodynamic model as a function of flow angles. During wind tunnel tests, the $\alpha $ and $\beta $ were varied from −5 to 50deg and −15 to 15deg, respectively, using the $\beta $ -mechanism. A six-component load balance was used to precisely measure aerodynamic forces and moments acting on a UAV. The obtained data is then processed to calculate the variation of nondimensional aerodynamic coefficients with flow angles, and the corresponding findings are presented in Fig. 1. From Fig. 1(a), it can be observed that the variation of nondimensional lift coefficient ( ${C_L}$ ) with $\alpha $ is linear and nonlinear from −5 to 11deg and 11 to 20deg, respectively. Maximum ${C_L}$ is observed at $\alpha = 21 {\rm{deg}}$ , which remains almost constant with $\alpha $ more than 21deg. From Fig. 1(b), it can be referred that the nondimensional side force coefficient ( ${C_Y}$ ), roll moment coefficient ( ${C_l}$ ) and yaw moment coefficient ( ${C_n}$ ) vary linearly from −5 to 5deg with $\beta $ . Based on the wind tunnel results, a detailed aerodynamics model is represented by Equations (11)–(23), and aerodynamic parameters that need to be estimated are given in Equations (24)–(27).

Figure 1. Wind tunnel results of CDRW UAV.

The aerodynamic model for the linear longitudinal (low angle-of-attack) flight regime is given by Equations (11)–(13) and the corresponding vector of aerodynamic parameters ${\Theta _{LG}}$ in Equation (24).

(11) \begin{align} {C_L} = {C_{{L_0}}} + {C_{{L_\alpha }}}\alpha + {C_{{L_q}}}\frac{{q\bar c}}{{2V}} + {C_{{L_{{\delta _e}}}}}{\delta _e} \end{align}

(12) \begin{align} {C_D} = {C_{{D_0}}} + kC_L^2 \end{align}

(13) \begin{align} {C_m} = {C_{{m_0}}} + {C_{{m_\alpha }}}\alpha + {C_{{m_q}}}\frac{{q\bar c}}{{2V}} + {C_{{m_{{\delta _e}}}}}{\delta _e} \end{align}

The aerodynamic model for the nonlinear longitudinal (moderately high angle-of-attack) flight regime is given by Equations (14)–(16) and the corresponding vector of aerodynamic parameters ${\Theta _{NL}}$ in Equation (25).

(14) \begin{align} {C_L} = {C_{{L_0}}} + {C_{{L_\alpha }}}\alpha + {C_{{L_{{\alpha ^2}}}}}{\alpha ^2} + {C_{{L_q}}}\frac{{q\bar c}}{{2V}} + {C_{{L_{{\delta _e}}}}}{\delta _e} \end{align}

(15) \begin{align} {C_D} = {C_{{D_0}}} + kC_L^2 \end{align}

(16) \begin{align} {C_m} = {C_{{m_0}}} + {C_{{m_\alpha }}}\alpha + {C_{{m_q}}}\frac{{q\bar c}}{{2V}} + {C_{{m_{{\delta _e}}}}}{\delta _e} \end{align}

The aerodynamic model for the quasi-steady stall [30] (near-stall angle-of-attack) is given by Equations (17)–(20) and the corresponding vector of aerodynamic parameters ${\Theta _{ST}}$ in Equation (26).

(17) \begin{align} X = \frac{1}{2}\left[ {1 - \tanh \{ a(\alpha - \tau \dot \alpha - {\alpha ^ \star })\} } \right] \end{align}

(18) \begin{align} {C_L} = {C_{{L_0}}} + {C_{{L_\alpha }}}\alpha {\left[ {\frac{{1 + \sqrt X }}{2}} \right]^2} + {C_{{L_q}}}\frac{{q\bar c}}{{2V}} + {C_{{L_{{\delta _e}}}}}{\delta _e} \end{align}

(19) \begin{align} {C_D} = {C_{{D_0}}} + kC_L^2 + {C_{{D_X}}}(1 - X) \end{align}

(20) \begin{align} {C_m} = {C_{{m_0}}} + {C_{{m_\alpha }}}\alpha + {C_{{m_q}}}\frac{{q\bar c}}{{2V}} + {C_{{m_{{\delta _e}}}}}{\delta _e} + {C_{{m_X}}}(1 - X) \end{align}

Equations (21)–(23) represent the linear (low angle-of-sideslip) lateral-direction aerodynamic model and the corresponding vector of aerodynamic parameters ${\Theta _{LD}}$ are given in Equation (27).

(21) \begin{align} {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} \end{align}

(22) \begin{align} {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} \end{align}

(23) \begin{align} {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}

(24) \begin{align} {\Theta _{LG}} = [{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}}}}}{]^T} \end{align}

(25) \begin{align} {\Theta _{NL}} = [{C_{{D_0}}},k,{C_{{L_0}}},{C_{{L_\alpha }}},{C_{{L_{{\alpha ^2}}}}},{C_{{L_q}}},{C_{{L_{{\delta _e}}}}},{C_{{m_0}}},{C_{{m_\alpha }}},{C_{{m_q}}},{C_{{m_{{\delta _e}}}}}{]^T} \end{align}

(26) \begin{align} {\Theta _{ST}} = [{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}}}}},a,\tau ,{\alpha ^ \star },{C_{{D_X}}},{C_{{m_X}}}{]^T} \end{align}

(27) \begin{align} {\Theta _{LD}} = [{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}}}}}{]^T} \end{align}

where ${\Theta _{LG}}$ , ${\Theta _{NL}}$ , ${\Theta _{ST}}$ and ${\Theta _{LD}}$ are the vectors of unknown parameters related to longitudinal linear, longitudinal nonlinear, longitudinal stall and lateral-directional flight regimes, respectively.

3.0 Formulation of ML-PSO

In general, nonlinear dynamics of an aircraft system without modeling uncertainties can be represented by Equations (28)–(30).

(28) \begin{align} \dot x(t) = f[x(t),u(t),\Theta ]{\rm{,}}\quad \quad x(0) = {x_0}, \end{align}

(29) \begin{align} y(t) = g[x(t),u(t),\Theta ], \end{align}

(30) \begin{align} z({t_k}) = y({t_k}) + Gv({t_k}) \end{align}

where

$x(t)$ the vector of state variables,

$u(t)$ the vector of control inputs independent of the system dynamics,

$\Theta $ the vector of unknown parameters,

$y({t_k})$ output vector at ${k^{{\rm{th}}}}$ discrete time,

$z({t_k})$ the vector of sensor outputs at ${k^{{\rm{th}}}}$ discrete time,

$f,g$ the nonlinear real-valued functions,

$G$ the measurement noise distribution matrix,

$\nu ({t_k})$ the measurement noise with zero mean and nonzero variance

Since sensor outputs are corrupted with measurement noise, $z$ can be considered as a vector-valued random variable of dimension ${n_y}$ . If $N$ measurements of $z$ are available, a likelihood function [31] using multivariate-Gaussian distribution can be defined as in Equation (31):

(31) \begin{align} p(z|\Theta ,{R_1}) = \{ (2\pi {)^{{n_z}}}|{R_1}{|\} ^{ - N/2}}\exp\!\left[ { - \frac{1}{2}\Sigma _{k = 1}^N{{[z({t_k}) - y({t_k})]}^T}R_1^{ - 1}[z({t_k}) - y({t_k})]} \right] \end{align}

where $p(z|\Theta ,{R_1})$ is the probability of $z$ with given $\Theta $ and ${R_1}$ . ${R_1}$ is the measurement noise covariance matrix. The maximisation likelihood function can be changed to the following minimisation problem by taking the negative logarithmic of Equation (31).

(32) \begin{align} L(z|\Theta ,{R_1}) = \frac{1}{2}\Sigma _{k = 1}^N{[z({t_k}) - y({t_k})]^T}R_1^{ - 1}[z({t_k}) - y({t_k})] + \frac{N}{2}{\rm{ln}}[\det ({R_1})] + \frac{{N{n_z}}}{2}{\rm{ln}}(2\pi ) \end{align}

If ${R_1}$ is unknown, it can be obtained by minimising $L(z|\Theta ,{R_1})$ w.r.t ${R_1}$ .

(33) \begin{align} {R_1} = \frac{1}{N}\Sigma _{k = 1}^N{[z({t_k}) - y({t_k})]^T}[z({t_k}) - y({t_k})] \end{align}

here ${R_1}$ changed to nothing but residual co-variance matrix and value of $L(z|\Theta ,{R_1})$ for above ${R_1}$ can be written as follows:

(34) \begin{align} L(z|\Theta ) = \frac{{{n_y}N}}{2} + \frac{N}{2}{\rm{ln}}[\det ({R_1})] + \frac{{N{n_z}}}{2}{\rm{ln}}(2\pi ) \end{align}

The maximum likelihood (ML) estimate of $\Theta $ can be obtained by minimising the following simplified cost function based $L(z|\Theta )$ :

(35) \begin{align} J(\Theta ) = \det (R) \end{align}

where $R = \frac{1}{N}\Sigma _{k = 1}^N{[z({t_k}) - y({t_k})]^T}[z({t_k}) - y({t_k})]$ is the residual co-variance matrix. Minimisation of $J(\Theta )$ is a nonlinear optimisation problem, which can be solved using the particle swarm optimisation (PSO) algorithm. The following steps are involved in implementing PSO:

STEP 1: (Initialisation) If $S$ is the search space for optimal solution ( ${\Theta _{gbp}}$ ), the possible solutions (positions of each particle) can be selected randomly as follows:

(36) \begin{align} {\Theta ^i}(0) \in S{\rm{,}}\forall i \in \{ 1,2,3,\ldots,{N_P}\} \end{align}

where ${\Theta ^i}(0)$ is ${i^{{\rm{th}}}}$ initial possible solution in the search space and ${N_P}$ is the number of possible solutions selected randomly. The initial personal best position ( $pbp$ ), the global best position ( $gbp$ ) and the velocity (change in position) of each particle are assigned as:

(37) \begin{align} \Delta {\Theta ^i}(0) = 0{\rm{,}}\forall i \in \{ 1,2,3,\ldots,{N_P}\} \end{align}

(38) \begin{align} \Theta _{pbp}^i(0) = {\Theta ^i}(0){\rm{,}}\forall i \in \{ 1,2,3,\ldots,{N_P}\} \end{align}

(39) \begin{align} {\Theta _{gbp}}(0) = \mathop {\min }\limits_{{\Theta ^i}(0)} \{ J({\Theta ^i}(0))\} {\rm{,}}\forall i \in \{ 1,2,3,\ldots,{N_P}\} \end{align}

where $\Delta {\Theta ^i}(0)$ , $\Theta _{pbp}^i(0)$ and ${\Theta _{gbp}}(0)$ are the initial changes in the position, personal best position and global best position of ${i^{{\rm{th}}}}$ particle.

STEP 2: (Position Update) Each particle’s velocity (change in position) is calculated using the previous personal best and the global best position of the population.

(40) \begin{align} \Delta {\Theta ^i}(k + 1) = w\Delta {\Theta ^i}(k) + {c_1}{R_1}(\Theta _{pbp}^i(k) - {\Theta ^i}(k)) + {c_2}{R_2}({\Theta _{gbp}}(k) - {\Theta ^i}(k)) \end{align}

(41) \begin{align} {\Theta ^i}(k + 1) = {\Theta ^i}(k) + \Delta {\Theta ^i}(k + 1) \end{align}

where $i$ and $k$ denote ${{\rm{i}}^{{\rm{th}}}}$ particle and ${{\rm{k}}^{{\rm{th}}}}$ iteration, respectively. $\Delta {\Theta ^i}(k)$ and ${\Theta ^i}(k)$ represent the change in the position of ${i^{{\rm{th}}}}$ particle at ${k^{{\rm{th}}}}$ iteration and the position of ${i^{{\rm{th}}}}$ particle at ${k^{{\rm{th}}}}$ iteration, respectively. $\Theta _{pbp}^i$ is the personal best position of ${i^{{\rm{th}}}}$ particle and ${\Theta _{gbp}}$ is the position of global best particle. $w$ , ${c_1}$ and ${c_2}$ are the inertial weight, the personal cognitive coefficient and the social cognitive coefficient, respectively. ${R_1}$ and ${R_2}$ are the diagonal matrices of random numbers from 0 to 1.

STEP 3: (Personal Best and Global Best Update) Personal best position ( ${\Theta _{pbp}}$ ) and global best position ( ${\Theta _{gbp}}$ ) are updated as follows:

(42) \begin{align} \Theta _{pbp}^i(k + 1) = \left\{ {\begin{array}{l@{\quad}l}{\Theta _{pbp}^i(k)} & {{\rm{if}}\,J({\Theta ^i}(k + 1)) > J(\Theta _{pbp}^i(k))}\\[5pt] {{\Theta ^i}(k + 1)} & {{\rm{if}}\,J({\Theta ^i}(k + 1)) \le J(\Theta _{pbp}^i(k))}\end{array}} \right. \end{align}

(43) \begin{align} \Theta _{gbp}^i(k + 1) = \left\{ {\begin{array}{l@{\quad}l}{\Theta _{gbp}^i(k)} & {{\rm{if}}\,J(\Theta _{pbp}^i(k + 1)) > J(\Theta _{gbp}^i(k))}\\[5pt] {\Theta _{pbp}^i(k + 1)} & {{\rm{if}}\,J(\Theta _{pbp}^i(k + 1)) \le J(\Theta _{gbp}^i(k))}\end{array}} \right. \end{align}

STEP 4: (Checking for Convergence) If the convergence criteria are met, terminate the iteration; else, go to STEP 2.

Brief convergence analysis of the proposed method and estimation of the confidence bound are given as follows [32]:

Definition 1: For a real-valued cost function $J(\Theta )$ , and ${\Theta ^ \star } \in S$ , if

(44) \begin{align} J({\Theta ^ \star }) \le J(\Theta ){\rm{,}}\forall \Theta \in S \end{align}

then ${\Theta ^ \star }$ is said to be a global optimal solution on $S$ , where $S$ is the search space for optimal solution.

Definition 2. $\forall \varepsilon > 0$ , let

(45) \begin{align} B_{\epsilon} =\left\{ \Theta \epsilon S \| J(\Theta)-J(\Theta^{\star})| < \epsilon \right\} \end{align}

Then $B_{\epsilon}$ is called a $\epsilon$ -optimal solution set, where ${\Theta}{\epsilon}B_{\epsilon}$ is said to be an -optimal solution.

According to PSO algorithm, the sequence $\{ J({\theta _{gbp}}(k))\} $ , $\forall k \in \{ 1,2,3,\ldots \} $ is a monotone-decreasing sequence. $p\{ {\Theta _{gbp}}(k) \in B_{\epsilon}\}$ and $p\{ {\Theta _{gbp}}(k) \in S\} $ are the probabilities of global best solution in $B_{\epsilon}$ and $S$ , respectively, at ${k^{{\rm{th}}}}$ iteration. Hence,

(46) \begin{align} \mathop {\lim }\limits_{k \to \infty } p\{ {\Theta _{gbp}}(k) \in B_{\epsilon}\} = 1 \end{align}

and

(47) \begin{align} \mathop {\lim }\limits_{k \to \infty } p\{ {\Theta _{gbp}}(k) \in S\} = 1 \end{align}

From the above two statements, it is evident that the probability of getting an optimal solution in the predefined search space is very high with a sufficient number of iterations.

Furthermore, the confidence in the estimates of unknown parameters using the ML-PSO method can be quantified in terms of the lower Cramer-Rao bound, and it can be found using diagonal elements of the Fisher information matrix ( $F$ ).

(48) \begin{align} F = \Sigma _{k = 1}^N{\left[ {\frac{{\partial y({t_k})}}{{\partial \Theta }}} \right]^T}{R^{ - 1}}\left[ {\frac{{\partial y({t_k})}}{{\partial \Theta }}} \right] \end{align}

(49) \begin{align} \sigma = \sqrt {{\rm{diag}}({F^{ - 1}})} \end{align}

Figure 2. Instrumented prototype of CDRW UAV [33].

Figure 3. Comparison of measured and estimated outputs of CDRW UAV in the low angle-of-attack flight regime.

Figure 4. Comparison of measured and estimated outputs of CDRW UAV in the moderately high angle-of-attack flight regime.

Figure 5. Comparison of measured and estimated outputs of CDRW UAV in the stall flight regime.

Figure 6. Comparison of measured and estimated outputs of CDRW UAV in the low angle-of-sideslip flight regimes.

Figure 7. Proof-of-match exercise for CDRW UAV.

4.0 Details of UAV and flight data acquisition system

A wing-alone cropped delta propeller-driven UAV, named CDRW, was used to generate various flight data sets. Its wing planform is designed with a root chord of 0.9m, a taper ratio of 0.167, and a wingspan of 1.5m. UAV’s pitch and roll motion are controlled using elevons of 0.125m chord and 0.45m span. In contrast, yaw motion is controlled by an all-movable dedicated vertical tail of 0.2m root chord, 0.08m tip chord and 0.42m span. The total takeoff weight of CDRW is 3.6kg, and a fully instrumented prototype of CDRW UAV is given in Fig. 2.

UAV was integrated with sensors to measure its response during flight tests, and a high-fidelity data acquisition system was used for data recording from various sensors and actuators. A high-accuracy 9 degree-of-freedom inertial measuring unit was used to measure components of the UAV’s linear acceleration, angular rates, and attitude angles. Initial bias error of accelerometer, gyro and magnetometers were $ \pm 0.002{\rm{g}}$ , $ \pm 0.25{\rm{deg/s}}$ and $ \pm 0.003\textrm{g}$ , respectively. Airflow angles and airspeed were measured with the help of vane-type in-house-fabricated sensors and pitot-static probe, respectively. Accuracy maintained by airflow angles, and airspeed sensor were $ \pm 0.1{\rm{deg}}$ and $ \pm 1.5{\rm{m/s}}$ , respectively. Pulse width modulated signals to the actuators of control surfaces and electronic speed controller for the motor to control the thrust produced were recorded and converted to their physical quantities using calibration relationships while post-processing.

5.0 Results and discussion

A total of eight compatible flight data sets generated using CDRW UAV pertaining to low angle-of-attack, moderately high angle-of-attack, stall angle-of-attack, and low angle-of-sideslip flight regimes are presented in Figs. 3–7 with legend ‘EXP’. These data sets are used for aerodynamic characterisation and validation of estimated parameters with the help of a proof-of-match exercise. Each flight data is represented with a specific name. The nomenclature of flight data is based on flight regimes. Abbreviations LG, NL, ST and LD represent flight data pertaining to the longitudinal linear, longitudinal nonlinear, longitudinal stall, and lateral-directional flight regimes, respectively. Simulated responses using OEM and ML-PSO have been denoted, in Figs. 3–7, with legends ‘OEM’ and ‘ML-PSO’, respectively.

Aerodynamic parameters of longitudinal linear flight regimes, given in Equation (24), are estimated with the help of the aerodynamic model presented in Equations (11)–(13) from the flight data sets LG1 and LG2 using the ML-PSO method. Estimated parameters are presented in Table 1, along with the measurements obtained from full-scale wind tunnel testing and OEM estimates. Though the control inputs are different for LG1 and LG2, the consistency in the estimates of ${C_{{L_\alpha }}}$ and ${C_{{m_\alpha }}}$ , can be observed from Table 1 for the CDRW configuration. Estimates of ${C_{{L_{{\delta _e}}}}}$ obtained from LG1 and LG2 using the proposed method are at the offset of 0.297 and 0.29 w.r.t wind tunnel results, whereas 0.099 and 0.02 offsets in ${C_{{L_{{\delta _e}}}}}$ estimates w.r.t wind tunnel results are observed with the OEM estimates. Relative errors of 0% & 0.2% and 2% & 0.5% are observed in estimates of ${C_{{m_{{\delta _e}}}}}$ w.r.t wind tunnel values from flight data LG1 & LG2 using OEM and ML-PSO, respectively. Although some of the estimated aerodynamic parameters have relative offsets with wind tunnel results, the simulated responses of the UAV using OEM and ML-PSO estimates match well with real flight data without any significant relative error that can be referred from Fig. 3.

Table 1. Longitudinal aerodynamic parameters of CDRW UAV at the low angle-of-attack

^*Values in the square brackets represent lower Cramer-Rao bounds.

Nonlinear aerodynamic parameters mentioned in Equation (25) and aerodynamic model represented by Equations (14)–(16) are used to characterise CDRW configuration in high angle-of-attack flight regime. A linear approximation is also made using the linear aerodynamic model given by Equations (11)–(13), to characterise UAV in nonlinear flight regime. Estimated aerodynamic parameters of CDRW configuration are presented in Table 2. It can be referred from the same table that estimates of pitching moment aerodynamic derivatives are not changing significantly with the linear and nonlinear aerodynamic model, which is also in good agreement with wind tunnel results given in Fig. 1 where ${C_m}$ vs. $\alpha $ plot is almost linear below 15deg angle-of-attack. Relative errors of 1.8 % and 4.7 % w.r.t wind tunnel values are observed in ${C_{{L_\alpha }}}$ estimates of CDRW obtained using a nonlinear aerodynamic model with OEM and ML-PSO method, respectively, whereas 5.7% and 5.1% relative errors are observed for the same UAV using the linear aerodynamic model with OEM and ML-PSO method, respectively. Estimated values of ${C_{{L_{{\alpha ^2}}}}}$ are −0.886 and −0.402 using OEM and ML-PSO method for CDRW, respectively. It can be referred that ${C_{{m_\alpha }}}$ estimates obtained for CDRW using OEM and ML-PSO method with the nonlinear aerodynamic model have relative offsets of 0.004 and 0.003 w.r.t wind tunnel values, respectively, whereas 0.001 and 0.001 offsets are observed for the same UAV using OEM and ML-PSO method with the linear aerodynamic model, respectively. Simulated outputs using estimated aerodynamic parameters obtained using OEM and ML-PSO method with the nonlinear aerodynamic model are compared with measured flight data, which are closely matched with measured flight data and that can be seen from Fig. 4.

Table 2. Longitudinal aerodynamic parameters of CDRW in moderately high angle-of-attack

^*Values in the square brackets represent lower Cramer-Rao bounds.

Aerodynamic parameters in stall flight regime and aerodynamic model given by Equations (26) and (17)–(20) are estimated using OEM and ML-PSO method, and estimated parameters of the UAV are presented in Table 3. Estimated values of $a$ , $\tau $ , ${\alpha ^ \star }$ , ${C_{{D_X}}}$ and ${C_{{m_X}}}$ for CDRW are 9.401, 14.187, 26.885deg, 0.09 and −0.055 using OEM, whereas, 2.368, 10.295, 25.61deg, 0.095 and −0.035 for the same UAV using ML-PSO method. From Fig. 5, it can be said that simulated ${C_D}$ , ${C_L}$ and ${C_m}$ using ML-PSO method estimates are closely matched with reconstructed data, whereas simulated outputs with OEM estimates have a relative mismatch with reconstructed data. A similar trend can also be seen in the stall hysteresis plot, where OEM predictions are relatively weak compared to the ML-PSO method.

Table 3. Longitudinal aerodynamic parameters of CDRW UAV in stall-flight regime

^*Values in the square brackets represent lower Cramer-Rao bounds.

The effectiveness of the proposed method is also analysed in low angle-of-sideslip flight regimes, where the vector of unknown parameters given in Equation (27) and aerodynamic model represented by Equations (21)–(23) are used for estimation. All the estimated aerodynamic parameters using the ML-PSO method are presented in Table 4 for CDRW UAV. From Table 4, it can be observed that estimated values of ${C_{{Y_\beta }}}$ using OEM and ML-PSO method have relative offsets of 0.059 and 0.005 and 0.051 and 0.072, from flight data LD1 & LD2, w.r.t wind tunnel values. Relative errors in ${C_{{l_\beta }}}$ of 10.9% and 10.9%, and 10.9% and 19.8% w.r.t wind tunnel values are observed using OEM and ML-PSO method from flight data LD1 and LD2. Estimated values of ${C_{{n_\beta }}}$ using OEM and ML-PSO method from data LD1 and LD2 have relative offsets of 0.001 and 0.001 and 0.001 & 0.002 w.r.t wind tunnel values. From Fig. 6, it can be observed that simulated outputs with OEM and ML-PSO method are in good match with measured flight data. A small apparent relative error w.r.t measured data can be seen in simulated acceleration along the body $y$ -axis.

Table 4. Lateral-directional aerodynamic parameters of CDRW UAV at low angle-of-sideslip

^*Values in the square brackets represent lower Cramer-Rao bounds.

Proof-of-match exercise is carried out with the average values of estimated aerodynamic parameters from two flight data sets of linear longitudinal and lateral-directional flight regimes. From Fig. 7(a), it can be seen that simulated outputs using OEM and ML-PSO method estimates are in very close agreement with measured flight data. Similarly, from Fig. 7(b), it can be noticed that simulated outputs using OEM and ML-PSO method estimates match well with measured flight data except the component of acceleration along the body $y$ -axis ( ${a_y}$ ). Simulated ${a_y}$ for CDRW using ML-PSO method estimates showing a maximum overshoot of $0.05{\rm{m/}}{{\rm{s}}^2}$ w.r.t measured data, while a maximum overshoot of $0.25{\rm{m/}}{{\rm{s}}^2}$ can be observed for CDRW using OEM estimates. Proof-of-match exercise is not performed in nonlinear and stall flight regimes due to limited flight data sets.

6.0 Conclusions

Current research demonstrates the applicability and effectiveness of the proposed method in estimating linear and nonlinear aerodynamic parameters of the cropped delta-wing unmanned aerial vehicle. A total of eight compatible flight data sets pertaining to linear, nonlinear and stall regimes are used for aerodynamic characterisation. A comparative study is performed for estimated aerodynamic derivatives belonging to different flight envelopes with wind tunnel testing and output error method. It is observed that longitudinal aerodynamic stability and control derivatives estimated using the proposed method have more relative offset than output error method w.r.t wind-tunnel values; however, simulated responses are consistent with measured flight data. A higher-order lift coefficient is used to model nonlinear aerodynamics before stall, and its estimate is satisfactorily justified by lower Cramer-Rao bound. While estimating stall characterising parameters, it is observed that the proposed estimation method is able to predict stall hysteresis with good accuracy than the output error method. It is noticed from lateral-directional aerodynamic estimates that static-stability derivatives obtained using the output error method are better consistent with wind-tunnel values than the proposed estimation method. A proof-of-match exercise is performed for linear longitudinal and lateral-directional flight regimes to ascertain more faith in the estimated parameters. One major advantage of using the proposed method over the output error method is that it is a non-gradient-based estimation method, which can reduce the computational burden. However, this can happen only when the search space of solutions, the number of search particles, and social and personal cognitive coefficients are optimal, which may lead to another research problem.