System identification of cropped delta UAVs from flight test methods using particle Swarm-Optimisation-based estimation

Abstract

In the era of Unmanned Aerial Systems (UAS), an onboard autopilot occupies a prominent place and is inevitable for many of their modern applications. The efficacy of autopilot heavily relies upon the accuracy of the sensors employed and the capability of the onboard flight controller. In general, aerodynamic behaviour and flight dynamic capabilities of Unmanned Aerial Vehicles (UAVs) govern the selection and the design of flight controllers. Precise modeling of linear aerodynamic characteristics from flight data can be achieved using many of the existing classical parameter estimation techniques such as Output Error Method (OEM), Equation Error Method (EEM), and Filter Error Method (FEM). However, all the classical methods may not be readily applicable for aerodynamic modeling in nonlinear flight envelopes. The current manuscript is an attempt to exploit the capabilities of the Artificial Intelligence (AI) technique, named Particle Swarm Optimisation (PSO), in combination with Least Squares (LS) cost function to perform linear as well as nonlinear aerodynamic parameter estimation. The aforementioned task is accomplished by considering flight data from manoeuvers pertaining to linear angles of attack, moderate and near stall flight envelopes of two different UAVs with cropped delta planform geometry. Parameters estimated using the proposed LS-PSO method are consistent with minimum standard deviation and are on a par with OEM estimates. The proposed LS-PSO method enhances the capabilities of LS-based EEM while estimating stall characteristic parameters, which was not possible with LS alone. The longitudinal and lateral-directional static parameters estimated from the full-scale wind tunnel testing of the two UAVs were also used to corroborate the results obtained from the flight data using the LS-PSO method.

Nomenclature

b

wing span in m

$\bar{c}$

mean aerodynamic chord in m

V

air speed in m/s

q

pitch rate in rad/s

p

roll rate in rad/s

r

yaw rate in rad/s

$I_X$

moment of inertial along body x-axis in $\mathrm{kg-m}^2$

$I_Y$

moment of inertial along body y-axis in $\mathrm{kg-m}^2$

$I_Z$

moment of inertial along body z-axis in $\mathrm{kg-m}^2$

$I_{XZ}$

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

S

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

g

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

$F_t$

thrust produced by engine in N

m

mass of UAV in kg

$C_L$

nondimensional lift force coefficient

$C_D$

nondimensional drag force coefficient

$C_Y$

nondimensional side force coefficient

$C_l$

nondimensional rolling moment coefficient

$C_m$

nondimensional pitching moment coefficient

$C_n$

nondimensional yawing moment 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_{D_0}$

drag force coefficient at zero lift

k

induced drag force correction factor

$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_{\delta_e}}$

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

$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

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

X

nondimensional distance of flow separation point

$a_1$

aerofoil static stall characteristics parameter

$C_{D_X}$

derivative of drag force coefficient w.r.t. X

$C_{m_X}$

derivative of pitching moment coefficient w.r.t. X

J

cost function

U

vector of independent states and control inputs

Greek symbol

$\alpha$

angle-of-attack in rad

$\beta$

sideslip angle in rad

$\rho$

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

$\phi$

roll angle in rad

$\theta$

pitch angle in rad

$\psi$

yaw angle in rad

$\tau_2$

hysteresis time constant

$\alpha^\star$

break point

$\delta_e$

elevator deflection angle in rad

$\delta_a$

aileron deflection angle in rad

$\delta_r$

rudder deflection angle in rad

$\Theta$

vector of unknown parameters

1. Introduction

Unmanned Aerial Vehicles (UAVs) are addressed with many names such as Remotely Piloted Vehicle (RPV), drone, robot plane, pilotless aircraft, to name a few, and are often employed to carry out various beyond the line-of-sight missions such as reconnaissance, command and control, deception, homeland security, combat surveillance, and arsenal delivery in military applications. With the advancement of communication technology, which is now at the common person’s reach, these UAVs also found their significance in industrial applications such as monitoring of intertidal reefs, drone journalism, fluvial remote sensing, assessment of existing structures, telemedicine, aerial imagery in natural resource management and survey of the landscape. Due to their wide applications and operations in civilian airspace, stringent safety measures are of paramount importance right from design-development to operations. According to the report on UAV accidents published by Defense Technical Information Center USA, 47% of UAV accidents during its flight occur due to human error. Autonomous capability complemented with a robust autopilot can significantly reduce human interference, which will enhance safe operations of UAVs even in civilian airspace. A stable autopilot with an efficient controller can be designed with the help of precise system dynamics modeling. Fixed-wing UAVs flight is heavily dominated by the associated aerodynamic forces and moments, which indeed are functions of flight conditions. Hence, detailed modeling of aerodynamic characteristics for various flight regimes is inevitable.

