New developments for experimental modal analysis of aircraft structures

This article presents an identification algorithm dedicated to the modal analysis of aircraft structures during flighttests. More specifically, this algorithm was designed to process short duration tests carried out with multi-input excitations. The identification problem is solved in the frequency domain and the limit effects are considered so as to avoid transient effects with short data sequences. To minimise the effects of the noise, a non-linear gradient-based optimisation method is used. Its performance is improved by the use of an appropriate over-parametrised matrix fraction descriptions. Because the cost function to be minimised is non-convex, this method is however sensitive to the initialisation. For this reason, an iterative instrumental variable method is used to find an initial estimate. This one gives a value of the cost-function sufficiently close to its global minimum so as to ensure a fast convergence of the optimisation. Thus, the algorithm presented in this article is a combination of two iterative methods that gives accurate mode estimations even with high level of noise, as shown on an illustrative example.


Introduction
The objective of a modal analysis is to identify the vibration modes of a structure.Because a structure can usually be modelled as a linear time invariant (LTI) system [1], a modal analysis can rely on LTI identification methods.Like any other identification application, a modal analysis is based on experimental data recorded during specific tests.The duration of these tests can have a major impact on the cost of the analysis.This is particularly true for the in-flight modal analysis of the aeroelastic modes of an aircraft structure.In order to reduce the test campaign durations, the modal analysis of civilian aircraft now mostly relies on short pulse tests.
Such tests imply several difficulties.First, the length of the input is too short [2] to use the conventional methods based on periodogram averaging for estimating transfer functions [3].Moreover, the energy injected into the structure during a short experiment is quite low.One solution to increase the energy of the excitation while preserving its brevity is to excite the structure at several points simultaneously.Thus, multi-input identification methods are relevant for processing short experiments.However, the noisy conditions still imply low signal to noise ratios.Although rather uncommon in modal analysis, an iterative approach is therefore well indicated so as to minimise the effects of the noise.
These specificities implied that existing modal analysis methods in the literature, usually developed for ground test conditions were not suitable [4].This motivated the development of a new identification algorithm.This one is based on a frequency domain formulation of the identification problem.Moreover, it relies on transfer a e-mail: jeremy.vayssettes@isae.frb e-mail: guillaume.mercere@univ-poitiers.fr functions modelled by over-parametrised matrix fraction descriptions.A gradient-based optimisation method is used to iteratively solve the formulated problem.The optimisation is initialised by using an iterative instrumental variable method.A first version of this algorithm was suggested in [5].The actual algorithm is the result of several developments carried on from this first version so as to improve its performance [5,6].
This article aims at presenting the final modal analysis algorithm and is organised as follows.Section 2 is dedicated to the identification problem formulation.The frequency domain problem is first recalled and a solution used to avoid transient effects with short data sequences is detailed.Then, a new parametrisation of matrix fraction descriptions is introduced in this section.This one has a maximum number of parameters and improves the performance of gradient-based optimisation methods.In Sect.3, the optimisation method and the iterative instrumental variable method used to find an initial estimate are detailed.Section 4 is dedicated to the evaluation of the final algorithm on an example representative of flight test conditions.In this section, the performance improvement thanks to the algorithm enhancements achieved since its first version introduced in [5] is discussed.Finally, Section 5 concludes this work.

Problem formulation
In the rest of the article, the integers n u , n y and n x designate the number of inputs and outputs of the system and its order, respectively.This section aims at defining the formulation of the system identification problem that must be solved in order to estimate the modal parameters.

Identification problem
The developed modal analysis algorithm is based on the use of an output-error identification method in the frequency domain.The goal of any output-error identification method is to find the parameters θ that minimise the energy of the error between the measured system outputs y(t) and the outputs y(t, θ) computed with a parametric model.Assuming that the considered system is causal, this energy is defined as where .denotes the 2-norm [7].The model output y(t, θ) is defined by the convolution between the impulse response h(t, θ) of the estimated model and the system inputs u(t) as From Parseval's theorem, this energy does not depend on the chosen domain and may be defined similarly in the frequency domain as where Y ( f ) and Y ( f, θ) are the Fourier transforms of y(t) and y(t, θ) respectively computed at the frequency f .From the properties of the Fourier transform, the values Y ( f, θ) are defined as where U ( f ) and H ( f, θ) are the Fourier transform of u(t) and h(t, θ), respectively.In practice, only the frequencies f in a set of interest, denoted F, are considered.Hence, using standard notations, the output-error identification problem consists in finding the values of the parameters θ that minimise the weighted quadratic cost function where W is a user-defined positive definite weighting matrix.

