The PolyMAX frequency-domain method : a new standard for modal parameter estimation ?

Recently, a new non-iterative frequency-domain parameter estimation method was proposed. It is based on a (weighted) least-squares approach and uses multiple-input-multiple-output frequency response functions as primary data. This so-called “PolyMAX” or polyreference least-squares complex frequency-domain method can be implemented in a very similar way as the industry standard polyreference (time-domain) least-squares complex exponential method: in a first step a stabilisation diagram is constructed containing frequency, damping and participation information. Next, the mode shapes are found in a second leastsquares step, based on the user selection of stable poles. One of the specific advantages of the technique lies in the very stable identification of the system poles and participation factors as a function of the specified system order, leading to easy-to-interpret stabilisation diagrams. This implies a potential for automating the method and to apply it to “difficult” estimation cases such as high-order and/or highly damped systems with large modal overlap. Some real-life automotive and aerospace case studies are discussed. PolyMAX is compared with classical methods concerning stability, accuracy of the estimated modal parameters and quality of the frequency response function synthesis.


Introduction
Experimental Modal Analysis (EMA) is currently one of the key technologies in structural dynamics analysis.Based on the academic fundaments of system identification, it has evolved to become a "standard" approach in mechanical product development.Essential in this evolution is that modal analysis research has, from the start, taken the point of view of industrial applicability, focusing on solving the specific problems related to testing and modelling of large industrial structures.The merit of each new method or new approach has always been checked against the added value it brought in terms of helping the application engineers to derive better models.The result is that EMA is now considered as a "commodity" tool, continuously expanding its application base.
In this context it is worthwhile reflecting on the current limitations or new challenges, as experienced by modal analysis users as these limitations constitute the background for further research.
One important consideration is that the role of Testing in the whole product design process is changing [1].Design cycles and costs are under pressure, while product quality demands continuously increase due to the competitive demands of the global market.As a consequence, Virtual Product Engineering using CAE methodologies is taking up a much larger part of the whole design cycle, trying to solve design problems as early as possible in the design cycle.Test results on critical components or on early product prototypes, however, remain essential, not only for model validation purposes, but also for direct use in "Hybrid" product models.Experimental models obtained on product variants or predecessors furthermore constitute a core knowledge base for starting of a new design.All these applications require EMA models with high accuracy, and which have to be obtained in the very short time available for prototype or component testing.Important problems faced in this process are: -The validity of the test results on which the analysis is performed, as there is often no chance for redoing the tests (test-data plausibility).Consistency between the measured test data is of the highest importance for applications such as Substructuring with other Test of FE models.-The complex task of selecting the correct model order and discriminating between spurious and structural system poles, in particular in the case of high-order and/or highly-damped structures.-The expertise of the users as the general acceptance of EMA has transferred it from the realm of the research experts to the product development workfloor.
A key requirement for EMA hence is that a reliable analysis of complex datasets should be possible with minimal user interaction.This is the context of the methodology developed in the present paper.
In Section 2, the PolyMAX method is reviewed.This method very successfully deals with "difficult" industrial cases.A particular strength, however, is that it is easy to use, making it very accessible to the non-expert.Section 3 positions the new method relative to the state-of-the-art methods.Section 4 discusses some difficult cases: the identification of trimmed car bodies and the analysis of flight flutter test data.

The "PolyMAX" or polyreference least-squares complex frequency-domain method
Originally, the least-squares complex frequency-domain (LSCF) estimation method was introduced to find initial values for the iterative maximum likelihood method [2].The method estimates a so-called common-denominator transfer function model [3].Quickly it was found that these "initial values" yielded already very accurate modal parameters with a very small computational effort [2,4,5].The most important advantage of the LSCF estimator over the available and widely applied parameter estimation techniques [6,7] is the fact that very clear stabilisation diagrams are obtained (Section 2.4).A thorough analysis of different variants of the common-denominator LSCF method can be found in [5].A complete background on frequency-domain system identification can be found in [8].
It was found that the identified common-denominator model closely fitted the measured frequency response function (FRF) data.However, when converting this model to a modal model Eq. ( 22) by reducing the residues to a rank-one matrix using the singular value decomposition (SVD), the quality of the fit decreased [4].
Another feature of the common-denominator implementation is that the stabilisation diagram can only be constructed using pole information (eigenfrequencies and damping ratios).Neither participation factors nor mode shapes are available at first instance [9].The theoretically associated drawback is that closely spaced poles will erroneously show up as a single pole.
These two reasons provided the motivation for developing "PolyMAX", which is a polyreference version of the LSCF method, using a so-called right matrix-fraction model.In this approach, also the participation factors are available when constructing the stabilisation diagram.The main benefits of the polyreference method are the facts that the SVD step to decompose the residues can be avoided and that closely spaced poles can be separated.The method was introduced in [9] together with a simulated example including a pair of closely spaced modes.The present paper reviews the theory and validates the method using industrial data sets.

