Physics-Based Aircraft Dynamics Identification Using Genetic Algorithms

: This research introduces a physics-based identification technique utilizing genetic algorithms. The primary objective is to derive a parametric matrix, denoted as A , describing the time-invariant linear model governing the longitudinal dynamics of an aircraft. This is achieved by proposing a fitness function based on the properties of the transition matrix and taking advantage of some of the capabilities of the genetic algorithm, mainly those of restricting the search ranges of the unknowns. In this case, such unknowns are related to the type of aircraft and flight conditions that are considered during the identification process. The proposed identification method is validated with a reliable nonlinear model that can be found in the literature, as well as with the calculation of the trim condition and linearization generally used in aircraft dynamics. In summary, this study suggests that the genetic algorithm provided with the adequate fitness function could be an appealing alternative for aircraft model identification, even when limited data are available. Furthermore, in some cases, linearization using a genetic algorithm can be more efficient than classical methods.


Introduction 1.Motivation
One of the fundamental pursuits in the scientific domain is the development of mathematical models that accurately represent physical systems.A fundamental criterion for evaluating such models is their usefulness.This may involve the ability to predict crucial aspects of the behavior of a physical system or to provide critical information to the design process through model parameters.In all cases, the mathematical model must achieve a balance between simplicity to be of practical use, and complexity to describe the most important characteristics of the system being modeled [1,2].
A long-standing need in aeronautics has been that of obtaining reliable and accurate models for aircraft dynamics.This goes beyond aircraft design and flight simulation; it involves the implementation and validation of control systems.In the current era, because of the impressive development of unmanned aerial vehicles (UAVs), the search for control systems has become more demanding, both in civil and military applications.The arrival of high-reliability, low-cost sensors and sophisticated control systems has broken down previous barriers.Small aircraft can now reach higher altitudes, higher speeds, and longer ranges.This provides the opportunity to implement cutting-edge algorithms, which are both smarter and more efficient [3].
However, the effort to identify an optimal model for advanced control system design objectives remains a costly and complex task.Given this reality, it is imperative to reevaluate the identification process, seeking an innovative approach that is simple, fast, and effective.It is crucial to avoid heuristic models and costly experiments, recognizing that the simplest option is often the most prudent.Several compelling and direct examples can be found in [4].Over the past two decades, identification in aeronautics, particularly for manned airplanes, has been vigorously explored by prominent civil military and government entities such as the Deutsches Zentrum für Luft-und Raumfahrt (DLR) and NASA [5].
A correct identification system incorporates experimental data and undergoes validation in flight.Nowadays, a huge amount of information on experimental identification and validation processes can be obtained from various technical reports and manuals [6], affording insights into many facets of authentic identification procedures.Recent research [5] offers an excellent overview of the state of the art in this field.Tools for the time domain are comprehensively examined in [7], with additional contributions from the German Aerospace Center featured in [8].Furthermore, system identification techniques for aircraft dynamics in the frequency domain have been provided in [9].
In [5], it is noted that advances in machine learning have yet to yield substantial benefits for aircraft system identification, primarily due to the absence of a physics-based model.This limitation becomes evident when employing neural networks to represent dynamics [10] or aircraft dynamics as a black box [11], or when resorting to genetic algorithms in pursuit of a heuristic model [12].There are also published works where two techniques are combined; for instance, in [13], genetic algorithms are part of a hybrid identification system originally based on optimization, and, in [14], the blend of neural networks and fuzzy systems is considered to obtain the longitudinal model of an aircraft.
When the integration of artificial intelligence is judiciously combined with an understanding of physics-based aircraft models, along with the aforementioned state-of-the-art knowledge, it enables the imposition of parametric constraints that facilitate the use of swift convergence algorithms.This approach circumvents the conventional issue of infinite solutions in computational algorithms and yields more cost-effective results compared to relying solely on artificial intelligence or physics-based algorithms [15].As demonstrated, this robust identification methodology proves invaluable not only for linear dynamics identification but also for the linearization under diverse trim conditions.It offers a simple, effective, and tailored solution to the formidable challenges encountered in today's aerospace industry.
It is generally recognized that the identification of nonlinear models poses a greater challenge compared to linear models.For example, successful identification through neural networks requires an appropriate determination of the number of layers, neurons per layer, and learning factor, all depending on the characteristics of the nonlinear system under consideration.Similarly, for fuzzy systems, an optimal approximation depends on the selection of an appropriate number of fuzzy rules, together with a judicious choice of membership functions.Even in the case of SINDY, customization is often imperative to address the unique complexities associated with the nonlinear system being identified [16].On the other hand, when the mentioned techniques are applied to identify linear systems, the problem reduces to parametric identification [17].
While a large number of approaches are available in the literature to successfully identify both linear and nonlinear systems, this paper focuses specifically on the utilization of a practical and effective technique based on genetic algorithms to linearize the longitudinal dynamics of an aircraft by introducing a new fitness function.

