Aircraft turbulence and gust identification using simulated in-flight data

Abstract Gust and turbulence events are of primary importance for the analysis of flight incidents, for the design of gust load alleviation systems and for the calculation of loads in the airframe. Gust and turbulence events cannot be measured directly but they can be obtained through direct or optimisation-based methods. In the direct method the discretisation of the Fredholm Integral equation is associated with an ill conditioned matrix. In this work the effects of regularisation methods including Tikhonov regularisation, Truncated Single Value Decomposition (TSVD), Damped Single Value Decomposition (DSVD) and a recently proposed method using cubic B-spline functions are evaluated for aeroelastic gust identification using in flight measured data. The gust identification methods are tested in the detailed aeroelastic model of FFAST and an equivalent low-fidelity aeroelastic model developed by the authors. In addition, the accuracy required in the model for a reliable identification is discussed. Finally, the identification method based on B-spline functions is tested by simultaneously using both low-fidelity and FFAST aeroelastic models so that the response from the FFAST model is used as measurement data and the equivalent low-fidelity model is used in the identification process.


Introduction
During flight, aircraft are expected to encounter atmospheric turbulence events with different levels of severity. Turbulence events not only reduce the comfort of the passengers but also produce additional loads in the airframe and in the case of severe turbulence events, the structure could be subjected to fatigue damage [1]. Turbulence can be considered as the movement of the air through which the aircraft passes. Any component of the air velocity normal to the flight path will change the effective incidence of the aerodynamic surface, the aerodynamic forces will change and consequently the dynamic response of the aircraft will also change [2]. The dynamics loads are significant and aircraft designers must ensure that the airworthiness requirement is met, i.e. aircraft is able to withstand vertical and lateral discrete gusts and turbulence.
In recent years with the development of lighter and more efficient aircraft, the effort to develop new active or passive gust load alleviation (GLA) systems has increased [3]. A large database of real gusts and turbulence events to which aircraft are subjected could lead to a more efficient design of the GLA system. Moreover, air-* Corresponding author.
E-mail address: 996702@swansea.ac.uk (D. Balatti). craft are equipped with digital flight data recorders, a device used to record specific aircraft performance parameters, which is designed to survive under extreme conditions that can result from an accident. In the analysis of flight incidents, gust and turbulence events are of primary importance in the estimation of limit loads [4]. Nowadays, the interest of the aerospace industry in digital twins is growing [5][6][7]. A possible benefit of a digital twin for aircraft is the possibility to calculate the loads in the airframe. To have a reliable calculation of the loads, it is necessary to measure the gusts and the turbulence events. A better estimation of the airframe loads is a major benefit for the aircraft operators. Indeed, knowing the loads at any location of the structure during the flight, or soon after the flight, can avoid unnecessary or expensive inspection of the structure and the inspections could be limited to specific parts of the aircraft [4,8].
The earliest approach for the study of the gust velocity profiles is known as the discrete-gust approach and dates back to the 40s and 50s. It is based on the analysis of peak vertical accelerations directly measured by the aircraft flying in gusts. These accelerations were assumed to originate from a series of isolated discrete gusts and were used to derive gust parameters as the distance parallel to the aeroplane's flight path for the gust to reach its peak velocity (called gradient distances) and the maximum gust velocities [9,10]. This approach was not able to return the real air turbulence, but was sufficient to evaluate normal accelerations on future aircraft designs [11]. In the 60s and 70s, further progress was made on the development of spectral techniques to design aircraft subject to gust encounters, but the complexity in the definition of the frequency response functions was a limitation for this technique [9]. In 1999, in the context of space and missile systems, a Monte-Carlo flight gust loads analysis approach was proposed [12]. The procedure uses forcing functions that were derived by extracting the short-duration, turbulent components of measured wind profiles. However, the method has limited applicability with a restriction to gust wavelengths greater than 500 ft (152 m). In 2009, a method consisting of a model-based approach with an observer for a non-linear aircraft model and a disturbance model for the estimation of gusts and structural loads was proposed [13]. This method uses on-board measured data and parameters available on commercial aircraft making the estimation of manoeuvre-induced structural loads easier. The only unknowns remain the gust velocities which were determined through a non-linear parameter optimisation that computed the gain matrix of the observer model. Recently, a neural network technique for wind gust identification using flight data recordings has been used [14]. In the last 80 years, different techniques for gust loads identification have been studied but mainly two approaches have been considered, i.e. the direct method and the optimisation method [15]. In the direct methods, the excitation f (t) is calculated directly from the measured responses y(t) by evaluating the inverse model. In structural dynamics the process of calculating the excitation force from the measured response is called an inverse problem. The optimisation methods, instead, use directly the model in an optimisation framework where the input is tuned until the model responses match the measured responses [4].
Many inverse problems of science and engineering have the form of a Fredholm integral equation of the first kind which is commonly ill-posed. Different regularisation techniques have been widely used for the numerical solution of this ill-posed problem, such as truncated singular value decomposition, Tikhonov regularisation, basis function expansion and collocation methods [16,17]. The Fredholm integral equation has been solved by expansion methods using different techniques, for example the Coifman wavelet method [18], the Rationalized Haar wavelet method [19], the Sinc-collocation method [20] and collocation methods based on cubic B-spline [21]. Regularisation techniques applied to force identification are mainly divided into two categories, namely frequency-domain and time-domain methods. Frequency-domain methods are based on fast-Fourier transformation (FFT) and they have been widely developed for stationary and pseudo-stationary conditions. However, these methods give poor results for nonstationary and transient responses [22]. On the other hand, timedomain methods show great promise for the study of transient and impulsive phenomena. Methods based on artificial neural networks [23], wavelet decomposition [24], Chebyshev polynomials [25] and cubic B-spline [26] have been used for impact force identification. In force identification, the basis functions used to regularise the problem could also be ill-posed. The inverse problem of force identification can become well-posed if the selected number of basis functions can accurately approximate the impact force [26]. In numerical solution of integral equations, more basis functions give better results. However in force identification, more basis functions could introduce oscillations in the solution [27]. Moreover, it is important to select proper basis functions that can reasonably represent the characteristic of the desired solution. Cubic B-splines have been used in a Newton form to reconstruct the impact force by solving two linear systems of equations [21] and as an efficient regularization method combined with the generalized crossvalidation criteria to identify impact forces [26].
Traditional research in force identification is based on deterministic assumptions. However, uncertainties are inevitable in practice due to material property variation, measurement imperfections, or other factors. Uncertainties can lead to a deterioration of the identification accuracy. The implementation of highprecision identification methods on stochastic structures with distributed dynamic loads is still an unsolved problem [28]. Researchers have considered dynamic loads with a defined probability density function [29], dynamic loads in the context of the probability theory [30] and unknown-but-bounded uncertainty in the structural system [28] [31].
The direct measurement of gusts and turbulence events is difficult, if not impossible, and for this reason, an indirect measurement is necessary. In this work, an extension of the theory for the identification of the impulse response of a dynamic system [26] based on cubic B-splines is applied to aircraft gust identification. In the case of noisy measurements, the comparison between gust reconstruction based on regularised matrices and on cubic B-spline functions is shown. Moreover, the precision required in the model used for the identification to reduce the reconstruction error is discussed. First of all, in Section 2, the general form of the aeroelastic equations used for the generation of the simulated data and for the identification, as well as the different types of gust and turbulence events used are introduced. In Section 3 the regularisation methods for the model inversion and the governing equations of force identification based on cubic B-spline functions are discussed. In Section 4 a numerical study on a low-fidelity system is used to confirm the propriety of the identification methods. After that, numerical studies are performed using a high-fidelity model, known as FFAST in the literature [32], and the identification results based on models with different accuracy are shown in Section 5, before the conclusions are given.

Aeroelastic model
Different fidelity computational models of aircraft may be used in the design process. In the preliminary design of an aircraft simplified models are used to study the general behaviour of the final aircraft and to help the designer select the main parameters. When the main parameters have been selected, detailed models can be created and validated using wind tunnel, ground vibration and flight tests. In this paper, to enable the performance of the identification methods to be assessed, the aeroelastic equations of motion are used not only as a model for the gust identification but also to create simulated data. Two models, with different levels of fidelity, representative of a civil jet aircraft with folding wingtips have been considered. These two models allow the measured data to be simulated with the high fidelity model and the gust to be identified with the low fidelity model, to determine the effects of modelling error. The details of both models are available in the literature, and hence only a brief summary is presented in this paper. The first model is a modified version of the FFAST aeroelastic model [32] whose structure has been modelled using a 'stick' model with lumped masses. Fig. 1 shows the aeroelastic model used for the analyses. The main objective was to investigate the benefits of using a flexible wing-fold device for load alleviation [33]. To avoid numerical errors the connection between the wing and the wingtip was modelled with a very low hinge spring stiffness. The doublet lattice panel method was used for the aerodynamic model [34]. The model is composed of a few thousand elements but exploiting modal reduction and considering only the modes with a frequency below 40 Hz it is possible to reduce the order of the model by retaining 55 modes. All the rigid body motions are constrainted except for the aircraft heave and pitch. Table 1 shows the natural frequencies below 10 Hz of the aeroelastic model at 200 m/s and sea level, see also Fig. 2. The modes related to the wingtip de-   flection and wing bending deflection are the modes that mainly contribute to the gust response as will be shown in the sequel. The second model represents a flexible aircraft and consists of a uniform, untapered, unswept flexible wing with a hinged wingtip, a rigid fuselage and a tailplane (more details of the model are given in [35] and the Appendix). This model has been used to evaluate the system performance and to perform parameter optimisation of an aircraft with an elastic wing and hinged wingtip subjected to gusts. The model has been defined in order to have the least number of degrees of freedom without losing the main features required. The aeroelastic equation of motion has been derived using the Lagrangian formulation and the aerodynamic model was based on strip theory. The model includes five degrees of freedom: the aircraft pitch and heave at the centre of mass, the wing bending and torsional modal coordinates and the wingtip deflection [35]. The total weight, the span dimension, the connection between the wing and the wingtip and the frequency of the elastic modes are the same as the FFAST model. Table 2a shows the elastic modes of the simplified aircraft model with no aerodynamic forces present. Table 2b Table 2a is also close to the wing bending dominant modes of the FFAST model at 2.22 Hz. The rigid body mode of the wing tip deflection of the simplified model is 0.0296 Hz, while the similar mode of the FFAST model is 0.0417 Hz. The extremely simplified version of the FFAST model lacks several details such as coupling between the modes, the effect of swept angle, etc., however, it can reasonably represent the gust response of the detailed FFAST model, as shown in Figs. 5 to 8. The proposed method used in this paper, will use flight measured data from a physical structure and utilize them along with the numerical model to identify the gust. In reality, the numerical models are not accurate and often include modelling errors. In this paper, one may assume that the 'measured' data is simulated using the FFAST model. In section 5.2, it is shown that the simulated measured data of the FFAST model may be identified using the simplified model, which highlights the robustness of the identification method to modelling errors. Regardless the model considered, the aeroelastic problem can be formulated in terms of the physical displacement y or modal coordinates q using the relation is the modal matrix. The aeroelastic equation of motion representative of an aircraft subjected to gust can be expressed as where A, D and E are the structural inertia, damping and stiffness matrices, B and C are the aerodynamic damping and stiffness matrices, q is the vector of the generalized coordinates and q and q are its first and second derivatives, ρ is the air density and V is the aircraft velocity [2]. On the right-hand side of Eq. (2) the forcing term f produces the initial trim and has different contributions (e.g. elevator deflection, incidence at zero lift and gravitational field), f W and f T are the force vectors associated with the gust on the wing and on the tailplane respectively and w g is the gust. The gust acts on the tailplane with a delay t * expressing the time required for the gust to travel from the wing to the tailplane. In this work, the contribution of the initial trim will not be considered; indeed due to the linearity of the system the introduction of a constant term only produces an offset in the results. Eq. (2) can be expressed in the frequency domain, as where ω is the frequency of interest. Eq. (3) can be used to calculate the transfer functions between the gust w g (ω) and the degrees of freedoms q(ω).

Gust and turbulence model
Atmospheric disturbance models are categorized into two idealized categories: discrete gusts and continuous turbulence [36]. In this paper, a typical '1 -cosine' gust disturbance is considered, and the profile is defined as where w g 0 is the maximum gust velocity and l g the gust wavelength. According to EASA regulations [36] for the case of a civil commercial aircraft the gust wavelength is varied between 18 m to 214 m and the gust velocity is calculated as show the wingtip deflection gust responses for the low-fidelity model and the high-fidelity model, respectively. The peak angular deflection of the wingtip is different in the two models because of slight differences in the hinge stiffness. The wingtip deflection has been defined such that a positive angle variation produces an upwards displacement. According, to EASA regulations [36], the power spectral density of atmospheric turbulence is described by the von Karman spectra as where is the spatial frequency, L is the scale of turbulence (commonly assumed to be 2500 ft). According to Hoblit [1] the turbulence velocity time history is obtained as the output of a shape filter with the input given by a stationary Gaussian 'white-noise' where η is the power spectral density of the white noise, τ is the ratio between L and the horizontal velocity of the aircraft and σ w is the component of the gust velocity. Fig. 9 shows the Bode diagram of the approximation of the von Karman turbulence model. The atmospheric turbulence is obtained in the time-domain as the   output of the state-space form of the transfer function Eq. (7) whose input is band-limited (from 0.01 Hz to 10 Hz) white noise.