Aerodynamic characterisation enables researchers to develop a mathematical description of a flight vehicle’s associated aerodynamics, which can either be a linear or nonlinear function of non-dimensional aerodynamic parameters. Quantification of these non-dimensional parameters for various flight envelopes is well explored using estimation methods such as Equation Error Method (EEM) [1-5], Output Error Method (OEM) [6-22], Filter Error Method (FEM) [23, 24] and Artificial Intelligence (AI) method [11, 25, 26, 27, 28] from the flight data of manned aircraft for various predefined manoeuvers. However in the case of UAVs, the availability of research related to the aerodynamic characterisation from flight data is minimal due to their classified applications. EEM is one of the simplest and computationally efficient flight test methods to estimate flight vehicle aerodynamic stability and control derivatives from flight test data [29]. However, the formulation of EEM restricts its application in estimating aerodynamic parameters from flight data pertaining to near stall and high angle-of-attack manoeuvers. In contrast, the application of OEM and FEM can be extended to estimate the aforementioned parameters from flight data of highly nonlinear flight regimes [30, 31]. It is well observed that both methods require priory information about initial conditions for better convergence and confidence of solutions. On the contrary, the AI estimation method based on Neural Networks does not require priory information about initial conditions and can be used to characterise the aerodynamic behaviour of a UAV in linear and nonlinear flight regimes [29, 31]. In general, the Neural Networks estimation method is based on OEM requires gradient computation to update the weights of networks, which makes it a high computational effort demanding method. Moreover, the confidence in estimates heavily depends upon the training criteria of Neural Networks, which will not be the same for different data sets. In the recent past, it has been observed that the aforementioned limitation can be resolved by another AI estimation technique based on Particle Swarm Optimisation (PSO). PSO is a data-driven bio-inspired search technique for an optimal solution. This technique can be used to find a local and global solution by providing proper bounds to the search space. PSO has been widely used in aerospace applications such as structural optimisation, controller optimisation, and flight route optimisation [32-36]. However, it started getting researchers’ attention as a tool of aerodynamic parameters estimation of flight vehicles recently. In recent publications, PSO based on Maximum Likelihood (ML) cost function formulation has been used to estimate the linear longitudinal aerodynamic parameters of the VTOL aircraft, symmetrical projectile and UAVs from flight data[37-40]. It can be observed that the ML cost function requires the numerical integration of equations of motion for each swarm particle and every iteration, which can make the aforementioned method computationally less efficient. PSO-based estimation method can be made computationally more efficient by augmenting with Least Squares (LS) cost function formulation, which does not require integration of equations of motion for each iteration. Indirectly, the proposed LS-PSO method can also enhance the EEM method’s capability to estimate the stall aerodynamic parameters.

Flight data pertaining to manoeuvers of various flight envelopes, generated using two UAVs, is used to exploit the aerodynamic parameters estimation capability of the proposed LS-PSO method. The aforementioned wing alone configurations share similar cropped delta planform and differ in cross-section area. One UAV with rectangular cross-section is named Cropped Delta Flat Plate (CDFP), and the other with reflex aerofoil as wing cross-section is named Cropped Delta Reflex Wing (CDRW) UAV. Flight tests were performed for various predefined manoeuvers, and respective flight data was recorded using an onboard dedicated data acquisition system. Aerodynamic parameters are estimated using the proposed method from linear, nonlinear and near stall regime flight data. Estimates obtained from the proposed method are compared with Maximum Likelihood Estimation (MLE) and wind tunnel estimates. The rest of the paper is organised as follows. Section 2 describes the detailed dynamic, kinematic and aerodynamic mathematical structure of UAVs. Details of the PSO algorithm and problem formulation are presented in section 3. Section 4 presents model description, instrumentation on board of UAVs and flight data generation. Section 5 presents the results obtained from the proposed aerodynamic parameter estimation technique and their comparison with the classical MLE method. Section 6 concludes the effectiveness of the PSO algorithm-based aerodynamic parameter estimation technique in all four flight regimes, limitations of the proposed method and future research opportunities.

2. Mathematical modeling CDFP and CDRW UAVs

The six-DOF simulation model has been developed using rigid body dynamics equations of motion. These equations, in general, are coupled in nature; a set of decoupled equations of motion are used in the current research of aerodynamic parameter estimation pertaining to various flight envelopes. Equations (1)-(4) and (5)-(8) represent the longitudinal and lateral-directional dynamics of UAVs, respectively. The assumed frame of reference while deriving the equations of motion is given in Fig. 1.

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

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

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

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