Contribution
A notable disadvantage of the application of genetic algorithms in obtaining mathematical models from real or simulated data is that, during the evaluation of the fitting function, the model obtained must be simulated in each iteration to determine the error between the data and the response of the intermediate model.This causes a very high com-putational cost that may be out of the reach of researchers with limited budgets.To reduce this computational cost, and considering that the linear system is sufficient to describe the longitudinal dynamics of an aircraft, in this work, the transition matrix [18,19] is considered instead of the complete simulation of the intermediate model in each iteration of the genetic algorithm.
Therefore, the main goal of the current work is that of presenting a fitness function based on the transition matrix, such that, when the fitness function is included in the genetic algorithm, the matrix, A, involved in the linear approximation of the aircraft, is promptly obtained, and relatively few data are needed.

Manuscript Organization
The rest of this paper is organized as follows.In Section 2, the nonlinear model of the aircraft is presented.Section 3 describes the trim condition considered.In Section 4, the traditional linearization of the nonlinear model is performed.Below, Section 5 presents a concise analysis of the parameters included in matrix A. The main result, i.e., the linearization achieved by the genetic algorithm, is discussed in detail in Section 6.Finally, in Section 7, some conclusions are presented.

Nonlinear 6-DOF Dynamic Model of a Research Aircraft
The primary function of a nonlinear model is to represent and capture complex relationships between variables that do not follow a linear pattern, allowing greater flexibility in modeling intricate patterns and dependencies.It is important to mention that, although in the current work only the longitudinal dynamics are excited, it is necessary to pay attention to the complete dynamics of the aircraft.Next, the nonlinear model used in the analysis and validation of the identification process of the genetic algorithm is presented synthetically.This complete and validated model was extracted from [20], which is an excellent summary of [21].
In accordance with rigid body dynamics and classical approximation from Newton-Euler, the following analysis can be carried out: Translation equation on body axes: where m is the mass, the sub-index cm/e means from the center of mass to Earth axis, and super-index b indicates body frame; this nomenclature comes from kinematics equations deduction (see [22]), with Vb cm/e = [ u v ẇ] T as the linear acceleration vector, V b cm/e = [u v w] T as the linear velocity vector, F b as the external forces, ωb b/e = [ ṗ q ṙ] T as the angular acceleration, and ω b b/e = [p q r] T as the angular rate.Rotational equations on body axes: where J b is the inertial matrix, and M b are the external moments.Kinematic Euler Equation: At this point, it is important to consider a transformation of angular rates from body axes to the general frame Euler system (inertial), as follows: where Φ = [ φ θ ψ] T is the vector of Euler angles rate and The full nonlinear model is considered nonautonomous and depends on the full state (even forces and moments are direct functions of the states) and control inputs, resulting in the following: with the state vector X = x 1 . . .x 12 T = u v w p q r ϕ θ ψ P N P E h T , where P N , P E , and P D are the x − y − z aircraft displacements on North-East-Down navigation frame system and h = −P D .
As mentioned above, external forces and moments are functions of states and control inputs, resulting in the following: which ultimately leads to In Equation ( 5), CS(ϕ, θ, ψ) is the body to navigation axis transformation: where c(•) and s(•) are cosine and sine function, respectively.