Identification theory
The direct measurement of many properties of real-world systems is not possible and this information can be deduced from other quantities which may be measured directly [37]. In general, an inverse problem consists of either reconstructing exciting signals acting on a system whose internal characteristics are known, or determining the characteristics of a system driven by controlled or known exciting signals. Fig. 10 shows the difference between direct and inverse problems. In the direct problem, the input, as well as the system, are known, so it is possible to calculate the output (see Fig. 10a). In the inverse problem, the input or the dynamics of the system is unknown and can be estimated by exploiting the known output (see Fig. 10b). Many problems in science and engineering have the form of a Fredholm integral equation of the first kind [37]. This class of problems are typically ill-posed or strongly ill-conditioned after discretization. The generic form of the Fredholm integral equation can be written as Eq. (8) links the unknown function f (t), w g (t) in the aeroelasticity application, over the interval [a, b] to the given kernel function k(x, t) and real-valued function y(x). The kernel function k(x, t) represents the mathematical characteristics of the system and y(x) is the observed data. The inverse problem takes the form of a deconvolution problem where the kernel function satisfies  (9) where τ is the time delayed operation satisfying t ≥ τ and the ker- where y and w g are n × 1 vectors composed of discrete values of the response y(t) and excitation gust w g (t), respectively. H is an n × n transfer matrix with Toeplitz structure composed of discrete values of the impulse response function h(t) as A possible strategy to determine w g when H is known and is non-singular and y is measured is to solve directly Eq. (10), to givẽ where w g is the identified gust. However, this method is very sensitive to the inversion of the transfer matrix due to its large condition number. Inverse problems are typically ill-posed problems, and thus small disturbances in the measured y(t) may result in a large error in the identified w g (t).

Regularisation methods
This section presents a brief overview of the regularisation methods used in this paper, and more information can be found in [17,38,39]. The linear system of equations is ill-posed if the singular values of H decay gradually to zero and the ratio between the largest and the smallest nonzero singular values is large [17]. When the matrix H is ill-conditioned, the problem of Eq. (13) is ill-posed in the sense that a small perturbation of y or H may lead to a large perturbation of the solution. Although different types of direct and iterative regularisation methods exist [17], this work considers three direct regularisation methods: Tikhonov regularisation, Truncated Singular Value Decomposition (TSVD) and Damped Singular Value Decomposition (DSVD).
The general version of Tikhonov's method takes the form where λ is the regularisation parameter defined as a positive constant chosen to control the norm of the solution vector and L can represent the first or second derivative operator but is often the identity matrix [39].
The TSVD defines a new well-posed problem, related to the illposed problem of Eq. (13) and has a solution which is less sensitive to perturbations. The method approximates the matrix H with a lower rank matrix H k . The matrix H k is defined as the rank-k matrix where k < n [17,38].
In the DSVD instead of neglecting n − k singular values, as in TSVD, a smoother cut-off is used by means of filter factors f i defined as These filter factors decay more slowly than the Tikhonov filter factors and thus, in a sense, introduce less filtering. The regularisation parameter λ is a positive constant and plays a similar role to the parameter in Eq. (14), although gust estimates from Eqs. (14) and (18) will be slightly different even if the same value of λ is used.
The selection of the regularisation parameter is a balance between the perturbation error and the regularisation error in the where H I is a matrix which produces the regularised solution w g when multiplied by y, i.e. w g = H I y [17,40].

