A small fixed-wing UAV system identification using metaheuristics

Abstract A novel method for system identification of small-scale fixed-wing Unmanned Aerial Vehicles (UAVs) using a metaheuristics (MHs) approach is proposed. This investigation splits the complex aerodynamic model of UAV into longitudinal and lateral dynamics sub-systems. The system identification optimisation problem is proposed to find the UAV aerodynamic and stability derivatives by minimizing the R-squared error between the measurement data and the flight dynamic model. Thirteen popular optimisation algorithms are applied for solving the proposed UAV system identification optimisation problem while each algorithm is tested for 10 independent optimisation runs. By performing the Freidman’s rank test, statistical analysis of the experiment work was carried out while, based on the fitness value, each algorithm is ranked. The outcomes demonstrate the dominance of the L-SHADE algorithm, with mean R-square errors of 0.5465 and 0.0487 for longitudinal and lateral dynamics, respectively. It is considered superior to the other algorithms for this system identification problem.


Introduction
Aviation has been being incorporated in many aspects of human life and has the potential to develop into an automotive Unmanned Aerial Vehicles (UAV) in recent future for various applications (Naser & Kodur, 2020), such as agriculture, mapping, military, surveying, etc. Since the UAVs do not require pilots, UAVs can be a better choice for hazardous conditions. The design of UAVs for a peculiar operation is often associated with numerous complicated functions that are hard to develop. For better and favorable UAV control system design, a designer needs to develop a strong mathematical model to incorporate the UAV's characteristics. The modeling approaches assist in finding the aerodynamic, stability and control derivatives of the UAVs, which can be accomplished by either of the following ways: wind tunnel test, computational fluid dynamics (CFD), and a system identification technique. The wind tunnel approach is prone to certain error which is typically associated with the design of the UAV structure and also results in higher investment cost. The CFD analysis requires a great level of expertise in UAV, furthermore, its computation is challenging along with its computational cost. On the other hand, the system identification technique is easy to implement and provides accurate results compared to the others. Therefore, the system identification technique is the most interesting choice among the 3 techniques of interest concerning mathematical modeling for fixed-wing UAVs. Moreover, the accuracy of a mechanical model depends on numerous factors, such as noise, operating condition, time delay, set of data, and processing procedure. The controlling of these factors is challenging and needs to be addressed carefully. Therefore, a system identification or parameter estimation technique was introduced in order to improve the precision and accuracy of the UAVs mathematic model (Gim et al., 2018;Korkmaz et al., 2013;Sanders et al., 2018;Wei et al., 2019).
The system identification can be determined in three ways: 1. Equation Error Methods (EEM; Morelli, 2006) 2. Output Error Methods (OEM; Kabaila, 1983) and 3. Filter Error Method (FEM; Grauer & Morelli, 2015). The EEM is easy to implement and offers sufficiently accurate models. However, the noise in data must be filtered out beforehand because using the data with high error yields a misleading system model. OEM compares the actual data from measurement with the achieved data from the system model, subsequently, it endeavors to reduce the difference between these two models to produce the most precise and accurate system model. FEM is the technique that reduces the error caused by the system noise that is coming from unknown signals and environment. The main tools for FME are Kalman Filter (KF; Durrant-Whyte, 2001), Extended Kalman Filter (EKF; Hoshiya & Saito, 1984), and Unscented Kalman Filter (UKF; Wan & Van Der Merwe, 2000). Unfortunately, all of the three methods require a prerequisite parameter which in certain cases it is difficult to identify. Nowadays, efficient system identification techniques based on artificial intelligence (AI) have been proposed such as Artificial Neural Network (ANN; Bagherzadeh, 2018;Fekih et al., 2007;Kirkpatrick et al., 2013) and Fuzzy set theory (Ghosh Roy & Peyada, 2017;Srivastava et al., 2019) while applying of optimisation tool such as meta-heuristic algorithms (MHs) in combination with AI has also been presented (Li & Yin, 2014). However, those techniques still required aircraft models. Therefore, to address this issue, the new approach of mathematical modeling of small fixed-wing UAVs is proposed based on using MH optimisers incorporating with the idea of Output Error Methods.
MH is a non-gradient-based optimisation method widely used in various engineering design optimisation (Abdullah-Al-Shafi et al., 2016;Buch et al., 2017;Digehsara et al., 2020;Edoziuno et al., 2020;Shuaibu Hassan et al., 2020). MHs are very popular nowadays and offer a highly flexible solution range for optimisation problems, due to its derivative independent nature and capability to solve sophisticated problems with the least computational efforts. The advantage of MH is that they do not require function derivatives, so they can deal with any kind of objective functions and design variables. However, a lack of search convergence and search consistency are inevitable in the MH due to some randomization in the search proposed. In this regard, numerous MH have been developed for various applications to overcome this issue. For the mathematical modeling of fixed-wing UAV, having MH techniques embracement is somewhat limited. Therefore, it becomes a potential research gap that can be addressed effectively. In this study, the 13 metaheuristics are used to optimize the error between the mathematical model and the actual value. Furthermore, comparative performance analysis for each algorithm is performed to optimize the parameter estimation. The 13 algorithms investigated in this experiment are the Ant Lion Algorithm (ALO; Mirjalili, 2015a)