Integration of the limit effects with short-data sequences
The Fourier transform computation is based on the integration of the signal over the entire (infinite) time scale.Obviously, in practical experiments, the signals are only measured on a limited time interval [t 0 , t 1 ].If the signals are not equal at the limit instant, Eq. ( 4) is not true any more.As first discussed in [8] for the case of SISO discrete-time transfer functions, it is necessary to integrate corrective terms in the Fourier analysis so as to avoid a biased estimation of the system parameters.Recently, we suggested a similar solution in [6] for the identification of MIMO continuous-time transfer functions.In the latter, it was shown that considering the limit effects explicitly in the identification problem implies to solving a similar identification problem as in Eq. ( 5).Indeed, this just involves the identification of a system with an additional fictitious input, noted V ( f ).This can be done by minimising the following identification criterion where gathers the system parameters θ plus the n x + n y additional parameters due to the limit effects integration, H( f, ) denotes the increased model of dimension Moreover, the numerical integration of the Fourier transform was discussed in [6].We showed that choosing provides an efficient numerical integration of the Fourier transform values.For more details, the interested reader can refer to [6,9].The integration of limit effects is of significant importance if one wants to consider short data sequences as illustrated in Fig. 1.This figure shows the response of a low-damped system excited with a short duration input.For this example, a pulse excitation of duration 0.2 seconds was used at the instant t 0 = 10 seconds.The time needed by the system to come back to its initial state can be long when the system exhibits very low damped modes.In such a situation, considering the limit effects is necessary in order to reduce the signal length considered for the identification without introducing parameter bias.Being able of decreasing the length of the considered data sequences can be suitable for two reasons.First, it enables to decrease the amount of data and consequently the system identification processing time, especially when high dimensional systems are studied.Second, this is an efficient way of improving the signal to noise ratio, and consequently the accuracy of the system identification [6].

New model parametrisation of transfer functions
Since the limit effects are considered, the problem ( 6) is well posed whatever the system states are at the limit instant considered.In order to solve this problem, a system description has to be chosen.In this article we consider 01001-p.2AVE2014 only left matrix fraction descriptions.Thus, H(s) is given as where D(s) and N(s) are polynomial matrices of dimension n y × n y and n y × (n u + 1) respectively.The vector ∈ R n then gathers the polynomial coefficients of D(s) and N(s).
We need however to derive a parametrisation of H(s), i.e., a degree structure for the numerator and for the denominator so as to satisfy suitable requirements for a system identification.First, the parametrisation should ensure that the order of the identified transfer is n x whatever the identified parameters.This requirement has to be satisfied when a fixed order system identification is carried out.Second, any real system has a bounded frequency response, i.e., its transfer function is proper [10].Thus, the parametrisation should ensure that the identified transfer function is always proper.Finally, in order to get a smooth functioning of the identification algorithms, every transfer function of order n x should be representable by the chosen parametrisation.
Finding such a parametrisation is quite easy and straightforward for SISO transfer functions.Indeed, because the denominator is scalar, the order is directly fixed by the degree of the denominator.Then, ensuring that the model is proper is achieved by setting the numerator degree equal to the order n x .However, things are much more complicated when polynomial matrices have to be dealt with, i.e., when MIMO systems are considered.A large number of studies were carried out, mostly before the 90's, to define appropriate parametrisations of MIMO systems [see e.g.[11][12][13][14][15][16].At this time, the authors pointed out that pseudo-canonical parametrisations are well-suited for a system identification.
So, from the system theory point of view, the pseudocanonical parametrisations seemed to be the best choice for system identification.However, optimisation methods are known to exhibit poor convergence properties when they are used to minimise an identification problem which relies on a pseudo-canonical parametrisation [17].More specifically, we showed in a recent study that the pseudo-canonical LMFD can sometimes lead to a numerical locking of the convergence.Such phenomena are caused by the constraint imposed on the highestdegree coefficients of the denominator [5].Relaxing these constraints was suggested in [5,6] in order to avoid this drawback, thus leading to a new parametrisation with n y over-parameters.
However, a more generic over-parametrised form can be found for LMFD.It can be shown that the parametrisation with the largest number of parameters, that still satisfies the three aforementioned essential requirements for the identification, is specified by the following degree structure [9] deg where ρ, µ 1 and µ 2 satisfy the relations Compared to the pseudo-canonical form used in [5], this new parametrisation has n 2 y additional parameters.As shown for state-space representation in [17], using a representation with a higher number of parameters improves the performance of gradient-based optimisation methods when this choice is combined with search dimension reductions at each iteration of the optimisation.As we will see in Sect.3, these dimension reductions can be carried out by using a singular value decomposition [7] of the Jacobian matrix.Hence, the new parametrisation given in Eq. ( 9) improves the performance of optimisation methods for the identification of matrix fraction descriptions.Compared to the parametrisations used in [5,6], the parametrisation (9) has also the great advantage to simplify the implementation work during the development phase of a new algorithm due to the simplified degree structure of the denominator.