(5) \begin{eqnarray}\dot{\beta} &=& -\frac{\rho S V}{2m}C_Y - \frac{F_t}{mV}\sin\beta + \frac{g}{V}\sin\phi - r \end{eqnarray}

(6) \begin{eqnarray}\dot{p} &=& \frac{1}{2} \rho S V^2 b [I_Z C_l + I_{XZ} C_n] \frac{1}{I_XI_Z - I_{XZ}^2}\end{eqnarray}

(7) \begin{eqnarray}\dot{r} &=& \frac{1}{2} \rho S V^2 b [I_{XZ}C_l + I_X C_n] \frac{1}{I_XI_Z - I_{XZ}^2}\end{eqnarray}

(8) \begin{eqnarray}\dot{\phi} &=& p \end{eqnarray}

where, V, $\alpha$ , q, $\theta$ , $\beta$ , p, r and $\phi$ are free stream airspeed, angle-of-attack, pitch angle, sideslip angle, roll rate, yaw rate and roll angle, respectively. Since these UAVs are propeller-driven with a single brushless motor and aligned to the body x-axis, the thrust force ( $F_t$ ) has only component along the X-axis of the body frame subsequently, and any external moment about the centre of gravity due to thrust force has been considered zero. S, $\bar{c}$ , m are wing planform area, mean aerodynamic chord and mass of UAV, respectively. $I_X$ , $I_Y$ and $I_Z$ are the mass moment of inertia about body X, Y and Z axes, respectively. $I_{XZ}$ is the product mass moment of inertia in the body XZ plane. $\rho$ and g are atmospheric air density and acceleration due to gravity at flight altitude, respectively.

Figure 1. Frame of reference.

Longitudinal aerodynamic parameters are estimated using the proposed LS-PSO method from the flight data pertaining to linear, nonlinear and near stall flight envelopes. From wind tunnel results Ref. [41], it is observed that two UAVs have a linear variation of aerodynamic coefficients with angle-of-attack from $-5$ to 10 deg, nonlinear variation from 10 to 16 deg and can be considered near stall (highly nonlinear) flight envelope if the angle-of-attack is more than 16 deg. Depending on the angle-of-attack achieved following aerodynamic models have been considered in order to estimate respective aerodynamic coefficients.

The following equations represent the linear low angle-of-attack ( $-5$ to 10 deg) aerodynamic model of UAV-

(9) \begin{eqnarray}C_L &=& C_{L_0} + C_{L_\alpha}\alpha + C_{L_q}\frac{q\bar{c}}{2V} + C_{L_{\delta_e}}\delta_e \end{eqnarray}

(10) \begin{eqnarray}C_D &=& C_{D_0} + kC_L^2 \end{eqnarray}

(11) \begin{eqnarray}C_m &=& C_{m_0} + C_{m_\alpha}\alpha + C_{m_q}\frac{q\bar{c}}{2V} + C_{m_{\delta_e}}\delta_e \end{eqnarray}

The non-linear high angle-of-attack (10–16 deg) longitudinal aerodynamics model of UAV can be represented as follows –

(12) \begin{eqnarray}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{eqnarray}

(13) \begin{eqnarray}C_D &=& C_{D_0} + kC_L^2 \end{eqnarray}

(14) \begin{eqnarray}C_m &=& C_{m_0} + C_{m_\alpha}\alpha + C_{m_q}\frac{q\bar{c}}{2V} + C_{m_{\delta_e}}\delta_e \end{eqnarray}

The following equations can give the quasi-steady stall model [42] –

(15) \begin{eqnarray}X &=& \frac{1}{2}\left[ 1 - \tanh\{a_1(\alpha - \tau_2\dot{\alpha} - \alpha^\star)\}\right] \end{eqnarray}

(16) \begin{eqnarray}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{eqnarray}

(17) \begin{eqnarray}C_D &=& C_{D_0} + kC_L^2 + C_{D_X}(1 - X) \end{eqnarray}

(18) \begin{eqnarray}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{eqnarray}

The following equations can give the lateral-directional model-

(19) \begin{eqnarray}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{eqnarray}

(20) \begin{eqnarray}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{eqnarray}

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

In aforementioned quasi-steady stall model, $a_1$ , $\tau_2$ , $\alpha^\star$ , $C_{D_X}$ and $C_{m_X}$ are considered as stall characteristic aerodynamic parameters, where $a_1$ represents aerofoil static stall characteristics parameter, $\tau_2$ is the hysteresis time constant, $\alpha^\star$ is break point and X is non-dimensional flow separation point. The effect of hysteresis on drag and pitching moment is modeled with $C_{D_X}$ and $C_{m_X}$ , respectively.

3. Parameter estimation methodology

3.1. Least square method

In general, mathematical model of a dynamic system in state space is given by following equations:

(22) \begin{eqnarray}\dot{x}(t) &=& f[x(t), u(t), \Theta] + Fw(t) \end{eqnarray}

(23) \begin{eqnarray}y(t) &=& g[x(t), u(t), \Theta] \end{eqnarray}

(24) \begin{eqnarray}z(t) &=& y(t) + Gv(t)\end{eqnarray}

where, f and g are assumed to be nonlinear, real-valued and differentiable functions, x(t) is state vector, u(t) is control input, y(t) is output vector and z(t) is measured output by sensors. In the above mathematical model, additive process and measurement noise are considered; however, these are considered zero for LS formulation. Further, the output is considered to be linearly dependent on states and control input, and can be given by the following equation:

(25) \begin{equation} Y = U\Theta + \epsilon\end{equation}

where, Y is measured response of size $N\times 1$ , U is vector of independent states and control input of size $N\times n$ , $\Theta$ is vector of unknown parameters of size $n\times 1$ and $\epsilon$ is modeling error.

LS cost function for above output equation can be given by following equation:

(26) \begin{equation} J (\Theta) = \frac{1}{2}\epsilon^T\epsilon = \frac{1}{2}[Y^T - \Theta^T U^T][Y - U\Theta]\end{equation}

Differentiation of above cost function w.r.t. to $\Theta$ leads to exact estimation of unknown parameters.

(27) \begin{equation} \hat{\Theta} = [U^TU]^{-1}U^TY\end{equation}

3.2. LS-PSO method

If the dependent variable (measured output) Y is a nonlinear function of U and $\Theta$ , exact estimation of unknown parameters is not possible using the above method. However, the current challenge of unknown parameter estimation can be addressed by an optimisation technique. In this paper, a solution to the aforementioned problem is proposed using PSO. Consider the following nonlinear output equation as follows:

(28) \begin{equation} Y = f(U,\Theta) + \epsilon\end{equation}

And the nonlinear LS cost function can be given as

(29) \begin{equation} J (\Theta) = \frac{1}{2}\epsilon^T\epsilon = \frac{1}{2}[Y^T - f(U,\Theta)^T][Y - f(U,\Theta)] \end{equation}

As PSO method search for global best solution of an optimisation problem based on given bounds of search space and can be used to optimize the above cast function w.r.t. $\Theta$ . The following steps are involved, while implementing PSO.

Step 1: (Problem definition) The population size (N), dimensions of particle position (D), inertial weight (w), personal cognitive coefficient ( $c_1$ ), social cognitive coefficient ( $c_2$ ), boundaries of search space and cost function need to be defined. Considered values of $c_1, c_2$ and N are 2, 2 and 50, respectively [43]. D is same as number of elements in $\Theta$ .

Step 2: (Initialisation) A swarm of particles, with problem definition, is generated using the randomisation method. Each particle holds a random position in a swarm. The velocity of particles also needs to be initialised.

Step 3: (Position and velocity update) The velocity of particles depends on their previous velocity, personal cognition and global cognition. The position of particles depends on their previous position and current velocity. The position and velocity of particles are updated as follows –

(30) \begin{eqnarray}w = 0.9^k \end{eqnarray}

(31) \begin{eqnarray}V_i(k+1) &=& wV_i(k) + c_1r_1\{P_{{best}_i}(k) - P_i(k)\} + c_2r_2\{G_{{best}_i} - P_i(k)\} \end{eqnarray}

(32) \begin{eqnarray}P_i(k+1) &=& P_i(k) + V_i(k+1)\end{eqnarray}

where, i denotes $i{\rm th}$ dimension of particle velocity and position. $P_i (k)$ and $V_i (k)$ are particle’s previous position and velocity, respectively. $P_i (k+1)$ and $V_i (k+1)$ are particle’s current position and velocity, respectively. w is inertial weight, $c_1$ is personal cognitive coefficient and $c_2$ is social cognitive coefficient. $P_{{best}_i}(k)$ and $G_{{best}_i}$ are previous best position of particle and best position of particle that have minimum cost in swarm.

Step 4: (Personal best and global best position update) In this step, the particle’s personal best position is updated with the current position if the cost function value associated with the current position is less than the cost function value associated with the previous best position. The global position is updated if any particle’s best cost value is less than the previous global best position’s cost value.

Step 5: (Termination) If the desired termination criteria satisfy stop the algorithm else, go to Step 3.

