Real time estimation of modal parameters and their quality assessment

In this paper the recursive method for modal parameters estimation is formulated and verified. Formulated algorithms are implemented in the FPGA electronic chip. As a result, the modal parameters and confidence bounds for the modal parameters are obtained in real time. The algorithms and their implementations are tested on laboratory test rig data and applied to – flight modal analysis of an airframe structure.


Introduction
Flutter can be a problem with any aircraft design and design engineers prefer to predict and modify the critical flutter speed of an aircraft before an incident occurs.The first step in a flutter analysis is performed by setting up a very detailed finite element (FE) model of the structure to find the eigenvalues and mode shapes.The FE model has to be validated by a Ground Vibration Test.The mode shapes from the FE are entered into the flutter analysis software to determine the critical flutter speeds.Balanced control surfaces, altitude and fuel load are some of the parameters that influence flutter speeds.The second and necessary step in the flutter test is confirmation of analytical results during a flight.The test flights are very expensive and time consuming.New methods of online flutter detection are intensively investigated to preserve test pilot safety (this is extremely important because flutter is an unstable phenomena) and to minimize test flight costs.Unstable vibrations of an airplane can be a reason of a catastrophic failure of the aircraft [1].In the literature [2][3][4] many cases of flutter phenomena are carefully studied.The procedure of in-flight flutter testing consists of measurements of structural vibration of an airplane and, based on these measurements, an estimation of modal model parameters [5,6].Each possible aircraft configuration should be tested separately.There are two cases: first, if artificial vibration excitation is used during a flight and second, if ambient excitation (none-measured) is a cause of structural vibration.Many papers are focused on development of operational modal analysis methods which can be directly applied for modal parameters estimation during a flight [7,8].Some are realized iteratively in real time, with less than 1 second interval between estimations [9].These methods applied for flight flutter testing shorten the testing procedure.Costs of the test procedure will then dramatically drop.During the proposed procedure, an aircraft is stimulated by operational excitation (turbulence, manoeuvres of airplane) or using dedicated excitation devices and its vibration responses are measured.Based on recorded signals, the modal parameters of the structure are estimated using any identification method.These parameters in the particular modal damping and natural frequencies are useful for flutter boundaries determination.The described procedure should be applied at each test point to determine Flight Clearance Envelope (FCE) [4].There are many different modal parameters identification methods that could be used for flight flutter testing.The methods can be categorized as experimental modal analysis methods and operational ones [8].The latter can be applied for an experiment scenario in which the excitation force is not measured.The methods can be implemented in the proposed flutterometer design.Both time domain and frequency domain methods match the requirements for on-line modal parameters identification [9].A wide range of methods can be considered for final implementation, many of them have been tested using simulated and real measurement data.The time-domain method, based on classical recursive RLS algorithm, is applied as the proposed solution and is efficient enough and relatively easy to implement.The major assumption applied in the procedure design is linearity of the system but the system can be treated as non-stationary due to changing of flight parameters and variation of system parameters due to different aerodynamic feedback at different flight points.Signal processing procedures should satisfy these assumptions.The wavelet based signal pre-processing is employed in the solution.The quality of modal parameters estimation methods depends on many different factors, but modal damping is a less accurate estimated parameter in the common modal analysis approach.Due to this, during the estimation process confidence bounds for modal damping should be computed to assess the quality of the estimation procedure.There are different methods for confidence bounds estimation in modal analysis techniques [11,15], but only a few can be realized on-line.The method implemented in the formulated procedure is based on the linearization of the relation between the ARMA model parameters variance and the standard deviation of the modal parameters.The method is applied for analysis of real flight vibration data.The results are carefully analysed to assess the quality and applicability of the approach to real time modal analysis of flying structures.