Solving the identification problem
In this section the iterative algorithm used to solve the identification problem ( 6) is detailed.The parameter vector calculated at iteration k is denoted k .

Gradient-based optimisation method
Gradient-based optimisation methods are known for their good convergence properties in the neighbourhood of the global minimum of the criterion function to be minimised [18].The Gauss-Newton method is one of the most popular gradient-based method.At the current iteration k, this method consists in replacing H( f, k ) in Eq. ( 6) by its first-order development about k−1 .For the model description (6), it can be shown that the first order term about k−1 is given by [6] ∂ H ∂ where k = k − k−1 .By using this expression, the Gauss-Newton criterion is with From the nonlinear identification criterion in Eq. ( 6), it can be noted that the Gauss-Newton method enables to get a linear cost function at each iteration with respect to k .Consequently, minimising Eq. ( 12) reduces to solve a least squares problem that can be written as 01001-p.3

MATEC Web of Conferences
where M gn is a non-square matrix constructed from the derivative terms (11) and Z is a vector that gathers the terms Y f for every frequencies.
Because the chosen description H( f, ) has n 2 y overparameters, the matrix M gn has always n 2 y singularities.Using the singular value decomposition (SVD) [7] to compute the Moore-Penrose pseudo-inverse of M gn enables to consider only the full rank part of the matrix M gn .In other word, the SVD enables to reduce the search space directions to the directions orthogonal to the kernel space of M gn only.With this approach, the new parameter vector calculated at each iteration is where α is a scalar value found by using a line-search technique [18].The vector k is the descent direction equal to where M # gn is the Moore-Penrose pseudo-inverse of M gn .The Gauss-Newton method is a powerful non-linear optimisation method that can provide accurate parameter estimations even in the presence of large amount of output noise.However, the minimised function is not convex.As a consequence, this method is sensitive to the initialisation and must be initialised in the neighbourhood of the minimum of the criterion so as to converge within a few iterations.

Initialisation with an iterative instrumental variable method
In order to find a good initialisation for the Gauss-Newton method, the use of an iterative instrumental variable method (IV) was suggested in [6].It was shown that the IV method converges within a few iterations to an estimation close enough to the optimal parameter value so as to ensure a fast optimisation by the Gauss-Newton method.To find an estimate that locally minimises J( ), the IV method consists in finding the parameters for which From the expresion of J( ) in Eq. ( 6), it can be shown that the IV method consists in solving [19]: where ∂ and (•) H denotes the Hermitian transpose.The left hand side of Eq. ( 17) being non-linear in , the IV method minimises iteratively the following linear problem with It can be seen that the values m gn ( f ) are the row vectors of the Gauss-Newton matrix M gn .Similarly, it can be shown that the row vectors m sk ( f ) are the rows of the matrix M sk constructed when the Sanathanan-Koerner algorithm is written [6].This algorithm indeed linearises the output error problem by approximating the denominator factor by the value estimated at the previous iteration.Thus, the linear problem in Eq. ( 18) can be written in a matrix form as As mentioned before, the matrix M gn is not a full rank matrix if the number of parameters is not minimal.The matrix M sk being always a full rank matrix, this implies that the product M gn M sk is a full rank matrix when the considered parametrisation is minimal.From Eq. ( 20), this implies to fix the n 2 y over-parameters of the parametrisation defined in Eq. ( 9) and to consider only the parts M gn (:, S) and M sk (:, S) calculated for the non-fixed parameters (S).The fixed coefficients of the parameter vector are denoted (F).By doing this, we have shown in [6] that the convergence scheme can be written in a similar way as the Gauss-Newton one as with

Application to the estimation of aeroelastic modes of aircraft structures
The algorithm detailed in this article has been developed in order to estimate aerolestic modes of aircrafts during flutter flight tests performed from short-duration input excitations.The main advantage is that the tests carried out are shorter than tests performed with sine-sweep or multi-step inputs more classically used in modal analysis.However, these inputs excite the structure with less energy.In order to provide more energy, multi-inputs excitations must be used.Without an appropriate MIMO modal analysis algorithm, such tests could not be carried out during the flight tests with a sufficient safety level.For this reason, no real data coming from such tests are available yet.However the performance of the developed algorithm was evaluated by mean of a simulation benchmark that was developed so as to be representative of the flutter flight-test conditions.

Description of the test case
The model used in this benchmark is derived from the aeroelastic model of a real aircraft.This model does not represent the entire aircraft but includes several difficulties encountered during flutter flight tests.The model has 3 inputs, 6 outputs, 5 aeroelastic and 1 one aperiodic modes.The model also has two perturbation inputs used to simulate the background non-white noise.To be more realistic, these inputs enable for the simulation of a noise 01001-p.4spectrum shaped by the aircraft model dynamics.During flight tests, the level of noise is not constant and depends on the flight conditions.Two level of noises were considered: a high level that represents unfavourable test conditions and a lower level that represents better conditions.Both noise levels were set according to noise levels measured on real dataset.For each level of noise, a Monte-Carlo simulation was performed with 100 noise sequences.The mean values of the signal to noise ratio (SNR) are 9.7 dB and 14.3 dB, respectively.The excitations that were used are 0.2 seconds pulses successively applied on the three system inputs.The data recorded from the start of the first pulse until 7.5 seconds after the end of the last pulse were used for the identification.

Identification results
The modes identified with the algorithm detailed in this article are shown in Fig. 2. For each aeroelastic mode, the identification rate (IR) is shown, i.e., the number of times the identified mode can be paired with the true one over the 100 runs of the Monte-Carlo simulation.This association is based on an enhanced formulation of the modal assurance criterion (MAC) [20].Moreover, the root mean square (RMS) error of the relative damping value and of the frequency value of each mode is also given in Fig. 2.These statistical results show that all the modes are almost always identified under normal flight test conditions.Even the two very close modes are well distinguished.Furthemore, the frequencies and dampings are estimated accurately.The results for the high level of noise shows that the algorithm still provides good results under very unfavourable test conditions.Hence, all the modes are identified in more than 50 percent of the cases with still accurate estimations of the damping and frequency values.

Achieved improvement for the estimation of the modes
The impact of the integration of limit effects, and the convergent schemes of the Gauss-Newton and IV methods have been studied in [6].An analysis of the influence of the over-parametrisation of LMFD can be found in [9].Table 1 and Table 2 show the performance improvement that is achieved with this algorithm compared to the more conventional version that does not include the developments detailed in this article, i.e., the limit effects consideration, the new over-parametrised LMFD and the new convergent schemes of the both Gauss-NEwton and IV methods.It can be seen that the identification rates of the modes are improved.Moreover, the relative damping errors and the frequency errors of the modes are decreased resulting in a more accurate estimation of the modes.More specifically, for the high level of noise, the identification rate as well as the damping and frequency accuracy of 4 modes is improved by more than 10%.Hence, it must be emphasised that the algorithm presented in this article enables to get good results with such a high level of noise whereas the more conventional versions of the Gauss-Newton and IV methods could not handle it.

Conclusion
This article presents a new modal analysis algorithm specifically developed in order to perform analysis from short duration tests.This algorithm is the result of several innovative developments.First, it is based on a frequency domain identification problem formulation that considers the limit effects.Thanks to this, the transient effects are avoided.Thus, the data sequences can be shorten which is crucial to increase the signal to noise ratio when short duration inputs are used.Second, the algorithm is based on a new over-parametrisation of left matrix fraction descriptions.This parametrisation improves the convergence of optimisation methods.Third, 01001-p.5 MATEC Web of Conferences a combination of an iterative instrumental variable method and of a gradient-based optimisation method is used to minimise the identification criterion.New convergent schemes have been suggested for these two methods that provide a fast and accurate convergence.The efficiency of the algorithm is shown on an example representative of the modal analysis of an aircraft structure during flight-tests.It was shown that accurate mode estimations are obtained even with a high level of noise.It was also shown that this efficiency is achieved thanks to the developments presented in this article.

Figure 1 .
Figure 1.Truncation of a low-damped system response to a short excitation.

Figure 2 .
Figure 2. Modal analysis results with the developed algorithm for the low (left) and high (right) level of noise.

Table 1 .
Improvement gains provided by the new developments for the lower level of noise.

Table 2 .
Improvement gains provided by the new developments for the high level of noise.