NEW HEURISTIC ALGORITHMS FOR ROLLING AIR FRAME AERODYNAMIC PARAMETERS ESTIMATION

In this research, based on heuristic optimization algorithms, three new strategies are developed for Aerodynamic Parameters Estimation (APE) of one pair ON-OFF actuator rolling airframe. In the 1st method namely EAM-PSO the aerodynamic parameters are directly estimated. While, the next two algorithms called EBM-PSO and SEBM-PSO are twostep strategies. In the 1st step the aerodynamic forces and moments are estimated, then after passing through a designed smoothing filter, in the 2nd step aerodynamic parameters are estimated. In EBM-PSO all the aerodynamic parameters are estimated at once by solving one optimization problem. In SEBM-PSO the APE is converted to solve four separate optimization problems. A modified particle swarm optimization algorithm is developed and used in estimation process. The performance of proposed algorithms is compared with that of state of the art algorithm EKF. The simulation results show that SEBM-PSO and EBM-PSO are better than EAM-PSO in term of accuracy and run time.


Introduction
An accurate aerodynamic model is crucial for developing high fidelity flight simulations and designing efficient control systems. Computer-based methods like theoretical-empirical relations and computational fluid dynamic (CFD) method as well as experimental approaches including wind tunnel and flight test are common ways to estimate the aerodynamic characteristics. The theoreticalempirical methods have a low accuracy relating to some simplifications done for their relations derivation. CFD methods are time consuming. Both of computer-based methods are sufficient for first design steps, but for more accurate results, we need to use experimental approaches. The wind tunnel test approach, itself, has uncertainties arising from scaled models and wind tunnel walls effect. It also has limitations in simulating the vehicle's whole flight envelope. Flight test method is the most accurate approach for all flight conditions. Therefore, during last decades, several approaches have been developed for APE based on flight test data. These approaches can be categorized as classic filter-based approaches, intelligent algorithms, and heuristics optimization methods.
There are various researches in literature utilizing classic filters for APE problem. The observer/Kalman filter identification method is applied to the problem of online system identification of an aircraft (Valasek & Chen, 2003). Based on real flight test data, three recursive parameter estimation algorithms including EKF, simplified Unscented Kalman Filter (UKF) and UKF are compared for APE of two aircrafts (Chowdhary & Jategaonkar, 2009). EKF based technique for the calibration of various parameters measured by air data system of an aircraft is studied in (Kamali et al., 2016). In Bayoğlu and Nalci (2016) an adaptive KF is applied to APE for a supersonic missile having rapid speed variations. A two-step methodology, studied in Moszczynski, Leung, and Grant (2019), in which an adaptive maximum aposteriori trajectory estimation scheme is adopted for accurate dynamic model identification. Two concepts namely Estimation After Modeling (EAM) and Estimation Before Modeling (EBM) are introduced for the first time by (Mohammadi & Massoumnia, 2000). According to EAM strategy, the aerodynamic parameters are directly estimated from the measured data. This is the most common approach have been used by researchers. While, in EBM method two extra sensors are needed and estimation is done in two steps. In that article, in first step, aerodynamic forces and moments were estimated using Frazier smoother, and in second step, APE is done using EKF. Classic filters are more convenient but they have some limitations: Kalman filter is optimal only for linear systems with Gaussian noise. EKF and UKF methods are used for nonlinear systems, but the problem associated with them is greatly attributed to their needed approximations in linearization process (Moszczynski, Leung, & Grant, 2019).
Intelligent algorithms like fuzzy and neural networks methods are knowledge-based techniques, can do tasks based on the data given for training or previous experiences. The learning process itself may take a very long time and there is usually no guarantee to success (Kumari, Sunita, & Smita, 2013). Intelligent algorithms like neural network and fuzzy control are also applied to the APE problem. Modeling and identification of spacecraft system using adaptive neuro-fuzzy inference systems is performed in (Hanafy, Al-Harthi, & Merabtine, 2014). Hatamleh et al. (2015) have focused on estimating unknown dynamics model parameters of an unmanned quadrotor under presence of noisy feedback signals, where iterative bi-section shooting method, artificial neural network method and a hybrid approach are used.
In current decade, nature inspired heuristic optimization algorithms methods are also implied to the parameter estimation problem. A review of intelligent algorithms including genetic algorithm, particle swarm optimization (PSO) and artificial bee colony algorithms is proposed in Wang et al. (2013) for aircraft APE. Taboo continuous ant colony system method (Nobahari & Sharifi, 2014), is used in Rezaei (2015) to estimate aerodynamic coefficients of a rolling airframe. Propeller z-force and pitching moment coefficients of a small unmanned aerial vehicle are estimated through a modified PSO in Tieying, Jie, and Kewei (2015), which in, accelerations, deflection angle, and pitch rate are taken as observations. A new optimization algorithm called adaptive chaotic mutation PSO is proposed to perform APE for a spinning symmetrical projectile (Guan et al., 2016), where, only aerodynamic drag and lift coefficients are estimated. Bian et al. (2016) presented an improved PSO method, which in every particle will do a local search before updating the global velocity and position. Then, the global best particle is created by a certain number of elitist particles. This algorithm is used to identify stability and control derivatives for an experimental small unmanned helicopter. The advantage of heuristics algorithms is their global search capability. But, they suffer from high computational time problem (Sone & Yadav, 2015).
The current research concentrates on developing new APE algorithms for a one pair ON-OFF actuator rolling airframe, based on heuristic optimization approaches. PSO algorithm is selected because of its capabilities in finding global optima in continuous optimization problems. Most APE algorithms are based on the conventional strategy namely EAM. In the current research, a new EAM based PSO algorithm, namely Estimation After Modeling Particle Swarm Optimization (EAM-PSO) is developed. In addition, two novel strategy namely Estimation Before Modeling Particle Swarm Optimization (EBM-PSO) and Separated Estimation Before Modeling Particle Swarm Optimization (SEBM-PSO) are developed and applied to the APE problem. EBM-PSO algorithm is a two-step strategy: In 1 st step, the aerodynamic forces and moments are estimated, then after passing through a designed smoothing filter, in 2 nd step, aerodynamic parameters are estimated. In EBM-PSO, the aerodynamic parameters are estimated via solving a single optimization algorithm based on a modified version of PSO. The derived formulation makes it possible to eliminate the numerical integration process and hence decreasing the computational time. Utilizing the symmetrical properties of rolling airframes, the aerodynamic parameters can be estimated by solving separate optimization algorithms. This resulted in SEBM-PSO algorithm. SEBM-PSO provides more exact results. The performance of newly developed algorithms and EKF, as a state of the art algorithm, are compared. It is shown that SEBM-PSO and EBM-PSO algorithms provide more exact results in comparison to those of EAM-PSO and EKF algorithms. EBM-PSO and SEB-PSO also have better run time than EAM-PSO because of eliminating the numerical integration during optimization algorithm iterations.
This paper is organized as follows: In section 1, the rolling airframe equations of motion are illustrated. Section 2 is devoted to the proposed modified particle swarm optimization (PSO) algorithm. In section 3, EAM-PSO method in addition to EBM-PSO method are explained in details. Simulation results are discussed in section 4, and finally conclusions are presented.