Since in the current problem, $C_L$ , $C_D$ , $C_m$ , $C_Y$ , $C_l$ and $C_n$ are dependent variables and can not be measured directly during flight test rather these are reconstructed using measured output data from various onboard sensors, which are given by following Equations (33)-(40). Now onward, these reconstructed outputs are considered as dependent variable (measured outputs) and ( $\alpha$ , $\beta$ , $V_\infty$ , p, q, r, $\delta_e$ , $\delta_a$ , $\delta_r$ ) as independent variables.

(33) \begin{eqnarray}C_X(i) &=& \frac{1}{\bar{q}(i)S}[m a_{X_{CG}}(i) - F_t(i)] \end{eqnarray}

(34) \begin{eqnarray}C_Z(i) &=& \frac{1}{\bar{q}(i)S} [m a_{X_{CG}}(i)] \end{eqnarray}

(35) \begin{eqnarray}C_L(i) &=& C_X(i)\sin\alpha(i) - C_Z(i)\cos\alpha(i) \end{eqnarray}

(36) \begin{eqnarray}C_D(i) &=& -C_X(i)\cos\alpha(i) - C_Z(i)\sin\alpha(i) \end{eqnarray}

(37) \begin{eqnarray}C_m(i) &=& \frac{1}{\bar{q}(i)S\bar{c}}[I_Y\dot{q}(i)] \end{eqnarray}

(38) \begin{eqnarray} C_Y(i) &=& \frac{ma_{Y_{CG}}(i)}{\bar{q}(i)S} \end{eqnarray}

(39) \begin{eqnarray}C_l(i) &=& \frac{1}{\bar{q}(i)Sb}[I_X\dot{p}(i) - I_{XZ}\dot{r}(i)]\end{eqnarray}

(40) \begin{eqnarray}C_n (i) &=& \frac{1}{\bar{q}(i)Sb}[I_Z\dot{r}(i) - I_{XZ}\dot{p}(i)] \end{eqnarray}

where, i denotes $i{\rm th}$ measurement. $a_{X_{CG} }$ , $a_{Y_{CG} }$ and $a_{Z_{CG} }$ are net acceleration along X-axis, Y-axis and Z-axis in body frame, respectively. $\bar{q}$ is dynamic pressure, $\dot{q}$ is pitch acceleration, $\dot{p}$ is roll acceleration, $\dot{r}$ is yaw acceleration.

In current research, vector of unknown aerodynamic parameters can be given as follows for different flight regimes.

(41) \begin{eqnarray}\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{eqnarray}

(42) \begin{eqnarray}\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{eqnarray}

(43) \begin{eqnarray}\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_1,\tau_2,\alpha^\star, C_{D_X}, C_{M_X}]^T\end{eqnarray}

(44) \begin{eqnarray}\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{eqnarray}

where, $\Theta_{LG}$ , $\Theta_{NL}$ and $\Theta_{ST}$ represents the vector of unknown parameters of longitudinal linear (low angle-of-attack), nonlinear (high angle-of-attack) and near stall flight regimes. While $\Theta_{LD}$ represents the vector of unknown parameters of lateral-directional flight regime.

The flow chart of the proposed LS-PSO method is given in Fig. 2. It can be referred from the flow chart that the proposed methodology takes control inputs, measured outputs, and a mathematical model to start. Based on inputs and outputs LS cost function is defined. A swarm of particles is generated using the cost function and initial conditions of the PSO algorithm in the following step. In the next step, iteration will start and end with termination criteria. These termination criteria can be a maximum number of iterations or minimum cost function value.

Figure 2. Flow chart of LS-PSO method.

4. Model description and flight data generation

As mentioned earlier, two cropped delta wing UAVs with similar planform geometry and differs in wing crossectional profile are used for flight data generation. Wing of CDFP and CDRW are designed with rectangular cross-section and NACA23110 reflex aerofoil, respectively. Directional control of both UAVs is achieved by all movable high aspect ratio dedicated vertical tail. These UAVs are controlled using ailevons, which are located at the trailing edge of the configuration. Pitch is controlled by symmetric deflection of control surfaces, and roll is controlled by differential deflection of control surfaces. Both UAVs have a wingspan of 1.5m, a root chord of 0.9m, a tip chord of 0.15m, a mean aerodynamic chord of 0.61m, an aspect ratio of 2.86 and the wing tapered about the wing trailing edge. Masses of CDFP and CDRW UAVs are 3.5 and 3.6kg, respectively.