Right matrix-fraction model
The polyreference LSCF method needs FRFs as primary data and identifies a right matrix-fraction model: where H(ω) ∈ C 1×m is the FRF matrix containing the FRFs between all m inputs and all l outputs; B(ω) ∈ C 1×m is the numerator matrix polynomial and A(ω) ∈ C m×m is the denominator matrix polynomial.Each row of the right matrix-fraction model can be written as: The numerator row-vector polynomial of output O and the denominator matrix polynomial are defined as: where Ω r (ω) are the polynomial basis functions and p is the polynomial order.In the LSCF method, a z-domain model is used (i.e. a frequency-domain model that is derived from a discrete-time model) and, by consequence, the basis functions are (with ∆t the sampling time): The polynomial coefficients β or ∈ R 1×m and α r ∈ R m×m are assembled in following matrices: The FRF model of Eq. ( 1) is now written as function of the coefficients H(ω k , θ).Note that in the LSCF formulation of this paper real-valued polynomial coefficients are assumed.It is possible to derive another variant of the method with complex coefficients.A discussion on these aspects for the common-denominator LSCF implementation can be found in [5].

Equation error formulation
The challenge is now to find all unknown model coefficients based on measurements (or better: non-parametric estimates) of the FRFs Ĥo (ω k ), where • denotes a measured quantity and ω k (k = 1, 2, . . ., N f ) are the discrete frequencies at which FRF measurements are available.The coefficients θ can be identified by minimising the following non-linear least-squares (NLS) equation errors where the scalar weighting function w o (ω k ) is introduced.This frequency-and output-dependent weighting function allows taking into account data quality differences that may exist between different outputs.The equation errors for all outputs and all frequency lines are combined in following scalar cost function: where • H is the complex conjugate transpose (Hermitian) of a matrix and tr{•} is the trace of a matrix (sum of the elements on the main diagonal).In fact, the 'trace' notation boils down to computing the squared 2-norm of the row vector ε NLS o (ω k , θ).The cost function is minimised by putting the derivatives of Eq. ( 7) with respect to the unknown model coefficients equal to zero.It is obvious that this leads to non-linear equations when expression Eq. ( 6) is used for the equation errors.This non-linear least-squares problem can be approximated by a (sub-optimal) linear least-squares one by right-multiplying Eq. ( 6) with the numerator matrix polynomial A, yielding equation errors ε LS o (ω k , θ) ∈ C 1×m that are linear in the parameters: The equation errors at all frequency lines are stacked in a matrix with: (11) where ⊗ denotes the Kronecker product.

Reduced normal equations
Similar as Eq. ( 7), the following cost function can be formulated based on the linearised equation errors: Minimising this cost function leads to a weighted linear least squares problem.Introducing the matrix notation of Eqs ( 9), ( 10) and ( 11), the cost function can be written as: is the so-called Jacobian matrix given by: In the case of real-valued coefficients θ, it can be shown that the expression J H J can be substituted by its real part Re (J H J) ∈ R (l+m)(p+1)×(l+m)(p+1) and the cost function Eq. ( 13) becomes: l LS (θ) = tr{θ T Re(J H J)θ} (15) with: in which following notation is used: The cost function Eq. ( 15) is minimised by putting its derivatives with respect to the unknown polynomial coefficients θ equal to zero: These equations are the so-called normal equations which can be written synthetically using Eq. ( 16) and Eq. ( 5): At first instance, we are only interested in the denominator polynomial α from which the poles and modal participation factors can be computed.This information is sufficient to construct a stabilisation diagram.By consequence, the least-squares problem can be reduced by eliminating the β o coefficients from Eq. ( 18): yielding the so-called reduced normal equations: where M ∈ R m(p+1)×m(p+1) is defined in above equation and can be computed from the measured FRF data.This equation can be solved for the denominator polynomial α in a least-squares sense.To avoid finding the trivial solution α = 0, a constraint is imposed on the parameters.Such a constraint also removes the parameter redundancy that exists in the right matrix-fraction model (multiplying numerator and denominator with the same matrix yields different numerator and denominator polynomials, but the same transfer function matrix).