Formulation of a system identification optimisation problem
An inverse problem-based optimisation for a small UAV system identification using MHs is presented in this study. Several established MHs are applied to solve the proposed system identification optimisation problem while their performance are investigated. To conduct the proposed scheme and investigate MHs performance, firstly, a flight test is conducted for flight characteristic data collection. Then, linear aircraft models for both longitudinal and lateral dynamics are selected. The model parameters are identified using MH optimisers by solving the proposed optimisation problem that posed to minimize the R-square error of the flight characteristic from the flight test and that obtained from the aircraft state space model.

Flight test
The fixed-wings UAV model mentioned in this paper is majorly made of thermocol foam sheets. The core of the UAV is made of lightweight wood to strengthen the UAV as shown in Figure 1.
Multiple parts of UAV are attached with the sensor for the quantification of aircraft body velocity, angular velocity, Euler angles (for state variables in accordance with air speed), roll rate, pitch rate, yaw rates, and barometric altitude information. Raspberry Pi 3 B+ and Arduino mega2560 microcontrollers programmed by MATLAB/Simulink are used to collect the flight test data (Figure 2). The parameters of fixed-wing UAV are shown in Table 1.
With this experiment, the aerodynamic stability state variable vectors that can be measured from the flight test are u p θ h ½ � T and q r ; ψ ½ � T while the inputs applied for flight data collecting are the multi-step input and the doublets input. The multi-step input is given to the elevator in order to identify the longitudinal modes of the aircraft whereas the doublets input is used in throttle, ailerons, and rudder. The data collection process took 15 seconds shown in Figure 3.

Aircraft model
In this study, the six degrees of freedom aircraft equations of motion can be linearized around a steady-level flight trim condition . The resulting equations have the well-known state space form: The matrices A and B contain the aerodynamic and stability derivatives and the actuator effectiveness coefficients, respectively. The mechanical model has been experimented under the assumption that there is a small external environmental factor that has no impact on system performance. For simplicity, state variables in the mechanical model are constructed by linear mathematical models. The longitudinal dynamics are decoupled from the lateral dynamics.

Aircraft longitudinal system
The measurement of aerodynamic stability of longitudinal dynamics has a state variable x long ¼ u wq θ h ½ � T and the control vector u long ¼ δ elevator δ throttle ½ � T . The longitudinal derivative can be used to construct the state space representation as:

Aircraft lateral system
The measurement of aerodynamic stability of lateral dynamics has a vector of state variables The lateral derivative can be used to construct the state space representation as: From Equation (2), we assume that the angle of attack (θ � ) and side-slip angle (ϕ � ) is small (θ � ,ϕ � ≈ 0, 0). Therefore, we can write Equation (2) in a simplified form as follows: It should be noted that the variables X, Y and Z are force components in x, y and z, directions while L, M and N are rolling, pitching, and yawing moments. The subscripts denote their derivatives with respect to the subscripted variables.

System identification optimisation problem
The aircraft system identification optimisation problem can be formulated as an inverse optimisation problem that is posed to find the dynamic stability derivative as shown in Equations (1) and (2) in order to minimise the R-square error between output data from the flight test and the mathematical model. Formulation of the optimisation problem can be expressed as follows.
where x m;i is the i th output response from the flight test x s;i is the i th output response from the mathematical model simulation � x m;i is the average value of x m;i t final is time operating = 15 seconds t is current time. Note that, the time step is set to be 5 milliseconds.
The output data used for calculating the objective function are u p θ h ½ � T and q r ; ψ ½ � T for longitudinal dynamics and lateral dynamics, respectively. The design variables to be found are all dynamic stability derivatives as shown in (1) and (2)