Flight data during flight tests have been measured and recorded using onboard sensors and a dedicated data acquisition system. A high accuracy 9-DOF Inertial Measuring Unit (IMU) has been mounted to measure body accelerations, body rates, and Euler angles. In-house fabricated and calibrated air data boom is used to measure airspeed, angle-of-attack and sideslip angle. Pulse Width Modulated (PWM) signals are used to control the speed of brushless DC motor and deflection of control surfaces using servos, and these signals are logged into the data acquisition system. Figure 3 shows instrumented CDFP and CDRW prototypes for flight testing. Eight sets of flight data pertaining to linear, non-linear, near stall and lateral-directional flight regimes are used in current research. The nomenclature adopted for each flight data set is as follows-each data set name starts with CDFP and CDRW, followed by an underscore and two alphabets. LG stands for linear longitudinal flight data, NL non-linear longitudinal flight data, ST stands for near stall longitudinal flight data and LD stands for lateral-directional flight data. E.g. CDFP_NL means non-linear flight data of CDFP configuration. Figures 4 and 5 represent linear, nonlinear, stall and lateral-directional flight data for CDFP and CDRW configurations, respectively.

Figure 3. Instrumented CDFP and CDRW.

Figure 4. Various flight data sets generated using CDFP.

Figure 5. Various flight data sets generated using CDRW.

From Figs 4(a) and 5(a), it can be seen that in flight data CDFP_LG a doublet kind of elevator control input is used to exit the longitudinal dynamics, while a 3-2-1-1 kind of control input is used in generating CDRW_LG flight data. Elevator deflections vary from 0 to $-8$ deg and 2 to $-2$ deg in flight data CDFP_LG and CDRW_LG, respectively. Slow time-varying control input is used to generate lateral-directional flight data sets. It can be observed from Figs 4(b) and 5(b) that aileron deflection are 5 to $-10$ deg and $-5$ to 5 deg for CDFP_LD and CDRW_LD, respectively. From Figs 4(c) and 5(c), it can be referred that to generate nonlinear longitudinal flight data CDFP_NL and CDRW_NL, maximum elevator deflections of $-18$ and $-9$ deg are given to achieve about 20 deg angle-of-attack. Precise elevator control inputs are given to generate quasi-steady stall data. Similarly, from Figs 4(c) and 5(c), it can be observed that angle-of-attack varies slowly from low angle-of-attack to the high angle-of-attack and drops rapidly from high angle-of-attack to low angle-of-attack, which makes flight data CDFP_ST and CDRW_ST suitable for stall hysteresis modeling.

5. Result and discussion

The LS-PSO parameter estimation technique is applied to estimate the longitudinal and lateral-directional aerodynamic derivatives of CDFP and CDRW UAVs from flight data sets belonging to different flight regimes. From flight data sets CDFP_LG and CDRW_LG, it can be observed that angle-of-attack (AOA) falls within the range of $-2$ to 7 deg and $-2$ to 6 deg, respectively, hence aerodynamic model presented in Equations (9)-(11) is used to estimate unknown aerodynamic parameters given by Equation (41). A nonlinear longitudinal aerodynamic model given by Equations (12)-(14) and flight data sets CDFP_NL and CDRW_NL are used to estimate unknown aerodynamic parameters given by Equation (42). It can be observed from flight data sets CDFP_ST and CDRW_ST that the maximum AOA achieved during flight tests is more than 16 deg; hence near stall aerodynamic model presented in Equations (15)-(18) is used to estimate unknown aerodynamic parameters given by Equation (43). Lateral-directional flight data sets CDFP_LD and CDRW_LD along with aerodynamic model given in Equations (19)-(21) are used to estimate the aerodynamic parameter given in Equation (44). All the estimated aerodynamic parameters of the above-mentioned flight regimes are corroborated with wind tunnel results and MLE method estimates. Confidence in estimates is quantified in terms of standard deviation, and these are represented with square brackets “[ ]” in the following tables. Furthermore, computed (simulated) outputs, using estimated aerodynamic parameters and measured control inputs, have been compared with measured flight data. Simulated outputs and measured outputs during flight tests have been labeled with LS-PSO and EXP in proceeding figures.

All the longitudinal aerodynamic parameters mentioned in Equation (41) have been estimated using the LS-PSO method from flight data sets CDFP_LG and CDRW_LG. These aerodynamic parameters have been tabulated in Table 1 along with wind tunnel and MLE method estimates. Computed outputs using these aerodynamic parameters and measured control inputs have been compared with measured outputs during flight tests. It can also be referred from Table 1 that estimated $C_{D_0}, C_{L_0}, C_{L_{\alpha}}, C_{L_{\delta_e}}, C_{m_0}, C_{m_{\alpha}}$ and $C_{m_{\delta_e}}$ using LS-PSO method have relative error of 0%, 0%, 0.03%, 0.38%, 0%, 0% and 0% w.r.t. wind tunnel results for CDFP UAV. Similarly, estimated $C_{D_0}, C_{L_0}, C_{L_{\alpha}}, C_{L_{\delta_e}}, C_{m_0}, C_{m_{\alpha}}$ and $C_{m_{\delta_e}}$ using LS-PSO method have relative error of 0%, 6%, 0.3%, 6.2%, 0%, 0.4% and 0% w.r.t. wind tunnel results for CDRW UAV. From the same table, it can be noticed that standard deviations of estimated aerodynamic parameters are very small and LS-PSO method estimates are close to MLE method estimates, which further strengthen the efficacy of estimates. Comparison of simulated and measured responses can be seen from Fig. 6.