Practical implementation
Theoretically, a complete right matrix-fraction model can be obtained by computing the numerator coefficients β o from Eq. ( 20).An alternative approach, to find a complete model is described in the following.This approach is very similar to the classical implementation of the least-squares complex exponential (LSCE) method [10].
In modal analysis applications, one is usually not only interested in a good model as such, but more important, in the accuracy of the estimated modal parameters.When trying to estimate the modal parameters from real data, it is generally a good idea to over-specify the model order considerably; i.e. to try to fit high-order models that contain much more modes than present in the data.Afterwards the true physical modes are separated from the spurious numerical ones by interpreting a so-called stabilisation diagram [11].The poles corresponding to a certain model order are compared to the poles of a one-order-lower model.If their differences are within pre-set limits, the pole is labelled as a stable one.The spurious numerical poles will not stabilise at all during this process and can be sorted out of the modal parameter data set more easily.The use of a stabilisation diagram solves a typical system identification problem, namely the determination of the model order.Examples of stabilisation diagrams will be shown in Section 4.
The poles and the modal participation factors are found from the LSCF denominator polynomial Eq. ( 21) by computing the eigenvalue decomposition of its companion matrix, which is rather classical [6].A p th order right matrix-fraction model yields pm poles.An efficient construction of the PolyMAX stabilisation diagram is possible by formulating the least squares problem Eq. ( 21) for the maximum model order p max .Considering submatrices of appropriate dimensions can then solve lower-order problems.
The interpretation of the stabilisation diagram yields a set of poles and corresponding participation factors.The mode shapes can then be found by considering the so-called pole-residue model (This model is the partial fraction description of the Laplace-domain version of Eq. ( 1)): where n is the number of modes; • * is the complex conjugate of a matrix; {v i } ∈ C l are the mode shapes; l T i ∈ C m are the modal participation factors and λ i are the poles, which are occurring in complex-conjugated pairs and are related to the eigenfrequencies ω i and damping ratios ξ i as follows: LR, U R ∈ R l×m have been introduced in Eq. ( 22) and are respectively the lower and upper residuals modelling the influence of the out-of-band modes in the considered frequency band.The numerators in Eq. ( 22) are the so-called residues {v i } l T i .They are rank-one matrices by definition.As the participation factors l T i and the poles λ i are estimated based on the α polynomial and the stabilisation diagram, the only unknowns in Eq. ( 22) are the mode shapes {v i } and the lower and upper residuals.They are readily obtained by solving Eq. ( 22) in a linear least-squares sense.This second step is commonly called least-squares frequency-domain (LSFD) method [6,10].

PolyMAX versus LSCE
As may be clear from previous section, the PolyMAX method proceeds along similar lines as the polyreference LSCE time-domain method: -Establishment of a set of linear equations for the maximum required model order, from which the matrix polynomial coefficients α can be computed in a least-squares sense.-Construction of a stabilisation diagram by solving an eigenvalue problem for increasing model orders.The information regarding eigenfrequencies, damping ratios and modal participation factors is contained in this diagram.-Based on the user-interpretation of the stabilisation diagram, computation of the mode shapes and the lower and upper residuals by solving Eq. ( 22) in a least-squares sense.
The difference between LSCE and PolyMAX lies in the first step.LSCE uses impulse responses to find the polynomial coefficients, whereas PolyMAX requires frequency response functions.
However, it turned out that the new PolyMAX estimator yields extremely clear stabilisation diagrams making it very simple to select the "physical" poles.In the LSCE method, the non-physical (and sometimes even the physical) poles tend to "wander" in the stabilisation diagram, especially at large model orders.The PolyMAX method has the interesting property that the non-physical poles are estimated with a negative damping ratio so that they are easily excluded from the set of poles before plotting them.Such a clear diagram does not mean that some of the poles are missing; on the contrary, more (physical) poles can be found than with the earlier state-of-the-art methods, as evidenced by the examples in next section.Other validation studies also revealed that the PolyMAX method has no problems in correctly estimating modes having a low damping ratio.It is sometimes stated that time-domain methods are preferred in case of low damping, and frequency-domain methods in case of high damping.The PolyMAX frequency-domain method excels in both cases.