Airframe equations of motion
The vehicle of concern is a roll-stabilized airframe with one pair ON-OFF actuator. For modeling this vehicle, it is assumed that the thrust vector is fixed through the center of gravity and coincide with the body frame x-axis, the airframe is rigid, Vehicle has an axis-symmetric and cruciform shape, so that the moments of inertia I yy (t) and I zz (t) are identical and products of inertia moments can be discarded.
Since, vehicle is modeled as symmetric in pitch and yaw planes (Aksu, 2013), the aerodynamic coefficients C y and C n are identical to C z and C m respectively. Therefore, we have; , , , Taken into consideration the above assumptions, the equations of motion in the body-fixed frame are as follows; cos sin sin tan cos tan d p q sin cos cos cos where, u, v, and ware velocity components in body-fixed coordinate system, p, q, and r are vehicle angular rates, ϕ, θ, and ψ are Euler angles, D is body diameter, S is a reference area and Q is the dynamic pressure. In presence of wind [u wind , v wind . w wind ], total velocity V T , angle of attack α, and sideslip angle β can be determined as follows (Vitale, 2013): Here, it is assumed that wind effect is negligible in comparison to vehicle velocity and one can say:

Modified particle swarm optimization algorithm
The original PSO, first developed by (Kennedy & Eberhart, 1995), is inspired by the social behaviors observed in flocks of birds and schools of fish. The main feature of PSO is that it works in continuous search spaces. Considering the algorithm proposed by (Karimi, Pourtakdoust, & Nobahari, 2011), which is a PSO-based constrained optimization algorithm, a non-constrained version of this algorithm is developed here. It starts with an initial population of randomly generated particles. For a search problem in a D-dimensional space, a particle represents a potential solution presented by its velocity and position. During a search process, each particle is attracted by its previous best particle (X best ) as well as its neighborhood best particle (X lbest ). The swarm uses two equations to update each particle's velocity and position.
( ) ( ) where, V i,j and X i,j are the velocity and position of particle i in j th dimension of search space, k = 1,…, N itr is number of iteration index, i = 1,…,N p is number of particles index, j = 1,…,D is number of dimension index, X besti,j is local best position, i.e. it represents the best experience of individual particle, X lbesti is the neighborhood best position or the best experience of i th particle neighbors, c 1 and c 2 are two positive constants which represent the confidence for the experience of individual particle and neighbors, here c 1 = c 2 = 1, both r 1 and r 2 are two uniformly distributed random values between [0,1] that avoid falling into the local optimal, w is an inertial weight which represents the confidence of the previous velocity, here w is set to be uniformly distributed random values between [0,1]. To prevent particles from leaving the boundaries of search space, the velocity of particles is also clamped with a predefined maximum value V max . In this way, each dimension of the search space is normalized to vary between [0, 1] and maximum velocity of particles in each dimension of search space is set to V max = 0.2. Therefore, the velocity vector of the saturated particles is randomized as:

The neighborhood structure
There are various neighborhood structures like star, ring, wheel and Von-Neumann (Engelbrecht, 2005) and Karimi, Pourtakdoust, and Nobahari (2011) have used a singlylinked ring structure for the neighborhoods. In singlylinked ring structure, as shown in Figure 1, a particle k has Figure 1. Singly-linked ring structure two neighbors, particle k-2 and particle k+1. In turn, particles k-2 and k+1 have not particle k as one of their neighbors. This topology keeps two neighbors for each particle, but breaks the mutual attraction between neighbors. In this way, the information is transmitted faster through the whole swarm than in the original ring topology.

Quadratic interpolation operator
Inspired by the mathematical fact that each three points produce a quadratic curve which it may have a minimum value in its valley, a quadratic interpolation operator was introduced by (Karimi, Pourtakdoust, & Nobahari, 2011). The quadratic interpolation operator utilizes the current experience of the intended particle, termed as main parent, and two other randomly chosen particles of the swarm. Subsequently, the fitness of the new particle is compared with the best experience of its main parent. The main parent's best position will be replaced by the new particle, if it is better. Note that, in this strategy, the position of the parent is not changed and only its best position may be changed. According to the quadratic interpolation, the position of the new particle, generated by particle j and two random particles nb1 and nb2, at t-th iteration and j-th dimension of the search space is calculated as: (20) where, x i (t,j), x nb1 (t,j), x nb2 (t,j) and x new (t,j) indicate the j-th component of i-th, n b1 , n b2 and new particle position vectors at t-th iteration, respectively.