Metaheuristics used in this study and numerical experiment set up
The MH were introduced before the 21 th century (Whitley, 1994). Prior to that, linear programming, nonlinear programming, and dynamic programming were three principal methods that were widely implemented for solving optimisation problems. However, there were still several problems that could not be solved by these methods as they are derivative dependent, computationally complex, and consecutively in the research study, they lost their position. Thus, to tackle numerous non-linear, derivative-free, and multi-modal challenging problems, novel optimisation methods were researched and developed. The optimisation methods were majorly inspired by the natural system behavior and have a systematic randomization system called Metaheuristics. Due to their simplicity, computational ease, and flexibility of such methods, they have seen a plethora of applications across various disciplines including UAVs and control system design. The process of MH is composed of two elements: design variables and objective function and the search process is executed by four steps to determine the optimal solution. The search process begins with an initial population or group of designed vectors initiated by randomization, followed by generating a child population (offspring) by applying the mechanism that imitates the behaviours of nature. Subsequently, the selection process begins to average the next generation's population. The probability of selection depends on the objective function or optimized solution of the member of the population shown in Figure 4.  In this study, 13 MH were tested by optimizing error between data from the flight test and the mathematical model of UAV, to determine the aerodynamic stability for both longitudinal and lateral dynamics. These algorithms parameters are set as follows: (1) The ant lion algorithm (ALO): The parameter settings are not required.
(2) Dragonfly algorithm (DA): The parameter settings are not required. (13) Adaptive differential evolution with optional external archive (JADE): The parameters are self-adapted during an optimisation process.
For each algorithm, the population size is taken as 100 while the maximum iteration is set as 1000 which results in 100,000 function evaluations for longitudinal dynamics. For the lateral dynamics, the population size is accounted as 100 while the maximum iteration is set as 400 which results in total 40,000 function evaluations. Each algorithm is tested by 10 independents optimisation runs.

Results and discussion
The performance of the considered algorithms is examined by observing error reduction capability between the mathematical model and the experimental data for both longitudinal and lateral dynamics of the UAV which can be verified by the proposed objective function of R-square error in (3). After performing 10 independent runs of all MHs, the mean values of objective function are used to measure the MHs search convergence while the standard deviation (STD) values are used to measure MHs search consistency. The lower mean and STD values is the better convergence and consistency, respectively. In addition, the Friedman test (Derrac et al., 2011) is used find the significant differences of the metaheuristic algorithms. The null hypothesis of the Friedman test states that medians among algorithms are equality. If the null hypothesis is rejected, the smaller Friedman score is the better efficient algorithm. In this work, statistical Friedman tool is employed for MHs ranking based on 95% significant level.  Figure 5 illustrates the best outcomes of the top three algorithms viz L-SHADE, SHADE, and CMAES concerning their fitness values for longitudinal dynamics analysis. It is the plot progress towards iteration value. The curve depicts a more descending fitness curve earlier, however, at a much slower rate than the early stage. To refine the results at large iteration value, the plot has been zoomed in and depicted with a small image. This illustrates the CMAES algorithm fitness value continuously descends at higher iteration while the SHADE and L-SHADE algorithm achieved their optimum value quite earlier. The performed assessment affirms that the L-SHADE outperforms others in terms of search convergence and consistency. However, a higher iteration or function evaluation could make the solution more accurate. Table 5 demonstrates the optimal solution obtained by the best run of the L-SHADE algorithm when exposed to the proposed fixed-wing UAV model for longitudinal dynamics. This solution has a value within the specified scope. The state space matrix and linear transfer function matrix obtained from those stability derivatives are shown in the Laplace domain s in Equation (4). This mechanical model is tested by 3-2-1-1 input for the elevator and doublets input for throttle at the same time for 15 seconds. The obtained solution of the L-SHADE algorithm from the modeling system has a minimal deviation against the actual flight. The modeling system can track the simulation system with precise every state variable shown in Figure 6.    Table 6 shows the results of each algorithm in solving the considered UAV problem compared to other algorithms for the lateral dynamics where the lower values of the mean, STD, and minimum values demonstrate the better solution. From the table, it is evident that the L-SHADE algorithm outperforms other algorithms in terms of Freidman ranking, minimum value, and mean value. The L-SHADE algorithm has the minimum fitness and mean value of 0.0475 and 0.0487, respectively, while the SHADE algorithm achieved the lowest STD value of 3.8 × 10 −4 . The L-SHADE algorithm is followed by the SHADE, CMAES, WCA, and MFO algorithms in terms of Freidman ranking. The SHADE, CMAES, WCA, and MFO algorithms are almost similar in every aspect except the minimum Figure 6. Time response match between flight data and modeling system signals. value of the SHADE algorithm is still higher than the other three. The rest of all other algorithms failed to track an acceptable solution and still have a high error. Similar to the longitudinal dynamic case, GOA is settled at the last rank compared to other considered algorithms. Figure 7 shows the lateral dynamics fitness value of the top four algorithms' best solutions. The fitness curves demonstrate a better convergence rate at initial iterations while continuing its trend at a slower rate at a high number of iterations. For a clear picture at high function evaluations, the obtained plot is zoomed in and the small Fig. illustrates the continuous descending fitness values of CMAES and SHADE algorithms at higher iteration. However, the L-SHADE algorithm successfully achieved the optimal solution at the initial level of iteration before CMAES, SHADE, and WCA. Moreover, from the small plot it can be perceived that at higher iteration or function evaluation, there is a likelihood of getting a more precise solution. The best solution achieved by the L-SHADE algorithm for the mechanical model of lateral dynamics is shown in Table 7. This solution has a value within the specified scope. The state space matrix and linear transfer function matrix obtained from those stability derivatives are shown in the Laplace domain s in Equation (5). This mechanical model is tested by doublets input for aileron and rudder command at the same time for 15 seconds. The solution of the L-SHADE algorithm has the lowest error between flight data and the modeling system signals. The modeling system can track the simulation system with very precise every state Variable shown in Figure 8