PolyMAX versus other frequency-domain methods
Typical for many frequency-domain parameter estimation methods is that they involve the inversion of a matrix containing powers of the frequency-axis of the data.Therefore one quickly runs into numerical conditioning problems and strict constraints apply to both the frequency range and the model order range of the analysis.In the past, it has been proposed to use an orthogonal polynomial basis for the frequency-domain model to solve the numerical problems.However, this increases significantly the computational time and memory requirements.
The PolyMAX method does not suffer from numerical problems as it is formulated in the z-domain (i.e. a frequency-domain model that is derived from a discrete-time model),whereas the existing frequency-domain methods use a Laplace-domain formulation (i.e. a frequency-domain model that is derived from a continuous-time model).In PolyMAX, the frequency axis that extends between f 0 and f end is shifted and mapped into half a unit circle in the complex plane Eq. ( 4): Similar to other frequency-domain methods, the PolyMAX method involves the inversion of a matrix containing powers of the frequency-axis of the data.The main advantage here is that taking powers of the z-variable does not increase the range of the values, as it boils down to a rotation in the complex plane.As a result, PolyMAX can deal with a large frequency range and very large model orders, speeding up the modal parameter estimation process considerably, as in many cases the complete frequency-band of interest can be processed at once.
There was some common belief that the numerical conditioning of frequency-domain methods is worse than time domain methods and that broadband analyses are preferably performed in time-domain [6].When using PolyMAX, these statements are definitely no longer true.

Computational efficiency
The advantages discussed here have no price in terms of computational time: PolyMAX is as fast as LSCE.LSCE became the industry-standard because of its high speed even for a very large number of measured outputs.A lot of research efforts [2,4,5,9] were spent to achieve this computational efficiency.On current PC platforms, calculation and display of the stabilisation diagram for a typical full car body model (like the trimmed car body examples discussed in next section) only takes a few seconds.

Highly-damped cases 4.1.1. Full vehicle
A typical example of a challenging modal analysis application is the identification of a trimmed car.The trim material turns a nicely resonating car body into a highly damped system with large modal overlap.In the present example, data from a Porsche 911 Targa Carrera 4 were used.The accelerations of the fully equipped car were measured at 154 degrees of freedom (DOF), while 4 shakers were simultaneously exciting the structure.More details about the test can be found in [12].Interesting about the data set are, next to the high damping, the large number of references.This allows investigating whether the property of very clear stabilisation diagrams present in the LSCF method [2,4,5] carries over to the polyreference case, which is the PolyMAX method.
The data were analysed in a frequency range from 3.5 to 30 Hz using three techniques: -The frequency domain direct parameter identification (FDPI) technique, which has traditionally been used to analyse data from highly-damped structures [6,10,13,14].-The least squares complex exponential (LSCE) method, which is the "industry-standard" [6,10,15].
-The new PolyMAX method; see [9] and Section 2 of this paper.
Figure 1 shows the stabilisation diagrams for the 3 methods.FDPI yields a clearer diagram than LSCE, confirming the common assumption that FDPI is the preferred method in case of high damping [6,10].However, the best diagram is obtained using PolyMAX: especially at lower frequencies it is much clearer than the FDPI diagram and also a larger number of stable poles is found.In the case of a large frequency band and high model orders, the z-domain formulation of the LSCF method is superior over the Laplace domain formulation of the FDPI method.As shown in [12], more stable FDPI poles can be found when splitting the considered frequency band in smaller sub-bands.For the subset of PolyMAX poles that have an FDPI counterpart, Table 1 shows that the estimates of frequency and damping are very close.The damping ratios are rather high, as can be expected from a trimmed car.Also the mode shapes are very similar, which is evidenced by the MAC values represented in Fig. 2. The modal assurance criterion (MAC) is defined as the squared correlation coefficient between two modal vectors v i , v j : The excellent identification results obtained with the PolyMAX method also show up when comparing the measured FRFs with the FRFs that can be synthesised from the modal parameters based on Eq. ( 22).Such a comparison is made in Fig. 3.The correspondence between two FRFs can also be quantified in single numbers [10].The correlation between a measured Ĥoi (ω k ) and synthesised FRF H oi (ω k , θ) is defined as: The normalised error between a measured Ĥoi (ω k ) and synthesised FRF H oi (ω k ) is defined as: As stated in the introduction of Section 2, one of the motivations to consider a polyreference implementation of the LSCF method was the fact that a better FRF fit is expected since there exist no need to reduce the residues to rank-one matrices using an SVD.This expectation is confirmed in this example.Both the common-denominator as the polyreference version of LSCF were compared and the quality of the FRF fit was assessed by averaging Eqs (26) and ( 27) over all 154 outputs 'o' and 4 inputs 'i'.The results of this assessment can be found in Table 2.The quality of the fit using the full-rank residues is comparable to the polyreference quality.However, when reducing the rank of the residues to find a true modal model Eq. ( 22), the quality significantly decreases.Another evidence of this fact can be found in [4].