Proposed heuristic estimation algorithms
As previously mentioned, most researches for APE are based on the conventional strategy namely EAM. The EAM is a single step estimation strategy, directly estimate aerodynamic parameters. The EBM algorithm first proposed in (Mohammadi & Massoumnia, 2000) for APE in a two-step estimation strategy. In the 1 st step the aerodynamic forces and moments are estimated, then after passing through a designed smoothing filter, in 2 nd step aerodynamic parameters are estimated. In current research, the modified PSO algorithm, previously described, is combined with both EAM and EBM strategies, and new algorithms are developed.

EAM-PSO algorithm
In practice, the proposed EAM-PSO algorithm requires inertial navigation system measurement data including linear accelerations and angular rates in rolling body frame, , and control surface deflection angle for estimating control surface contribution in aerodynamic forces and moments. Therefore, in this research, after preparing the measurement data from 6DOF simulation model explained in section 4, an optimization problem is formed. The objective of this problem is to match measured values of [ , , , , , ] m xm ym zm m m m a a a p q r = y with their estimated values [ , , , , , ] x y z a a a p q r By defining estimation state vector as , and estimated aerody namic para meters vector as x z z zq l lp m m mq , rolling airframe equations of motion, Eq. (2)-Eq. (10), can be stated in the following form: As can be seen, here, the estimation model is stated as a function of estimated aerodynamic coefficients. The estimated accelerations are calculated as: The EAM-PSO algorithm flow chart and pseudo code are shown in Figure 2 and Table 1, respectively. According to this algorithm, in initialization step, Aero X  is randomly selected and the estimation model of Eq. (24) is numerically integrated by 4 th order Runge Kutta method (Tan et al., 2012) and the estimated data vector is calculated. Note that numerical integration must be done each time that the objective function J EAM is calculated. This makes EAM-PSO to an algorithm having high computational time. During each iteration of estimation process, velocity and position vectors are updated. After that, particles personal and local best positions are calculated. Stopping criterion is set to be a predefined number of iterations.

EBM-PSO algorithm
In EBM-PSO, in addition to INS, an extra accelerometer set is also needed. The first accelerometers set is fixed at the position [R x1 , 0, 0] T from the C.G, and supplies the measured data Substituting dp dt , dq dt and dr dt from Eq. (5)-Eq. (7) into Eq. (28) for both accelerometers yields (Larsson, 2013): (29) Or in compact form: Two steps of EBM-PSO are described in next subsections.

EBM-PSO First Step
In first step of EBM-PSO, the aerodynamic forces and should be estimated using Eq. (29). But the problem is that matrix A is not a full rank matrix (rank of A≠6). Rolling airframe symmetry characteristics causes this problem. Hence, by adding the relation of second row with that of fifth row, a new set of full rank equations will be found:

EBM-PSO smoothing step
The aerodynamic forces and moments calculated in first step are smoothed using Savitzky-Golay filter (Nirmal et al., 2016). In this filter a set of integers (A -n , A -(n-1) , …, A n-1 , A n ) are weighting coefficients to carry out the smoothing operation. The use of these weighting coefficients, known as convolution integers, turns out to be exactly equivalent to fitting the data to a polynomial, and it is computationally more effective and much faster. The smoothed data point by the Savitzky-Golay algorithm is given by the following equation: where, for a filter width (2n + 1) = 7 weighing coefficients are as follows: A -3 = -2, A -2 = 3, A -1 = 6, A 0 = 7, A 1 = 6, A 2 = 3, A 3 = -2.

EBM-PSO second step
In second step, the aerodynamic parameters are estimated by forming an optimization problem.
x z z zq l lp m m mq as aerodynamic parameters to be estimated (or optimization design variable), the estimation model can be summarized in the following vector form: The objective of this problem is to match calculated forces and moments with their estimated values, during the intended interval time, i = 1 … N, i.e.: Solving the optimization problem, all nine aerodynamic parameters Aero X  can be estimated. The EBM-PSO optimization algorithm showed in Figure 3. The pseudo code is also presented in Table 2. Noting that in first step F A and  M are estimated before aerodynamic modeling in second step, so this strategy is named EBM. Advantage of EBM-PSO algorithm is that numerical integration is removed in it, and this provides a low computational run time.