Conclusions
A new approach of system identification of the small fixed-wings Unmanned Aerial Vehicle (UAV), equipped with the low-cost sensor, is successfully demonstrated. They are successfully applied to solve the proposed system identification optimisation problem. The 13 distinguished MH are employed to solve the proposed system identification optimisation problem in order to identify the aerodynamic stability of the UAV. Overall, the best algorithm for this case study is the L-SHADE algorithm compared to the others for both longitudinal and lateral dynamics ranked based on the Freidman's statistical test where the L-SHADE algorithm has a minimum R-square error of 0.5465 and 0.0487 for longitudinal and lateral dynamics, respectively. Moreover, the fitness curve also demonstrates the dominance of the L-SHADE algorithm over the other contrasted algorithms. The best obtained aerodynamic and stability derivatives from the L-SHADE are adequately accurate  and precise to be applied in reality. The trend of output response from the mathematical modeling can precisely track the simulation system for every state variable. The improvement and modification of the algorithm could increase solution consistency. Moreover, noise and disturbance will be further explored in the future research work. X u ; X w and X q = stability derivative of force component in X-direction derivative with respect to u; w; and q; respectivly X δe and X δt = control derivative of force component in X-direction derivative with respect to δ e and δ t ; respectively Y v ; Y p ; and Y r = stability derivative of force component in Y-direction derivative with respect to v; p; and r; respectively Y δa ; and Y δr = control derivative of force component in Y-direction derivative with respect to δ a and δ r , respectively Z u Z w and Z q = stability derivative of force component in Z-direction derivative with respect to u; w; and q; respectively Z δe ; and Z δt = control derivative of force component in Z-direction derivative with respect to δ e and δ t ; respectively M u , M w ; and M q = stability derivative of pitching moment derivative with respect to u, w, and q, respectively M δe and M δt = control derivative of pitching moment derivative with respect to δ e and δ t , respectively L v ; L p ; and; L r = stability derivative of rolling moment derivative with respect to v; p; and r; respectively L δa , and L δr = control derivative of rolling moment derivative with respect to δ a and δ r , respectively N v , N p , and N r = stability derivative of yawing moments derivative with respect to v; p; and r; respectively N δa and N δr = control derivative of yawing moments derivative with respect to δ a and δ r , respectively 1 Sustainable and Infrastructure Research and