Partially trimmed car
Although the FDPI method performed satisfactorily in the Porsche example discussed in previous subsection, this is not always the case.In present subsection, data from a partially trimmed car are used.The modal damping ratios are a bit lower (2-3%) than in the fully equipped car.The considered frequency range is 35-75 Hz.There were more than 500 output DOFs and 2 inputs.Figure 4 compares the FDPI and the PolyMAX stabilisation diagram.The FDPI stabilisation diagram is hardly to be interpreted at all and, again, PolyMAX is very successful in finding the system poles.

Noisy case
Cases exist where the FRF data are highly contaminated by noise.Flight flutter testing is such a case.In the example considered here, both wing tips of a plane are excited during the flight with so-called rotating vanes.These vanes generate a sine sweep through the frequency range of interest.The forces are measured by strain gauges.Next to these measurable forces, also turbulences are exciting the plane resulting in rather noisy FRFs.
Figure 5 shows some typical multiple coherences and the corresponding FRFs.They clearly show the noisy character of the data.During the flight, accelerations were measured at 9 locations, while both wing tips were excited (2 inputs).The data were analysed using both the LSCE and the PolyMAX method.Figure 6 shows both stabilisation diagrams.Also in this example, PolyMAX has an advantage over the LSCE method, in the sense that it is much easier to select the true system poles.
When looking at the synthesised FRFs (Fig. 7), it appears that the PolyMAX method is able to accurately model the data.Obviously this becomes more difficult as the FRF becomes noisier (Fig. 7, Right).
In this context, it is interesting to mention some recent developments in which the turbulence excitation of the airplane is not just considered as noise, but as useful input to the system.Very promising results are obtained with these so-called combined deterministic-stochastic identification methods [16,17].

Conclusions
The proposed PolyMAX or polyreference least-squares complex frequency-domain (LSCF) modal parameter estimation method offers very good prospects in view of the current EMA challenges.The method in general yields very clean stabilisation diagrams, easing dramatically the problem of selecting the model order and the best structural system poles.The analysis results compare favourably with current best-of-class commercially available methods, without increasing the computational effort.Its very good identification behaviour for noisy data sets as well as for high order, highly damped structures allows to position PolyMAX as a potential new standard general purpose method [18].

Fig. 2 .
Fig. 2. MAC values assessing the mode shape correlation between FDPI and corresponding PolyMAX mode shapes.

Fig. 3 .
Fig. 3. Comparison of the measured FRFs (grey) with FRFs synthesised from the identified modal model (black).The FRFs between the 4 inputs and 1 typical output are shown.

Fig. 4 .
Fig. 4. Stabilisation diagrams obtained by applying different parameter estimation methods to data from a partially trimmed car: (Top) FDPI; (Bottom) the new PolyMAX method.

5 .
Flight flutter test data.(Left) multiple coherences of a sensor at the wing tip close to the excitation (black) and a sensor at the back of the plane (grey); (Right) corresponding FRFs.The frequency axis is blind for confidentiality reasons.

Fig. 6 .
Fig. 6.Stabilisation diagrams obtained by applying different parameter estimation methods to the flight flutter test data: (Top) LSCE; (Bottom) the new PolyMAX method.

7 .
Comparison of the measured FRFs (grey) with FRFs synthesised from the identified modal model (black).(Left) Sensor at the wing tip; (Right) Sensor at the back of the plane.

Table 1
Eigenfrequencies and damping ratios obtained by applying the FDPI and the PolyMAX method

Table 2
Averaged correlations and errors between measured and synthesised FRFs.A comparison is made between the polyreference (i.e.PolyMAX) and the common-denominator implementation of LSCF