The Aircraft
The original airplane is a twin-engine research civil aircraft with a weight mass of 120,000 kg and inertial moment I y of 7.68 × 10 6 kg • m 2 , a mean aerodynamic chord of 6.6 m, and a wing platform area of 260 m 2 .The dynamics of the control input and the allowable limits on the elevator are specified as −15 ≤ δe ≤ 15 degrees per second and −25 ≤ δ e ≤ 10 degrees, respectively.For the two motors, the rate and end limits are −1.6 ≤ δ th ≤ 1.6 degrees per second and 0.5 ≤ δ th ≤ 10 degrees.

Trim Condition
In the context of identification, the equilibrium condition will be a steady-state flight with wings level and no sideslip [22].Determining the equilibrium point in a coupled, non-linear, multivariable system, such as Equation ( 9), is nontrivial and requires a specialized algorithm [23].Before linearization, a fitting condition for the system must be selected in the general form Ẋ = f (X, U): f (X e , U e ) = 0.
In accordance with Equation ( 5), navigation variables X e ∈ R 9 and U e ∈ R 5 are neglected, i.e., u = v = ẇ = ṗ = q = ṙ = 0 for steady-state flight.At this point, the methodology presented in [22] is considered, whose problem statement is the following: where H is a weight diagonal matrix, and C ∈ R 9+1 is the restriction vector.The process is then performed with the Nelder and Mead Simplex minimization search algorithm, using the command fminsearch() MATLAB, which is a nonlinear programming solver to find the minimum of an unconstrained multivariable function using the derivative-free method.
The cost function and initial lookup values are required.In this case, all functions in C that are equal to zero are in Equation ( 9), V e is the required operation speed in trim condition, |V b cm/e | = x 2 1 + x 2 2 + x 2 3 , and γ = x 8 − atan(x 3 /x 1 ) is the flight path angle.The weight matrix was defined as H = diag(3) after a trial and error process that provided valid results.The final function cost is f (X) = C T HC.
Using fminsearch() with the cost function defined above, after 9 iterations and with a minimization error of 1 × 10 −19 , the clipping condition X e obtained is x 1 = 110 m/s, x 3 = −2.0882m/s, x 8 = −0.0634rad, and all remaining states are zero.

Linearization
For the coefficient identification process, it is necessary to obtain the time-invariant linear system around the equilibrium point X e .To this end, one may consider Equation ( 4) with ẋ ∈ R n , x ∈ R n , u ∈ R m , and the mapping function f : R n+m −→ R n .This expression can be rewritten in an implicit form as follows: where f : R 2n+m −→ 0 n .By linearizing Equation (10) through a Taylor expansion series around the equilibrium point and considering only low-order terms (linear), the following results are obtained: Now, the equilibrium condition, must be evaluated.Then, from Equations ( 10) and ( 12), ẋ = ẋ0 , x = x 0 , u = u 0 , and define the partial derivatives as and, by simplifying Equation ( 11), one has Thus, solving for δ ẋ yields the following: where It is important to mention that the LTI system ẋ = Ax + Bu with matrices Equations ( 18) and ( 19) will provide values of perturbations (δx) from an initial equilibrium condition X e , as in the case of the control inputs.
The partial derivatives appearing in Equations ( 13)-( 15) can be obtained by using a numerical algorithm that employs a finite differential approximation, a feasible approach in simulation where one has control over the reliability and integrability of the data.While the process is relatively straightforward, careful consideration must be afforded to numerical errors, which can be mitigated by avoiding excessively small values of h.In the literature, it is widely acknowledged that central differentiation is more accurate compared to one-sided difference formulas [24,25].Then, for any partial derivatives, one has where h = x(t k+1 ) − x(t k ).This method is considered to determine Ê, Â, and B column by column, i.e., To obtain data and calculate the matrices in Equation ( 21), the initial condition and the input signals are used as perturbations.For instance, to produce a response in all states of the nonlinear model with no input signal, x(0) 1 = 10 m/s and x(0) 2 = 5 m/s are chosen as initial conditions; these values are additional to the trim condition.Then, after an increment of the state variable x j from h j = 1 × 10 −12 , the function F i is computed from each state x j to finally calculate the partial derivatives that make up matrices as in Equation (21).According to Equations ( 18) and ( 19), the state-space matrices for the 6-degree-of-freedom linear system are determined.Such a system is described by the differential equation ẋ = Ax + Bu, where x ∈ R 9 and u ∈ R 5 .However, the longitudinal dynamics can be separated as in the following reduced system: The matrix A is directly related to Equation (24) and is considered parametric, implying a physical interpretation.