Identification using cubic B-spline collocation method
Splines are piece-wise polynomials of degree k and continuity C k−1 [41] and they can be used to solve the Fredholm integral equation of the first kind [26,42]. Cubic B-splines are piecewise polynomials of degree three with C 2 continuity at the junction points between adjacent segments. This section only presents a general overview of cubic B-splines, and more detail of the method used can be found in reference [26]. In this work, cubic B-splines with uniformly distributed knots are considered. Let interval can be defined in a piece-wise polynomial function form [42,26]. Fig. 11 illustrates one cubic B-spline, which is composed of four cubic polynomials. The cubic B-spline functions S(t) is constructed as a weighted sum of m + 2 cubic B-spline basis functions B i (t), namely where the matrix A + is the Moore-Penrose generalized inverse matrix. Finally, the calculation of the unknown vector w g can be reformulated in matrix vector notation as The number of collocation points governs the degree of expansion of the cubic B-spline functions employed for to approximate the unknown force. The residual is used to determine the optimal number of collocation points, by minimising r = ||y − Hw g || 2 (26)

Identification results for the low-fidelity model
In this section, the identification results based on the lowfidelity model are presented. The measurement data are simulated by using Eq. (2) and the identification is based on the corresponding model in the frequency domain, Eq. (3). The response is obtained as a time history of 1000 equally spaced samples representative of 20 seconds. This section analysis compares the results of gust identification based on the regularisation with the one obtained by approximating the gust as a summation of cubic B-spline functions. In addition, for the identification based on cubic B-spline functions, the effects of collocation points location and the ability to identify atmospheric turbulence are investigated. The results assume that the degrees of freedom of the model in modal coordinate may be measured directly (e.g. the bending or the torsional modal coordinate of the wing). In practice, the modal coordinates would be estimated from the physical coordinates using Eq. (1). Nevertheless, this is not a limitation of the method and similar results can be obtained using different measurements.

Gust identification considering measurement noise
An important feature of a reliable gust identification technique is the ability to reconstruct the gust in the presence of noisy data, since in the measurement process, measurement noise cannot be avoided. In this work, the noisy measurements data (ŷ) are created by summing the measurements and the noise aŝ  where the MATLAB script function std(•) denotes the standard deviation of the vector, the MATLAB script function rand(n,1) returns an n × 1 vector containing uniformly distributed random numbers on the interval (0, 1) and l noise is the current noise level of the simulated response. In the results 10%, 20% and 30% of measurement noise correspond to l noise = 0.1, 0.2 and 0.3, respectively. Fig. 12 shows the gust reconstruction considering the model inversion technique of Eq. (12) and the measurement of the bending mode in the case of clean data and 10% measurement noise. The reconstruction in the case of noisy data has a large error, while the reconstruction from the clean data gives good results. In this case, the order of magnitude of the condition number of the matrix H is 10 19 ; so the problem is ill conditioned and a regularisation method is required. Figs. 13, 14 and 15 show the GCV function considering the noisy data, in the case of Tikhonov (TIKH) regularisation, TSVD regularisation and DSVD regularisation, respectively. The GCV function considering clean measurement, in the case of the DSVD and Tikhonov regularisation methods, is a monotonically increasing function and in the case of the TSVD method is a monotonically decreasing function. On the contrary, in the case of noisy measurement the GCV function has a minimum and the lowest GCV achievable increases with the increase of the noise. Figs. 16 and 17 show the gust reconstruction from the model inversion method considering different regularisation techniques, in the case of 10% and 30% of measurement noise, respectively. The results have been obtained by setting the regularisation parameter to the value associated with the minimum GCV. The comparison between the gust identification from noisy data of Figs. 12 and 16 show that the DSVD regularisation method is able to reduce the error in the   reconstruction of the gust. Moreover, Tikhonov regularisation and TSVD regularisation can further reduce the reconstruction error. Figs. 18 and 19 show the results of the gust identification based on cubic B-spline functions. Fig. 18 shows the residual in the cases of clean and noisy measurement from the bending mode. It shows that for low numbers of cubic B-splines and low noise levels, the noise does not affect the reconstruction of the gust. Moreover, increasing the noise level the lowest residual value achievable increases. Fig. 19 shows the gust reconstruction based on cubic B-spline in the case of 10% and 30% of measurement noise. The reconstruction based on 73 collocation points gives good results for both 10% and 30% of measurement noise. The reconstructions for 30% measurement noise in Fig. 19 confirm the trend of the  residual in Fig. 18. Indeed, the reconstruction based on 200 collocation points is less accurate than the one based on 73 collocation points. Thus, Fig. 19 shows that the cubic B-spline identification acts as a filter and decreasing the distance between collocation points the identification never converges to a smooth result. These gust identification techniques are also tested using different levels of coloured noise using 'pinknoise' command in MATLAB and similar results were obtained.
Figs. 16, 17 and 19 show that the matrix regularisation through TSVD and the approximation of the solution through cubic B-spline functions give similar results. Indeed, in the cubic B-spline method the number of cubic B-splines plays a similar role to the regularisation parameter [26]. The B-spline method has been found to perform better than the regularisation methods and therefore the effects of the location of the collocation points will be shown in the next section. Moreover, turbulence identification using cubic Bspline functions will be used for gust identification for the FFAST model.