A method of flutter margin estimation
The time domain method is formulated for linear systems with time-varying parameters.The model parameters are estimated for each time step of the iteration based on values from a previous calculation.The recursive identification methods are widely used in adaptive control systems [16].The flowchart of the recursive procedure for damping estimation is depicted in Fig. 3.
Parameters estimators are computed at each sampling interval according to the minimum of objective function which is defined as follows: where: α is the forgetting factor, which describes the influence of the j-th sample on the i-th estimator, 0 < α 1.
If the model parameters are not changing, α should be set as 1.The proper choice of the α value gives better parameters tracking properties of the algorithm.The following heuristic procedure is proposed for proper choice of the forgetting factor: where; α r and α d oscillate between 0.95 and 0.99, e d is chosen experimentally for a given case.The basic assumption for the algorithm is that the matrix P is positive definite at each iteration step: The matrix P −1 can be singular if: -system is not stimulated well enough, -direct influence of input to output is observed, -over-fitting of polynomial A and C, -sampling frequency is to large, -improper value of α was chosen.The ARMA model's parameters in the presented method are estimated from the iterative formula: The procedure was tested based on the simulation of two degrees of freedom model with varying damping parameters.The same test was carried out with the time-varying forgetting factor.In order to find the damping ratio for given vibration natural mode, poles of the system should be found.The poles can be found by finding the roots of the polynomial A(z −1 ) = 0. Order of the polynomial depends on the number of modes, which should be considered.In a case of flutter, monitoring two or three vibration modes should be tested on-line during a flight [13,14].The roots are poles of the sampled (discrete) system.The relationship between continues and discrete roots are as follows [17]: where; λ i is the characteristic polynomial root for discrete case, δ i , Ω i are the damping ratio and natural frequency for continuous case, T is the sampling period.Based on the above formulas, the modal damping and natural frequencies have been estimated on-line.The second order model is assumed for each analysed mode and many independent models can be investigated concurrently.It allows for identifying modal parameters for many modes simultaneously.The number of modes that can be estimated depends on the applied hardware solution.The procedure should give results for all modes in timeless rather than the sampling period to avoid loss of any signal samples.The procedure can be applied to any number of measured signals but implementation is limited to two acceleration signals measured on an airplane structure.
To achieve a better performance of the proposed algorithms, the wavelets based signal filtering can be applied.Wavelets allows for filtering in both the time and frequency domain, thus for separation of vibration modes in time-varying systems.The Morlet type wavelet is used in the proposed algorithm as shown in Fig. 2.
Proposed flutterometer design is based on the application of a new technology that employs the system-on-aprogrammable-chip.This choice makes it possible to optimize the computational efficiency of implementation for used algorithms and makes the solution flexible enough to modify its parameters or even architecture.

Confidence bounds for flutter margin real time estimation results
During recursive identification, covariance matrix elements P (i) are estimated at each step of the procedure.Diagonal parameters of the covariance matrix are variances of the particular ARMA model parameters.These parameters are inputs of confidence bounds estimation for modal parameters for all modes under investigation.
The relationship between the ARMA model parameters and the modal model parameters is strongly nonlinear because the ARMA model is defined in the discrete time domain (sampled signal) but the modal parameters are identified in the continuous time domain.To find the confidence bound based on known ARMA parameters the variance linearization of relation between the ARMA and the modal parameters has to be linearized.The Taylor expansion method is proposed for this purpose.The covariance matrix for modal parameters can be estimated at each step of the procedure by applying the following formula: where: κn is the vector of the modal parameters at sample n, J() is the Jacobian matrix, θN is a vector of the ARMA model parameters, θ 0 is the vector of exact value of the ARMA model parameters, P(θ n ) is the covariance matrix of the ARMA model parameters, P k (κ n ) is the covariance matrix of the modal parameters.
The Jacobian matrix can be obtained from the numerical differentiation using the central difference theorem.The ith column of matrix J (θ N ) can be numerically estimated by: where; dθ i is a small ARMA parameter perturbation of the value one order of the magnitude smaller than the required accuracy of standard deviation estimation; f is a function expressing the relationship between the modal parameters and parameters of the ARMA model (Eqs ( 5) and ( 6)).
The standard deviation of modal parameters can be numerically calculated, according to Eqs ( 7) and ( 8), from the diagonal elements of matrix P κ (κ n ).