SEBM-PSO algorithm
In separated EBM-PSO or SEBM-PSO algorithm, the estimation model of Eq. (33) is separated into the following set of equations: Considering the above set of equations, the nine aerodynamic coefficients can be estimated by solving the following four separate optimization problems: First is an optimization problem with one-dimensional design vector [ ] x C  , and its objective function is formulated as: Second optimization problem has three design variables [ , , ] z z zq C C C α δ    , with the objective function: Third, has a two dimensional design vector l lp C C   , with the objective function: And finally, fourth problem has a three dimensional design vector [ , , ] m m mq C C C α δ    , having the objective function: SEBM-PSO is similar to EBM-PSO, but instead of one optimization problem, four separated optimization problems are formed.

Simulation results
A 6DOF simulation model, based on 6DOF equations of motion in Eq. (2)-Eq. (7), is built for generating the needed flight test data. In this model, thrust, mass, inertia and vehicle's dimensions are obtained from experimental tests, and the aerodynamic model is obtained from CFD and Missile DATCOM softwares.
As stated before, the control actuation system is a one pair ON-OFF. Therefore, an actuation signal δ(t) = sin(2pt) + Asin(pt+φ) is applied to the vehicle. Where, Asin(pt+φ) is input signal and sin(2pt) is the linearization signal (Nobahari & Mohammad Karimi, 2011).
The model outputs are: accelerations  Because of noisy nature of test sensore, in order to simulating real sensor outputs, a rational noise is also added to simulated data. Therefore, here, a uniformly distributed random noise is added to the original signal (sig m ), obtained from simulation (Tieying, Jie, & Kewei, 2015); (1 (2 1)) m m where, ξ determines noise to signal ratio, and rand stand for a uniformly distributed real number between [0, 1]. The noise effects are added according to Eq. (43), and adding noise to the data obtained from simulation will make the simulation data more close and equivalent to the real experimental data. Noise effect is considered by setting ξ = 0.05 for acceleration and angular rate data.
Sampling rate of measurement signal is set to be 500 Hz, and the simulation data are saved during interval time 5-6 seconds from flight regime. The selected flight regime is a high speed regime with limited Mach number variations (1.68-1.7 Mach). So, the aerodynamic coefficients are supposed to be constant. As we know, in high Reynolds numbers the viscosity effect are small. In other words, during this nearly constant speed flight condition, effect of aerodynamic phenomenon like viscosity, Mach numbers and Reynolds number are negligible.
As the current research is focused on the estimation algorithms, using this model will not affect the efficiency of the proposed methods. Even it can be used as a benchmark to validate the proposed algorithms. At last, simulation outputs is saved and used as the measurement data. All routines run on a PC with a 2.0 GHZ CPU and 4.0 GB of RAM.

APE with EAM-PSO
Estimation model of Eq. (23) is used in EAM-PSO. Figure  4 shows the EAM-PSO objective function convergence behavior during 100 iterations, for a single run.

