Next Article in Journal
Deep Learning-Based Infrared Image Segmentation for Aircraft Honeycomb Water Ingress Detection
Previous Article in Journal
Longitudinal Motion System Identification of a Fixed-Wing Unmanned Aerial Vehicle Using Limited Unplanned Flight Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Bayesian Identification of High-Performance Aircraft Aerodynamic Behaviour

by
Muhammad Fawad Mazhar
1,
Syed Manzar Abbas
2,
Muhammad Wasim
1 and
Zeashan Hameed Khan
3,*
1
Department of Aeronautics and Astronautics, Institute of Space Technology, Islamabad 44000, Pakistan
2
Department of Avionics Engineering, Air University Islamabad, Aerospace and Aviation Campus Kamra, Attock City 43570, Pakistan
3
Interdisciplinary Research Center for Intelligent Manufacturing & Robotics (IRC-IMR), King Fahd University of Petroleum and Minerals, Dhahran 31261, Saudi Arabia
*
Author to whom correspondence should be addressed.
Aerospace 2024, 11(12), 960; https://doi.org/10.3390/aerospace11120960
Submission received: 20 July 2024 / Revised: 26 October 2024 / Accepted: 18 November 2024 / Published: 21 November 2024
(This article belongs to the Section Aeronautics)

Abstract

:
In this paper, nonlinear system identification using Bayesian network has been implemented to discover open-loop lateral-directional aerodynamic model parameters of an agile aircraft using a grey box modelling structure. Our novel technique has been demonstrated on simulated flight data from an F-16 nonlinear simulation of its Flight Dynamic Model (FDM). A mathematical model has been obtained using time series analysis of a Box–Jenkins (BJ) model structure, and parameter refinement has been performed using Bayesian mechanics. The aircraft nonlinear Flight Dynamic Model is adequately excited with doublet inputs, as per the dictates of its natural frequency, in accordance with non-parametric modelling (Finite Impulse Response) estimates. Time histories of optimized doublet inputs in the form of aileron and rudder deflections, and outputs in the form of roll and yaw rates are recorded. Dataset is pre-processed by implementing de-trending, smoothing, and filtering techniques. Blend of System Identification time-domain grey box modelling structures to include Output Error (OE) and Box–Jenkins (BJ) Models are stage-wise implemented in multiple flight conditions under varied stochastic models. Furthermore, a reduced order parsimonious model is obtained using Akaike information Criteria (AIC). Parameter error minimization activity is conducted using the Levenberg–Marquardt (L-M) Algorithm, and parameter refinement is performed using the Bayesian Algorithm due to its natural connection with grey box modelling. Comparative analysis of different nonlinear estimators is performed to obtain best estimates for the lateral–directional aerodynamic model of supersonic aircraft. Model Quality Assessment is conducted through statistical techniques namely: Residual Analysis, Best Fit Percentage, Fit Percentage Error, Mean Squared Error, and Model order. Results have shown promising one-step model predictions with an accuracy of 96.25%. Being a sequel to our previous research work for postulating longitudinal aerodynamic model of supersonic aircraft, this work completes the overall aerodynamic model, further leading towards insight to its flight control laws and subsequent simulator design.

1. Introduction