Implementation of the damping monitoring algorithm
Performance of calculations depends on how fast and how precisely individual parts of the algorithms are realized.Increased precision means longer time needed for calculations.Imposed tight execution time restrictions can be satisfied using a fast, expensive processor, a multiprocessor system or by mixed hardware-software implementation.The latter option was selected so as to allow for the greatest flexibility, future improvements and modifications.
All calculations are shared between software (dark gray blocks) and hardware (light gray blocks), as shown in Figs 2a and 3, on a system created on a Stratix FPGA chip.Software means that parts of the algorithm are written in C and then compiled for a Nios II soft-processor.Hardware fragments are realized in the logic of the FPGA and the Avalon Bus is used for data exchange.
The flutter monitoring algorithm consists of five steps (Fig. 2a): data read from the analog-to-digital converters (ADC), convolution of a stored wavelet with signals from sensors, signals reconstruction, on-line recursive least square (RLS) routine for estimation of parameters of the ARMA models, and determination of damping coefficients for flutter detection.The wavelets are generated during system initialization by the processor and stored in buffers in custom hardware.The floating-point operations of the software part are accelerated by custom instructions created in the Nios II ALU hardware, namely, floating-point addition, multiplication, reciprocal and square root.Data acquisition from the DACs is performed with a programmable frequency in the range of 10-200 Hz.
Normally, the most time-consuming part of the algorithm -convolution and signal reconstruction -is accelerated by a custom hardware accelerator and executed in fixed-point arithmetic.It is possible because of restricted resolution of input data and fixed size of both the wavelet table and the input buffer.Signal reconstruction is carried out simultaneously with convolution.The hardware accelerator (Fig. 2b) is controlled by a single master control unit  The processor and hardware accelerator work in a master-slave arrangement with the processor as master.Completion of calculations by the hardware accelerator is registered in status bits and generates a DMA request to transfer the results to software.The RLS algorithm and damping coefficients calculations are performed for each wavelet-filtered signal in the floating-point arithmetic by the processor being supported by custom instructions.The flutter appearance can then be determined using the damping coefficients thresholds table indexed by the actual flight conditions.Figure 3 shows overall hardware architecture of the flutterometer.

Verification of the designed flutterometer
The performance of the flutterometer has been verified on acceleration signals recorded with a sampling frequency of 160 Hz during flight test of a TS-11 Iskra military jet aircraft that performed at variable speeds within the 250-750 km/h range (Fig. 4).During the flutterometer verification experiment, two signals were used: one measured a fin and the second the tail of the plane.At both points, two modes of vibration dominate: the first with a frequency of 27 Hz and the second 47 Hz.They were decoupled by convolution of 512 signal samples with Morlet wavelets around 100 points long.Results (Fig. 5) were compared with those obtained in Matlab with direct RLS algorithm and ARMA models of a higher order.
Convolution of such lengths of signals and wavelets buffers with signal reconstruction takes 2.5 ms per sample, irrespective of the number of modes analyzed.Recursive RLS algorithm damping ratio calculations consume 0.5 ms of processor time per sample per mode.All calculations for the four modes take 4.5 ms per sample and results in over 200 Hz of maximum sampling frequency.Thus, processing 20 modes would take 12.5 ms, allowing for signals sampling with 80 Hz.
Assuming the confidence level to be 99.7%, uncertainty bounds for the modal parameters as estimated using the procedure described in Section 3, are shown in Fig. 6.
As it can be noticed in Fig. 6, the confidence bounds for the system under analysis are relatively small, but only at the start of the recursion process variability of parameters have an unacceptable level.

Conclusions
A formulated algorithm allows computing modal parameters of complex structures on-line during structure operation.Implementation of the flutter monitoring algorithm is proposed with the Hardware-Software Co-design approach, i.e. a part realized by hardware and the remaining part by software running on a Nios II soft-processor contained in the FPGA.The flutterometer is an example of the System-on-Chip, which allows for high level of integration and flexibility -it can be altered, e.g. to optimize for different algorithms, or to add some functionality, by reprogramming the FPGA without modifications of the PCB.For the tested structure, 12 first modes parameters have been identified simultaneously during 0.001 second.Confidence intervals for all parameters are relatively small and the method can be applied for flight flutter testing based on in-flight measurements.

Fig. 2 .
Fig. 2. a) Hardware-software partitioning of the flutter monitoring algorithm; b) architecture of the custom hardware accelerator for multi-channel wavelet transformation.

Fig. 5 .
Fig. 5. Modal damping of first vibration mode obtained in experiment, compared to calculated in Matlab without wavelet filtration.