Effect of the location of the collocation points
In the identification through B-spline functions, a fundamental aspect is related to the choice of the number of collocation points. Gusts defined by EASA can have different frequencies and amplitudes and the optimal choice of the number of collocation points for one gust could lead to some error in the identification of a different gust [36]. Fig. 20 shows the trend of the residual in the interval 1 to 200 cubic B-spline considering as measurement the bending modal coordinate. In the considered interval the residual is not a monotonic decreasing function but has local minima at

Gust and turbulence event identification
EASA regulation requires to consider the response of the aircraft to discrete gusts and continuous turbulence event separately [36]. In this section, the identification based on cubic B-spline functions in the case that the gust and the turbulence event are combined together is considered. For all the time histories the turbulence is present and the gust acts after 5 seconds of the simulation.    of the maximum peak is 0.9% and the error in the identification of the turbulence field is in the interval ±0.25 m/s.

Identification using measurements from the FFAST model
In this section, the results of the identification based on cubic B-spline functions of the FFAST model are presented. The simulation is performed for 100 seconds and is composed of 4000 equally spaced points. In all the cases, the measurement data are generated considering the 55 modes model. The time response from 5 seconds to 100 seconds is considered for identification and the rotation of the centre of gravity of the aircraft in pitch is used as a measurement because this information is typically available on the aircraft. The first 5 seconds of the measurement are not con-  sidered for the identification because, due to the complexity of the model and the numerical errors without any gusts, the response of the system is not exactly zero initially. Moreover, for the identification process, the impulse response function assumes zero response and zero first derivative of the response at the initial time. In the B-spline function, it would be possible to introduce controlled end conditions in order to reduce the error at the extremities of the identification [43]. In this work this has not been done because it is always possible to change the initial and final time of the identification. In any case controlled end conditions will alleviate but not cancel the identification errors at the extremities. The models used for the identification are based on two approaches; in the first case, the model is obtained through a modal reduction and in the second case a low-fidelity model is used.

Identification model: modal reduction
The definition of the model through modal coordinates allows us to introduce or exclude modes in the model in order to have different levels of accuracy. The identification is performed considering that the model used for the identification has an increasing number of modes. Fig. 26 shows the residual when the model used for the identification has 4 modes, 5 modes and 55 modes. The gust identification based on the 4 modes model is not able to converge to the correct result, but the identification based on the 5 modes model gives results similar to the identification based on 55 modes. Fig. 27 shows the identification results considering the 5 modes model, the 8 modes model and the 55 modes model for 600 collocation points. Fig. 28 shows the difference between the real gust and the reconstructed gusts of Fig. 27. As expected,  decreasing the accuracy in the model increases the error in the reconstruction. The error in the first seconds of identification is related to the aforementioned problem of the initial conditions and disappears after two seconds of identification. This error can be mitigated by increasing the time of the simulation. Moreover, after the initial condition, the identification error using 55 modes is in the interval ±0.2 m/s and it can be reduced further by increasing the number of B-splines. The identification based on the 5 and 8 modes models has the maximum error during and soon after the peak of the gust where the higher frequency modes are exited. The error in the identification of the peak of the gust is 25% for the 5 modes model and 4% for the 8 modes model. Moreover, Fig. 28 shows that in the second part of the identification where the contribution is only from the turbulence event, the three models give similar results. Fig. 29 reports the gust and turbulence event response in the frequency domain obtained using the 5, 6, 7, 8 and 55 modes models. The response of the 55 modes model has a higher amplitude at low frequencies and then the amplitude decreases for higher frequencies. The good estimation of the gust and turbulence event with the low order model can be explained considering that they are obtained through a modal reduction, so considering only a subset of the low-frequency modes. Indeed, in Fig. 29 at low frequency, the response of the reduced order models is similar to that of the complete model.