Aerodynamic coefficients of aircrafts are formulated using Computational Fluid Dynamics (CFD) and Wind Tunnel testing, which require a large amount of empirical data and costly experimental setup for its numerical analysis. Alternatively, system identification (SI) depicts a mathematical model of a dynamical system using an observational dataset of a system using a parametric modelling technique categorized as white, grey, and black box modelling techniques. Grey box modelling uses observation data along with a priori knowledge of the system to model its dynamics. Bayesian estimation being closely linked with grey box modelling, thus, becomes a natural selection for parameter identification [1]. OE and BJ model structures have been adopted in this paper for obtaining the lateral–directional mathematical form of an aerodynamic model under varying flight and noise conditions. An aerodynamic model of a supersonic aircraft must capture practical dynamics related to its nonlinear behaviour as well as have sufficient complexity to accurately predict its dynamic response over a given flight regime. Figure 1 presents an orientation of this research, which mainly encompasses Model postulation, Reduced Order Modelling, and Parameter optimization.
Aircraft parameter estimation is an optimization problem [1]; therefore, a technique proposed in this research is the optimization of parameters using Bayesian framework using a priori knowledge of aircraft constraints and parameters obtained from OE and BJ polynomial models. A reduced order model is obtained using AIC, and error minimization is acquired through L-M algorithm. The Bayesian framework is better than other classical estimators as it gives interval estimates instead of point estimates calculated by Maximum Likelihood and Least Square estimators, which is much simplistic as compared to Bayesian estimates [1].
To start with the most important, sensitive, and tedious part of SI process; input design was tackled. An optimized input as per dictates of non-parametric (FIR) modelling optimally excites the aircraft dynamic modes [2]. Subsequently, MATLAB (https://matlab.mathsworks.com, accessed on 20 March 2024) Simulink-based Flight Dynamic Model (FDM) of F-16 was developed as shown in Figure 2. Aircraft trim points were obtained under two flight conditions, i.e., Translational Steady flight and coordinated pull up flight [3]. After attaining a stable response from FDM, a time-domain Multiple Input Multiple Output (MIMO) database was established based upon set of aileron (δa)-rudder inputs (δr) and roll (P)-yaw rate (R) outputs. Dataset was smoothed using data pre-processing techniques and further segregated for testing and validation. Initial parameters have been obtained using the Output Error (OE) model and a further BJ model was postulated using initial values of the OE model. Reduced order modelling was performed using AIC. Parameter error minimization was performed using L-M algorithm and optimization of parameters using the Bayesian technique was performed.
Overall gains accrued from this research includes development and obtaining deep insight into nonlinear flight control laws of an agile aircraft, understanding pilot–aircraft interactions in complex lateral–directional maneuvers, knowing flying qualities, and carrying out future modifications for a flight training simulator.

1.1. Related Works

The roots of SI goes back to the eighteenth and nineteenth century by Gauss and Bayes [1], since then, different methodologies have been explored, which mainly revolved around two main streams, i.e., model structure determination and parameters’ estimation. Herbert et al. [4] explained Bayesian identification with grey box modelling to estimate model parameters and its explicit advantages over known physical modelling and black box modelling. Prolific research on Aircraft SI [5] has been performed by the NASA Langley Research Centre (NaLC) for testing and application of state-of-the-art SI techniques on different aerial vehicles; however, SI of supersonic aircrafts still needs further exploration with respect to reduced order modelling and parameter estimation. SI on the F-16 aircraft nonlinear model was performed by Morelli et al. [6], who presented a comparison of EE and OE models. Dynamic Mode Decomposition (DMDC) was performed by Oznurlu et al. [7] in which aerodynamic model of fighter aircraft was discovered. The estimated mean errors of 0.25 and 0.17 were achieved using the DMDC technique. Standard grey box modelling structures of OE, ARX, and ARMAX were implemented by E. Belge et al. [8], on UAV, using Least Square estimation. The BJ method, being more stable than the OE method, was implemented by Omar et al. [9], for SI of quadcopters; however, a more parsimonious model would have been obtained through model order reduction algorithms. Pillonetto et al. [10] in his research of SI, explained regularization techniques to optimize model order and achieved an accuracy of more than 82% model fit. A M Kwad et al. [11] analyzed online and offline techniques of SI using L-M algorithms and obtained good fit results coupled with fast computations. Hierarchical least squares and a multi-innovation least squares algorithms were proposed by Yan Ji et al. [12] for parameter estimation. Millidere [13] performed reduced order modelling in SI using Binary Particle Swarm Optimization (BPSO) method.
Eugene A. Morelli [14] analyzed the Equation Error (EE) method for parameter estimation on simulated data of F-16 aircraft and compared it with OE method. Mathias et al. [15] introduced a wavelet transform technique with OE Method using a compact time–frequency representation for resolving SI problem. This enabled multi-axis inputs for complex maneuvers with fewer parameters. A joint input–output SI technique was implemented by Lopez et al. [16] for parameter estimation of highly redundant control surfaces of the Bell V-280 aircraft. The OE Method with Vortex Lattice technique was presented by Simmons et al. [17] for aircraft SI and improved the model. The OE Method being most popular technique to resolve an SI problem was also implemented by Dimas Abreu Archanjo Dutra [18], who discussed its limitations on aircraft SI from flight-test data using collocation formulation technique, which resulted in better convergence and robustness of the model. It is important to highlight the research contribution by Mukhopadhaya et al. [19] utilizing Gaussian Process Regression, close to Bayesian estimation, for improvement of the uncertainty of aerodynamic database. The technique combines data from multiple sources with variable uncertainty levels, and then quantify uncertainty of one dataset. However, this is a non-parametric regression technique, which comes under the data-driven SI domain. Bayesian estimation using parametric modelling utilizing grey structure needs exploration.
A substantial research contribution [20,21,22,23,24] has been made for parameter estimation and its optimization using varying nonlinear optimization algorithms. A comparative analysis of traditional SI model structures on Unmanned Aerial Vehicles has been conducted by Fatima et al. [25]. Best fit results were obtained using the ARMAX model. The Single Decomposition Method with ARX structure was applied by Bagherzadeh [26] for SI of elastic aircrafts, who further used the L-M algorithm for parameter optimization. A reduced order modelling technique on the SI of a simulated F-16 fighter jet was performed by Michele et al. [27] by employing Particle–Bernstein polynomials with Principle Component Analysis (PCA). Ram et al. [28] dealt with SI in Bayesian framework to handle parameter and structural model uncertainties. Despite noteworthy research contribution on SI of aircrafts using different techniques, a wide range of research exploration is still needed in discovering aerodynamic model of supersonic aircrafts with a model configuration that amply covers both deterministic and stochastic parts. Limited use of Bayesian estimation in aircraft SI has been clearly explained in [1,2] due to the requirement of knowing prior probability distribution function (p.d.f); therefore, this research gap needs exploration to investigate pros and cons of identifying parameters of supersonic aircraft using this unique probabilistic framework. Refs [29,30,31] also presented different accounts of SI based structures.
A summary of the most recent research work using Bayesian approaches for aircraft design and modelling has been presented in Table 1, comparing them with Bayesian framework adopted in this paper. It was concluded that our research work uses the Bayesian approach with a unique perspective and application as compared to earlier work conducted.

1.2. Research Motivation

Discovering a mathematical model for the highly intricate, complex, and unstable dynamics of supersonic aircrafts is the most challenging and time-consuming task during the design and simulation phase. SI methods made it easy, through wide-ranging techniques used by white, grey, and black box modelling. SI opened the gateway towards efficient investigation of the dynamic behaviour of complicated systems. All mathematical models are built on observational data collected from the FDM of aircraft and subsequently its mathematical form is postulated. However, the tough part is to design a model that is simple, compact, yet precise. Therefore, the major motivation in this research is to postulate a grey box lateral–directional mathematical model, which is parsimonious, reliable, and robust and helps to design a controller with minimum computational cost. As per the recent research work mentioned above, methods adopted in this paper have been rarely applied on the SI of an agile aircraft, to contribute to the research work so far conducted in this field.

1.3. Research Contribution

The research presents the novelty of the technique as well as the novelty of application. The novelty of the technique by implementing SI embedded with Bayesian theory and the novelty of application stem from its implementation on a high-speed, agile, and high-performance aircraft. Further research has been focused on both contours of SI, i.e., Model structure determination and parameter estimation using time-domain grey modelling techniques. In comparison to previous research work [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,37], this paper emphasized determination of the open-loop lateral–directional aerodynamics of supersonic aircraft using the BJ model with optimization through the L-M technique. For this, simulated F-16 FDM was developed, and appropriate inputs were designed, and against them time-based outputs were recorded. The BJ structure is selected, being the most appropriate model structure in time-domain grey box modelling for complex plants. Model order reduction was performed using AIC. Parameters of reduced order BJ structure were refined using L-M optimization and Bayesian estimation. Parameter estimation accuracy of 96.25% has been achieved. In this way, a parsimonious and accurate model has been designed, opening the potential and gateway to propose a suitable adaptive controller for the aircraft lateral–directional dynamics.
This paper is divided into five sections. Section 2 gives a detailed account on mathematical foundations and aerodynamics parameters of the aircraft modelling. Section 3 gives description of methodologies containing algorithms and techniques applied to discover the lateral–directional dynamical model of supersonic aircraft. Section 4 presents results and discussion. Lastly, conclusion has been presented to give wholesome review of the paper, leading towards possible future research work.

2. Aircraft Modelling

2.1. Mathematical Formulations

6DOF nonlinear ODEs derived to demonstrate aircraft dynamics [3,37,38] were, flat earth, rigid-body mechanics, with constant mass and thrust force along its longitudinal-axis shown in Figure 3. Subscripts ‘B’, ‘s’, and ‘w’ written with force components x, y, z represent body, stability, and wind axes, respectively. Equations (1)–(10) lay the foundation of our research, which are used to develop the ‘Flight Dynamic Model’ of the aircraft in MATLAB/Simulink as shown in Figure 2. Once the Flight Dynamic Model is ready, then simulated inputs are given to record input–output time histories. This input–output dataset is utilized for estimation of Aerodynamic behaviour of the aircraft using grey box modelling and Bayesian estimation. Angle of attack, slide Slip Angle, velocity, Forces, kinematics and Moment equations have been implemented in MATLAB/Simulink to develop Flight Dynamic Model of F-16. Aerodynamic database and equations have been acquired from [3,38].
α and β are angle of attack and sideslip angles, respectively. V t is the velocity vector. u, v, and w are velocity components in the x, y, and z axis.
α = t a n 1 w u β = s i n 1 v u V t = u 2 + v 2 + w 2
The full derivation has not been presented due to brevity purpose. Forces, Kinematics, and Moments are elaborated one by one in Equations (2)–(10).
Translation forces can be presented as:
U ˙ = R V Q W g D sin θ + X A + X T / m
V ˙ = R U + P W + g D sin ϕ cos θ + Y A + Y T / m
W ˙ = Q U P V + g D cos ϕ cos θ + Z A + Z T / m
where U ˙ , V ˙ , and W ˙ are the components of linear velocities along the three body axes. P, Q, and R are the angular velocities (roll, pitch, and yaw rates) along body axes, and m is mass. ( X A + X T ), Y A + Y T , and ( Z A + Z T ) are resultant forces in the x, y, and z axes. g D is the gravitational constant.
The basic relationship for rotational kinematics in an inertial frame of reference is presented by ϕ ˙ , θ ˙ , and ψ ˙ , which are euler angle rates:
ϕ ˙ = P + tan θ ( Q sin ϕ + R cos ϕ )
θ ˙ = Q cos ϕ R sin ϕ
ψ ˙ = ( Q sin ϕ + R cos ϕ ) / cos θ
The implementation of rotational moment equations when applied in the body frame of reference is transformed as:
P ˙ = J x z J x J y + J z Γ P Q J z J z J y + J x z 2 Γ Q R + J x z n Γ + J z l Γ
Q ˙ = J z J x J y P R J x z J y P 2 R 2 + m J y
R ˙ = J x J y J x + J x z 2 Γ P Q J x z J x J y + J z Γ Q R + J x z l Γ + J x n Γ
Γ = J x J z J x z 2
where Jx, Jy, and Jz are the moments of inertia of the aircraft and l, m, and n are roll, pitch, and yaw moments in x, y, z directions, respectively. J x z is the product of inertia in x and z axis, and Γ is the inertia matrix component.
Non-dimensional Force and moment coefficients of the aircraft were calculated using the following equations [2]:
For the force coefficients X -axis, Y axis, and Z axis:
C x = 1 q ¯ S ( m a x T )
C Y = m a y q ¯ S
C Z = m a z q ¯ S
where a x , a y , and a z are the acceleration components in the x, y, and z axes, and T is the engine thrust along x axis.
Combining Equations (5a) and (5c), the coefficient of Lift ( C L ) and coefficient of drag ( C D ) are presented as:
C D = C X cos α + C Z sin α
C L = C Z cos α + C X sin α
C L = 1 q ¯ S b J x p ˙ J x z p q + r ˙ + J z J y q r
where terms used in Equation (7a) have already been explained in Equation (4).
Non-dimensional moment coefficients (roll, pitch, and yaw) of the aircraft were calculated using following equations [2]:
C l = 1 q ¯ S c ¯ I y q ˙ + ( I x I z ) p r + I x z p 2 r 2 I p p r
C m = 1 q ¯ S c ¯ I y q ˙ + ( I x I z ) p r + I x z p 2 r 2 I p p r
C n = 1 q ¯ S B I z r ˙ I x z p ˙ q r + I y I x p q + I p p q
where C l , C m , C n are roll, pitch, and yaw coefficients of moments, respectively. I p and p are the inertia of the rotating mass and angular velocity, respectively.
Equations (1)–(10) lay mathematical formulation for development of FDM and subsequently plant model for observational data recording for SI of the lateral–directional aerodynamic model.

2.2. Development of FDM

In accordance with [3,38], 6DOF simulation environment has been developed using scripted MATLAB (2021b) and Simulink (Figure 2) and implementing nonlinear Equations (1)–(10). Aerodynamic data obtained from NASA-Langley wind tunnel tests [38] on a subscale model of F-16 airplane, having the trimmed speed range (between 400 ft/s to 800 ft/s), angle of attack (−10 to 45 degrees), and of sideslip angle (−30 to 30 degrees). Within this speed range, variations in leading edge flap (lef) and speed brake (sb) factors have been eliminated. In this research, evaluation of the aerodynamic behaviour of F-16 has been conducted at 502 ft/s, which is the trimmed flight speed at nominal sea level flight conditions [3].

2.3. Analysis of Aircraft Lateral-Directional Model

In order to proceed with SI of the aircraft, one of the mandatory steps is to find the frequency response of the aircraft for designing optimized inputs [2]. Following steps (1 and 2) have been implemented to achieve the frequency spectrum of the aircraft.

2.3.1. Step 1: Numerical Linearization to Obtain Trim Points

A trim point is a point where all state derivatives of the system become zero. Aircraft nonlinear state equations were linearized to two flight conditions, steady translational flight, and coordinated turn flight. Equations (2)–(10) were linearized under constraints of wings level, non-slide slipping, and steady state flight. Numerical techniques [3] were implemented to determine trim points for each flight condition by changing the values of independent variables and to achieve a cost function value, which is based upon state derivatives computed from aircraft FDM, subsequently minimizing the cost function through a minimization algorithm. The numerical solution obtained, thus, provided the required state and control information for linearization. Algorithm 1 formulated and achieved a complete trim model of the aircraft. There are three categories of variables; input and state vectors are used to find the state derivatives and then a scalar function is determined, which is optimized to obtain a trim point information.
Algorithm 1: Aircraft Steady State Points
Input:
1. Specify inputs and states ( x T , z T ,   γ ); inputs x T = α , β , δ t , δ e , δ a , δ r ;
     States z T = [ V t , h ,   γ , ϕ ˙ , θ ˙ , ψ ,   ˙ C G ]
Iterations:
2. Set tolerance value = 1 × 10−8
  2.1. Compute z ˙ (state derivatives) from inputs x T and states z T .
  2.2. Compute cost function: J = a1∗ z ˙ 1 ( 2 ) , a2∗ z ˙ 2 ( 2 ) an z ˙ n ( 2 ) .
  2.3. Apply minimization algorithm (Nelder Mead Algorithm) on cost function.
  2.4. Stopping criteria: Tolerance value achieved—terminate the iterations.
Output:
3. Display Trim data: Vt, α ,   β ,   δ t ,   δ e ,   δ a ,   δ r ,   δ e ,   P ,   Q ,   R   for   flight   conditions   Y T .

2.3.2. Step 2: Non-Parametric (FIR) Modelling

Non-parametric modelling is mandatory, as it is aimed to define the input, which excites dynamical modes of the plant (aircraft). The FIR model determines the natural frequency, bandwidth, and discrete time input–output delays (nk) of the aircraft. In addition to above elements, insights to aircraft design characteristics, such as gains and time constant are also obtained. The FIR model is obtained as per Equations (11a) and (11b):
y k = l = 0 M g l u k l
G e j ω = n = 0     g n e j ω n
where g is the coefficient of impulse, k is sampling instant for discrete-time, and n for continuous time. l denotes the delay. Frequency parameters obtained, such as highest frequency, sampling rate, bandwidth, time delays, and time constant, lay a foundation for specialized input design for aircraft lateral–directional SI. MATLAB ‘impulsest’ has been implemented to achieve this step.

3. Methodology

In accordance with the research framework formulated in Figure 1, proposed methodologies in this paper along with their algorithms are elaborated in this section.

3.1. Input Design

Inputs comprise two main categories [40], i.e., heuristic and non-heuristic. Non-heuristic inputs (doublet, 3-2-1-1, 2-1-1) have been used in this research, as they pertain to time-domain modelling of aircraft dynamics, and a priori aircraft information in terms of frequency responses has also been obtained [2]. An iterative process was followed to design an optimal multiple input design such that changes in P and R remain in the linear range as per Equations (12) and (13) and flow chart Figure 4.
C r 2 = m a x t U 2 ( t ) 1 / N t = 1     N U 2 ( t )
Pulse   Width = 0.7 2 f n
where Cr is the crest factor and u(t) is the elevator input ( δ e ) over time period t from 1 to N. The peak values (numerator) and Root Mean Square values (denominator) are used to calculate crest factor of the doublet inputs, which gives us a refinement of the peak of doublet. A smaller crest factor defines a good wave. The pulse width of doublet input is determined using Equation (13), where f n represents the natural frequency of the aircraft.
For lateral–directional dynamics of a high-performance aircraft, doublet input 2-1-1 was considered more suitable than the 3-2-1-1 input [2], as the 3 pulse is quite long that it tends to drive the aircraft away from trim flight conditions. Supersonic aircrafts are open-loop unstable; therefore, the 3-2-1-1 phenomenon becomes more pronounced, as it tends to cross the limitations in terms of aircraft natural frequency as well as aerodynamic and structural constraints. Subsequently Algorithm 2 presents model postulation using BJ Structure.
Algorithm 2: Determine BJ Structure
Input:
  • Prepare Dataset: Aileron ( δ a ) and rudder ( δ r ) inputs as u(t) over time and record output (P) roll and yaw (R) rates.
  • Estimate input–output delay through the FIR Model.
  • Pre-processing data: Detrend, fill missing data, removing drifts and outliers
  • Segregate training and validation dataset

Iterations:
5.
Initial estimate of polynomial orders nb, nc, nd, and nf. Start with a guess of [1, 1, 1, 1]
6.
BJ Structure: B ( q 1 ) F q 1   u k + C ( q 1 ) D q 1   e k
7.
Stopping criteria: Manually increase model orders until max best fit % is achieved.

Output:
8.
Model order, polynomials, coefficients, best fit %, FPE, and MSE

3.2. Model Postulation—OE and BJ Structures

In grey box modelling, the OE model forms the basis for obtaining initial parameters of the model. However, the model lacks accuracy due to the non-availability of the stochastic component, which is provided by the BJ model structure. Built in function of MATLAB have been used for this purpose. Inputs as per Figure 4 were given to the FDM (Figure 2), and outputs were recorded. The initial parameters obtained from the deterministic part of the OE model were subsequently used in constructing a BJ polynomial model. The BJ model, being a superset structure to model the plant and noise characteristic, gives an edge over other structures (OE, ARX, ARMAX), which exhibit limitations in the modelling plant and noise dynamics, distinctly, as well as lack accuracy if the system exhibits nonlinearity [4,9]. For input u(k), output y(k), process noise e(k), and measurement noise v k , basic OE and BJ models are given by Equations (14) and (15), respectively.
y k = B ( q 1 ) F q 1   u k + e k
y k = B ( q 1 ) F q 1   u k + C ( q 1 ) D q 1   e k
where B and F are deterministic polynomials and C and D are stochastic polynomials.

3.3. Reduced Order Model—(AIC)

Reduced order modelling is performed by following a rule of parsimony, which states that given two or more models with same accuracy, the model with lesser terms is selected as a suitable model [2,4]. An algorithm with trade-offs between model complexities, model fitness, loss function value, and information theoretic metrics is selected for attaining a reduced order model. For aircraft SI, AIC has been selected, which does not depend upon sample size of a dataset or distribution of errors and only caters for likelihood function and number of parameters of the model. AIC metrics use Kullback–Leibler’s Information (KLI) theoretic measure [41,42], which examines distance between probability densities of two functions using Equations (16) and (17):
I   ( f ,   g ) = f ( x )   log   [ f   ( x ) / g ( x ) ]   d x
AIC j = 2 l j θ M L E ( j ) 2 k j
Equation (16) is known as the KLI equation, where f and g are the p.d.f of two different approximations. Equation (17) proposes the use of expected KLI distance and establishes its connections with log-likelihood. Where θ M L E ( j ) is the log-likelihood function, l j = dim ( θ ), and the second term k is the penalty for number of terms to achieve a parsimonious model. In this paper, after obtaining the basic BJ model, AIC metrics have been applied to BJ polynomials and reduced model has been achieved after losing a nominal accuracy of the model. The AIC scheme was implemented after obtaining a set of BJ models. Iterations were carried out using a MATLAB snippet that compares AIC values of each BJ model and selects the model with lowest AIC value.

3.4. Error Minimization—L-M Algorithm

An L-M algorithm has been implemented in this research work for achieving parameter error minimization of the BJ model structure. As the BJ model uses a nonlinear Least Squares technique, its parameter error is curtailed using the L-M algorithm. L-M Algorithm 3 implements a hybrid approach and utilizes Gauss–Newton and Steepest Descent to converge to an optimal solution. A mathematical presentation [40] is given Equation (18).
H c i = H i + λ I p p
H c i represents a modified hessian matrix with H i as the initial hessian matrix obtained from Newton–Raphson technique. λ denotes a regularization parameter (initially set to 0.001), I is identity matrix, and p is the number of parameters to be optimized. The convergence criteria set for error minimization is equated in Equation (19). J is the cost function for θ ^ k estimated parameters.
J θ ^ k J θ ^ k 1 J θ ^ k 1   < 0.001

3.5. Model Optimization—Bayesian Approach

The only framework that fully articulates and explores prior uncertainties in parameters of a grey box model is the Bayesian approach. While dealing with grey box modelling parameters of aircraft SI, Bayesian probabilistic theory appears to be one of the natural approaches as the user has portion of knowledge of θ as well as understanding of aircraft constraints, boundaries, and flight conditions. In this section, model parameters obtained from the grey box BJ structure are refined using Bayesian estimation. A major difference between Bayesian and other traditional approaches is that it assumes parameters to be random. Precisely, it is not the randomness of the parameter, rather it is the knowledge of the uncertainty that makes it random [1]. The aim is to calculate posterior p.d.f P( θ | Z ) using Bayes’ rule [2]:
P ( θ   | Z )   P Z θ P θ =   P Z θ P θ
P ( θ   | Z )   = P Z θ   P ( θ ) P Z θ   P ( θ )
where
P ( θ ) = ( 2 π ) n p Σ p 1 2 e x p 1 2 θ θ p T Σ p 1 θ θ p
and
P ( Z θ ) = ( 2 π ) N | R | 1 2 e x p 1 2 ( Z H θ ) T R 1 ( z H θ )
θ ˆ = m a x θ   p ( θ z )
Objective function has been calculated as:
J ( θ ) = 1 2 ( z H θ ) T R 1 ( z H θ ) + 1 2 θ θ p T Σ p 1 θ θ p
where P( θ | Z ) represents posterior probability distribution of parameter θ , given observed data Z (input–output dataset). p ( Z θ ) is the likelihood function representing the probability of observing data given that the specific value of θ . θ p represents the prior mean of parameter vector obtained from BJ model. N and n p are dimensions of data Z (input–output dataset) and parameter vector θ p , respectively. R is the covariance matrix of noise term (e(k)), assumed to follow normal distribution with zero mean. Σ p represents the covariance matrix of posterior distribution of θ . H represents observation matrix between θ and data Z .   θ ˆ is th maximum a-posterior estimate, and J is the objective function.
Algorithm 3: L-M Algorithm
Input:
1. Initial parameters obtained from Nonlinear Least square estimation of BJ structure.
2. Set regularization value to 0.001.
Iterations:
3. Set tolerance value.
4. Initial hessian matrix (∇2) using Newton-Raphson Technique: xn + 1 = xn − (xn)/f’(xn)
5. Compute Jacobean matrix
6. Compute Modified Hessian Matrix
7. Compute cost function: J θ ^ k J θ ^ k 1     J θ ^ k 1 < 0.001
8. Update θ ^ vector.
9. Stopping criteria: Tolerance value achieved
Output:
10. Optimized parameters θ ^
Figure 5 describes a simplistic flow chart for implementation of the Bayesian algorithm for parameter refinement of the BJ polynomial model. Equations (20)–(24) have been stage-wise implemented and further expanded in Algorithm 4 to explain the structured way to for Bayesian implementation. An initial guess of parameters has been obtained from the BJ polynomial model. Subsequently, convergence criteria and maximum iterations have been set. Likelihood and prior distribution have been calculated. Posterior distribution has been obtained followed by objective function. Update of parameter θ and the loop continues till convergence criteria or maximum iterations are achieved.
Algorithm 4: Bayesian Estimation
Input:
1. Initial guess of parameters obtained from BJ model θ 0 , , θ k
Iterations:
2. Set Convergence criteria: ϵ = 1 × 10−5
3. Iterate for k = 1 ,   k m a x
4. Compute P Z θ and P ( θ )
5. Compute P ( θ Z )
6. Compute θ ^ = m a x θ   p ( θ z ) as the solution of J ( θ ) = 1 2 ( z H θ ) T R 1 ( z H θ ) + 1 2 θ θ p T Σ p 1 θ θ p
7. Stopping criteria: until k = k m a x or θ ( k ) = θ ( k 1 )
Output:
8. Estimated parameters θ ^

Bayesian Sensitivity Analysis

An essential step to perform sensitivity analysis of Bayesian posterior distribution with respect to its prior estimates of parameters obtained from the Box–Jenkins polynomial model has been implemented in this section. In accordance with [42,43,44], Relative Entropy or Kullback–Leibler (K-L) Divergence between two posterior probabilities density functions P ( θ Z ) and θ p have been determined. Using Equations (20)–(23), K-L Divergence has been further derived in Equations (25) and (26). Minor steps have been excluded to maintain brevity.
D   ( P ( θ Z ) ǁ θ p ) = ( θ 0 , , θ k ) P ( θ Z ) l o g P ( θ Z ) θ p
D   ( P ( θ Z ) ǁ θ p ) = E ( p ) log p ( x ) q ( x )
where P ( θ Z ) is the posterior probability, and θ p is the prior belief log of P ( θ Z ) , and θ p determines the dependence of posterior over prior over each value of θ . E(p) is the expectation value of the probability distribution. The value of divergence (D) specifies how much posterior is dependent on prior or how much divergence has occurred after posterior ( P ( θ Z ) has observed data Z . The larger value of D is desirable as it quantifies that data have more influence over posterior, than the prior.
An additional term has been implemented known as ‘Bayes Factor (BF)’ in accordance with [43] to ascertain that posterior is data-driven and not prior-dependent. BF is the ratio of two hypothesis H 1 and H 2 . In accordance with [43], a higher value of BF depicts that posterior estimates are data-driven, and lower values depict that posterior is dependent on priors. Mathematical formulation of a discrete form is given in Equation (27)
BF   [ H 1 : H 2 ] = P ( Z H 1 ) P ( Z H 2 )
where H 1 is the hypothesis that prior has least influence on posterior, and H 2 represents the hypothesis that prior has more influence on posterior. Table 2 quantifies the BF values as under.

4. Results and Discussion

4.1. Aircraft FDM

F-16 FDM has been used for implementation of set of algorithms containing L-M and Bayesian techniques. Aircraft geometry and nominal flight conditions [3] are shown in Table 3 and Figure 6. In accordance with [3], the FDM of F-16, as shown in Figure 2, was developed in MATLAB Simulink using Equations (1)–(10) to obtain a simulated MIMO input–output dataset of aircraft at different noise levels and trim flight conditions. The main objective of this research is to obtain a discrete time dynamical polynomial model of coupled lateral–directional dynamics; therefore, Aileron Deflection (δa) and Rudder deflection (δr) have been given to the aircraft FDM and output in terms of roll (P) and yaw (R) rates have been recorded for development of the mathematical BJ model.

4.2. Aircraft Trim Conditions

Aircraft trim conditions were obtained by implementing Algorithm 1 for two flight conditions: Steady Straight and Level flight at sea level and Coordinated Turn flight in accordance with [3,38], in which nominal flight conditions have been defined for F-16 simulation examples and subsequent controller design. Table 4 represents F-16 trim flight conditions.

4.3. Optimum Input Design

As discussed in Section 3, an optimum input design strategy has been chosen for this research as it offers number of benefits in terms of flight test efficiency, enhanced excitation in lateral–directional maneuvers and identification of roll and yaw control effects [6]. An iterative procedure as per flow chart (Figure 4) and Equations (11a) and (11b) have been implemented in this research. Non-parametric modelling (FIR) and Bode plots have been implemented to obtain frequency modes of the aircraft in terms of natural frequency, input–output delays, and resonant frequency. Figure 7 and Figure 8 and Table 3 shows the FIR model, Bode plot of the aircraft, and impulse response parameters. Apparently, the input–output delay (nk) seems to be zero; however, to implement a practical aircraft flight scenario, a delay of 1 sample, i.e., nk = 1 has been assumed and a resonant input frequency of 5 Hz has been implemented.
Table 5 shows that maximum resonant frequency is 4.76 Hz, which implies frequency of input signal to be greater than maximum frequency and pulse width to be around 0.1 s to excite lateral–directional dynamic modes of the aircraft.
2-1-1 multiple time-skewed doublet inputs (δa and δr) with amplitude of ±2 degrees, pulse width of 0.1 s, and frequency of 5 Hz has been excited at 1 s and 6 s, respectively. 2-1-1 doublet is considered more suitable than 3-2-1-1 doublet input for supersonic fighter aircraft as the aircraft may leave its steady state/trim condition [2]. Simulated multiple input measurements time histories are shown in Figure 9. This research focuses on MIMO coupled lateral–directional dynamics; therefore, both inputs have been given and output has been analyzed for model postulation. Distinct input–output datasets for testing and validation have been formulated based upon time-skewed-uncorrelated 2-1-1 doublet inputs. Output in the form of Roll and Yaw angles and rates are depicted in Figure 10 and Figure 11. Responses in roll rate (P) and yaw rate (R) till 3 decimals have been rounded off for simplicity purpose.
Figure 10 and Figure 11 demonstrate time histories of responses acquired by giving simulated doublet inputs to the FDM of supersonic aircrafts. It is clearly seen that 2-1-1 doublet inputs at 1 s and 6 s are producing roll and yaw rate responses as well as roll angle changes. However, pitch angle remains the same as inputs are aileron and rudder-based, which merely affects the pitch angle of the aircraft during trimmed maneuver. It can also be seen that the aileron doublet input is producing significant changes to roll rate and rudder doublet input is producing changes to yaw rates. The MIMO lateral–directional dynamics of supersonic aircraft are not decoupled [3]; therefore, aileron doublet input also brings slight changes to yaw rates, and rudder doublet input brings changes to roll rates of the aircraft.

4.4. Model Identification and Parameter Refinement

Optimum input design further leads towards attaining an accurate and parsimonious mathematical grey box model for lateral–directional aerodynamic model of an agile fighter aircraft. In accordance with theoretical concepts of aircraft SI [1,2], the OE model was obtained to construct initial parameter foundation for development of a discrete time BJ Model. A parsimonious BJ model was attained using AIC, and its parameters were optimized using the Bayesian approach. The complete mechanism in the form of flow chart is graphically explained in Figure 12.
Each model has been validated for two flight conditions i.e., Case 1: Steady Straight and Level flight and Case 2: Steady Coordinated Turn Flight Condition. Both Case 1 and Case 2 have been implemented and further validated using “Gaussian noise”. For Case 1, Figure 13 shows the time-domain OE and BJ model structures implemented using Equations (14) and (15), and their gradual improvement with respect to achieving a parsimonious model. Similar simulations were conducted for Case 2, and results are shown in Figure 14.
In accordance with [45], noise patterns to be incorporated in the inputs and outputs are considered to be from stationary random processes with variances σ u 2 and σ y 2 for input and output, respectively. Noise inputs depend upon the type of the aircraft and may be determined using previous experiments as well. Finally, inputs and outputs were corrupted with Gaussian white noise with standard deviation of 0.1 for inputs: aileron deflection (δa) and rudder deflection (δr) and 0.3 for outputs: (P) and (R). Model order reduction has been performed using Equation (17) of AIC. Each result has been presented with its evaluation statistics in terms of its number of coefficients, model order, fit %, FPE and MSE values. A mathematical model has been presented in Equations (27) and (28). A Comparative Analysis matrix has also been drawn in Table 4, Table 5, Table 6 and Table 7, which shows a comparison of postulated model parameters using NLS, MLE, and Bayesian estimators.
Definition 1. Model accuracy criteria is based upon definition presented in [2]: “An estimator is the best linear unbiased estimator if it has minimum MSE among the class of unbiased estimators”.
Case 1: Steady Straight and Level Flight
(P = Q = R = 0, remaining flight parameters as per Table 1).
Figure 13a presents the initial OE model with MSE of 8.1 × 10−3 and accuracy of 95.18%, Figure 13b shows the Reduced Order OE Model with compromise on accuracy (94.09%) using Algorithm 3, Figure 13c shows the initial BJ model, which has used initial parametric values of the OE model, Figure 13d shows the optimized BJ model after incorporating Bayesian estimates with a maximum accuracy of 96.25% and an MSE of 3.9 × 10−3. It was observed that Bayesian estimation improved the model accuracy as well as the MSE value of the model while staying within same model order and number of coefficients. Figure 13e shows residual correlation of optimized BJ Model within 95% confidence interval (blue region), and Figure 13f shows the pdf of each parameter attained during Bayesian estimation process. Figure 13g shows posterior sensitivity analysis (posterior dependence on prior) of Bayesian estimates, where a larger value of 31.3017 shows that Bayesian posterior has diverged from prior after observing the dataset and is not much dependent on prior estimates. The BF values calculated from two hypothesis comes as 348.54, which confirms that posterior is data-driven and not prior-dependent. Equation (28) shows the mathematical model form of discrete time BJ Model of order 5, obtained for lateral dynamics of a supersonic aircraft. B and F represent plant dynamics, and C and D represent noise dynamics. z is the shift operator.
B(z) = −2.479 − 2.942 z−1 + 0.7024 z−2 + 1.285 z−3
C(z) = 1 − 0.9291 z−1
D(z) = 1 − 1.138 z−1 + 0.2772 z−2 + 0.2227 z−3 − 0.4422 z−4
F(z) = 1 + 0.4055 z−1 − 1.211 z−2 − 0.3027 z−3 + 0.3932 z−4 − 0.009715 z−5
Case 2: Steady Coordinated Turn Flight
(P = −0.01, Q = 0.29, R = 0.06, remaining flight parameters as per Table 1).
Figure 14a presents the initial OE model with 83.99% accuracy, Figure 14b shows Reduced Order OE Model with compromise on accuracy (72.96%) using the L-M algorithm, Figure 14c shows initial BJ model, and Figure 14d shows the optimized BJ model after incorporating Bayesian estimates with a maximum accuracy of 80.07% and a minimum MSE of 0.1955. It was concluded that Bayesian estimation improved the model accuracy while staying within same model order and number of coefficients. Figure 14e shows residual correlation of the optimized BJ Model within a 95% confidence interval (blue region). Figure 14f shows the pdf obtained during Bayesian process. Figure 14g shows the posterior sensitivity analysis (posterior dependence on prior) of Bayesian estimates, where larger value of 37.4666 show that the Bayesian posterior has diverged from prior after observing the dataset and is not much dependent on prior estimates. BF values calculated from two hypothesis comes as 299.31, which confirms that posterior is data-driven and not prior-dependent. Equation (29) shows the mathematical model form of grey box BJ Model obtained for lateral dynamics of the aircraft. B and F represent plant dynamics, and C and D represent noise dynamics. ‘z’ is the shift operator. For brevity purpose, similar plots obtained for yaw rates (R) have not been included; however, parameter estimates have been included in Table 6 and Table 7.
B(z) = −2.487 − 0.4589 z−1
C(z) = 1 − 0.9973 z−1 + 0.134 z−2 + 0.1181 z−3 − 1.096 z−4 + 0.8415 z−5
D(z) = 1 − 1.626 z−1 + 0.6488 z−2 + 0.1839 z−3 − 0.3126 z−4 + 0.1056 z−5
F(z) = 1 − 0.6974 z−1 + 0.1451 z−2 + 0.07067 z−3 − 0.2726 z−4 + 0.2265 z−5
Estimated parameter obtained from Bayesian technique for both flight cases i.e., Straight and level flight and Coordinated Turn Flight have been enlisted in Table 6, Table 7, Table 8 and Table 9. Straight and level flight having simplistic flight conditions have lesser number of coefficients as compared to coordinated turn flight. Furthermore, a comparison matrix has been developed in which parameters obtained from Bayesian estimation were compared with Nonlinear Least Square (NLS) estimation and Maximum Likelihood Estimation (MLE). Obviously, Bayesian being linked with grey box modelling due to its a priori probability knowledge of parameters produces better results in terms of MSE than other two techniques. Otherwise, calculation of pdf restricts its use in aircraft SI [2]. The first column represents the number of parameters for both flight conditions and the second and third columns represent estimators for respective flight conditions. Parameters estimated under different flight conditions have different values due to obvious change in roll and yaw rates. MSE and fit % for each estimation technique is written under respective heading in blue colour and uncertainties of each parameter is mentioned with symbol (±) in a purple colour. Results indicate that estimators were able to capture straight and level flight better than Coordinated Turn flight because of the obvious simplicity of flight conditions.

5. Conclusions

The open-loop, MIMO, coupled, lateral–directional dynamics of an agile aircraft in a straight and level flight condition and a coordinated turn flight condition were identified using discrete time grey box modelling SI and Bayesian techniques. Aircraft dynamical model structure has been identified using a BJ polynomial model, and parameter refinement has been achieved using a Bayesian probabilistic method. An Aircraft Flight Dynamics Model was designed and excited using multiple customized time-skewed inputs of aileron and rudder deflection as per the dictates of non-parametric model (FIR model). Input–output time histories were recorded and segregated as test and validation datasets. Initially, the OE model structure was developed to obtain initial parameters for the BJ model structure. Reduced order modelling was achieved using AIC and parameter error minimization was attained using an L-M algorithm. After obtaining the optimum model structure, parameter fine-tuning was conducted using the Bayesian approach on both flight conditions. The accuracy of the initial the OE model is less than the BJ model as it lacks capturing of the system if it exhibits nonlinearity. Different noise levels were also injected to supersonic aircraft MIMO system to postulate a model near to practical conditions. Controller design requirements are better fulfilled using mathematical models obtained from grey box modelling as they give precise and robust behaviour characteristics in comparison to constant values of aerodynamic stability and control derivatives obtained through first principle modelling.
A unique combination of Bayesian estimation with grey box modelling provides an opportunity to achieve a parsimonious and robust framework for parameter estimation of a dynamical system. The same has been seen in this research in which aircraft prior knowledge of parameters has been linked with MIMO observational data and improved MSE and percentage fit has been attained as shown in Table 6, Table 7, Table 8 and Table 9, wherein Bayesian estimation (with percentage fit of 96.25% and MSE of 3.9 × 10−3) has proved to be more accurate than other techniques such as NLS and MLE. However, this approach is mainly avoided in aircraft parameter estimation due to difficulties in finding pdf of random variables. In this research, Bayesian has been applied with lesser difficulty due to a lesser amount of dataset and simple flight conditions. The comparison matrix of NLS, MLE, and Bayesian also showed better refinement of parameters using Bayesian estimators. The same approach can be further extended to investigate the uncertainties in estimation of parameters leading towards online SI and robust controller design.

Author Contributions

All authors contributed to the research conception, design and results. M.F.M., contributed to the conceptualization, data formulation, analysis, methodology, software, writing—original draft and editing. S.M.A. contributed to the data correction and improvement of results. M.W. contributed to the conceptualization, methodology, and validation. Z.H.K. contributed to the methodology, editing and reviewing the results. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data presented in this study is available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Nomenclature

Roman
CD, CY, CL: Aerodynamic force coefficients
Cl, Cm, Cn: Aerodynamic moment coefficients
Throttle position: δt [0 1]
Inertia in x, y and z axis: Jx, Jy, Jz
Roll, Pitch and Yaw Rates: p, q, and r
Roll, Pitch and Yaw acceleration: p ˙ , q ˙ , and r ˙
F r e e   s t r e a m   d y n a m i c   p r e s s u r e : q [ N / m 2 1 b / f t 2 ]
I X , I Y , I Z moments of inertia about X , Y , and Z body axes, k g m 2 (slug- f t 2 )
I x z = p r o d u c t   o f   i n e r t i a   w i t h   r e s p e c t   t o   x   a n d   z body axes, kg m 2 (slug- f t 2 )
C X : X -axis force coefficient along positive x body axis, A e r o d y n a m i c   x a x i s   f o r c e q S ¯
C X , t total X -axis force coefficient
C Y : X -axis force coefficient along positive x body axis, A e r o d y n a m i c   x a x i s   f o r c e q S ¯
C Y , t total X -axis force coefficient
C z         Z-axis force coefficient along positive Z body axis, A e r o d y n a m i c   z a x i s   f o r c e q S
C Z , t total z -axis force coefficient
C m pitching-moment coefficient about Y body axis, A e r o d y n a m i c   p i t c h i n g   m o m e n t q S c
C m , t : total pitching-moment coefficient
c : wing mean aerodynamic chord, m ( f t )
Probability Density Function: p.d.f

Greek
Angle of attack: α (AOA)
Sideslip angle: β
Roll, pitch and yaw angles: Phi(φ), theta   ( θ ) , pssi(ψ)
Deflections in elevator, rudder and aileron: δe, δr, δa
Flight path angle: γ
Air Density: ρ
System Identification: SI
Finite Impulse Response: FIR
Flight Dynamics Model: FDM
Centre of Gravity: CG
Quadratic Cost Function: J
Throttle: THTL
Elevator: ELEV
Aileron: AIL
Rudder: RDR
Mean Square Error: MSE
Fir Percentage Error: FPE
Degree of Freedom: DOF
Mach Number: M
Aircraft mass: m [kg (slugs)]
Wing Span: b [ m , (ft)]
Wing Area: S [ m 2 ( f t 2 ) ]
Airplane body axes: x, y, z
C e n t e r o f   g r a v i t y   l o c a t i o n :   X c g
Single input Single Output: SISO
Multiple Input Multiple Output: MIMO
Output Error: OE
Box–Jenkins: BJ
Auto-Regressive Exogenous: ARX
Auto-Regressive Exogenous Moving Average: ARMAX
Weight: w
Feet: ft
Pound square feet: psf
Velocity: Vt
Height: h
Radians: rad

Superscripts
.: Time derivative
T: Transpose
^: Estimate

References

  1. Tangirala, A.K. Principles of System Identification Theory and Practice; CRC Press: Boca Raton, FL, USA, 2015. [Google Scholar]
  2. Klein, V.; Morelli, E. Aircraft System Identification Theory and Practice; American Institute of Aeronautics & Astronautics: Reston, VA, USA, 2006. [Google Scholar]
  3. Stevens, B.L.; Lewis, F.L. Aircraft Control and Simulation, 2nd ed.; John Wiley & Sons: Hoboken, NJ, USA, 2003. [Google Scholar]
  4. Tulleken, H.J.A.F. Grey-box Modelling and Identification Using Physical Knowledge and Bayesian Techniques. Automatica 1993, 29, 285–308. [Google Scholar] [CrossRef]
  5. Morelli, E.A.; Grauer, J.A. Advances in Aircraft System Identification at NASA Langley Research Center. J. Aircr. 2023, 60, 1354–1370. [Google Scholar] [CrossRef]
  6. Morelli, E.A.; Grauer, J.A. Practical Aspects of Frequency-Domain Approaches for Aircraft System Identification. J. Aircr. 2020, 57, 268–291. [Google Scholar] [CrossRef]
  7. Oznurlu, C. Data-Driven Model Discovery and Control of Lateral-Directional Fighter Aircraft Dynamics. Master’s Thesis, Middle East Technical University, Ankara, Turkey, 2022. [Google Scholar]
  8. Belge, E.; Kaba, H.K.; Parlak, A.; Altan, A.; Hacioğlu, R. Estimation of small unmanned aerial vehicle lateral dynamic model with system identification approaches. Balk. J. Electr. Comput. Eng. 2020, 8, 121–126. [Google Scholar] [CrossRef]
  9. Bnhamdoon, O.A.A.; Hanif, N.H.H.M.; Akmeliawati, R. Identification of a quadcopter autopilot system via Box–Jenkins structure. Int. J. Dyn. Control 2020, 8, 835–850. [Google Scholar] [CrossRef]
  10. Pillonetto, G.; Dinuzzo, B.F.; Chenc, T.; Nicolao, G.D.; Ljung, L. Kernel methods in system identification, machine learning and function estimation: A survey. Automatica 2014, 50, 657–682. [Google Scholar] [CrossRef]
  11. Kwad, A.M.; Hanafi, D.; Omar, R.; Rahman, H.A. Development of system identification from traditional concepts to real-time soft computing based. IOP Conf. Ser. Mater. Sci. Eng. 2020, 767, 012050. [Google Scholar] [CrossRef]
  12. Ji, Y.; Jiang, X.K.; Wan, L.J. Hierarchical least squares parameter estimation algorithm for two-input Hammerstein finite impulse response systems. J. Frankl. Inst. 2020, 357, 5019–5032. [Google Scholar] [CrossRef]
  13. Millidere, M. Optimal Input Design and System Identification for an Agile Aircraft. Ph.D. Thesis, Middle East Technical University, Ankara, Turkey, 2021. [Google Scholar]
  14. Morelli, E.A. Practical Aspects of the Equation-Error Method for Aircraft Parameter Estimation. In Proceedings of the AIAA Atmospheric Flight Mechanics Conference, Keystone, CO, USA, 21–24 August 2006. [Google Scholar]
  15. Roeser, M.S.; Fezans, N. Method for designing multi-input system identification signals using a compact time-frequency representation. CEAS Aeronaut. J. 2021, 12, 291–306. [Google Scholar] [CrossRef]
  16. Lopez, M.J.S.; Ruckel, P.; Berrigan, C.S. Bell V-280 System Identification and Model Validation with Flight Test Data Using the Joint Input-Output Method. In Proceedings of the Vertical Flight Society’s 76th Annual Forum & Technology Display, Virginia Beach, VA, USA, 5–8 October 2020. [Google Scholar]
  17. Simmons, B.M.; McClelland, H.G.; Woolsey, C.A. Nonlinear Model Identification Methodology for Small, Fixed-Wing, Unmanned Aircraft. J. Aircr. 2019, 56, 1056–1067. [Google Scholar] [CrossRef]
  18. Dutra, D.A.A. Collocation-Based Output-Error Method for Aircraft System Identification. In Proceedings of the AIAA Aviation Forum, Dallas, TX, USA, 17–21 June 2019. [Google Scholar]
  19. Mukhopadhaya, J.; Whitehead, B.T.; Quindlen, J.F.; Alonso, J.J. Multi-fidelity modeling of probabilistic aerodynamic databases for use in aerospace engineering. Int. J. Uncertain. Quantif. 2020, 10, 425–447. [Google Scholar] [CrossRef]
  20. Wang, Q.; Zheng, F.; Qian, W.; Ding, D. A practical filter error method for aerodynamic parameter estimation of aircraft in turbulence. Chin. J. Aeronaut. 2022, 36, 17–28. [Google Scholar] [CrossRef]
  21. Srivastava, A.; Kumar, A.; Ghosh, A.K. Determination of Parameters during Quasi-Steady Stall Maneuver Using Genetic Algorithm. Int. J. Aviat. Aeronaut. Space 2019, 6, 4. [Google Scholar]
  22. Zainuddin, F.A.; Samad, M.F.A. Crossover of genetic algorithm for linear and nonlinear system identification. In Proceedings of the Innovative Research and Industrial Dialogue, Melaka, Malaysia, 17 December 2020. [Google Scholar]
  23. Ordóñez, C.; Lasheras, F.S.; Roca-Pardiñas, J.; Juez, F.J.d.C. A hybrid ARIMA–SVM model for the study of the remaining useful life of aircraft engines. J. Comput. Appl. Math. 2019, 346, 184–191. [Google Scholar] [CrossRef]
  24. Gen, M.; Cheng, R. Genetic Algorithm and Engineering Optimization; Institute of Technology: Ashikaga, Japan, 2000. [Google Scholar]
  25. Fatima, S.K.; Abbas, S.M.; Mir, I.; Gul, F.; Forestiero, A. Flight Dynamics Modeling with Multi-Model Estimation Techniques: A Consolidated Framework. J. Electr. Eng. Technol. 2023, 18, 2371–2381. [Google Scholar] [CrossRef]
  26. Bagherzadeh, S.A. Flight dynamics modeling of elastic aircraft using signal decomposition methods. J. Aerosp. Eng. 2019, 233, 4380–4395. [Google Scholar] [CrossRef]
  27. Alessandrini, M.; Falaschetti, L.; Biagetti, G.; Crippa, P.; Turchetti, C. Nonlinear Dynamic System Identification in the Spectral Domain Using Particle-Bernstein Polynomials. Electronics 2022, 11, 3100. [Google Scholar] [CrossRef]
  28. Ram, P.R.M.; Römer, U.; Semaan, R. Bayesian Dynamical System Identification with Unified Sparsity Priors And Model Uncertainty. arXiv 2021, arXiv:2103.05090. [Google Scholar]
  29. Piga, D.; Breschi, V.; Bemporad, A. Estimation of Jump Box–Junction Models. Automatica 2020, 120, 109126. [Google Scholar] [CrossRef]
  30. Carr, J. 2014. An Introduction to Genetic Algorithms [Online]. Available online: https://www.whitman.edu.USA (accessed on 25 May 2024).
  31. Hamilton, J.; Wall, D.G.; Saddington, A.J.; Economou, J.T. Hypersurface normalised gain-scheduled controller for a non-linear 6-DOF fast jet. Aerosp. Sci. Technol. 2020, 106, 106155. [Google Scholar] [CrossRef]
  32. Rémy, P.; Hugo, G.; Ian, C.; Stephane, D.; Youssef, D.; Nathalie, B. An efficient application of Bayesian optimization to an industrial MDO framework for aircraft design. In Proceedings of the AIAA Aviation Forum, Virtual, 15–19 June 2020. [Google Scholar]
  33. Botero, E.M. Generative Bayesian Networks for Conceptual Aircraft Design. Ph.D. Thesis, Stanford University, Stanford, CA, USA, 2019. [Google Scholar]
  34. Kim, D.; Oh, H.S.; Moon, I.C. Black-box Modeling for Aircraft Maneuver Control with Bayesian Optimization. Int. J. Control Autom. Syst. 2019, 17, 1558–1568. [Google Scholar] [CrossRef]
  35. Saves, P.; Bartoli, N.; Diouane, Y.; Lefebvre, T.; Morlier, J. Constrained Bayesian Optimization over Mixed Categorical Variables, with Application to Aircraft Design. In Proceedings of the AeroBest, Lisbonne, Portugal, 21–23 July 2021. [Google Scholar]
  36. Scoggins, J.B.; Wignall, T.J.; Nakamura-Zimmerer, T.; Bibb, K. Multihierarchy Gaussian Process Models for Probabilistic Aerodynamic Databases using Uncertain Nominal and Off-Nominal Configuration Data. In Proceedings of the AIAA SCITECH 2023 Forum, National Harbor, MD, USA, 23–27 January 2023; NASA Langley Research Centre: Hampton, VA, USA, 2023. [Google Scholar]
  37. Napolitano, M.R. Aircraft Dynamics: From Modelling to Simulation; Scitus Academics LLC: Wilmington, DE, USA, 2012. [Google Scholar]
  38. Nguyen, L.T.E.A. Simulator Study of Stall/Post-Stall Characteristics of a Fighter Airplane with Relaxed Longitudinal Static Stability [F16]; NASA: Washington, DC, USA, 1979.
  39. F 16 Jet Design, Collimator. 2024. Designing an F 16 Fighter Jet [Online]. Available online: https://www.collimator.ai/tutorials/simulating-a-nonlinear-f16-model (accessed on 25 May 2024).
  40. Mohajerani, M.H. Frequency-Domain System Identification for Unmanned Helicopters from Flight Data. Master’s Thesis, Concordia University, Montreal, QC, Canada, 2014. [Google Scholar]
  41. Kullback, S.; Leibler, R.A. On Information and Sufficiency. Ann. Math. Stat. 1951, 22, 79–86. [Google Scholar] [CrossRef]
  42. Muhammad, W.; Ahsan, A. Airship aerodynamic model estimation using unscented Kalman filter. J. Syst. Eng. Electron. 2020, 31, 1318–1329. [Google Scholar] [CrossRef]
  43. Rundel, M.C.M.Ç.-R.C.; Banks, D.; Chai, C.; Huang, L. An Introduction to Bayesian Thinking. A Companion to the Statistics with R Course; Taylor and Francis: Abingdon, UK, 2020. [Google Scholar]
  44. Cover, T.M.; Thomas, J.A. Elements of Information Theory, 2nd ed.; Wiley: Hoboken, NJ, USA, 2006. [Google Scholar]
  45. Grauer, J.; Morelli, E. Method for real-time frequency response and uncertainty estimation. J. Guid. Control Dyn. 2014, 37, 336–344. [Google Scholar] [CrossRef]
Figure 1. Research Framework.
Figure 1. Research Framework.
Aerospace 11 00960 g001
Figure 2. Top Level Simulink Model of Aircraft Flight Dynamic Model (MATLAB-2021b).
Figure 2. Top Level Simulink Model of Aircraft Flight Dynamic Model (MATLAB-2021b).
Aerospace 11 00960 g002
Figure 3. F-16 6-DOF Dynamics [39].
Figure 3. F-16 6-DOF Dynamics [39].
Aerospace 11 00960 g003
Figure 4. Optimal Input Design Flowchart.
Figure 4. Optimal Input Design Flowchart.
Aerospace 11 00960 g004
Figure 5. Bayesian Implementation Flowchart.
Figure 5. Bayesian Implementation Flowchart.
Aerospace 11 00960 g005
Figure 6. F-16 Kinematics Variables [39].
Figure 6. F-16 Kinematics Variables [39].
Aerospace 11 00960 g006
Figure 7. Non-Parametric (FIR) Model of Aircraft.
Figure 7. Non-Parametric (FIR) Model of Aircraft.
Aerospace 11 00960 g007
Figure 8. Bode Plot Aircraft Lateral Dynamics.
Figure 8. Bode Plot Aircraft Lateral Dynamics.
Aerospace 11 00960 g008
Figure 9. Simulated Time-Skewed 2-1-1 Doublet Inputs—Aileron (δa) and Rudder (δr).
Figure 9. Simulated Time-Skewed 2-1-1 Doublet Inputs—Aileron (δa) and Rudder (δr).
Aerospace 11 00960 g009
Figure 10. Roll and Yaw Rate Time histories in repose to 2-1-1 Doublet Inputs.
Figure 10. Roll and Yaw Rate Time histories in repose to 2-1-1 Doublet Inputs.
Aerospace 11 00960 g010
Figure 11. Roll and Pitch Angle time histories to 2-1-1 Doublet Inputs.
Figure 11. Roll and Pitch Angle time histories to 2-1-1 Doublet Inputs.
Aerospace 11 00960 g011
Figure 12. Aircraft Parameter Refinement Flow chart.
Figure 12. Aircraft Parameter Refinement Flow chart.
Aerospace 11 00960 g012
Figure 13. (a) Initial OE Model; (b) Reduced Order OE Model; (c) Initial BJ Model; (d) Optimized BJ Model; (e) Residual Correlation; (f) pdf of Model Parameters; (g) Posterior Sensitivity Analysis (K-L Divergence)—Straight and Level Flight.
Figure 13. (a) Initial OE Model; (b) Reduced Order OE Model; (c) Initial BJ Model; (d) Optimized BJ Model; (e) Residual Correlation; (f) pdf of Model Parameters; (g) Posterior Sensitivity Analysis (K-L Divergence)—Straight and Level Flight.
Aerospace 11 00960 g013aAerospace 11 00960 g013b
Figure 14. (a) Initial OE Model; (b) Reduced Order OE Model; (c) Initial BJ Model; (d) Optimized BJ Model; (e) Residual Correlation; (f) pdf of Model Parameters; (g) Posterior Sensitivity Analysis (K-L Divergence)—Coordinated Turn Flight.
Figure 14. (a) Initial OE Model; (b) Reduced Order OE Model; (c) Initial BJ Model; (d) Optimized BJ Model; (e) Residual Correlation; (f) pdf of Model Parameters; (g) Posterior Sensitivity Analysis (K-L Divergence)—Coordinated Turn Flight.
Aerospace 11 00960 g014
Table 1. Summary of Bayesian Approaches used for Aircraft Design and Modelling.
Table 1. Summary of Bayesian Approaches used for Aircraft Design and Modelling.
Ref.Bayesian MethodologyApplication
Herbert et al. [4]Bayesian framework to estimate aerodynamic model parametersAircraft parameter estimation
Mukhopadhaya et al. [19]Gaussian Process Regression, close to Bayesian estimation, for improvement of uncertainty of aerodynamic databaseAerodynamic Modelling using database
Rémy Priem et al. [32]Using Bayesian approach for optimization in Aircraft design configurationsAircraft Design Configurations
Emilio M. Botero [33]Using Generative Bayesian Network for conceptual design of aircraftAircraft Design
Kim et al. [34]Bayesian network optimization for data-driven controller of the aircraft using black box modellingData-driven controller design for aircraft maneuver
Paul Saves et al. [35]Aircraft design optimization using Bayesian approach and Gaussian Process.Aircraft Design
James et al. [36]Mathematical framework for modelling probabilistic aerodynamic datasets using conditional coupled with Gaussian ProcessesAerodynamic Modelling
Table 2. Interpreting Bayes Factor [43].
Table 2. Interpreting Bayes Factor [43].
BF   [ H 1 : H 2 ] Evidence   Against   H 2
1 to 3Not worth a bare mention
3 to 20Positive
20 to 150Strong
>150Very strong
Table 3. Aircraft Geometry and Nominal Flight Conditions.
Table 3. Aircraft Geometry and Nominal Flight Conditions.
SymbolValueUnit
b30ft
c 11.32ft
S300 f t 2
W20,500Lbs
gd32.17 ft / s 2
Vt502ft/s
h30 ft
q 300psf
Xcg 0.3 c ft
Ixx9496 slug-f t 2
Iyy55,814 slug-f t 2
Izz63,100 slug-f t 2
Ixz982 slug-f t 2
Table 4. Linearized Flight Conditions.
Table 4. Linearized Flight Conditions.
SymbolSteady Straight and Level FlightCoordinated Turn FlightUnit
Vt502502ft/s
h30 30 ft
q 300300psf
Xcg 0.35 c 0.35 c ft
α 0.030.24rad
β 00rad
φ01.3rad
θ 0.150.05rad
P0−0.01rad/s
Q00.29rad/s
R00.06rad/s
ψ ˙ 00.15rad/s
Table 5. Frequency Parameters Obtained from Impulse Response.
Table 5. Frequency Parameters Obtained from Impulse Response.
PolesDampingFrequency
(rad/s)
Time Constant
(s)
−4.78 × 10−31.004.78 × 10−32.09 × 102
−2.221.002.224.50 × 10−1
−9.24 × 10−1 + 4.67i1.94 × 10−14.761.08
−9.24 × 10−1 − 4.67i1.94 × 10−14.761.08
Table 6. Lateral Dynamics (Straight and Level).
Table 6. Lateral Dynamics (Straight and Level).
Parameter Estimates and Errors
ParametersNLS
MSE: 9.3 × 10−2
% fit—87.99%
MLE
MSE: 5.2 × 10−3
% fit—95.75%
Bayesian
MSE: 3.9 × 10−3
% fit—96.25%
b11.2319 ± 0.52241.1129 ± 0.01131.1036 ± 0.0690
b2−1.5237 ± 1.4569−0.6976 ± 0.2309−0.9870 ± 0.6327
b30.6885 ± 0.26850.5832 ± 0.43120.6799 ± 0.5177
b4−1.2097 ± 0.2557−1.1414 ± 0.6314−1.1567 ± 0.4477
c1−0.5983 ± 0.3882−0.2188 ± 0.2424−0.3291 ± 0.3066
d11.1822 ± 0.28400.5998 ± 0.329840.7939 ± 0.2365
d2−0.3579 ± 0.2333−0.4501 ± 0.2069−0.3026 ± 0.1699
d30.0206 ± 0.08410.0566 ± 0.09880.0716 ± 0.1455
d4−0.7206 ± 0.1603−0.0619 ± 0.1422−0.0752 ± 0.1399
f1−2.3801 ± 0.2184−1.444 ± 0.1977−1.3332 ± 0.2503
f21.6348 ± 0.19541.5924 ± 0.31661.7412 ± 0.3417
f30.2456 ± 0.11341.1868 ± 0.29591.1912 ± 0.2789
f4−0.7353 ± 0.1399−0.8423 ± 0.0978−0.7926 ± 0.1255
f50.2426 ± 0.2953−0.3576 ± 0.0463−0.3525 ± 0.0132
Table 7. Lateral Dynamics (Coordinated Turn).
Table 7. Lateral Dynamics (Coordinated Turn).
Parameter Estimates and Errors
ParametersNLS
MSE: 0.2838
% fit—77.51%
MLE
MSE: 0.1549
% fit—72.27%
Bayesian
MSE: 0.1076
% fit—80.07%
b1−2.7703 ± 0.5224−2.6948 ± 0.6282−2.8416 ± 0.3352
b22.0599 ± 1.45691.8639 ± 1.05051.2786 ± 0.8046
c1−0.2802 ± 0.26850.0832 ± 0.41410.0058 ± 0.9224
c20.1900 ± 0.25570.6910 ± 0.25990.7222 ± 0.2309
c3−0.0267 ± 0.38820.7470 ± 0.33670.8252 ± 0.2990
c4−0.8072 ± 0.2840−0.6734 ± 0.3529−0.6921 ± 0.2244
c50.7626 ± 0.23330.7639 ± 0.18130.4692 ± 0.9028
d1−1.3011 ± 0.0841−0.5260 ± 0.1892−0.7591 ± 0.8219
d21.3436 ± 0.16031.3588 ± 0.29411.0851 ± 0.9780
d3−1.0987 ± 0.2184−1.0939 ± 0.2717−0.9738 ± 0.4596
d40.7316 ± 0.19540.8825 ± 0.24870.6657 ± 0.4727
d5−0.1844 ± 0.1134−0.4063 ± 0.1648−0.0999 ± 0.2929
f10.4115 ± 0.1399−1.1599 ± 0.3476−1.5465 ± 0.2038
f2−1.8137 ± 0.29530.5206 ± 0.18620.4772 ± 0.1453
f30.0845 ± 0.17300.0990 ± 0.08240.6454 ± 0.0428
f4−0.2993 ± 0.0574−0.3032 ± 0.0737−0.2414 ± 0.0276
f50.2120 ± 0.05770.2268 ± 0.09620.4398 ± 0.0652
Table 8. Directional Dynamics (Straight and Level).
Table 8. Directional Dynamics (Straight and Level).
Parameter Estimates and Errors
ParametersNLS
MSE: 6.8 × 10−1
% fit—78.66%
MLE
MSE: 1.3 × 10−1
% fit—80.75%
Bayesian
MSE: 1.4 × 10−2
% fit—95.49%
b1−0.0399 ± 0.2376−0.2424 ± 0.2980−0.3906 ± 0.0531
b20.3101 ± 0.13750.3612 ± 0.15190.1206 ± 0.2529
c1−1.9760 ± 1.0797−1.9883 ± 0.7384−0.6182 ± 0.854
c20.9760 ± 0.92340.9895 ± 1.29430.8582 ± 0.7253
d1−0.8397 ± 0.2274−0.9174 ± 0.2912−2.5870 ± 1.212
d2−0.4616 ± 0.5261−0.6852 ± 0.47272.3138 ± 2.3003
d30.5737 ± 0.56130.8832 ± 0.3967−0.6969 ± 1.215
f1−0.6638 ± 0.5958−1.5629 ± 0.3434−1.4591 ± 0.728
f2−0.7063 ± 0.60850.3123 ± 0.90050.4699 ± 0.7319
f30.4186 ± 0.45581.2015 ± 0.95781.1084 ± 0.1574
f40.0585 ± 0.31240.1958 ± 0.83340.0446 ± 0.1438
f50.2305 ± 0.31420.4038 ± 0.39820.0814 ± 0.1458
Table 9. Directional Dynamics (Coordinated Turn).
Table 9. Directional Dynamics (Coordinated Turn).
Parameter Estimates and Errors
ParametersNLS
MSE: 4.4 × 10−2
% fit—82.93%
MLE
MSE: 4.3 × 10−2
% fit—83.75%
Bayesian
MSE: 3.9 × 10−2
% fit—85.49%
b10.6851 ± 0.0904−0.3410 ± 0.0980−0.3465 ± 0.0987
b2−2.6887 ± 0.45300.3602 ± 0.12860.3642 ± 0.1184
c12.4674 ± 0.6864−0.4502 ± 0.7700−0.6261 ± 0.863
c2−0.7786 ± 0.27890.2385 ± 0.73960.3914 ± 0.8208
d1−2.2853 ± 0.1387−1.6688 ± 0.6378−1.6878 ± 0.614
d22.3718 ± 0.28711.4156 ± 1.07891.4251 ± 1.0266
d3−1.3191 ± 0.2246−0.7082 ± 0.6029−0.7017 ± 0.5427
f10.3695 ± 0.0686−1.9609 ± 0.3877−1.9085 ± 0.360
f2−2.1439 ± 0.11891.2114 ± 0.54841.1118 ± 0.5213
f32.1391 ± 0.2527−0.1710 ± 0.4755−0.1088 ± 0.484
f4−1.0981 ± 0.23410.0248 ± 0.49610.0107 ± 0.502
f50.1920 ± 0.1106−0.0220 ± 0.3423−0.0126 ± 0.344
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Mazhar, M.F.; Abbas, S.M.; Wasim, M.; Khan, Z.H. Bayesian Identification of High-Performance Aircraft Aerodynamic Behaviour. Aerospace 2024, 11, 960. https://doi.org/10.3390/aerospace11120960

AMA Style

Mazhar MF, Abbas SM, Wasim M, Khan ZH. Bayesian Identification of High-Performance Aircraft Aerodynamic Behaviour. Aerospace. 2024; 11(12):960. https://doi.org/10.3390/aerospace11120960

Chicago/Turabian Style

Mazhar, Muhammad Fawad, Syed Manzar Abbas, Muhammad Wasim, and Zeashan Hameed Khan. 2024. "Bayesian Identification of High-Performance Aircraft Aerodynamic Behaviour" Aerospace 11, no. 12: 960. https://doi.org/10.3390/aerospace11120960

APA Style

Mazhar, M. F., Abbas, S. M., Wasim, M., & Khan, Z. H. (2024). Bayesian Identification of High-Performance Aircraft Aerodynamic Behaviour. Aerospace, 11(12), 960. https://doi.org/10.3390/aerospace11120960

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop