A time domain method for the identification of dynamic parameters of structures

Mechanical Systems and Signal Processing (1993) 7(1), 45-56

A TIME D O M A I N M E T H O D F O R THE I D E N T I F I C A T I O N OF D Y N A M I C P A R A M E T E R S OF S T R U C T U R E S C. GONTIER, M. SMAIL Ecole Centrale de Paris, Chatenay-Malabry, France AND P. E. GAUTIER SNCF D~partement des Essais et des Laboratoires, Paris, France (Received 8 May 1991, accepted 17 November 1991) When using classic frequency methods, the identification of the modal parameters of a structure frequently fails in the event of high damping ratios and close modal frequencies, especially if associated with noisy measurements. This paper shows how, in such cases, time domain methods result in a correct and unbiased identification for a large range of parameter values, generally with a reduced sampling set of measures. The first part gives a short review of the mathematical framework. Formulations of the Laplace transfer function and z-transfer function are developed, resulting in a procedure for retrieving mechanical matrices of mass, damping stiffness. In the second part the ztransfer function is used to define the ARMA model of a structure. It is well-known in the field of automatics, that due to measurement noise, the least squares time domain estimates of the ARMA parameters are usually biased. An adaptive procedure for statistical estimation of the ARMA parameters with elimination of that bias is presented. Finally, the performances and limitations of the above procedure are determined by means of a large sequence of numerical simulations. The regions for correct identification in terms of damping ratio, sampling frequency and noise level are determined.

1. BASIC FORMULATION The experimental identification of the dynamic parameters of a structure is a problem whose importance from the engineering point of view is rapidly growing. This stems from the difficulty of determining a reliable dynamic model of the structure using classic finite element analysis. Experimental analysis therefore often remains unavoidable, especially when damping characteristics of the system (for instance in vibro-acoustics problems) are required. Most classic methods available in the field belong to the family of frequency methods. The performances and limits of these methods are well known. Some examples of classic misleading situations are: the presence of coupled modes with high damping ratios; the presence of coupled modes with close frequencies and measurements polluted with high level noise. In these cases, frequent methods can result in a dramatic lack of accuracy. Over the past decade, new methods derived from automatics, such as time domain methods, have been successfully transferred [1-8]. A R M A methods belong to this family. The main characteristic of time domain methods is that the determination of eigenmodes and eigenfrequencies is not considered to be the main purpose. A time domain procedure tries to 45

46

C. G O N T I E R E T AL.

fit the mathematical equation of a system for which an input sequence (excitation) and an output sequence (observation) are provided [8]. Here the modal parameters appear only as intermediate values towards the determination of mechanical parameters. The following demonstrates the ability of A R M A models to deal with the difficult situations mentioned above for a large range of parameter values. More precisely it will be shown that the bias resulting from the use of a simple least-squares method can be avoided through the use of an appropriate algorithm. I.I. DYNAMIC EQUATION Let us consider a mechanical system whose dynamic equation is

M2 + CYc+ Kx =f(t). Using the following notation

this equation may be written as

B ) ( - QX= F.

(1)

The inverse of matrix B being

B-t=IM-t_M-~CM-t] 0

m -I

premultiplying both sides of dynamic equation (1) by B -~ results in the classic state equation:

J(=AX+B-IF withA=B 1.2.

_,[0

Q = _M_~ K

,]

_M_~

ORTHOGONALITY PROPERTIES

From equation (1), the dynamic equation of the free system

B f ( - QX=O results in a standard eigenvalue form so that a set of 2n eigenvalues ;I.pand 2n eigenvectors ~,p can be found which verify

(;~pB- Q)~,,=O. From spectral theorems, the matrix Vt of the column eigenvectors is such that

vrBvt=b,

~trQv=q

(2)

where b and q are diagonal matrices. Normalisation of the column vectors of matrix ~, is arbitrary so that we can choose

Vr Q~t = I

(2a)

It is easy to prove that the diagonal matrices b and q are related with the diagonal matrix A of the eigenvalues 2p of A. Assuming that q = I it is found that:

b=A -~

(3)

TIME DOMAIN 1.3. N A T U R A L

IDENTIFICATION

47

SYSTEM RESPONSE

Due to the nature of the vector X, it is a classical result that eigenvector ~/'p can be written in the form

leading to the natural dynamic response:

x= E '~,e~"4'.+ E a* e~;%*. p=l,n

p=l,n

Finally, the eigenvalue and eigenvector matrices have the form

1.4. T R A N S F E R

FUNCTION

WITH LAPLACE TRANSFORM

Taking the Laplace transforms of the two sides, dynamic equation (1) becomes ( s B - Q)X(s) = F(s). Let the excitation equation be

and the observation equation be

the transfer function results in: H(s) = C ( s B - Q)-'D. 1.5. T R A N S F E R F U N C T I O N 1N T H E M O D A L S P A C E In the modal basis, the forced vibration solution is defined as

x(t)= Z ~At)wp p= I,.2n

From this expression and the dynamic equation (1), taking the Laplace transforms of both sides, the solution is found in the form:

v,.~P(s)

X(s)= Z x p - - ~ , . =

E x~

v,.v,. ~

#(s).

p= l ,2n

Splitting this equation in its two conjugate halves, the transfer function is:

p=~..

s---~ ]'