Parameters
The matrix coefficient in a state-space form has a well-known analytic expression in the field of flight dynamics.In this particular context, only longitudinal dynamics are considered.The matrix A is square and corresponds to the vector state x = [u w q θ] T .Each coefficient is assigned a specific nomenclature, delineating moments and forces as functions of the state variables, i.e., In Equation ( 24), X θ and Z θ represent the linearized force projections due to gravity on body axes and are direct functions of the trim condition defined by θ e ; M refers to moments on body axes, specifically controlling pitch.Data for the common coefficient range parameters and their physical interpretation can be found in various sources, including [22,[26][27][28].These significant values are contingent upon the specific aircraft and prevailing flight conditions.

Main Result: Linearization by Means of the Genetic Algorithm (GA)
In this section, the nonlinear model introduced in Section 2 is employed as the benchmark for the evaluation of the performance of the linearization methods outlined in Section 4, including the genetic algorithm-based linearization approach.The maneuver considered corresponds to a cruise flight condition wherein the aircraft is traveling at a speed of 122.8357 m/s at sea level.The specified state variables include an angle of attack of 1.1354 degrees, a pitch angle of −3.7701 degrees, and all remaining states set to zero.The trim elevator input is −5.9505 degrees, and both motors operate at the 67.5% of total power.Utilizing the nonlinear model defined by Equations ( 1)-( 6), the longitudinal behavior of the aircraft is simulated, deviating from the trim condition with perturbations of ∆u nl 0 = 10 m/s, ∆α nl 0 = 2.3309 degrees, ∆q nl 0 = 12 degrees per second, and ∆θ nl 0 = 0 degrees over a duration of t = 180 s with a sampling time of T = 0.05 s.
In order to establish a reference for comparing the genetic algorithm (GA) linearization, the linear system with matrix A, derived through the linearization procedure detailed in Section 4 and represented by Equation (22), is simulated under identical conditions.The numerical results are illustrated in Figure 1, considering only perturbation data from the trimmed flight condition above.The mean squared error (MSE) between the linear states and their corresponding nonlinear counterparts is computed as 0.0034.Next, the nonlinear system is linearized using a genetic algorithm.Before presenting the results, a concise description of the genetic algorithm is provided.

Genetic Algorithm (GA)
A genetic algorithm (GA) is an optimization technique that mimics the processes of evolution and natural selection.Usually, it is applied to solve complex optimization problems.
In a genetic algorithm, a potential solution to the problem is represented as a chromosome, typically as a binary string, but other representations are also possible, depending on the problem at hand.A population is a collection of such chromosomes.The size of the population is a user-defined parameter.The fitness function evaluates how well each chromosome in the population solves the problem.It assigns a fitness score to each chromosome based on the goal of the problem.The fitness function is problem-specific and provided by the user.
Chromosomes with higher fitness scores have a greater chance of being selected for reproduction.The most common selection methods include roulette wheel selection and tournament selection.Crossover is the process of combining two parent chromosomes to create one or more offspring.Various crossover operators are used, such as one-point, two-point, and uniform crossover.The crossover rate (probability of crossover) is a userdefined parameter.Mutation involves randomly changing some bits in a chromosome with a low probability.It introduces genetic diversity into the population.The mutation rate (probability of mutation) is another user-defined parameter.Genetic algorithms terminate when a stopping criterion is met.Common termination conditions include a maximum number of generations, reaching a satisfactory fitness level, or a time limit.
A typical genetic algorithm operates in a loop, which includes selection, crossover, mutation, and replacement of the old population with the new one.This process continues for a predefined number of generations [29][30][31][32][33].It is important to mention that the GA has been successfully used to solve a number of problems in aeronautics.For instance, a methodology involving a multiobjective genetic algorithm for exploring parametrized microstructures is presented in [34].In [35], the authors investigate the aerodynamics of desert locust tandem wings through computational fluid dynamics (CFD) simulations.Using 2D and 3D Navier-Stokes equations, they explore wing interactions and corrugation effects.Optimization with genetic algorithms and Nash game theory improves gliding performance by at least 77%.In [36], an automated framework for developing interpolatable aeroelastic reduced-order models (AE ROMs) across diverse flight conditions is presented.By combining system identification, state-consistency enforcement (SCE), and a genetic algorithm (GA), the approach addresses the issue of state inconsistency in ROMs.
In MATLAB, optimization via genetic algorithm (GA) can be performed using the function ga, which is part of the Global Optimization Toolbox.The syntax is as follows: x = ga( f un, nvars, A, b, Aeq, beq, lb, ub, nlcon, opts), (25) where f un is the fitness function to be optimized; x is the vector of variables included in the fitness function; nvars is the dimension of the vector x; the pair (A, b) defines the inequality restrictions in the form of A • x ≤ b; the pair (Aeq, beq) defines the equality restrictions in the form of Aeq • x = beq; the vectors lb and ub, both of dimension nvars, define the lower and upper limits for the search of x, respectively; nlcon defines the nonlinear constraints; opts allows the user to modify parameters such as tolerance, the type of plot to be depicted, the maximal number of generations, etc.For more details, please refer to [37].
An essential challenge when working with GA is the need to define an appropriate fitness function.This is a crucial factor as the fitness function plays a critical role in evaluating how effectively a solution addresses the specific problem.The selection of the fitness function can significantly influence the overall performance of the algorithm and the quality of the solutions it produces.
With this in mind, a simple yet practical fitness function is introduced in the following sections.The function is designed to derive a parametric matrix A that characterizes the longitudinal dynamics of an aircraft.It utilizes data obtained from the simulation of the highly reliable nonlinear aircraft model described in Section 2.

Linearization via GA
In this section, the objective is to find a matrix A in the form of Equation ( 24) for the linearization of the longitudinal dynamics of the aircraft described in Equations ( 1)-( 6) using genetic algorithms (GA).To accomplish this, a fitness function is derived in a practical manner.The unknown x ′ s values to be determined by the GA are organized in matrix A GA as follows: x(7) −9.7903 x(2) x( 5) 107.9602 x(9) x(3) x (6) x( 8) where the constants appearing in Equation ( 26) are fixed values dependent on the characteristics of the modeled aircraft.In other words, only nine of the sixteen elements of the matrix regarding Equation ( 26) need to be determined.It is proposed that the fitness function evaluates the matrix described in Equation ( 26) in each iteration, comparing the linear states with their corresponding nonlinear states at arbitrary sample instants.To minimize computational load, the inclusion of the transition matrix is considered in this process.

Transition Matrix
The transition matrix, denoted as Φ(t − t 0 ) = e A(t−t 0 ) , describes the evolution of a time-invariant linear system over time, considering a given initial condition x(t 0 ).It shows how the state of the system at an initial time t 0 is transformed to the state at a later instant t.For continuous-time systems described by ẋ(t) = Ax(t), the transition matrix can be used to obtain x(t) from the initial condition x(t 0 ) (with t 0 = 0) as follows: Please refer to [18,19] for further details.In MATLAB, the transition matrix can be computed using the function expm(), and its application in this work is explained below.

Fitness Function
For optimized computational resource utilization, the proposal is to apply the fitness function to a reduced set of samples.Rather than utilizing the entire set of simulated data, consisting of 3601 data points collected over a duration of 180 s with a sampling interval of T = 0.05 s, a smaller subset must be considered.Such a subset is suggested to comprise the region where the linear approximation is desired to closely match the nonlinear data.
Taking this into account, from the simulation of the nonlinear model Equations ( 1)-( 6) and based on the transition matrix mentioned above, the following algorithm describes the fitness function used to linearize the nonlinear dynamics of the aircraft. Explanation: • Line 1: The name of the function is f itness, and the argument x ∈ R 9 is the vector of the variables to be found by the GA in order to minimize this function.• Line 2: Global variables ndata, x NL , x 0 , and t are used in this function, where ndata is the number of samples considered during the linearization, x NL ∈ R ndata×4 is a matrix containing the values of the nonlinear states u, w, q, θ, from the numerical simulation at desired instants, x 0 represents the vector of initial conditions, and t ∈ R ndata is a vector with the corresponding time instants at which matrix x NL has been obtained.• Line 3: Matrix A GA is constructed from elements of vector x from a previous iteration.• Line 4: Matrix x L ∈ R ndata×4 is initialized as a matrix of zeros with ndata rows (number of samples) and four columns (states u, w, q, θ).

•
Lines 5-7: A loop iterates over the elements of t, using the exponential matrix operation to compute the states of x L at each one of the time instants included in vector t.• Lines 8-11: Squared errors sqe i for each state, with i = 1, . . ., 4. • Line 12: The total error f is computed as the square root of the sum of squared errors.• Line 13: The function returns f as the result of the fitness evaluation.

Numerical Simulations
In this section, the results of two linearizations using the genetic algorithm (GA) are presented.These experiments were performed on a personal computer with the following specifications: Intel(R) Core(TM) i5-10500H CPU @ 2.50 GHz, 16 GB of RAM, running MATLAB 2019b.

Simulation 1: GA Linearization Using 31 Samples
In this case, only the samples obtained during the first 3 s of the nonlinear simulation with a sampling interval of T = 0.1 s are considered.This involves selecting every other sample from the nonlinear simulation, resulting in a total of 31 samples.In this example, according to the fitness function outlined in Algorithm 1, one has ndata = 31, x NL ∈ R 31×4 , x 0 = [1050.20940]T , and t ∈ R 31 .

Algorithm 1 Function f itness
global ndata, x NL , x 0 , t 3: x(7) −9.7903 x(2) x(5) 107.9602 x(9) x(3) x (6) x( 8) x L ← zeros(size(t), for k = 1 to ndata do 6: x end for 8: 10: f ← √ sqe 1 + . . .+ sqe 4 13: return f 14: end procedure During the execution of the GA in MATLAB, neither inequality nor equality restrictions nor nonlinear constraints are considered.However, due to the limitations of the available computer resources, the decision was made to search for the solution around the values obtained by the linearization procedure explained in Section 4.
Therefore, the GA in MATLAB was executed using the following three lines of code: with lb = [−0.06,−0.3, −0.005, −0.002, −1, −0.05, 5.5, −2, 0] ub = [0, 0, 0, 0, 0, 0, 7, 0, 0.7], It is worth mentioning that the options chosen in Equation ( 28) are primarily intended for monitoring the GA's behavior during execution, which may lead to a higher consumption of computer resources.A more refined selection of these parameters could potentially yield more efficient results.For a comprehensive description of MATLAB's function ga, interested readers are referred to [37].
The elapsed time is 355.863055s, and the matrix obtained by the GA is Then, matrix Equation ( 30) is used to simulate the linear system for 180 s with a sampling time of T = 0.05 s in order to compare the states of the linear identified system with the nonlinear ones.The results are presented in Figure 2. The overall mean square error (MSE) between the nonlinear states and their linear approximations is 0.0167.Notice that the MSE obtained in this experiment is greater than the one obtained by the linearization of Section 4, but the GA has been executed considering only 31 samples.In the next experiment, a larger set of samples is considered.6.3.2.Simulation 2: GA Linearization Using 66 Samples Now, the samples that occurred during the initial 3 s of the nonlinear simulation with a sampling time of T = 0.1 s, along with samples taken from the remainder of the signal with T = 5 s, are considered, resulting in a total of 66 samples.Therefore, one has ndata = 66, x NL ∈ R 66×4 , x 0 = [1050.20940]T , and t ∈ R 66 .The execution of the GA in MATLAB is exactly the same as presented in Equations ( 28) and (29).The elapsed time is 1304.814475s, and the matrix obtained by the GA is As before, the linear system with state matrix Equation ( 31) is simulated for 180 s with a sampling time of T = 0.05 s, and the comparisons with the nonlinear states are depicted in Figure 3.The overall mean square error (MSE) between the nonlinear states and their linear approximations is notably smaller: 4.3984 × 10 −4 , but the execution time is larger than in the previous case. .Simulated data with initial condition located away from the equilibrium point.Linearization achieved using a GA using the initial three seconds of nonlinear simulated data, sampled at T = 0.1 s, followed by T = 5 s for the rest of the nonlinear data.

Linearization via Least Squares (LS) Method
In contrasting the potential of the genetic algorithm (GA), a presentation of a simple linear regression utilizing the least squares method is provided.It is widely recognized that perfect identification results when applying this method to data generated by a linear time-invariant (LTI) system without noise or external perturbations.However, complications in identification may arise when applying this method to data from a nonlinear system [7].The method was, nonetheless, applied using the same nonlinear simulation data employed in the GA algorithm.Initially, the analysis was conducted using a limited dataset comprising 31 samples, with a final time of 3 s and samples taken at 0.1 s intervals.The matrix obtained is presented in Equation (32), and the results of the numerical simulation are depicted in Figure 4.The overall MSE between the nonlinear states and their linear approximations is 6.2725 × 10 6 .
Finally, the entire dataset produced by the nonlinear simulation was considered, spanning 180 s with a sample interval of 0.05 s.In this case, the matrix obtained is shown in Equation (33), while the numerical results are illustrated in Figure 5.In this case, the overall MSE is 235.2556.
This analysis is summarized in Table 1.

Conclusions
The traditional linearization algorithm is based on the idealized assumption of perfectly clean and complete data, often considering h = 1 × 10 −12 , rendering it impractical for real-world experimental data.Nonetheless, it has been demonstrated in this study that the genetic algorithm (GA) not only generates matrices with clear physical interpretations but also exhibits the essential robustness necessary for reliable experimental identification.Furthermore, the automation of the linearization process is achievable.
In situations characterized by limited computational resources, where applying GA from scratch to linearize nonlinear dynamics is impeded by unknown search ranges, the GA approach based on the proposed transition matrices proves valuable as a post-processing tool.This is particularly advantageous when measurements of the nonlinear system are available, as demonstrated in this study, showcasing its capacity to enhance results obtained through traditional linearization methods.
Most identification methodologies, including the one presented here, can demonstrate robustness against noise when appropriately designed.With this proposal, assurance is provided that matrix A possesses a physical interpretation, and the presence of nonlinearity does not impede the linear identification process.
In contrast to other methodologies that face challenges in deriving a matrix A with a clear physical interpretation, this approach excels in this regard.Furthermore, the model obtained through this methodology proves functional across a broader operational zone compared to alternative methods.
Considering these aspects, this approach emerges as an alternative for identifying linear time-invariant (LTI) equations from experimental data.The exploration of this potential, particularly in broader experimental scenarios, constitutes an essential avenue for future research.

Figure 1 .
Figure1.Simulated data with an initial condition located away from the equilibrium point.The linearization process, as described in Section 4, is applied to obtain the corresponding linearized representation.

Figure 2 .
Figure 2. Simulated data with initial condition located away from the equilibrium point.Linearization achieved via GA using the first three seconds of nonlinear simulated data, sampled at T = 0.1 s.

Figure 3
Figure3.Simulated data with initial condition located away from the equilibrium point.Linearization achieved using a GA using the initial three seconds of nonlinear simulated data, sampled at T = 0.1 s, followed by T = 5 s for the rest of the nonlinear data.

Figure 4 .
Figure 4. Linearization achieved using linear regression with least squares method using the initial three seconds of nonlinear simulated data, sampled at T = 0.1 s.

Figure 5 .
Figure 5. Linearization accomplished through linear regression employing the least squares method, utilizing the complete dataset of the nonlinear simulated data sampled at T = 0.05 s.

Table 1 .
MSE and execution time.
Matrix for output definition in state form E Jacobian matrix of nonlinear equations with respect to the first-order state variables F b Force vector in body frame F Jacobian matrix of nonlinear equations with respect to the first-order state variables h Step size for numerical derivative, h = x(t k+1 ) − x(t k ) Aircraft velocity from the center of mass to Earth axis, super-index b indicates body frame V e Flight velocity in trim condition v Lateral component for acceleration from equilibrium point, in body frame