Longitudinal Motion System Identification of a Fixed-Wing Unmanned Aerial Vehicle Using Limited Unplanned Flight Data

Abstract

Acquiring knowledge of aircraft flight dynamics is crucial for simulation, control, mission performance and safety assurance analysis. In the fast-paced UAV market, long flight testing campaigns are hard to achieve, leaving limited controlled flight data and a significant amount of unplanned flight data. This work delves into the application of system identification techniques on unplanned flight data when faced with a shortage of dedicated flight test data. Based on a medium-sized, fixed-wing UAV, it focuses on the system identification of longitudinal dynamics using structural routine flight test data of pitch down and pitch up manoeuvres with no specific guidelines for the control inputs given. The proposed solution uses first- and second-order parameter-based models to build a non-linear dynamic model which, using a least square error optimisation algorithm in a time domain formulation, has its parameters tuned to converge the model behaviour with the real aircraft dynamics. The optimisation uses a combination of pitch, altitude, airspeed and pitch rate responses as a measure of model accuracy. Very significant improvements regarding the UAV model response are found when trimmed flight manoeuvres are used, resulting in proper estimation of important aerodynamic and control derivatives. Pitching moment and control derivatives are shown to be the crucial parameters. However, difficulties in estimation are shown for untrimmed flight manoeuvres. Better results were obtained when using multiple manoeuvres simultaneously in the optimisation error metric, as opposed to single manoeuvres that led to system bias. The proposed system identification procedure can be applied to any fixed-wing UAV without the need for specific flight testing campaigns.

1. Introduction

In a broad context, system identification is defined by the usage of techniques to characterise a system from which a mathematical model can be created that represents the same system. In aircraft design, system identification can be performed using simulation tools, real data analysis and processing, or a mixture of both. It is used to build a virtual model capable of predicting the aircraft flight dynamics based on the control inputs provided to the system to allow for behavioural analyses of an aircraft, and it is relevant in regard to simulating any flight conditions, evaluating the overall static and dynamic stability of the aircraft, reducing certification flight testing needs, and creating feedback loops to reinvent and correct design decisions.

System identification procedures in the aviation field are usually associated with flight data processing. However, system identification encompasses all these processing studies, including all the simulated environment techniques used to gather knowledge on the aircraft behaviour. These can also be used to either develop an initial model for the system or to perform trade-off studies regarding the design of an aircraft before flight test data are available.

Examples of such research include a comparative study between F16 and F22 aircraft using high-fidelity computational aerodynamics [1], a preliminary aerodynamic model development based on wind tunnel test data and similar aircraft databases [2], different aerodynamic analysis methods spanning from empirical to high-fidelity for the determination of the lift and drag profiles of an aircraft [3], and geometrics-only-based identification [4].

Meaningful and validated data, however, can only be verified and obtained using flight data, where an existing model can be compared directly to its real counterpart. Several different work approaches have been used to perform system identifications based on collected flight data. This includes a classical planned flight testing analysis comparing model and real data [5,6,7,8]; usage of big data and artificial neural networks [9,10,11]; machine learning techniques used to train a model to behave as the real aircraft [12]; aerostructural system identification assuming a flexible aircraft scheme [13]; remote uncorrelated pilot input excitation [14]; and aero-propulsive modelling [15]. Other approaches have tried to continuously innovate system identification techniques using new methods, such as the usage of a weighted recursive least squares algorithm for parameter identification [16] or a wavelet-based approach [17].

System identification validation or usage requires real-world data; thus, the data collection techniques and processing prove to be significantly connected with the accuracy of the results obtained. Extensive studies have been created on the best practices in flight testing experiments for the construction of classical aircraft dynamic models [18,19,20]. Similarly, data acquisition plays an important role in creating a validated system model. Poor data availability is shown to affect model estimation [18,19,21].

Longitudinal motion identification is of the utmost importance in terms of reliably assessing the design capabilities and addressing the differences between design and real aircraft longitudinal motion, bridging the gap between what is simulated and what is real behaviour. Assessing the real longitudinal motion of the aircraft can result in better performance estimation [22,23] and/or minimisation of flight test requirements for certification purposes [24]. Additionally, it provides important data about the system that allows for reliably assessing the design capabilities and addressing the differences between design and real aircraft models, allowing for better judgement regarding future aircraft generations.

In this work, the longitudinal dynamics of a medium-sized fixed-wing UAV are estimated using a parameter estimation error with a least square error optimisation scheme. Due to the fast-paced UAV market, access to planned flight test data is incompatible with product demands, so a solution is presented to cover the gaps in system identification methodologies. This work does not aim to present a completely new idea to the system identification world but rather to prove that applying established system identification methods can work even if the only flight data available are unplanned or scarce.

Contrary to the popular implementations of system identification methods [5,6,7,8,9,10,12,13,14,15], the presented formulation aims to prove that similar identification techniques can be used even if dedicated data availability is difficult to obtain. This study shows that using routine flight data is sufficient for acquiring important knowledge of the longitudinal motion of the system. This study proves that system identification techniques can be used in unplanned/unprepared flight data, even if these data have been collected prior to the development of the system identification scheme and with no control input guidance whatsoever.

In traditional system identification techniques used for aircraft, flight experiments feature meticulously planned and standardised inputs—such as frequency sweeps, single impulses, and doublets—that are applied to ensure accuracy in identifying aircraft dynamics. For instance, Dorobantu et al. [2] used frequency mapping techniques to generate reliable data for model validation, and even online estimation approaches, like those used by Valasek and Chen [25], rely on controlled manoeuvrers to achieve accurate, real-time estimates. Such methods typically demand dedicated flight tests with well-defined control inputs to ensure data quality and repeatability.

In contrast, this work demonstrates that system identification does not inherently require controlled manoeuvrers or dedicated flight campaigns. By applying standard identification techniques to unplanned, operational flight data, this study shows that meaningful longitudinal dynamics can still be accurately captured. This approach reveals the potential to derive essential aerodynamic characteristics of the aircraft from routine flight data, proving that the constraints of planned testing are not mandatory. Furthermore, it suggests a more flexible and cost-effective method of system identification that utilises historical or operational data, reducing reliance on costly and time-intensive dedicated flight tests, from plan to execution.

Using the unplanned flight data of consecutive pitch up and pitch down manoeuvres (used in structural flight validation) and easily available dynamic data of the system, a dynamic model is built and optimised to minimise the error between collected and simulated data.

Furthermore, key aspects of system identification are identified as error sources, and an overview of best practices and common pitfalls are also presented, with considerations being given to the best guidelines for proper system identification when there is a lack of suitable dedicated flight testing data.

The studies hereby presented were based on the Tekever AR5 aircraft [26]. Thus, due to the sensitivity of the data, the results are presented in non-dimensional form. However, this does not affect the reliability of the analysis and its reproducibility in similar aircraft with unplanned flight data.

This work is organised as follows: Section 2 and Section 3 overview the background theory used in this work, including the equations of motion, system identification techniques, and aircraft system modelling; Section 4 explains the aircraft data gathering and the simulation environment; Section 5 determines the parameter estimation formulation used, along with the validation scheme used and the formulation of the optimisation problem; Section 6 presents the optimisation results for both single and multi-manoeuvre cases; and Section 7 presents the conclusions found in this work.

2. Background

A brief overview of the most relevant subjects required to understand this work is presented.

2.1. Equations of Motion

The aircraft dynamic response can be described by the equations of motion, assuming a rigid body with 6 degrees-of-freedom (DOF), expressed as $\dot{\overrightarrow{A}}=\overrightarrow{f}(\overrightarrow{A},\overrightarrow{C})$, where $\overrightarrow{A}$, $\overrightarrow{C}$ and $\overrightarrow{f}$ are, respectively, the state vector, control vector, and the vector-valued non-linear function of the individual states and controls [27].

For the overall understanding, in this work, the Earth-fixed frame (^E) is considered as a global reference system, with its origin being in the center of the Earth and its axes (north, east and down) being fixed in relation to Earth surface; additionally, the body-fixed frame (^B) is considered as a local reference, with its origin being in the center of gravity of the aircraft, its longitudinal axis pointing towards the nose, its vertical axis pointing down, and its lateral axis being to the side.

The Earth-fixed frame (ECEF) is a global reference system that has its origin at the center of the Earth, with axes that are fixed relative to the Earth’s surface. In contrast, the body-fixed frame is centered on a specific body (like a planet or satellite) and moves along with that body. The Earth-fixed frame is typically used for global positioning and navigation, while the body-fixed frame is used for local analysis of bodies in motion.

For notation, the position, linear and angular velocities, forces, and moments are, respectively,

$$R^E=\left[\begin{array}{c}\hfill x^E\\ \hfill y^E\\ \hfill z^E\end{array}\right],U^B=\left[\begin{array}{c}\hfill u\\ \hfill v\\ \hfill w\end{array}\right],{\mathsf{\Omega }}^B=\left[\begin{array}{c}\hfill p\\ \hfill q\\ \hfill r\end{array}\right],F^B=\left[\begin{array}{c}\hfill X\\ \hfill Y\\ \hfill Z\end{array}\right],M^B=\left[\begin{array}{c}\hfill L\\ \hfill M\\ \hfill N\end{array}\right]$$

where the superscripts E and B refer to the Earth-fixed or body-fixed frame, respectively. Following the standard definitions, the body frame is centered at the aircraft center of gravity (CG), having its $x^B$-axis pointing to its nose and the $z^B$-axis pointing downward.

For the equations of motion, Newton’s second law is used [28,29]. It states, for translational motion, that the sum of external forces $\overrightarrow{F}$ is equal to the change in linear momentum and, for rotational motion, that the sum of external moments $\overrightarrow{M}$ regarding the CG is equal to the rate of change in angular momentum,

$$\overrightarrow{F}={\displaystyle \frac{d}{dt}}{\left(m\overrightarrow{U}\right)}_E={\displaystyle \frac{d}{dt}}{\left(m\overrightarrow{U}\right)}_B+\overrightarrow{\mathsf{\Omega }}\times {\left(m\overrightarrow{U}\right)}_B,$$

(1)

$$\overrightarrow{M}={\displaystyle \frac{d}{dt}}{\left(\left[I\right]\overrightarrow{\mathsf{\Omega }}\right)}_E={\displaystyle \frac{d}{dt}}{\left(\left[I\right]\overrightarrow{\mathsf{\Omega }}\right)}_B+\overrightarrow{\mathsf{\Omega }}\times {\left(\left[I\right]\overrightarrow{\mathsf{\Omega }}\right)}_B,$$

(2)

where m and $\left[I\right]$ represent the mass and inertia tensor, respectively.

If an Earth-fixed frame is used for the position and a body-fixed frame is used for the attitude of the aircraft, neglecting the structural dynamic behaviours, the final set of equations for the dynamic system are the rotation rate equations,

$$\left\{\begin{array}{c}\dot{\varphi }=p+qsin\varphi tan\theta +rcos\varphi tan\theta \hfill \\ \dot{\theta }=qcos\varphi -rsin\varphi \hfill \\ \dot{\psi }=q{\displaystyle \frac{sin\varphi }{cos\theta }}+r{\displaystyle \frac{cos\varphi }{cos\theta }}\hfill \end{array}\right.,$$

(3)

where $\psi$, $\theta$ and $\varphi$ are the Euler angles that correspond to yaw, pitch, and roll; the force equations, which describe the translational motion, neglecting wind, appear as

$$\left\{\begin{array}{c}X-mgsin\theta =m(\dot{u}+qw-rv)\hfill \\ Y+mgsin\varphi cos\theta =m(\dot{v}+ru-pw)\hfill \\ Z+mgcos\varphi cos\theta =m(\dot{w}+pv-qu)\hfill \end{array}\right.;$$

(4)

and the moment equations, which describe the rotational motion assuming a symmetric aircraft ($xz$-plane), and include the gyroscopic effects from the propulsion system with angular momentum $(h_x,h_y,h_z$), present as

$$\left\{\begin{array}{c}L=I_{xx}\dot{p}+I_{xz}\dot{r}+(I_{zz}-I_{yy})qr+I_{xz}pq+h_{z}q-h_{y}r\hfill \\ M=I_{yy}\dot{q}+(I_{xx}-I_{zz})rp+I_{xz}(r^2-p^2)+h_{x}r-h_{z}p\hfill \\ N=I_{xz}\dot{p}+I_{zz}\dot{r}+(I_{zz}-I_{xx})pq-I_{xz}rq+h_{y}p-h_{x}q\hfill \end{array}\right..$$

(5)

Note that the derived equations of motion assume, as mentioned, a rigid-body motion, which neglects the aerostructural behaviour of an aerodynamically loaded aircraft during flight. Evidently, the impact of structural flexibility is important to completely identify the dynamics of an aircraft [13]. However, for small- to medium-sized UAVs with robust composite material structures, such as the one seen in this work, the structural flexibility is assumed to have a lower impact on the overall longitudinal qualities of the system.

The forces and moments are generally described using non-dimensional coefficients as

$$\begin{array}{c}{C}_X={\displaystyle \frac{X}{q_{d}S}},C_Y={\displaystyle \frac{Y}{q_{d}S}},C_Z={\displaystyle \frac{Z}{q_{d}S}}\\ C_l={\displaystyle \frac{L}{q_{d}Sb}},C_m={\displaystyle \frac{M}{q_{d}Sc}},C_n={\displaystyle \frac{N}{q_{d}Sb}}\end{array},$$

(6)

where $q_d={\displaystyle \frac{1}{2}}\rho {\left|U^B\right|}^2$, S, c, and b represent the dynamic pressure, reference area, reference chord, and reference span, respectively. For convenience, the aerodynamic coefficients are retrieved in the stability frame, in which the x axis points towards the longitudinal velocity, as the stability derivatives exhibit minimal variations with the trim state [30] where the body frame is rotated by the angle of attack $\alpha$ about $y^B$. Using these aerodynamic coefficients (lift and drag), the transformation from the body-fixed frame to the stability frame becomes

$${\left[\begin{array}{c}\hfill C_X\\ \hfill C_Y\\ \hfill C_Z\end{array}\right]}_B=\left[\begin{array}{ccc}-cos\left(\alpha \right)& 0& -sin\left(\alpha \right)\\ 0& 1& 0\\ sin\left(\alpha \right)& 0& -cos\left(\alpha \right)\end{array}\right]\left[\begin{array}{c}\hfill C_D\\ \hfill C_Y\\ \hfill C_L\end{array}\right].$$

(7)

This work assumes a state and control function dependency of each force and moment coefficient. For a generic coefficient in the longitudinal motion, its formulation becomes

$$C_A=f_0\left(\alpha \right)+f_1\left(\dot{\alpha }\right)+f_2\left(q\right)+f_3\left({\delta }_e\right)+f_4\left({\delta }_t\right),$$

(8)

where $f_i$ represents the state and control dependency function, with the angle-of-attack $\alpha$, its rate $\dot{\alpha }$, and the pitch rate q being the relevant state variables and the elevator deflection ${\delta }_e$ and throttle command ${\delta }_t$ elements acting as the relevant control variables. The longitudinal velocity u is not considered due to its effect being highly linked with a propulsion system that is assumed externally validated. The purpose of this work is to find the set of linear or non-linear functions $f_i$ that minimise the differences between simulated and real aircraft-dynamic behaviour.

2.2. System Identification

System identification is the process of determining the mathematical model of a real system that best predicts the outputs of a system that has been given certain inputs [18]. This process involves four steps [31]: (1) data collection; (2) selection of model structure; (3) selection of system identification scheme; and (4) model optimisation.

Data collection is related to the data acquisition techniques used to retrieve the real system behaviour, namely the input and output quantities of interest. It not only depends on the data acquisition systems, i.e., sensors, but also on the methodology used to obtain the data. As explained by Milliken [32], the quality of system identification depends on the information stored in registered outputs, which, in turn, depends on the inputs used to obtain this information. Typically, the data retrieved can be processed in the time domain [5,6,7,8], frequency domain [33,34], or a combination of both [35]. Depending on the formulation intended, knowing the a priori behaviour, using a simplified model of the system can help excite the needed dynamic response, thus enabling a better and more accurate system identification.

The choice of the model structure is one of the most important decisions in system identification methods [31]. These are the mathematical principles used to model the digital twin of the real system. Choosing an incorrect model can lead to a partially described system [31]. Different approaches have been used to create the model structure of fixed-wing models: linear force and moment coefficient models with non-linear dynamic models [5,6,7,8,15,21], usage of structural deformation factors due to the presence of flexible wing phenomena [36], and a physics-based neural network approach [10]. This work focuses on polynomial-based equations to describe the influence of the state and control variables on the forces and moments acting on the aircraft, which are used to feed a non-linear dynamic model equation set.

The next most prominent aspect of system identification is the scheme selection. The most common methods used for aircraft identification are output error methods and equation error methods [18]. While an equation-error method focuses on the dynamic equation description of the system where the force and moment errors between model and real system are minimised, an output error focuses on the observable data of the system, i.e, the state variables error between simulated and real systems [18,19]. Since state variable output is more easily accessible and force and moments calculations depend on the gathering of these quantities from the real system (and feature a high amount of noise), the output error method is used here, where a free choice on the observable data and weight given to each state is implemented.

The final step is the minimisation of the differences between real and simulated models, hence the usage of the optimisation term. In this work, an automatic correction of the digital twin model is used where sequential dynamic simulations of the system are computed until the changes conducted on the model are sufficient enough to minimise the difference between the observable data of the real system and the simulation data. Figure 1 depicts an overview of the system identification methodology used in this work.

3. Preliminary Aircraft Model

To enhance the quality of the system identification and allow for a better a priori understanding, an initial digital model of the aircraft is created using low-fidelity aerodynamic tools, propulsion data gathered from flights, and point mass inertia descriptions. This leads to the collection of the initial aerodynamic coefficients of the model.

The aircraft model used in this work is a physics-based model where each force and moment is described as a series of functions that depend on aerodynamic coefficients given by the aerodynamics of the aircraft. Other methods, such as the usage of machine learning [12] or deep learning techniques [9,10,11], are also a possibility, but these present a high dependency on the availability of data, usually requiring flight data points to the order of thousands or tens of thousands [11]. This work aims to present a model that can be identified using a shortage of data—in the order of tens or even less. As such, the used physics-based model is deemed to be more appropriate. The description of the model construction is explained here.

3.1. Aerodynamic Modelling

The Tekever AR5 has a conventional wing-fuselage configuration with a H-tail. The first aerodynamic assessment for its dynamic behaviour was collected through the usage of the Athena Vortex Lattice (AVL v3.40) [37], a flight-dynamic analysis software for rigid aircraft of arbitrary configuration. The computational model is based on the vortex lattice method (VLM), a potential model with the assumptions of inviscid, incompressible, and irrotational flow. However, viscous corrections can be provided, and these were collected using XFoil v6.99 [38]. The aircraft model used for the analysis is presented in Figure 2, where a gap wing scheme is used to overestimate the lift reduction region due to the fuselage. The aerodynamic model is assumed to be close enough to reality to allow for the exploration of the aerodynamic coefficients.

3.2. Propulsion Modelling

The Tekever AR5 has a twin-engine configuration with clockwise rotating propellers when seen from the aircraft tail. These were modelled using the internal engine and propeller modules in JSBSim, along with an experimental correction for the initial estimated curves for power and thrust of the propeller. The dynamic model includes the roll moment and angular momentum induced by the propulsion system caused by the rotating engines.

The lack of wind tunnel data to describe the propeller in flight creates a gap in our understanding of the propulsion system behaviour. To mitigate this and gain knowledge beyond the static bench thrust, initial flight test data for different flight regimes (cruise, climb and descent) and different engine working points were used.

Using the initial aerodynamic model of the aircraft, a correlation was built so that, for similar flight conditions, the same output (power, thrust, RPM and fuel consumption) was achieved. Assuming trimmed flight conditions, the equilibrium of forces yields

$$\left\{\begin{array}{c}{C}_{X}sin\left(\theta \right)-C_{Z}cos\left(\theta \right)+C_{T}sin\left(\theta \right)-C_W=0\hfill \\ C_{X}cos\left(\theta \right)-C_{Z}sin\left(\theta \right)+C_{T}cos\left(\theta \right)=0\hfill \end{array}\right.,$$

(9)

where $C_T={\displaystyle \frac{T}{q_{d}S}}$ and $C_W={\displaystyle \frac{W}{q_{d}S}}$ are the thrust and weight coefficients, respectively. Knowing the X and Z forces behaviour of the aircraft based on the aerodynamic model and the real flight data of trimmed flight conditions, the thrust component can be retrieved and corrected in the propulsion model. The power, RPM, and fuel consumption are then calibrated to match the real system output.

This model correction creates a coupling between the aerodynamic model and the propulsive system behaviour, as it becomes very dependent on the drag of the aircraft. Thus, the propulsion model is only as accurate as the aerodynamic model used to correct it. Nonetheless, this approach assures that, for the same range of flight data points, the same equilibrium conditions are found.

The propulsion model ultimately converges to the collection of the power and thrust coefficient curves ($C_P$ and $C_T$) used by JSBSim to evaluate both the thrust and power in the simulated environment.

3.3. Inertia Modelling

Point masses are used to describe the inertia of the aircraft, assuming that the mass body inertia of each component about its center of mass is negligible when compared to its contribution to the inertia about the aircraft CG. The fuel point mass is allowed to change position according to the fuel volume used in each condition, emulating the real fuel center of mass. Thus, the resulting moments and products of inertia are, using the parallel axis theorem, respectively,

$$I_{ii}=I_{C_{ii}}+m{d}_i^2\approx m{d}_i^2$$

(10)

$$I_{ij}=I_{C_{ij}}-m{d}_i{d}_j\approx -m{d}_i{d}_j,$$

(11)

where m is the point mass and $d_i$ and $d_j$ represent the distance of the point mass to the i and j axis, respectively. A generic view of the aircraft point mass system can be observed in Figure 3.

3.4. Dynamic Model

This work uses the open-source JSBSim v1.2.0 [39] package where the description of forces and moments for an arbitrary rigid body with 6 DOFs can be described in a wide range of formats and schemes, allowing for creativity freedom and accurate reliability in the description of any model, from static contributions to n-variable functions, such as the selected formulation expressed in Equation (8).

As mentioned in Section 2.1, each longitudinal moment and force contribution is described, in the dynamic model, as a set of functions aggregated to each state and control variable. This, coupled with the inertia definition of the aircraft, allows for the generation of the dynamic model of the system. The lateral dynamic behaviour, although modeled into the design, is not used in the simulation beyond maintaining the flight with longitudinal motion only, as this helps match the real and simulation environments.

The dynamic model solves the non-linear and differential equations of motion presented in Section 2.1.

As illustrated in Figure 4, from a current time $t_i$ and given state, each force and moment contribution is calculated and fed into the system of equations, then solved to reach a new state at time step $t_{i+1}$. This cycle is repeated until the desired final simulation is reached. This dynamic simulator assumes, by design, a user input initial condition that allows for the inclusion of the $t_0$ state of the system.

To solve the equations of motion, JSBSim [39] uses two methods: the second-order Runge–Kutta (RK) method to solve the kinematic Equation (3) and the second-order Adams–Bashforth–Moulton (ABM) method to solve both the translation and rotation rate Equations (4) and (5).

The time step ${\Delta }_t$ in the simulations presented in this work was $0.0083$ s, corresponding to a frequency of simulation data acquisition and processing of 120 Hz. This is believed to capture sufficient detail regarding the dynamic behaviour of the UAV.

4. Flight Data and Simulation

The flight data obtained and the simulation environment used to gather data for the optimisation conducted are of the utmost importance for the success or failure of the procedures shown here. This work uses routine structural test flight data with large longitudinal motion excitation to allow for capturing of the longitudinal characteristics of the system. The test starts from a trimmed flight condition, with subsequent pull-up and pull-down manoeuvres being used to evaluate the structural reliability of the aircraft.

An example test is shown in Figure 5, where six plots are presented: pitch $\theta$, pitch rate q, estimated angle of attack $\alpha$, elevator deflection delta from the initial trim state $\Delta {\delta }_e$, calibrated airspeed $V_{CAS}$, and altitude h.

4.1. Aircraft Data Gathering

The onboard systems of the aircraft have a data collection module to determine the angular position, rates, accelerations, and speed, as well as its position. The data gathering system also includes a Kalman filter to better estimate the state variables of the aircraft by reducing both the noise and uncertainty coming from affected readings. For the longitudinal motion of the aircraft, the most important observable data are the longitudinal velocity, pitch angle and both pitch rate and acceleration. These are retrieved using gyroscopic sensors and calibrated Pitot tubes. The control data are provided by the outputs given from the auto-pilot system and are linked to both the throttle and elevator.

In terms of data rates and frequencies, the auto-pilot system works at 400 Hz; however, the control actuation is limited to 50 Hz. For this reason, the control output is only evaluated at the control actuation frequency. Since the simulated environment acts at 120 Hz, there is a sufficient frequency margin to capture the control input commands provided to the real aircraft system. A linear interpolation is used to have a continuous function linked to the provided control inputs.

It is important to note that, although real-state flight data estimation is prone to error, the observable data are considered to be exact. Thus, the error observed between the simulation and real aircraft is assumed as is. The analysis of the degree of uncertainty in the sensor acquisition data or their post-processing is out of the scope of this work.

Nonetheless, having proven the method, including new techniques for a more robust and accurate data gathering system would only benefit the final system identification scheme.

4.2. Simulation Environment

The purpose of the simulation environment is to model the longitudinal dynamic behaviour of the aircraft by mimicking the same control inputs provided to the system. The simulation environment set for this work is based on two phases: initial conditions and control commands.

For the first phase, the simulation assumes that the aircraft starts from a trimmed leveled flight condition. With the help of a simulated auto-pilot system, the modeled UAV is instantiated in the simulation environment and its state is converged towards the one found in leveled flight conditions by the aircraft based on the equivalent log data. The convergence criteria are based on the evaluation of the steadiness of carefully chosen variables: calibrated airspeed $V_{CAS}$, vertical velocity $\dot{h}$, pitch rate q, roll rate p, and yaw rate r. This phase sets the reference trimmed-flight state and control conditions as the ones found in the beginning of the real manoeuvre.

As will be demonstrated, the real flight manoeuvre should have a portion of trimmed flight to ensure high-quality conditioning of the simulation environment.

The second phase is simply an injection of the control inputs into the simulated model, allowing the equations of motion to drive the state of the system. The simulated auto-pilot is turned off for the longitudinal motion and the longitudinal controls are replaced by mimicking the real flight data. To ensure equal dynamic excitation on both the real and simulated models, and to avoid trimmed flight offset errors, delta ($\Delta$) control inputs are used instead of absolute control data. Thus, the elevator deflection command fed into the simulated environment is not the actual elevator deflection but the variation from the reference trimmed-flight condition. A comprehensive overview of the simulation environment used in this work is shown in Figure 6.

The simulation environment is a consequence of the dynamic equations solved through time. The accuracy of its results is highly dependent on the time step used and on the physical description of the model. The simulation environment is capable of predicting complex non-linear phenomena as long as these non-linear acting forces and moments are described in the model. In this work, simple linear or polynomial functions are used to determine the dynamic behaviour of the system, which are valid for small regimes of attitude of the aircraft where the behaviour is assumed to be linear.

For higher attitude amplitudes, research on this non-linear behaviour modelling can be studied, and the same simulation environment can be used, as it is a robust non-linear dynamic solver.

5. Parameter Estimation Formulation

Before addressing the optimisation and system identification scheme, the construction and usage of the parameter estimation formulation is explored. First, the problem is formulated along with the error function and the algorithm is validated using a sample scenario. Then, the design variables chosen for the optimisation are studied.

Both single-manoeuvre and multi-manoeuvre algorithms are hereby presented.

5.1. Problem Formulation

The purpose is to identify the error between simulated and real flight manoeuvre data and provide an algorithm that minimises this error by determining the tuned aerodynamic coefficient values of the simulated aircraft. To achieve this, a python-based framework was built to handle both the simulation model runs and all data analysis and processing between the simulation model and test flight data. The framework has the added benefit that the real and simulated models can be chosen by the user, meaning two simulated aircraft can be used to allow the convergence of one to another. Additionally, different models can be implemented without requiring the creation of a new framework to handle the same algorithm procedure.

As stated, the formulation is an error-output-based algorithm, where the error between observable data of both simulated and real aircraft is measured and used. The error is evaluated in regard to four different variables: calibrated airspeed ($V_{CAS}$), altitude (h), pitch angle ($\theta$) and pitch rate (q). Through the usage of maximum likelihood methods, the similarity between the real and modeled flight is evaluated.

To calculate the error between observable and simulated data, different approaches can be conducted based on the type of error calculation assumed. These can be categorised in different terminologies. The most popular are the mean absolute error (MAE), mean square error (MSE), root mean squared error (RMSE), and mean absolute percentage error (MAPE) [40]. The choice of error method is dependent on the problem at hand and its inherited properties [40]. Thus, due to simplicity in the error calculation, reliability, and the usage of different amplitude error variables, in this work, a normalised version of the root mean squared error (RMSE) is used.

For the single-manoeuvre case, all errors are evaluated and aggregated as

$$E_s=\sum _i{\displaystyle \frac{RMS{E}_{y_i}}{RMS{E}_{y_{i0}}}}=\sum _i{\displaystyle \frac{\sqrt[]{{\sum }_{t=1}^M{\left({\delta }_{y_{i_t}}\right)}^2}}{\sqrt[]{{\sum }_{t=1}^M{\left({\delta }_{y_{i{0}_t}}\right)}^2}}},$$

(12)

where the index i represents the variable, M is the number of temporal data points in the manoeuvre, and ${\delta }_{y_{i_t}}$ is the error between the test flight and simulated model values for a given variable $y_i$, where the subscript ₀ represents the initial error. This work uses two different error formulations: one formulation uses the absolute value of the error variable,

$${\delta }_{y_{i_t}}=y_{i_t}-{\widehat{y}}_{i_t},$$

(13)

where y refers to the flight data and $\widehat{y}$ is the simulated model data; the other uses a step difference between each state,

$${\delta }_{y_{i_t}}=(y_{i_{t+1}}-y_{i_t})-({\widehat{y}}_{i_{t+1}}-{\widehat{y}}_{i_t}),$$

(14)

where i and t represent the variable index and the time step index, respectively. The first formulation evaluates the absolute similarity between dynamic behaviours and the latter evaluates the similarity between state-variable curve shapes. Both formulations can be used for any given longitudinal error variable and the framework allows the choice to be made by the user.

For the multi-manoeuvre case, the error is the normalised sum of all single-manoeuvre errors,

$$E_m=\sum _j{\displaystyle \frac{E_{s_j}}{E_{s_{j0}}}},$$

(15)

where the index j represents the manoeuvre. Note that the error-output variables chosen ($y_i$) are user-based and can range from one of the aforementioned variables to the full set.

All design variables can be seen in

Table 1. Design variables considered for the optimisation problem.

Variable	Definition
$C_{m_{\alpha }}$	Derivative of pitching moment coefficient with respect to angle of attack.
$C_{m_q}$	Derivative of pitching moment coefficient with respect to pitch rate.
$C_{m_{\dot{\alpha }}}$	Derivative of pitching moment coefficient with respect to angle of attack rate.
$C_{m_{{\delta }_e}}$	Derivative of pitching moment coefficient with respect to elevator deflection. Described as linear or quadratic function. Positive and negative deflections have different terms.
$C_{L_{{\delta }_e}}$	Derivative of lift coefficient with respect to elevator deflection. Described as linear or quadratic function. Positive and negative deflections have different terms.
$C_{L_{\alpha }}$	Derivative of lift coefficient with respect to angle of attack.
$C_{L_0}$	Lift coefficient at zero angle of attack.
$C_{L_{\dot{\alpha }}}$	Derivative of lift coefficient with respect to angle of attack rate.
$C_{L_q}$	Derivative of lift coefficient with respect to pitch rate.

, but some were not used in the optimisation procedure (see Section 5.3). In addition, two variables are neglected from the analysis, namely the pitching moment coefficient at zero lift ($C_{m_0}$) and all drag contributions.

The former is due to the simulation being handled via delta control inputs and not absolute values. The latter is due to the deeply connected drag-propulsion behaviour, where changing the drag behaviour of the aircraft would induce, by consequence, an invalid propulsion system behaviour.

To validate the multi-manoeuvre optimisation scheme, the combined manoeuvre optimisation needs two types of data subsets for the proper validation of the model. Using as reference the surrogate modelling theory [41], from all available manoeuvres, $75\%$ are used as training data and $25\%$ are used as validation data.

The built framework uses the tool openMDAO v3.34.2 [42] to efficiently handle the optimisation process. Among the available optimisers in the pyOptsparse v2.12.0 package [43] used by openMDAO, the Sequential Least Squares Programming (SLSQP) method stands out. Following the works of Lyu et al. [44], the SLSQP method was shown to be the best performing method in the Rosenbrock optimisation problem. For this reason, the SLSQP algorithm was chosen, coupled with a first-order finite-difference method to calculate the gradient vector of the objective function with respect to the design variables.

The optimisation follows a standard optimisation scheme, as shown in Figure 7, where $\tau$ represents the vector of the design variables.

The optimisation problem is then defined as

$$\begin{array}{c}min{\displaystyle \sum _j}{\displaystyle \frac{E_{s_j}}{E_{s_{j0}}}}\\ w.r.t.f_y\end{array}$$

(16)

where $f_y$ are the functions that describe each force and moment coefficients considered (as presented in Equation (8)). The optimisation stopping criteria included a convergence tolerance of ${10}^{-3}$ and a maximum number of iterations of 200.

5.2. Error Formulation Selection and Procedure Validation

To both select the best error formulation and validate the optimisation algorithm, four different combinations of variables were tested with a validation JSBSim model example. This validation was performed using two JSBSim dynamic models: an unmodified initial Tekever AR5 model and a modified version of it. The purpose was to confirm that the optimisation algorithm is successful in matching the dynamic description of the aircraft by closing the response output gap between the unmodified and modified models.

The error study was conducted incrementally, starting with the pitch step error ${\delta }_{\theta }$ and then including altitude h, pitch rate q, and calibrated airspeed $V_{CAS}$ errors. The results are compiled in Figure 8, where, for each optimisation, the results for all single variable errors are shown (Figure 8a) and the error for each variable is also shown (Figure 8b). The plotted improvement value is determined, for each variable, as one minus the error given by the RMSE ratio from Equation (12). This work opted for the usage of step error formulation for the pitch variable ($\delta \theta$) and absolute formulation for all other variables (altitude h, pitch rate q, and calibrated airspeed $V_{CAS}$). This was shown to better capture the dynamics and more reliably provide a superior system identification.

The results indicate that increasing the variables in the error formulations helps by improving the overall score in all variables whilst having better single score in some cases. Additionally, the similarity between aerodynamic coefficients also increases with the increased number of error source variables. The usage of quadratic terms for the elevator deflection function seems to provide a small influence in the overall longitudinal dynamics. All other variables provide a higher influence on the error minimisation.

As for the validation, Figure 9 provides a graphical overview of the conducted optimisation obtained using the defined complete error combination of pitch, pitch rate, calibrated airspeed, and altitude errors. The optimisation results show a clear convergence of the dynamic behaviour, thus validating the usage of this algorithm and the associated complete error formulation, at least when data are clean and not prone to noise and other sources of error.

5.3. Sensitivity Analysis

To assess the impact of the considered design variables on the formulated optimisation problem and the created system identification procedure, a sensitivity analysis was conducted by performing optimisation routines for each single variable error in both absolute and step error formulations. The compiled results are shown in Figure 10, where the accuracy improvements for each design variable are demonstrated using absolute (Figure 10a) and step (Figure 10b) errors.

Results have shown not only the importance of each design variable on the overall error but also the significance of the error method chosen for each optimisation. In both error formulations, the design variables with biggest impact are $C_{m_{\alpha }}$, $C_{m_q}$, $C_{L_{{\delta }_e}}$, $C_{m_{{\delta }_e}}$, $C_{m_{\dot{\alpha }}}$, $C_{L_{\alpha }}$, and $C_{L_0}$. In contrast, $C_{m_{\dot{\alpha }}}$ and $C_{L_q}$ present a minimal impact, so these were removed from the design variable poll.

This study also encompassed the difference obtained by using a linear or a quadratic function for the elevator $f_3\left({\delta }_e\right)$. Using a quadratic formulation is shown to reduce the error between the model and real aircraft, so the quadratic term coefficients for the elevator were used for the remainder of the work.

5.4. Applicability and Usability

Before addressing the successful case study developed in this work, some of the shortcomings of the applied method are revealed and explained. The overall method presented was heavily dependent on two significant factors: (1) data collection and reliability and (2) flight manoeuvre initial conditions.

The former is explained by the assumption that all collected data were assumed to reflect the exact behaviour of the aircraft during the flight manoeuvre. However, this is not the real case, as there are a multitude of sensor and post-processing data errors. A more robust model could have been developed by accounting for these effects, but it would require knowledge of the internal Kalman filter and the modelling of all sensors, thus increasing complexity and reducing computational efficiency.

The latter involves the starting simulation conditions for the developed work. As mentioned, the system identification technique assumes that the aircraft starts from a levelled, trim-state condition from which the variations in control are accounted for and imported into the simulated model. Since this is conducted with unplanned flight data, this behaviour does not occur in every flight manoeuvre, which can affect the data quality and optimisation accuracy.

The aforementioned effects were studied on the developed optimisation scheme. Two distinct aircraft were submitted to the system identification scheme of this work, where non-optimal manoeuvres in the standard provided before were used to study the effectiveness of the method. The overall optimisation results are shown in Figure 11, where the validation data case for both of these aircraft is shown.

There was a clear overall gain in the system identification, with improved stability and similarity between the model and real aircraft. However, as will be shown, these results were not in sync with the best possible outcome of this method. Extensive analysis of the worst results led to the conclusion that the initial trim points of each manoeuvre were not real trim points in real aircraft behaviour, which affected the optimiser interpretation of the real control commands and the overall dynamic behaviour of the system. Nonetheless, the optimiser was still capable of tuning the stability and control derivatives in the aircraft model for improved digital twin response. This shortcoming presents itself as an issue that can be resolved in future work, leading to a more robust and applicable method for complete unplanned flight behaviour.

6. Longitudinal Model Identification

With the problem formulated and the early preliminary model created, an optimisation routine for system identification could be performed in which both single- and multi-manoeuvre cases were studied.

6.1. Single-Manoeuvre Optimisation

Using the presented problem and error formulation, as well as the design variables choice from the sensitivity analysis, the initially built model was corrected using, at first, a single flight test manoeuvre. Figure 12 shows the solution (Figure 12a) and the error improvement (Figure 12b).

Clearly, the optimised model response more accurately follows the real aircraft behaviour, since the optimiser minimised the error between the flight data manoeuvre and the simulated environment outcome, successfully identifying key values of the dynamic model representation. Furthermore, the initial model exhibits divergent behaviour that developed into instability as opposed to what was found in the real flight data (high-pitch oscillations, very low angle of attack, and high altitude loss). Even when a single cycle of the manoeuvre is used in the optimisation, the optimised aircraft model exhibits considerable improved accuracy, within an acceptable margin of error. The most significant improvements were found for the pitch and altitude values, which exhibited the highest initial error.

Although the system identification results were good, single-manoeuvre optimisations are prone to error due to a lack of variability. The optimiser finds dynamic coefficients that minimise the error for one particular manoeuvre but lack the broadness of acquiring general aircraft behaviour dynamics, being heavily single-targeted. To address this issue, the solution found was tested against two other similar flight manoeuvres of the same aircraft. Figure 13 yield the results for the manoeuvres tested (Figure 13a,b), where the initial solution results are also presented.

As expected, the simulated dynamic response of the optimised aircraft model is considerably closer to the real flight than the initial model, with the latter quickly diverging from the real aircraft dynamics wand the former being capable of maintaining an equivalent behaviour. In both manoeuvres, the error is minimised, but the difference is more prominent in the first manoeuvre shown, where the altitude h and calibrated airspeed $V_{CAS}$ responses still significantly diverge from the real response. Nonetheless, the single-manoeuvre optimisation study presented significant improvements in terms of minimising the difference between simulated environment and real flight data. Although prone to single-target optimisation errors, the optimiser was able to capture some of the real aircraft dynamic coefficients, reducing the gap between modelled and real systems, thus effectively generating a first approach to the system identification scheme that is presented here.

To gather some statistics, the single-manoeuvre optimisation case was tested with eight additional unique similar manoeuvres, from which the mean and standard deviation of the error were calculated. Figure 14 yields these results for the observable error data: pitch Figure 14a), pitch rate (Figure 14b), altitude (Figure 14c), and calibrated airspeed (Figure 14d).

The statistical data show that the optimised solution has both lower mean and standard deviation values for all error variables, effectively demonstrating the improved accuracy of the new tuned model. Interestingly, some of the results proved to be better against the manoeuvre used in the optimisation, namely comparing manoeuvre #3 with manoeuvre #1 used in optimisation, which highlights the capability of the proposed system identification procedure in regard to capturing the longitudinal dynamic behaviour of the aircraft.

6.2. Multi-Manoeuvre Optimisation

In order to improve the target bias of the single-manoeuvre optimisation results and augment the system identification capability of the procedure used, the multi-manoeuvre approach was implemented. To this end, multiple flight test datasets performed by the same specific Tekever AR5 aircraft were collected and equivalent flight patterns were identified. A total of nine manoeuvres were selected, with six (75%) being used for training and the remaining three (25%) being used for validation, creating the arbitrary split.

Using the error formulation in Equation (15) in the optimisation problem (Equation (16)), the same preliminarily assessed initial model was optimised using the training dataset of manoeuvres.

Figure 15 exhibits the aircraft response for two manoeuvres used in the optimisation (#4 in Figure 15a) and #7 in (Figure 15b)), with both the initial and optimised models being presented.

The optimised model is found to have increased dynamic behaviour equivalence and, simultaneously, more stability and capability in regard to maintaining proper flight dynamics. It should be noted that only the first cycle in the pull-up/pull-down manoeuvre was fed to the optimiser by limiting the number of time points M in Equation (12). The remaining dynamics were not used during the system identification procedure and serve only to evaluate the consequent aircraft behaviour. However, sufficient dynamic behaviour was captured to allow the simulated aircraft to follow the flight pattern with multiple cycles with the control input injected. Altitude and pitch divergence still occurs after two cycles for the second manoeuvre shown.

When compared to the single-manoeuvre optimisation results, the increased number of manoeuvres fed into the optimiser provided additional dynamic behaviour data that allowed for a significantly more robust and accurate dynamic response from the simulated model.

The final model obtained was capable of predicting the real aircraft dynamics within the manoeuvres used in the optimisation and during the considered time scope, while it also revealed sufficient stability to follow the real aircraft flight pattern with minor divergence when maintaining the real control inputs.

To further demonstrate the final model and validate the proposed system identification process, the response to the validation manoeuvres, unused in the optimisation, was analyzed. Firstly, the mean and standard deviation of the error of all solutions—initial model, single- and multi-manoeuvre optimised models—for all observable variables are shown in Figure 16, namely pitch (Figure 16a), pitch rate (Figure 16b), altitude (Figure 16c), and calibrated airspeed (Figure 16d). The highlighted manoeuvres, #2, #5 and #6, are the validation manoeuvres, and the remaining manoeuvres were used during the optimisation process.

The data show that there is a decrease in both mean and standard deviation errors in the multi-manoeuvre optimisation procedure, hence an increase in similarity between the model and real aircraft. The only exception to this involves manoeuvre #9, where both error metrics presented are maximum, potentially due to poor flight log data or poor initial point conditioning (see Section 5.4).

Looking at the validation manoeuvres (#2, #5 and #6), the optimised model was able to significantly reduce the error between flight and simulation, thus accurately capturing the real dynamics of the aircraft using only single-cycle portions of flight data. Furthermore, data shows that the error metrics found for the validation manoeuvres are within the ones found for the optimisation manoeuvres. Hence, the system identification procedure was able to capture system specific quantities that allowed the construction of a representative model for the longitudinal motion.

For a deeper look, Figure 17 yields the response of each validation manoeuvre using the real, initial, and final model solutions for comparison.

A clear improvement from the initial model to the final optimised model can be seen. The final model solution has minimal divergence in sequential pitching up and down cycles and accurately follows the real aircraft dynamics. The available flight tests clearly provide data that resulted in a substantial improvement of the initial model and assures equivalent dynamic behaviour, especially during the first cycle.

7. Conclusions

This work focused on the development of a system identification procedure for the longitudinal motion of a fixed-wing aircraft. It used a six degrees-of-freedom, non-linear, dynamic model, with forces and moments described by aerodynamic coefficients composed of linear and quadratic polynomial functions of chosen state and control variables to accurately predict real aircraft longitudinal dynamics. Through the usage of gradient-based optimisation, a small set of unplanned manoeuvres with pitching up and down motion were demonstrated to be sufficient for capturing the relevant aircraft dynamic characteristics, significantly reducing the difference between modelled and real aircraft behaviour, such that validation manoeuvres, not used during the optimisation procedure, maintained similar mean and standard deviation errors.

This study has shown that system identification techniques can be used for unplanned or unprepared flight data, where neither specific flight testing campaigns or specific control inputs are needed. Simple known system identification methods based on error-output methodologies, proven to work in planned and guided flight testing campaigns, are shown to be successful in scenarios with a shortage of flight data coupled with unguided or unplanned flight patterns.

The shortcomings of poor simulation conditioning were also discussed, highlighting the importance of correctly choosing flight manoeuvres. Key points for future-work global system identification were also identified.

Interestingly, the single cycle portions of flight data used were shown to be sufficient for modelling the UAV but also for greatly minimising the divergence between model and real aircraft for multiple pitch up and pitch down cycle manoeuvres.

Furthermore, the relevance of using a multi-manoeuvre approach, opposed to a single-manoeuvre approached, was demonstrated. The former allowed for a more robust and broad identification of the system. The multi-manoeuvre approach was shown to be an effective way of mitigating single-target bias error. Additionally, the split between training and validation flight data subsets was proven to be an important step in the system identification process as it not only allowed for validation of the final model but also provided clear data on the effectiveness and applicability of the system identification process.