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

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. [Reference Saderla41], 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 [Reference Fischenberg and Jategaonkar42] –

(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 [Reference Kennedy and Eberhart43]. 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.