APE with EBM-PSO and SEBM-PSO
Eq. (28) Table 3 summarizes the EAM-PSO results. It shows that the average value of estimation error is about 5.1%. The best estimated coefficient is C lp with 1.2% error, and the worst is C x with 9.3% error. On average, each run of EAM-PSO with 100 iterations takes 92 minutes run time. As can be seen, applying smoothing filter has deleted high frequency noises. Then second step of EBM-PSO have been done for both EBM-PSO and SEBM-PSO. The estimated and nominal values of aerodynamic coefficients C x , C zq , C lp and C mα , for a single run, are compared in Figure 13 to Figure 16. Summarized results of a single run are presented in Table 4. For EBM-PSO, the average value of estimation error is 2.4%, the best estimation result is for C lp with 0.5%, and the worst result is the estimation of C zα with 4.6% error. SEBM-PSO results shows that the average value of estimation error is 2.2%, the best estimation is for C l0 with 0.2% error, and the worst result is the estimation of C mq with 5.5% error. As shown, SEBM-PSO gives more accurate results than EBM-PSO because of solving separated optimization problems. Here, in addition to three newly developed algorithms, an EKF algorithm is also implemented for the problem at hand. The EKF state vector is X = [u, v, w, p, q, r, φ, θ, ψ], its unknown parameters vector to be estimated is Θ= [C x ,C zα ,C zδ ,C zq , C l0 ,C lp , C mα ,C mδ ,C mq ] and measurement vector is Y=[a X1 , a Y1 , a Z1 , a X2 , a Y2 , a Z2 , p , q , r]. The results of these four algorithms, for 100 successive runs, are compared and referred to nominal values calculated by CFD and Missile Datcom simulations and proposed in Table 7. It can be seen that the most accurate estimations are obtained by SEBM-PSO. The average estimation error for SEBM-PSO is 2.2%, for EBM-PSO is 2.4%, for EAM-PSO is 5.1% and for EKF is 8.4%. Three proposed algorithms give better estimations rather than EKF. However EKF gives a better estimation for C lp . In term of run time, EKF has the best performance with average 180 seconds. After that, EBM-PSO and SEBM-PSO have a comparable computational run time, while EAM-PSO has a high run time (92 min). Both EBM-PSO and SEBM-PSO algorithms have considerably better run-time in comparison with EBM-PSO, because the numerical integration process is deleted in EBM strategy. Best value for each aerodynamic parameter is highlighted in Table 5. SEBM-PSO gives best accuracy for C x ,C zα , C zδ , C l0 , C mα and C mδ .The best value for C zq is given by EAM-PSO. The best estimation for C mq is given by EBM-PSO. The best estimation error for C lp is obtained by EKF algorithm.

Measurement noise effect on EBM-PSO performance
The noise amplitude effect on EBM-PSO performance is studied here by performing 100 successive runs with three different noise magnitudes. The results are shown in Table 6. As was expected, increasing the noise amplitude causes the estimation accuracy to decreases. In another simulation, the effect of noise magnitude when smoothing filter is deleted from EBM-PSO is evaluated. The average values of estimated aerodynamic coefficients resulted by 100 successive runs are presented in Table 7. By comparing results of Table 6 it can be seen that Savitzky-Golay smoothing filter help the algorithm to decrease noise effects. For example when ξ = 0.05, using smoothing filter enhances the average estimated error from 6.2% to 3.5%.

Conclusions
In this article, based on a modified version of particle swarm optimization algorithm, three heuristic estimation algorithms are proposed to perform aerodynamic parameter estimation of a typical rolling airframe. It is shown that EBM-PSO uses an extra accelerometer in practice and the measurement unit needs to be fixed far from C.G. position. It was shown that EBM-PSO is more rapid in run time due to canceling the time consuming numerical integration procedure which is needed in EAM-PSO algorithms. SEBM-PSO algorithm provides more exact results by separating the estimation problem to a set of low-dimensional optimization problems. All aerodynamic coefficients may be estimated by proposed algorithms at once. Comparing the proposed algorithms performance with that of EKF shows its more exact results, while having comparable run time in EBM-PSO and SEBM-PSO algorithms. The evaluation studies show that aerodynamic parameter estimation accuracy is affected by measurement noise, and applying of Savitzky-Golay smoothing filter is very useful for eliminating the noise effect and enhancing the estimation accuracy. The simulation results revealed that the proposed methods can be used in practical applications. The performed work may be extended to involve heuristic real-time estimation algorithms that consider aerodynamic coefficients variations with flight parameters like Mach number.