Identification model: realistic example, simulating measured data using FFAST and using low fidelity model for identification
In this section, the low-fidelity model is used for identification purposes while the higher fidelity model is utilized to simulate measurement data. Fig. 30 shows the residual considering the measurement of the pitch of the aircraft. Fig. 31 shows the identification based on 450 and 600 collocation points. The comparison between Figs. 27 and 31 show that the identification based on the low-fidelity model gives similar results to the identification based on a reduced order model. It is worth highlighting that the two models have been defined using different techniques and the aerodynamic models are different; thus the identification method is robust with respect to modelling errors. Fig. 32 shows the reconstruction error in the case of identification based on 450 and 600 collocation points. The error in the identification of the peak is 8% in the case of 450 collocation points and 3% in the case of 600 collocation points. The positive and negative peaks of the gust identification error are related to a time delay in the identification of the peak of the gust. Indeed, Fig. 31 shows in the identification of the peak a time delay of 0.1 seconds in the case of identification considering 450 collocation points and 0.07 seconds in the case of identification considering 600 collocation points. This is due to fact that the low-fidelity model reacts faster than the high-fidelity model. However, the faster reaction of the low-fidelity model does not affect the turbulence event reconstruction that has an error in the interval of ±0.5 m/s.

Conclusion
The aim of this work is to demonstrate a robust technique for aircraft gust identification based on cubic B-splines. Two aeroelastic models have been used, a low-fidelity model and a higher fidelity model. The results found in the case of the low-fidelity model show that the approximation of the gust through cubic Bspline functions gives similar results to those obtained through the matrix regularisation. Moreover, the identification with cubic Bspline functions is able to identify the gust with a small error even when selecting a non-optimal number of cubic B-splines or considering noisy measurements. The identification in a realistic scenario, where there is a small turbulence event and a gust, gives promising results. Indeed, the turbulence is well captured except for the components at high frequency, because the selection of the distance between consecutive collocation points limits the maximum identifiable frequency. This is not a limitation since the contribution to the response given by the higher frequencies is minimal due to the filtering action played by the structure.
The high-fidelity model has been used to generate data and the identification by cubic B-splines has been performed considering reduced-order models of the high-fidelity model and the low-fidelity model. The results have shown that it is possible to reconstruct the gust using models with a small number of degrees of freedom. The low-fidelity model is able to identify the main peak of the gust with better accuracy than the reduced-order model with the same number of degrees of freedom. However, in the reduced-order model, the introduction of additional degrees of freedom can increase the accuracy in the identification. The identification of the turbulence event gives similar results considering different reduced-order models and, better results can be obtained considering reduced-order models than the low-fidelity model.
To validate the proposed method, it is essential to have a direct measurement of the gust. Hence, a wind tunnel test campaign is planned in the future; after the gust generator is characterised, the response of an elastic wing subjected to known gusts will be measured. Furthermore, the identification method proposed is based on the assumption that the system is linear. However, the method proposed can be extended to nonlinear systems and the identification based on the linear system can be used as the starting point of an iterative process.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A.1. Low-fidelity model
In this section the low-fidelity aeroelastic model used in Sections 4 and 5 is defined. This model is based on the model published in [35]. The structural and aerodynamic models have been selected to obtain a model with the least number of degrees of freedom without losing the main aeroelastic characteristics. The only deformable parts of the model are the wing, in bending and in torsion, and the wingtip. Additional degrees of freedom could have secondary effects in the study of gust load alleviation. The aerodynamic model is based on strip theory. Although the panel method could provide more accurate results it would lead to further aerodynamic degrees of freedom [44]. In the aeroelastic model the following degrees of freedoms are considered: displacement z c (downwards) and pitch α (nose up) at the centre of mass, the torsional q t (nose up) and bending q b (downwards) modes of the wing [2] and the wingtip rotation θ . θ is defined such that a positive angle variation produces a upward displacement. Fig. 33 shows the model. The wingtip's span is assumed to be 20% of the half span and there is a relative angle (γ ) between the hinge axis and the free stream velocity. It is worth noting that in the region of s 2 the chord is not constant, but is a function of the flare angle γ , of the elastic wing span s and of the span position y. The geometric relation is The model of the symmetric aircraft has been obtained through the Lagrangian formulation. The simplified flexible aircraft consists of a uniform, untapered, unswept flexible wing of chord c and semi-span s, plus a rigid fuselage and tailplane, as shown in Fig. 33. The wing is assumed to have a uniform mass distribution and the wing mass axis (WM) lies at distance l W M ahead of the aircraft centre of mass.
The mass and pitch moment of inertia of the aircraft fuselage will be represented by discretization into three 'lumps' of mass m F , m C and m T . These discrete masses are located, respectively, at the front fuselage (at distance l F forward of the CM), at the whole aircraft CM and at the tailplane aerodynamic centre (at distance l T aft of the CM). The wingtip is considered as a rigid body of mass m wt with the centre of mass at = ( x , y ) defined in a reference system with the origin at the leading edge of the elastic wing tip and with the x-axis parallel to the hinge axis, as shown in Fig. 33.
The wing elastic axis (WE) is assumed to lie at distance l E ahead of the WM axis. The wing aerodynamic axis (WA) is at the wing quarter chord and is at distance l W ahead of the centre of mass and at distance l A ahead of the elastic axis. In order to minimize any coupling between the rigid body modes and the flexible mode equations, the mean axis reference frame has been used [2].
The displacement z W A (y, t) (downwards positive) of the elastic wing aerodynamic axis is (29) where A and B are unknown constants defining the amount of bending and twist present along the wing span and k e 0 and γ e 0 are constants defining the wing root displacement and twist deformation. Considering small rotation of the wingtip and assuming the aerodynamic centre of the wingtip is at the quarter chord and halfway along the wingspan its vertical displacement is The wingtip centre of mass is at position = ( x , y ), and so its vertical displacement is (31) where x f is the longitudinal position of the elastic axis measured from the wing leading edge. The displacement z T (t) (downwards positive) of the tailplane aerodynamic centre is The aerodynamic terms due to the wing and the tailplane have to be determined. To this end, the tailplane and the wingtip are considered as rigid, while the elastic wing contribution involves integration using a strip dy because of the flexibility. The lift of a strip dy at the position y along the wing span is where α 0 is the incidence for zero wing lift and a w is the sectional wing lift curve slope. There is also a zero lift pitching moment for the wing The wingtip lift is given by where the contribution θ sin(γ ) is the component of the rotation around the hinge perpendicular to the free air-stream. The tailplane lift considering the contribution of the downwash k , the effective incidence due to the nose up pitch rate and the increment of lift due to a rigid vertical displacement, is (36) where a E is the tailplane curve slope defined with respect to the elevator angle and η is the elevator angle and has been included to provide trim. The effect of the vertical gust on the aerodynamics is a change of angle of attack. Thus on the elastic wing the increment of lift on a strip dy is given by while on the wingtip is On the tailplane the gust will act with a delay given by the ratio between the distance between the wing aerodynamic centre and the tailplane aerodynamic centre and the free stream velocity (t * = l W +l T V ). Hence The kinetic energy due to the rigid motion and the dynamic motion is T = where m is the total mass of the aircraft, I y is the aircraft pitching moment at the centre of mass, and m b and m t are respectively the bending and torsional modal masses. The wingtip gives a negligible contribution to the total inertia, although an important contribution to the local inertia of the wingtip The elastic potential energy corresponds to the strain energy in bending, torsion, and the spring at the hinge between the elastic wing and the wingtip, such that Finally, the virtual work done by lift forces and moment and the gravitational field force is Exploiting the Lagrange formulation the full aeroelastic equation in the classical second order form is obtained as (43) where A, D and E are the structural inertia, damping and stiffness matrices, B and C are the aerodynamic damping and stiffness matrices, q is the vector of generalised coordinates and the force vector f on the right hand side of the equation has contributions due to the elevator (f η ), zero incidence (f 0 ), gravitational field (f g ), gust on the wing (f W g ) and gust on the tailplane (f T g ).

A.2. Numerical data used
The dimensions and total weight have been estimated from the FFAST aeroelastic model [32] of a representative civil jet aircraft whose structure was modelled using a 'stick' model with lumped masses. The weight distribution, the dimensions and the main parameters are reported in Tables 3 and 4. The effect of engine mass is also considered in this model. The engine is modelled as a lumped mass m M located at the longitudinal distance x M from the elastic axis (positive aft) and y M from the symmetric axis as reported in Table 5. The coefficients A, B, γ e0 , k e0 , and the bending and the torsion modal masses have been obtained through a modified minimization process in order to consider the effects of the engine [2]. The bending and torsional modal stiffnesses are determined to give a bending modal frequency of 2.5 Hz and a torsion modal frequency of 4.5 Hz.  Table 5 Civil commercial aircraft engine parameters.

Mass
1680 kg x position from EA 0.0 m y position from fuselage 9.344 m