-

1.6. TRANSFER FUNC'rION WITH z TRANSFORMS The solution of the dynamic equation (1) is well known X(t) = ~ ( t , O)X(O) +

fo

~ ( t , r)B-~F(r) dr

48

c. GONTIER

ET AL.

¢, being the transform matrix given as • (t, to) = e At'-'°). This equation being discretised at time intervals T, the system's evolution is described by the equation Xk+~ = ~ X k +

q ~ ( t - r ) B -~ d r Fk÷l.

(4)

where a shortened notation was used ~ ( t ) = ~ ( t , 0), • = ~ ( T ) . Letting .(2 =

fo

~ ( t - r) d r =

fo

e m'-~ d r

then taking the z-transform of (4), the system's response results in

if= z(z:- ~ )-' n s - ' [~ where the symbol A holds for the z-transform of the function A. Introducing the same excitation and observation equations as defined in section 1.4, the z-transfer function is found: H(z) = C z ( z I - ~ )-l l2B-t D. 1.7. z-TRANSFER FUNCTION IN THE MODAL SPACE Matrix .(2 can be expressed in modal co-ordinates .(2 =

~ e A t ' - ~ -I d r = 7/

e A('-~) d r g - i .

Letting E = fo T e A ( ' - r ) d r = - A - I [ea~'- ~)]r = - A - l ( l - - e at) so that = g E ~ -I and also ( z I - q~)-I = ~t(zi_ea,)-l~t-I the z-transfer function becomes H(z) = C ~ z ( z l - e A' )-i iit -if2 B-J D. Replacing the expression of 12, the following expression is obtained : H(z) = Cv/z ( z I - e a' )-I E ~ -t B-I D. Then from equations (2) and (3) H(z) = C ~ ( l - z -I e A ' ) - J E A ~ r D

(5)

TIME DOMAIN IDENTIFICATION

49

and splitting the eigenvectors in their two conjugate halves, the transfer function becomes

(

d~pckr t-2,* e~, - - H(z)=p=E,.,, ~pep 1-z-'zp l-z-'z*] in which zp holds for zp=e x"r

(6)

and ep is an element of the diagonal matrix E. 1.8. RETRIEVING M, C, K MATRICES Conversely, if the whole transfer function of the structure has been determined, retrieving the dynamic matrices M, C, K can be performed from the same equations as above. From either of the two methods, Laplace transform or z-transform, retrieving M, C, K matrices from the transfer matrix can be performed by the following procedure: (a) Determine transfer function H as a rational fraction whose numerator is a matrix. The zeros of the denominator lead directly to the eigenfrequencies. In the case of the ztransfer function, matrix E is then determined from these poles through equation (5). (b) Decompose the transfer function into a sum of elementary fractions. It is worth noticing that for each mode a single column of the numerator matrix is sufficient. For instance, from the z-transfer matrix, the first columns of the numerator matrices for all pmodes Rp,, = ,a.pepqSpCp, i ;

p = t , rl

determine the whole transfer matrix. This property permits working with a single excitation point, if observation is performed at n points of the structure. (c) From these column vectors, the eigenmodes are easily derived; f f o m R p . , , t h e value Of~p.i

(¢..,)=._ R..,, ;L.e leads to the modal vector p--

Rp'l

(d) Then from the eigenvectors 4, retrieve M, C, K matrices using (2), (2a) and (3)• 1.9. R E M A R K : D E T E R M I N A T I O N OF E I G E N F R E Q U E N C I E S

The determination of eigenvalues ~.p from the poles Zp of H(z) through relation (6) calls into question the uniqueness since the argument ap of ~.p is not determined. ap=~pp+ 2kJr ,

- l r ~<~pp~
(7)

It will be verified that this argument leads directly to the frequency response cop ap a-~ k f" - 2--/- 277"- 2-7T + r so that relation (7) is equivalent to the Shannon's condition :fp
50

c. GONTIER

ET

AL.

2. IDENTIFICATION OF ARMA PARAMETERS 2.1. BASIC NOTATION

From the transfer function, each observation may be related to the input function through a 2n order recursive equation. Assuming f is any persistently exciting function [5], this relation has the form 2n

2n

(8)

y ( k ) = - ~.. diy(k - i)) + ~ b i f ( k - i) i= I

i= I

so that a set of 2n observation equations may be arranged in the matrix form Y = MO

with Y= (y(2n) . . . . . y ( N ) ) r, N = 4 n -

I

-y(?.n_l)

M

=

L-y(N-1)

...

-y(2n-2) .

.

l,

.

.

.

-y(N-2)

...

-y(0) .

.

.

.

...

f(En-1) .

-y(N-2n)

.

.

.

.

f(N-1)

.

.

...

f(0) .

.

f(N-2n)

and 0 = (dj, d2 . . . . . d2., bl, bz . . . . , b2n) (parameter vector). If there is no noise added to the observation y, the parameter vector may be directly computed with O=M-Iy.

Actually observation y is always corrupted by some added measurement noise so that a stochastic approach is necessary. 2.2. LEAST SQUARES METHOD

The optimal parameter vector 0 is now determined as the minimising vector for some error function which is the quadratic error function: J( O) = ( r ( k ) - MO) r( Y(k) - M ( O)

so that the best predictor is given by

0J(0)

- - = 0

00

leading to = (MrM)-'MTy.

(9)

In the standard least squares estimation it is assumed that M is a deterministic matrix and that the output y is corrupted by some independent measurement noise v: Y= MO + v so that the mean value E ( O ) = E { ( M r M ) - I M T ( M O + v) results in 0. Unfortunately this argument cannot be defended here since matrix M includes the nondeterministic terms y(k). It follows that equation (9) results in a biased estimate. To eliminate this distortion various methods in the field of adaptive algorithms may be used. We will concentrate on recursive methods, which appear best adapted to real-time applications and parallel processing.

TIME D O M A I N I D E N T I F I C A T I O N

51

2.3. P A R A M E T R I C ADAPTIVE A L G O R I T H M

With the notation ~pr(k) = (y(k), y ( k - 1) . . . . .

y ( k - 2n + 1), f ( k ) , f ( k - 1) . . . .

, f ( k - 2n + 1))

the best predictor for the observation y at time k + 1 is r(k + 1) = ~r(k) 0(k). Minimisation of the quadratic error J ( O ) = Y, lkl ( y ( i ) - - ~ o r ( i ) O ( k ) ) 2 leads to the best estimate of the parameter vector: O(k) =

q~(i- l)q~r(i-1) i

I

2 Y(i)q)(i-l)=G(k) i=1

2 y(i)rp(i-1). i=1

A recursive formulation of this equation is classical: let ~°(k) = y ( k ) - q~r(k - 1) O(k - I) the estimate of 0 is defined by O(k) = O ( k - 1 ) + G ( k ) ~ p ( k - 1)e°(k) with G - J ( k ) = G - t ( k - 1) + rp(k- 1)tpr(k - !). Using the

matrix inversion lemma: G(k) = G(k-

I)

G ( k - 1)~p(k- 1)qgr(k - 1)G(k- 1)

1 + q~r(k- 1)G(k-l)¢p(k- 1)

we find O(k)=O(k-

1)+ G ( k - 1 ) ~ ( k - 1 )

e°(k) 1 +cpr(k - l ) G ( k - 1)tp(k- 1)"

Up until now, the procedure is equivalent to the direct least-squares method, and is termed the "recursive least-squares method". The advantage is that this formulation allows adaptive filtering through various methods. The method proposed by Landau and M'Sirdi [2] to eliminate the measurement noise in M, tp and G consists of the replacement of the measured values y in tp by their best predictor rio. r°(k) = q~r(k-l)(3(k-l). From this choice, the values in the observation matrix do not depend directly on the measurements so that the noise can be filtered. Slightly different methods, such as the extended least-squares method, maximum likelihood methods and especially the instrumental variable method can be used. The latter method is often said to be best adapted with respect to noise elimination. However the choice of the instrumental variable often leads to loss of cross-correlation between this instrumental variable and the input signal, resulting in non-accurate identification.

3. IDENTIFICATION OF ARMA PARAMETERS 3.1. SYSTEM S I M U L A T I O N To our knowledge the performances of ARMA methods in the presence of noise have not been investigated systematically. We have tested the above adaptive method of ARMA identification on simulations of mechanical systems defined by their impulse responses. Up to now we have concentrated our efforts on test case systems with one or two degrees of freedom (dof).

52

c. GONTIER

ET AL.

The first test case corresponds to a system with a single dof and the impulse response h ( t ) = A e-"' sin cot,

a=og0q,

o9 = 2rrf

o9 = og0xfi-s~.

Its z-transfer function is n(z) =A

e -aT sin coT z -

1 - 2 e - " r COS OgTz -~ + e ~ r z -2"

The second test case corresponds to a system with 2 dof with the impulse response h ( t ) = A e-"'r sin col t + A 2 e-"~r sin r02t

Its z-transfer function has the form blz -I + b2z-2 + b3 z-3 H(z) = 1 + dlz -I + d2z-2 + d3z-3 + d4z -4"

The coefficients di and bi are functions of a,, a2, o9~, o92, T which for the sake of brevity will not be detailed here. In each case, the transfer function leads to a recursive equation of type (8). Used with any persistently exciting input function, for instance a gaussian white noise, this equation generates a system response x(k). This response is then corrupted with some uncorrelated noise simulating a measurement error. The resulting signal is submitted to the identification procedure. 3.2. SYSTEM WITH A SINGLE DEGREE OF FREEDOM

The identification procedure is tested on the case of a system with a unique modal frequency f = 27r/o9 = 10 Hz. Figure 1 compares the performances of the proposed method with those of the classic least squares method. The noise ratio is fixed to a value of 20%, the sampling frequency

200

150

n

J

n

m

_jj J

100 O

8. ~

501 /

O

"

"

,.--0---"01 •

-50 0"0(

•

•

•

0"20

•

0'40

II

0"60

I

0j

0"80

1-00

Domping rotio

Figure 1. Comparison of the adaptive filtering method ( t ) with the least squares method for 20% noise.

TIME DOMAIN IDENTIFICATION

53

to a value of 50 Hz. As mentioned above, the least squares estimates are systematically biased, and that bias increases with the damping ratio. On the contrary the tested procedure provides correct values up to an 80% damping ratio. This is the most important result of these investigations, it means that with this method, in the presence of an important measurement noise, A R M A identification can be performed correctly for any value of the damping ratio. Figure 2 shows the convergence rate of the procedure, in terms of the required sampling size, for different values of the damping ratio. Each curve holds for a particular value of the sampling frequency. The procedure is said to be convergent if modal parameters (co, r/) are identified within a maximum error of 5%.

ilh

o

so

0 5 0 Hz +75Hz Z~ 1 O0 Hz 0 1 5 0 Hz o 2 0 0 Hz I

0"00

0"20

~

I

0-40 Domping rotio

J

I

0"60

,

0"80

Figure 2. Required sampling size vs. damping ratio, for different sampling frequency values (noise level=

2o%).

The best sampling frequencies to use are clearly an outcome of the figure. A good choice in the range of five to ten times the modal frequency (50-100 Hz) allows an accurate identification with a reduced sampling set. Moreover the drawback of oversampling appears from this figure. Keeping in mind that the convergence of the least squares method is exponential [5], in the case of a persistently exciting input signal it is clear that the first samples are much more important for the identification than the last ones. In the case of oversampling, the measurement noise affects dramatically this initial estimation which cannot be corrected further. Figure 3 shows convergence rates of the procedure for various combinations of the measurement noise and damping ratio. Convergence is assumed if the modal parameters are identified within a tolerance of 5% for less than 20 000 time sample measurements (equivalent to 20 blocks of 1024 time samples in frequency analysis). What is interesting here is that a good result can be expected at noise levels greater than 100%, even up to 250% if the optimal sampling rate has been chosen for the mode. Figure 4 is a 3-D plot-graph for the regions of successful identification in terms of damping ratio, sampling frequency and noise level.

54

ET AL.

C. GONTIER

oooL o25Hz

iii

-~ :p

• 75Hz

0 100Hz • 150Hz a 200Hz

A 300Hz

,

0-0

I

,

I

0-2

,

A

0.4

0 6

,

-¢~ 08

~ 10

Damping rotio Figure 3. Allowable noise level vs. damping ratio, For different sampling frequency values.

I o o ' ~ ~ ~

1

(O,o)/

~oo~ | ~ o

Noise level •

/

Damping ratio

>'ou

300

~

om

.

).80 0-70

o.~o

u.-u

~,~ Sampling frequency (Hz)

Figure 4. 3-D graphoplot or allowable noise level vs. damping ratio for different sampling frequency values.

3.3. SYSTEM wrrH TWO DEGREES OF FREEI)OM In this ease, the procedure is tested for its ability to separate very close modal frequencies. The reference frequency is 10 Hz and the sampling flequency 50 Hz. Figure 5 shows the rate of convergence of the procedure as a function of the difference of modal frequencies in the range from 0 to 1 Hz (0 to 10%). The convergence criterion is the same as above. The damping ratio of the reference mode has a value of 10%. The different curves hold for different values or the second damping factor.

TIME D O M A I N I D E N T I F I C A T I O N

55

100

80

4w

60 Domping r o t i o : o0.01 00.25 40

2( 0"00

,

I 020

,

I 0"40

,

I 060

,

I 0"80

, 1"00

( f 2 - f l ) HZ

Figure 5. Convergence rate vs. Frequency difference, for different second damping ratio values (noise level =

IO%). Note that the frequential discrimination succeeds for frequency differences down to 0.1 Hz, for most values of the second damping ratio. However no curve is drawn for coincident damping ratios 0.1 for the two modes. The reason is that in this case of close modal frequencies and simultaneous Close damping factors the procedure's performance breaks down because the two modes combine in a single mode and must be identified as such.

4. CONCLUSION An A R M A formulation of the modal identification problem has been presented. After recalling that the classic identification procedure based on the least-squares algorithm provides biased estimates, an unbiased adaptive recursive procedure, suggested by Landau, was adapted to the problem. Its performance in such cases as close coupled modes with high damping ratios and high level measurement noise was assessed. A systematic determination of the optimal sampling frequency and the required sampling size was performed. Further investigations will involve a multivariate implementation of the algorithm and its application to an industrial example.

REFERENCES I. P. E. GAUTmR 1986 June Ecole Centrale de Paris Chatenay, Malabry, France, Doctoral thesis. Comportement dynamique des v6hicules ferroviaires. 2. I. D. LaNoau and K. N. M'SIRDI 1986 Traitement du signal 3/4 5, 183-204. Techniques de mod~lisation r+cursive pour I'analyse spectrale param&rique adaptative. 3. N. C. MICKLEBOROUGH and Y. L. Pt 1989 hzternational Journal for Numerical Methods in Engineering 28, 2307-2321. Modal parameter identification using z-transforms. 4. N. C. MXCKELBOROUGHand Y. L. PI 1989 Journal of Engineering Mechanics: 115(10), 22322250. Modal identification of vibrating structures using ARMA model.

56

C. GONTIER ET AL.

5. K N. M'SmDI 1988 Universit~ de Grenoble France, Doctoral thesis. Mod~lisation param~trique adaptative et application ~ l'analyse spectrale. 6. D. BONNECASE, M. PREVOSTO and A. BENVENISTE 1990 IMAC Proceedings, 382-385. Application of a multidimensional ARMA model to modal analysis. 7. F. R. VIGNEnON 1986 Journal of Applied Mechanics 53, 33-38. A natural modes model and modal identities for damped linear structures. 8. J. M. LEURIOAN,D. L. BROWN and R. J. ALLEMANG 1986 Journal of Vibration, Acoustics, Stress, and Reliability in Design 108, I-8. Time domain parameter identification methods for linear modal analysis: a unifying approach. 9. D. J. EwiNs 1986 Modal Testing, Theoo, and Practice. Letchworth : Research Studies Press.