Table 1. Linear flight regime longitudinal aerodynamic parameters

Figure 6. Measured and estimated states of CDRW and CDFP in longitudinal linear flight regime.

Lateral-directional aerodynamic parameters have been estimated from flight data using the proposed method. All the estimated parameters have been tabulated in Table 2 along with wind tunnel and MLE estimates. It can be seen that estimated $C_{Y_0}, C_{Y_\beta}, C_{Y_{\delta_r}}, C_{l_0}, C_{l_\beta}, C_{l_{\delta_a}}$ , $C_{l_{\delta_r}}, $ $C_{n_0}, C_{n_\beta}$ and $C_{n_{\delta_r}}$ using LS-PSO method have relative absolute offset of 0, 0.041, 0.08, 0, 0, 0, 0.001, 0, 0.001 and 0.001 w.r.t. wind tunnel results for CDFP UAV. Similarly, estimated $C_{Y_0}, C_{Y_\beta}, C_{Y_{\delta_r}}, C_{l_0}$ , $ C_{l_\beta}, C_{l_{\delta_a}}, C_{l_{\delta_r}},C_{n_0}, C_{n_\beta}$ and $C_{n_{\delta_r}}$ using LS-PSO method have relative absolute offset of 0, 0.052, 0.03, 0, 0.011, 0.006, 0, 0., 0.001 and 0.002 w.r.t. wind tunnel results for CDRW UAV. Standard deviations of estimated parameters using proposed method are small which signifies good confidence in estimates. From Fig. 7, it can be referred that computed outputs, using LS-PSO method estimated aerodynamic parameters and measured control inputs, are having small relative error w.r.t. to experimental outputs. From Table 2, It can also be referred that LS-PSO method estimates are comparable to MLE.

Figure 7. Measured and estimated states of CDRW and CDFP in lateral-directional linear flight regime.

Table 2. Linear flight regime lateral-directional aerodynamic parameters

Nonlinear longitudinal aerodynamic parameters are estimated from high angle-of-attack flight data using the proposed method. A second-order aerodynamic derivative $C_{L_{\alpha^2}}$ has been used to model the nonlinear variation of $C_L$ w.r.t. angle-of-attack, whereas a linear aerodynamic model is used to estimate aerodynamic derivatives with MLE in Ref. [31]. All the aerodynamic parameters mentioned in Equation (42) have been estimated using LS-PSO method and tabulated in Table 3 along with MLE and Wind tunnel results. It can be referred from same table that estimated $C_{D_0}, C_{L_0}, C_{L_{\alpha}}, C_{L_{\delta_e}}$ , $C_{m_0}, C_{m_{\alpha}}$ and $C_{m_{\delta_e}}$ using LS-PSO method have relative absolute offset of 0, 0.01, 0.059, 0.101, 0, 0 and 0 w.r.t. wind tunnel results for CDFP UAV. Similarly, estimated $C_{D_0}, C_{L_0}, C_{L_{\alpha}}, C_{L_{\delta_e}}, C_{m_0}, C_{m_{\alpha}}$ and $C_{m_{\delta_e}}$ using LS-PSO method have relative error of 0, 0.007, 0.024, 0.144, 0, 0.001 and 0 w.r.t. wind tunnel results for CDRW UAV. From Fig. 8, it can be seen that computed outputs using LS-PSO method estimated aerodynamic parameters and measured control inputs are consistent with experimental outputs. From Table 3, it can also be referred that the standard deviations of estimated parameters are low, which indicates high confidence in estimates. Higher-order aerodynamic term $C_{L_{\alpha^2}}$ also has a standard deviation order of $10^{-3}$ , which provides extra confirmation about the efficacy of the proposed nonlinear longitudinal aerodynamic model.

Figure 8. Measured and estimated states of CDRW and CDFP in longitudinal nonlinear flight regime.

Table 3. Nonlinear flight regime longitudinal aerodynamic parameters

Quasi-steady stall aerodynamic parameters given in Equation (43) are estimated from flight data sets CDFP_ST and CDRW_ST using the proposed method. Linear flight regime $C_{L_\alpha}$ is considered while modeling quasi-steady stall. However, nonlinearity in the stall model is addressed by a nondimensional distance of flow separation point. All the estimated parameters are tabulated in Table 4 along with wind tunnel and MLE estimates. Since nondimensional aerodynamic coefficients are not measured directly, the reconstructed nondimensional aerodynamic coefficients from measured outputs are considered experimental outputs. From Figs 9(a) and 10(a), it can be observed that simulated outputs, using LS-PSO method estimated aerodynamic parameters and measured control inputs, are following the trend of measured flight data. The lowest value of the nondimensional distance of flow separation point achieved for CDFP and CDRW UAVs are 1 and 0.5, respectively, which indicates that flow separation happens at the leading edge of the wing in the case of CDFP UAV and at half of the mean chord in case of CDRW UAV. From Table 4, break points estimated using the LS-PSO method are consistent with wind tunnel and MLE estimates, while other stall characterizing parameters have offset. It can be seen from Figs 9(b) and 10(b) that the proposed method is able to predict the stall hysteresis. A higher-order standard deviation of estimated parameters is observed due to the presence of about 5% white Gaussian noise in data, which indicates low confidence in estimates.

Figure 9. Measured and estimated states of CDFP in near stall flight regime.

Figure 10. Measured and estimated states of CDRW in near stall flight regime.

Table 4. Stall flight regime longitudinal aerodynamic parameters

The nonlinear cost function mentioned in Equation (29) can also be solved with nonlinear optimisation techniques other than PSO such as Gauss-Newton Optimisation (GNO), Trust Region Optimisation (TRO), Newton-Raphson Ridge Optimisation (NRRO), Neural Networks (NN) and Evolutionary Algorithm (EA). The techniques GNO, TRO and NRRO rely on the gradient and Hessian matrix, which can be more computationally demanding than the search-based methods. The NN method requires rigorous training, which can be time-consuming. On the other hand, PSO is a search-based technique that does not require gradient or Hessian matrix calculation. In order to understand the computational burden using the LS-PSO method, the convergence of cost function values is plotted w.r.t. CPU time. A total of 50,000 measured data samples are considered to evaluate the performance of the proposed method on a computer running with Intel(R) CORE(TM) i5-10300H CPU @ 2.50GHz processor. It takes 4.8 seconds to estimate all the aerodynamic parameters mentioned in Equation (41) without parallel computing enabled and the same can be referred from Fig. 11.

Figure 11. Cost function values vs CPU time.

Two new flight data sets belonging to the longitudinal linear flight regime, which is not considered for parameter estimation, are used to test the prediction capability of the LS-PSO method. All the simulated responses of CDFP and CDRW UAVs are produced using estimated aerodynamic parameters mentioned in Table 1 and compared with measured flight responses. From Fig. 12, it can be observed that simulated $\alpha$ , $\theta$ , q, $V_\infty$ , $a_x$ and $a_z$ show good consistency with experimental outputs of CDFP UAV, while simulated $V_\infty$ is following the trend of the measured value with a slight offset in case of CDRW UAV. This analysis validates the acceptable prediction capabilities of the proposed method.

Figure 12. Proof-of-match performed in longitudinal linear flight regime.

6. Conclusion

Current research work is aimed to formulate and implement AI-based LS-PSO aerodynamic parameter estimation technique for aircraft. During its formulation, least-squares of the error cost function, augmented with reconstructed nondimensional aerodynamic coefficients as the dependent variable, is used, making it one of the computationally efficient method. In this paper, the LS-PSO method is proposed to estimate the linear and nonlinear aerodynamic parameters of UAVs. The aforementioned capabilities of the proposed method are demonstrated with flight data sets of two mini-cropped delta wing UAVs pertaining to the low angle-of-attack, moderate angle-of-attack, near stall and lateral-directional flight regimes. It is observed that most of the estimated aerodynamic parameters belonging to longitudinal linear flight regimes using the proposed method are close to wind tunnel results and are par with MLE method estimates, whereas nonlinear aerodynamic parameters estimates are comparable to MLE method estimates. Implementation of the LS-PSO method also covers the estimation of stall characterising parameters, which was a limitation of the ordinary EEM. The proposed method is able to model stall hysteresis of both UAVs satisfactorily. The efficacy of estimated aerodynamic parameters is quantified in terms of standard deviation. Confidence bounds (standard deviations) of estimates using the LS-PSO method are low, indicating reasonable belief in estimates. Further study is required to investigate its applicability in identifying coupled aerodynamic derivatives of unstable aircraft.