ON THE PROBLEM OF STOCHASTIC EXPERIMENTAL MODAL ANALYSIS BASED ON MULTIPLE-EXCITAION MULTIPLE-RESPONSE DATA - PART II: THE MODAL ANALYSIS APPROACH S.D. \Fassois and J.E. Lee2 Report UM-MEAM-91-15 @1991 by S.D. Fassois and J.E. Lee All rights reserved. 1 Assistant Professor, Mechanical Engineering and Applied Mechanics 2 Graduate Research Assistant, Mechanical Engineering and Applied Mechanics

ON THE PROBLEM OF STOCHASTIC EXPERIMENTAL MODAL ANALYSIS BASED ON MULTIPLE-EXCITATION MULTIPLE-RESPONSE DATA - PART II: THE MODAL ANALYSIS APPROACHt S.D. Fassois* and J.E. Lee** Department of Mechanical Engineering and Applied Mechanics The University of Michigan Ann Arbor, Michigan 48109-2121 ABSTRACT In this part of the paper the stochastic multiple-excitation multiple-response experimental modal analysis problem is considered. The relationship between the actual structural and noise dynamics and their discrete special-form ARMAX-type representation is studied for each one of the vibration displacement, velocity, and acceleration data cases, and a novel and effective modal analysis approach, that, unlike previous schemes, is capable of operating on any one of these types of data records, is introduced. By accounting for issues such as the required excitation signal type and stochastic model form, algorithmic instability occurrence and other well-known estimation difficulties, model structure estimation and model validation, as well as model reduction and analysis based on the Dispersion Analysis methodology introduced in the first part of the paper [1], the proposed approach not only overcomes the limitations and drawbacks of current schemes but also constitutes the first comprehensive procedure for stochastic multiple-excitation multiple-response experimental modal analysis. The effectiveness of the approach is demonstrated through numerical experiments with structural systems characterized by well-separated and closely-spaced modes, and data records of various lengths and noise-to-signal ratios. Comparisons with the classical Frequency Domain Method and the deterministic Eigensystem Realization Algorithm are also made, and the approach is finally used for the experimental modal analysis of a three-span beam froni laboratory data. t OCopyright 1991 by S.D. Fassois and J.E. Lee. All rights reserved. * Assistant Professor, Mechanical Engineering and Applied Mechanics. (TO WHOM ALL CORRESPONDENCE SHOULD BE ADDRESSED) ** Research Assistant, Mechanical Engineering and Applied Mechanics.

IIA 1-4 -.1 f,'- S slews' < - lz~

1. INTRODUCTION In the first part of this paper [1] the problem of quantitatively assessing the relative importance of a structural system's vibrational modes was examined, and an appropriate and physically meaningful Dispersion Analysis methodology introduced. According to this methodology the significance of a particular mode is determined by its contribution (dispersion) to the total vibration signal energy. Physical interpretations of modal dispersion were provided in both the correlation and spectral domains, and the phenomenon of modes characterized by negative dispersion was investigated and properly interpreted. In this second part of the paper the complete stochastic experimental modal analysis problem is considered, and a comprehensive and effective multiple-excitation multiple-response [also referred to as Multiple-Input Multiple-Output (MIMO)] approach, of which the Dispersion Analysis methodology is an essential part, developed. As is well-known, and due to a number of difficulties associated with extending stochastic single-response into multiple-response methods [1], very few such approaches are currently available [2,3]. The proposed one is not only the first that is capable of operating on either displacement, or velocity, or acceleration vibration data and overcoming their major weaknesses and limitations (see the discussion in the introduction of part I), but also the first comprehensive scheme for the solution of the stochastic multiple-excitation multiple-response experimental modal analysis problem, as it considers all of its important aspects and takes into account all major issues ranging from the selection of an appropriate type of force excitation signal through the assessment of the relative importance of the estimated vibrational modes. Indeed, the proposed approach accounts for and attempts to minimize all major types of errors that are typically associated with modal parameter estimates, namely: * Type I Errors due to the use of an inappropriate discrete model form (for instance a model with inappropriate orders or of a structure incompatible with the nature of the response measurements or the excitation signal type). * Type II Errors due to the parameter estimation procedure. * Type III Errors due to the modal parameter extraction procedure (for instance errors due to the use of a discrete-to-continuous model transformation that is incompatible with the excitation signal type and the discrete model structure). In addition, the proposed approach is characterized by a number of important features that enable it to overcome well-known weaknesses of general modal analysis methods that are due to the estimation stage (algorithmic instability occurrence, wrong convergence, the need for initial guess

parameter values, high computational complexity), inappropriate selection of the triple consisting of the discrete model structure, excitation signal type, and discrete-to-continuous transformation, as well as the lack of an effective model structure estimation scheme and a methodology for assessing the relative importance of the estimated vibrational modes (as necessary for model reduction, distinction of "structural" from "extraneous" modes, and subsequent structural analysis). The proposed stochastic modal analysis approach specifically consists of: * Appropriately selected and mutually compatible force excitation signal type and special-form discrete-time stochastic multivariate models for the combined representation of the structural and noise dynamics in each one of the vibration displacement, velocity, and acceleration measurement case. This selection is instrumental in overcoming commonly encountered Type I errors. * A realistic and effective Linear Multi-Stage (LMS) parameter estimation algorithm that overcomes difficulties such as algorithmic instability occurrence, wrong convergence, the need for initial guess parameter values, and excessive computational complexity, while also attempting to minimize Type II errors. * Effective model structure estimation and model validation procedures that are also important in overcoming Type I errors. * An appropriate discrete-to-continuous model transformation that forms a compatible triple with the force excitation signal type and the stochastic model form and minimizes Type III errors, and, an effective and physically meaningful procedure for model analysis and reduction (including the distinction between "structural" and "extraneous" modes and the determination of the former's global and local characteristics) based on the Dispersion Analysis methodology. The performance characteristics of the proposed experimental modal analysis approach are evaluated via a number of numerical experiments. For this purpose the analysis of structural systems characterized by both closely-spaced and well-separated modes is considered, and attention is paid to issues such as the achievable estimation accuracy and resolution at various N/S ratios, the effects of data record length on estimation accuracy, the effectiveness of the model structure estimation and model validation schemes, as well as that of the Dispersion-Analysis-based model reduction procedure. In addition, comparisons with the classical Frequency Domain Method (FDM) and the deterministic Eigensystem Realization Algorithm (ERA) [4] are made, and experimental results concerning the application of the proposed approach to the modal analysis of a three-span beam from laboratory data are presented. The rest of this paper is organized as follows: The proposed experimental modal analysis

approach is presented in Section 2, and its performance characteristics studied in Section 3 by using simulated structural systems. In Section 4 the approach is used for the analysis of a threespan beam from experimental data, and the conclusions are finally summarized in Section 5. 2. THE MODAL ANALYSIS APPROACH Consider a linear, viscously-damped, n degree-of-freedom structural system described by the vector differential equation: M i;(t) + C. r(t) + K v(t) = f(t) (1) where M, C, K represent the real and symmetric nx n mass, viscous damping, and stiffness matrices, respectively, {f(t)} the n-dimensional force excitation signal, and {v(t)} the resulting n-dimensional vibration displacement signal. The problem examined in this paper may be stated as follows: "Select an appropriate type of vector force signal {f(t)} with which the structural system under study is to be excited. Based on uniformly-sampled measured excitation {f[t]}N1=, and corresponding noise-corrupted displacement {y[t]}N=, or velocity {y,[t]}N=l, or acceleration {ya[t]}N=1 vibration test data, estimate an accurate modal representation (consisting of the complete set of natural frequencies, damping ratios, and mode shapes) of the structural system under study." The measured vibration displacement {y[tl}, or velocity {y,[t]}, or acceleration {ya[t]} data are of the form: y[t] = v[t] +n[t] (2a) Yv[t] = vv[t] + n[t] (2b) Ya[t] = Va[t] + n[t] (2c) respectively, with {v,[t]}, {va,[t]} representing the sampled theoretical velocity and acceleration data, and {n[t]} a corrupting-noise sequence, which is, for the sake of simplicity, denoted by the same symbol in all three cases. This noise sequence is subject to the following standard assumptions: Al. {n[t]} is a zero-mean second-order stationary and ergodic stochastic process characterized by rational, Hermitian, and positive definite spectra. A2. The scalar components of the noise process {n[te} may be cross-correlated, so that its covariance matrix may, in general, be non-diagonal. The approach introduced for the solution of this problem is presented in the sequel as follows: The required force excitation signal type, along with the discrete-time representation of the

structural and noise dynamics via special-form multivariate AutoRegressive Moving-Average with eXogeneous input (ARMAX) models, are discussed in subsection 2.1. The multivariate ARMAX parameter estimation algorithm is summarized in 2.2, and the model structure estimation scheme presented in subsection 2.3. Model validation is discussed in 2.4, and the model transformation, reduction, and analysis procedures are finally outlined in subsection 2.5. 2.1 The Excitation Signal Type and the Required Stochastic Model Forms for the Combined Representation of the Structural and Noise Dynamics In this subsection the form of the required multivariate ARMAX representation of the stochastic structural dynamics is examined under specified force excitation conditions, and separately for each one of the vibration displacement, velocity, and acceleration measurement cases. 2.1.1 The Displacement Data Measurement Case By taking the Laplace Transform of Eq.(1) the following expression is obtained: V(s) = [M s2 + C - s + K]-1 - F(s) G(s)~ F(s) (3) in which s represents the Laplace Transform variable and G0(s) the indicated receptance transfer matrix. By combining the discrete-time "equivalent" of this with Eq.(2a) the following expression is obtained: y[t] = GO(B) ~ f[t] + n[t] (4) in which the signals {f[t]}, {y[t]}, and {n[t]} are of dimensions m, s, and s, respectively (m, s < n). In Eq.(4) G0(B) is the discrete transfer matrix corresponding to G0(s), and B the backshift operator defined such that B. y[t] y[t - 1]. Based on Al, and by invoking the spectral factorization theorem [5], the system equation (4) may be put into the multivariate ARMAX(na, nb, nc) form: A~(B)~ y[t] = B~(B)~f[t] + CO(B)~ w[t] (5) in which {w[t]} represents an unobservable s-dimensional zero-mean white noise (innovations) process which is independent of {f[t]} and has a nonsingular (and, in general, non-diagonal) covariance matrix E~. The autoregressive (AR), exogeneous (X), and moving-average (MA) polynomial matrices A~(B), B~(B), and C~(B), respectively, are of the forms: A~(B) I, + A~(1). B + + A~(na) ~ Bn [s x s] (6) B~(B) - B~(0)+B~(1).B +. +B~(nb) - Bnb [s x m] (7) C~(B) A Is+ CO(1). B +..+C~(nc). BnC [sxs] (8) 6~~~~~~~~~~~~~8

with I, standing for the identity matrix of order s and the quantities in the brackets indicating matrix dimensions. The stability of the original structural system (l) imposes the following condition on the ARMAX system (5): Cl. All zeros of the determinant of A0(B) lie outside the unit circle ( "stability condition"). In addition the following standard assumptions are imposed on (5): A3. All zeros of the determinant of C0(B) lie outside the unit circle ( "invertibility assumption"). A4. The Smith form of [A0(B),B0(B), C0(B)] is [I, 0], or, equivalently, the greatest common left factor of A0(B), B0(B) and C0(B) is a unimodular matrix (i.e. a matrix with nonzero constant determinant). Since the mapping between the continuous-time transfer matrix G0(s) and its discrete counterpart G"(B) is not one-to-one but depends on the selected excitation signal type [6], the force excitation signals are selected to satisfy the following condition: C2. The force excitation signals consist of trains of "impulsive" functions occurring at the sampling instants. Based on C2, and by following a procedure similar to that in Ben Mrad and Fassois [6], the transfer matrices G0(s) and G0(B) may be shown to be related through the transformation expression: G~(s) (Rk + R;). - -(Rk. lI nA +. In AX). T-1 s= + 2(kwnkS + wnk ~ Go(B) E - (Rk + A*)- B +.k + )' B2 k= 1 (A k +A)-. k.. in which n represents the number of degrees of freedom, (Rk, R*) a pair of complex conjugate residue matrices of the discrete-time transfer matrix, (Ak, A*) a pair of complex conjugate eigenvalues of the discrete transfer matrix, Wnk, 4k, the k-th mode natural frequency and damping factor, respectively, and T the selected sampling period. The required special-form ARMAX -type representation of the structural and noise dynamics may be now obtained by comparing the discrete transfer matrix G~(B) in (9) to that implied by (5), on one hand, and, its continuous counterparts in (9) and (3), on the other. Indeed, the first part of (9) implies that the discrete-time transfer matrix can be rewritten as: G(B D(O)+ D(1). B +. + D(2n - 1) B2n-1 A D(B) () ) ld(l). B + + d(2n). B2n d(B)

with: n D(0)- E(Rk + R~) (11) k=1 where d(B) represents the discrete equivalent of the structural system's characteristic polynomial. By comparing the ARMAX representation (5) to the transfer matrix expression (10), the following relationship is obtained: G~(B) = [AO(B)]-1. B(B) = adj[A~(B)] B~(B) (12) detlA~(B)I in which detllI indicates determinant and adj[-] the adjoint of the indicated matrix. By further comparing (10) to (12) we obtain: na s = 2n (13) na.(s-1)+nb = 2n-1 (14) from which: nb = na- 1 (15) In addition, the continuous-time receptance transfer matrix G0(s) may be expressed as: G~(s) = [Ms2.+C s+K]-1 -o * s2n-2 +... + A2n-2 =o - S (16) 32n + a,1. 2n-1 +. +.. (+ 2n Based on the second part of (9) G~(s) may be also written as: -' 2n-1 +''' + d2n-2 Go (s) = (17) G~(s)=s + 2+ S2 s2-1+ -. + a2n with: n A-1 - (Rk + Rk) (18) k=l From Eqs.(16) and (17) it is obvious that: -,:- A -0 (19) and therefore [compare (11) to (18)]: D(0) _ 0 (20) The monic nature of the AR polynomial (6), the fact that D(0) = o0, and (13), (15) impose the following conditions on (5):

C3. The required multivariate ARMAX representation (5) is characterized by unit time-delay, that is B~(O) = 0. C4. The AR order na of the required s-variate ARMAX representation is such that na. s = 2n. C5. The X order nb of the required s-variate ARMAX representation is equal to the AR order reduced by one, that is nb = na - 1. For the proper representation of the structural and noise dynamics the special-form of the multivariate ARMAX representation (5) that satisfies the assumptions A3 and A4, while also conforming to conditions C1-C5, is therefore required in this case. Remark: For the special case in which an n-variate (n being the number of degrees of freedom) ARMAX form is selected (which means that the number of vibration measurement locations is equal to n), Eqs.(13) and (15) imply that na = 2 and nb = 1, so that an ARMAX(2,l,nc) model is appropriate for the structural and noise dynamics representation. For reasons that will be discussed in later sections, however, slightly higher AR and X orders may be sometimes necessary. 0 2.1.2 The Velocity Data Measurement Case Consider, now, the case of vibration velocity measurements. The mobility transfer matrix G'(s), which represents the relationship between the force excitation F(s) and the velocity vibration response Vv(s), may be expressed as [see (3) and (16)]: Gv(s) = s G~(s) 0do s2n-1 +... + A2n-2. S s2n + q. s2n-1 +.. + 2n (Rk'k + R,') - (R,'k ln A +-R' 1lnAk) T-' (21) 2E 2 (21) k=1 S2 + 2~kwnk S +,,nk The transformation expression (9) then implies that the equivalent discrete-time transfer matrix G~(B) may be obtained as: Go(B) (R'k + R')- (R'k. A + R*. Ak). B k=l 1 k B2 D'(O) + D'(1). B +..- + D'(2n - 1) B2n-1 l+d(l) B+ ~ +d(2n)-B2n = [AO(B)]-1 BO(B) (22) with A~(B) and B~(B) representing the AR and X polynomial matrices, respectively, of the appropriate ARMAX(na, nb, nc) representation: A~(B) ~ yv[t] = B~(B) f[t] + CO(B) w[t] (23)

that relates the sampled force signal {f[t]} and the innovations {w[t]} to the sampled and noisecorrupted velocity signal {yv[t]}. The polynomial matrices A~(B), B~(B), and C~(B) are of the forms (6)-(8), respectively, and also satisfy assumptions A3 and A4, and conditions C1, C2, C4, C5. D'(O) is, however, not null any more, and hence condition C3 needs to modified to: C3'. The multivariate ARMAX representation (23) is zero-delay, that is B~(O) y# 0. 2.1.3 The Acceleration Data Measurement Case For the case of vibration acceleration measurements, the inertance (or accelerance) transfer matrix G(s), which represents the relationship between the force excitation F(s) and the acceleration vibration response Va(s), may be expressed as [see (3) and (16)]: G,(s) = s2 Go(s) Ado S2n + _ * +,2n-2' S2 2n +,1.s2n-1+ +2n A'. 2n-1 + + n- o +4 0n s2n+ _a ~ s2n-1 +-* + a2n =o + d(R + R"*) - (RR"I ln + R"*.ln Ak) T-' kc=1 ~ s2 + 2(kWnks + Wnk (24) Based on C2 and the transformation expression (9), the equivalent discrete-time transfer matrix GO(B) may be obtained as: G"'(B) = a+,(k+Rk (R"k +R") -(R'k- A + R"; Ak) B a I -k 1(Ak + A*) B + Ak. A*B2 k==l D"(0) + D"(1) - B + *.. + D"(2n - 1). B2no 1l+d(1).B+..- +d(2n).B2n D"'(O) + D"'(1) * B + -.. + D"'(2n). B2n 1+d(1) B+... +d(2n).B2n = [A~(B)]-1. B~(B) (25) with A~(B) and B~(B) representing the AR and X polynomial matrices, respectively, of the appropriate ARMAX(na, nb, nc) representation: A~(B) y,[t] = B~(B) f[t] + CO(B)~ w[t] (26) that relates the sampled force excitation signal {f[t]} and the innovations {w[t]} to the sampled and noise-corrupted acceleration signal {y,[t]}. 10

The polynomial matrices A"(B), BW(B), and CG(B) are of the forms (6)-(8), respectively, while also satisfying assumptions A3 and A4, and conditions C1, C2, C3', and C4. Based on (25) we, however, now have: na*s = 2.n (27) na.(s-l)+nb = 2 n (28) from which: nb = na (29) This implies that condition C5 should be altered to: C5'. The X order nb of the required s-variate ARMAX representation is equal to the AR order, that is nb = na. The following special case is of particular importance as it is often encountered in practice: Special Case: Diagonal Mass Matrix and Transfer Inertance Functions In this case the transfer receptance G.(s), i $ j, of G0(s) in Eq.(16) may be expressed as: Gob(s) _ Y(s) Aijl s2n-3 + + Aij(2n-2) i # j (30) =S Fj(s) s2n+1a s2n- + + a2n with,ijl representing the ij-th element of the matrix Al. The transfer inertance G~,j(s), i 0 j, then becomes: G\ Yai(S) Aijl, 52n-1 +... + zLij(2n_2) S2 aij S- Fj(s) $2n + a1. S2n-1 + + at2n (R"'ijk + R"'j* ) s - (R"'ijk In A* + R *k ln Ak) T-31) -=1 Z2+2(kWnk n+W ik Based on C2 and the transformation expression (9) the equivalent discrete-time transfer function G~ (B) then is: Gn (R"'1 ij + R"''*) (Rijk A + Rt*k Ak) B aGj() = + 1(Ak+A)B+Ak) AB* B2 k=1 Di'"'(0) + Dii'"'(1) B +... + Dii"1(2n - 1) B2n-1 l+d(l).B+.- +d(2n)B2n [A~(B)]-1 B~(B)Ij (32) This expression implies that the ij-th element of the polynomial matrix B~(B) is of degree nb- 1 = na- 1 for i $ j. 11

2.2 Model Parameter Estimation Based on the previous discussion, once the numbers of excitation and vibration measurement locations (m and s, respectively) have been selected and the data collected, an s-variate ARMAX model needs to be estimated within the model structure [compare with (5)]: M: A(B, 0) y[t] = B(B, 0) f[t] + C(B, 0) ett/0] E([t/0] eT[tl/]} = E(8) (33) In this expression the 4-tuple 0 (A,B,C,E) denotes the complete set of parameters to be estimated, {y[t]} any one of the vibration displacement, velocity, or acceleration signals, {e[t/8]} the one-step-ahead prediction error, while the autoregressive A(B, 0), exogeneous B(B,8), and moving-average C(B, 8) polynomial matrices are of the form (6), (7), and (8), respectively. Despite their phenomenological resemblance to their univariate counterparts, multivariate ARMAX models have a much richer structure and are much more difficult to estimate. Indeed, multivariate ARMAX estimation methods are not only much more prone to difficulties related to local extrema/wrong convergence, algorithmic instabilities, excessive computational complexity, and the need for good initial guess parameter values, but also the estimation of a model within the model structure M is itself an ill-posed problem because of the lack of identifiability [8,9]. This is due to the fact that an ARMAX representation of the form (5), or (23), or (26), and subject to the corresponding assumptions and conditions, is not uniquely defined, in the sense that many different such models representing the same dynamics exist, and in order to obtain a unique representation additional constraints have to be imposed. These constraints are related to the internal structure of the AR, X, and MA polynomial matrices and the innovations covariance matrix, and lead to various types of identifiable parametrizations [9]. Within the context of the proposed approach the fully-parametrized ARMAX pseudo-canonical form is used, which implies that: C6. The AR, X, and MA are full polynomial matrices and the innovations covariance is a full matrix. and identifiability is then warranted under the additional assumption [10]: A5. Rank[A~(na) B~(nb) C~(nc)] = s. This selection (and the imposed assumption A5) does not represent any real practical limitation, whereas it significantly simplifies the structure estimation problem in the sense that only three structural indices, namely the orders na, nb, and nc, need to be determined [7,8].

The Estimation Algorithm The parameters of a model within the fully parametrized ARMAX pseudo-canonical form are estimated by the Linear Multi-Stage algorithm recently introduced by Fassois and Lee [7,8]. This algorithm is a novel prediction-error-type multi-stage scheme that overcomes many of the-drawbacks that render alternative techniques unrealistic in many practical applications, and is based on the replacement of the original nonquadratic Prediction Error (PE) estimation problem by an appropriate sequence of PE and/or linear subproblems for which uniquely-determined closed-form solutions may be obtained. A few brief remarks on the special features and advantages of the LMS algorithm over alternative schemes are presented at the end of this subsection; for a more detailed discussion the interested reader is however referred to [7,8]. A brief summary of the stages of the LMS algorithm, including all the equations needed for its implementation, is given in the sequel: Stage 1: Estimation of a truncated ARX representation. A multivariate AutoRegressive with eXogeneous input (ARX) model of the form: p p Ml E> Hy(j, t). y[t - j] = E Hf(j, 7). f[t - j] + e[t/lt] (34) j=o j=o in which H (O) I,, E[t/li] represents the one-step-ahead prediction error associated with this model, and 7 E R[sx(p.s+p.m+m)] the matrix:; -_- [ Hy(1)Hy(2) -. Hv(p)' Hf(O) Hf(1)... Hf(P) (35) is estimated. This is achieved by estimating each one of the rows of 7- separately, through the following linear least-squares estimator: h(i) = u[t U] y[t] (i = 1,2, *..'s) (36) In this expression the (ps + pm + m)-dimensional vector hT(i) represents the i-th row of'I, yi[t] the i-th element of the output vector y[t] = [yl[t].* yi[t]... ys[t]]T, and u[t] -yT[t — 1] yT[t-2].yT[t _- p ] fT[t] fT[t - 1].. - fT[t- p] (37) Stage 2: Initial estimation of the MA polynomial matrix C~(B). An initial estimate of the MA polynomial matrix is obtained by solvineg the following block Toeplitz and Hermitian linear system of equations:

RH(01)RH(-l)... RH(1-n) T(1) RH(1) RH(1) RH(O) - RH(2- nc) | T(2) RH(2) ( (38) RH(nc- 1) RH(n - 2)... RH(O) dIL dT(nc RI(nc) ( in which RH(i) is calculated as: EZPi=I-) Hy(i) Hy (i + kj) RH( - j) (k > j) RH(k,j) = (39),p-(j-k) HIy(i + j - k). (i) = RH(k - j) -i=l Y (k < j) with ({H(i)} representing the impulse response matrix estimates obtained in Stage 1. 0 Stage 3: Estimation of the AR and X matrices. The AR and X matrices are obtained by estimating a model within the model structure M2 defined as: nb na M2' YF[t] = ZFF[t - j]' col[B(j, 02)]- YF[t- j]- col[A(j, 02)] + e'[t/02] (40) j=O j=1 with the col[.] operator transforming a matrix into a vector by stacking its columns one underneath the other, E'[t/02] denoting the associated one-step-ahead prediction error, and: yF[t] - [yT[t] 0 C'-(B)] col[Is] [s x 1] (41) YF[t - j] yT[t - j] C(B) [S x 2] (42) FF[t - j] _ fT[t- -j] -i(B) [s X m s] (43) 02 - col[A(1)... A(na)' B(0).. B(nb)] e 2[(2'n'+m's'(nb+l))xl] (44) with 0 indicating Kronecker product [11] and C(B) the available estimate of the MA polynomial matrix. Since the model structure M2 is linear in 02, minimization of the Prediction Error criterion defined by the trace of the sample covariance matrix of the one-step-ahead prediction errors leads to the linear least-squares estimator: 02 = [ UF[t] ] *[ ~ UTF[t] F] (A45) in which the matrix UF[t] E sx(s2.na+m.S."nb+m.s) is defined as: UF[t] _ [-YF[t-1]... YF[t-na] FF[t]... FF[t-nb]] (46) UF~~~~~t]~1

In implementing this estimator the filtering operations of Eqs.(41)-(43) may be performed by decomposing the matrix YF[t] as: YF[t] [ Yl [t] Y2[t] ]- Y[tl ] (47) with Yi[t] E RSXS being the i-th block element of YF[t], and then rewritting (42) as: EYi[t - j] C(j) = yi[t] (i = 1,2,...s) (48) j=o with C(O) Is. From this expression the sequences {Yi[t]}N (i = l,...,s) can be obtained by assuming that Yi[t]- 0 for t < 0 (i = 1,2,...,s). {yF[t]}N l thus becomes also available, whereas the sequence {FF[t]}Y=1 is similarly obtained by assuming that FF[t]- 0 for t < 0. Stage 4: Estimation of the MA matrices and the innovations covariance. After the AR and X matrices have been estimated, the MA matrix estimates are updated by using the following estimator expression:, (j) y(i - j) = A(i) (i = 0,1,2,...,nc) (49) j=o in which A(i) _ O for i > na and the sequence {Hy(i)} is available from Stage 1. The prediction error sequence {e[t/01]} with: 1 col[ A(1) A(2) A(2) A(na)' B(O) B(1)... B(nb) C(1) C(2)... C(nc)] (50) is then estimated by using the model expression (33) in which A(B, 0), B(B, 0), and C(B,0) have been replaced by their respective estimates A(B), B(B), and C(B). The innovations covariance is finally estimated as: 1 N = N Ee[t/6] eT[t] (51) t=l Remarks: (a) Unlike alternative schemes, this algorithm is characterized by mathematically guaranteed stability! This is due to the fact that the initial MA estimator (38) is guaranteed to yield a strictly minimum phase Co(B) polynomial matrix [8]. (b) The LMS estimation algorithm is characterized by modest computational complexity, requires no initial guess parameter values, and offers a uniquely-determined estimate for a given data set.

(c) Although statistical decision theory criteria may be used for the determination of an appropriate order p for the ARX model (34), it has been found that, at least for modal analysis problems, a simple rule of thumb is to select p in the range p' (2.5 e 5) x max(na, nb, nc). (d) Stages 3 and 4 of the estimation algorithm may be iterated until the resulting sequence of estimates of ~O converges. 2.3 Model Structure Estimation As it has been already indicated, the use of the fully-parametrized ARMAX pseudo-canonical form simplifies the structure estimation problem considerably and reduces the problem to that of estimating the AR, X, and MA orders na, nb, and nc, respectively. In experimental modal analysis, in particular, the problem is further simplified by the special form of the multivariate ARMAX models used, which additionally satisfy conditions C3 through C5, or their alternative forms. Indeed, in the special case in which the number of degrees of freedom is a-priori known and the number of output measurement locations is selected equal to that (s = n), an n-variate ARMAX(2, 1, nc) (for displacement or velocity vibration data) or ARMAX(2,2, nc) (for acceleration vibration data) with appropriate time-delays (see subsection 2.1) is theoretically adequate. If, on the other hand, the number of vibration measurement locations is different from the number of degrees of freedom, an s-variate ARMAX(na, na - 1, nc) (for displacement or velocity vibration data), or ARMAX(na, na, nc) (for acceleration vibration data), model with na selected such that na x s = 2 x n, is theoretically required. In practice, however, the number of degrees of freedom is often unknown, and, in addition, slight overdetermination may be necessary due to a number of reasons, including the fact that the dynamics of imperfect instruments are present in actual experimental data, and, therefore, have to be properly accounted for. A systematic and effective strategy for the estimation of the required orders of an adequate s-variate ARMAX representation is therefore necessary. Within the context of the proposed approach the following procedure is followed: A sequence of ARMAX(k, k-l, k) models (vibration displacement or velocity measurement case) or ARMAX(k, k, k) models (vibration acceleration measurement case) with unit (vibration displacement measurement case) or zero (vibration velocity or acceleration measurement case) time-delay are successively fitted by increasing the value of k, until a statistically adequate representation is achieved for, say, k = k*. In this process the values of the AR order k are selected such that the product k x s is even (see C4). A sequence of ARMAX(k*, k* - 1,1) (vibration displacement or velocity measurement

case) or ARMAX(k*, k*, I) (vibration acceleration measurement case) models is subsequently fitted for I = k* - I,k* - 2,.., and, also, I = k* + 1,k* + 2,..-, and the "best" representation finally selected for, say, I = l*. It is remarked that in this procedure the MA order is initially equal to the AR since the resulting noise model then has the flexibility of representing a number of stochastic processes, including "white noise". The exact required value of the MA order needs to be, however, accurately determined, and this is subsequently done by examining MA orders that are lower (initially) and higher (subsequently) than the AR. This procedure is motivated by experimental evidence with structural systems which has indicated that for low noise levels the required MA order is often smaller or equal to the AR one. Statistical adequacy is judged by means of appropriate statistical decision theory criteria, and in this context the Bayesian Information Criterion (BIC) criterion [12] according to which the ARMAX(na, nb, nc) model that minimizes the index: BIC(8) = log detlE(9)l + d x logN (52) is selected as adequate, is used. In this expression 0 represents the estimated parameter vector, d the number of estimated scalar parameters, and N the length of the data records used in the estimation. Tests of the form (54) (see subsection 2.5) that use the notion of dispersion percentages are also used in diagnosing order overdetermination, and therefore indirectly judging model adequacy. 2.4 Model Validation The finally selected multivariate ARMAX model is accepted as an accurate discrete-time representation of the structural and noise dynamics after successfully passing the following two-phase validation procedure: Phase 1: Since an accurate representation of the system under test should be characterized by an uncorrelated innovations sequence, the uncorrelatedness of the estimated model's innovations is examined through a standard statistical test [13]. 0 Phase 2: The predictive ability of the model is evaluated within a data set that is referred to as the validation set, and which is different from the estimation set used for model estimation. For this purpose a standard predictor corresponding to the estimated model is built, and the measured vibration signal {y[t]} is compared to the model-predicted signal {$r[t/9]}. In addition the statistic Q defined as: q - ~'~ 1 (53)

with Eyy representing the sample covariance matrix of the measured vibration signal {y[t]} and Lyi the sample cross-covariance between {y[t]} and the model-predicted vibration signal {y[t/9]}, is constructed. For an accurate model Q may be shown to tend to the identity matrix as the noise-to-signal ratio tends to zero. In practice, the closer the Q to the identity matrix the better the predictive ability of the structural model, and the validity of a given model is judged based on the magnitude of the discrepancy between Q and the identity matrix. 0 2.5 Model Transformation, Reduction, and Analysis After the estimated discrete-time model has been validated, its part representing the structural dynamics ([A(B,0)]-1 - B(B, 0)) is used for the estimation of the actual number of structural degrees of freedom and the extraction of the modal parameters. Towards this end this part of the model is first transformed from the discrete back into the continuous-time domain by using the transformation expression (9). The transformed model may, in certain cases, be overdetermined, which means that its number of degrees of freedom may be higher than that of actual structural modes. This may be due to a number of reasons, including the fact that the dynamics of the various instruments used are inevitably present in the experimental data and may be thus "captured" by the modeling procedure (see section 4). As a consequence the actual structural modes in the model need to be identified and separated from any "extraneous" modes which have to be discarded, and the model thus reduced. This is achieved by using the Dispersion Analysis methodology introduced in the first part of this paper [1]. Indeed, the relative importance (quantified in terms of dispersion percentages) of every estimated mode is evaluated, and only the modes that are found to be nonnegligible contributors to the total vibration energy are considered to be important and kept in the model. This is quantitatively achieved through a test of the form: Ilakll < (%) (54) which, once affirmative, indicates that the k-th estimated mode has negligible contribution in the energy of the measured vibration signals and can be considered as "extraneous". In the above expression lIAkll represents an appropriate norm of the k-th dispersion percentage matrix with ml-th element computed as described in [1], and e a small positive threshold value. After the actual structural modes have been thus determined complete modal information, in the form of natural frequencies, damping factors, and mode shapes, may be readily obtained. Indeed, the global structural characteristics (natural frequencies and damping factors) can be ob

tained as functions of the discrete eigenvalues (Ak, Ak) (k = 1, 2,.., n) [14], whereas the k-th mode shape is given as [see Eq.(9)]: Ri2k Rink T k-[ Rlk'= ~ii - (k 1, 2,...,n) (55) with R1k,..., Rink representing the estimated elements of the i-th row of the k-th residue matrix of the transfer function relating the excitation force to the vibration displacement. In addition, the relative importance of the system's vibrational modes may be quantitatively assessed based on the dispersion percentage matrices, and this information subsequently used for further structural model reduction and analysis. 3. PERFORMANCE EVALUATION WITH SIMULATED STRUCTURAL SYSTENMS In this section the performance of the proposed MIMO Experimental Modal Analysis approach is evaluated with numerically simulated structural systems. The excitation force signals used consist of pseudo-random vector Gaussian processes characterized by zero means and approximately flat spectra. The resulting vibration displacement signals are obtained by integrating the equations of motion by using the Wilson-O method, and subsequent sampling with the indicated periods. Noisecorrupted versions of these signals are then obtained by adding pseudo-random Gaussian noise that is zero-mean, uncross-correlated with the excitations, and characterized by a general non-diagonal covariance matrix. The N/S ratio is defined as the standard deviation of the scalar corrupting noise over that of the noise-free vibration signal, and is kept the same in all output measurement channels. Unless otherwise stated, data sets consisting of 1000 samples each and divided into estimation (each one consisting of 900 samples) and validation (each one consisting of 100 samples) subsets are used in each case. In order to investigate the statistical characteristics of the proposed approach, Monte Carlo simulations are also performed. For each Monte Carlo analysis twenty data records, corresponding to the same structural system and N/S ratio, but obtained with different pseudo-random force excitation signals, are used. Each such set is analyzed by the proposed approach, and the sample mean values and standard deviations of the modal parameter estimates are subsequently computed. In all test cases the discrete parameter estimation algorithm is allowed to iterate, with the result corresponding to the minimum estimated innovations covariance trace selected as best. In addition, comparisons with the classical Frequency Domain Method and the deterministic time-domain Eigensystem Realization Algorithm [4] are presented. The purpose of these comparisons is to illustrate some of the limitations and difficulties of the Frequency Domain

Method and deterministic schemes. It is, however, noted that the comparisons with the ERA are not formal, as the latter operates on either free or impulse response data while the proposed approach uses forced vibration data. The main objective of the simulation work of this section is the evaluation of the performance characteristics of the proposed approach, with special emphasis placed on issues such as the: (a) achievable estimation accuracy at different N/S ratios, (b) achievable resolution (ability to distinguish neighboring vibrational modes) at different N/S ratios, (c) effects of data record length on the achievable accuracy, (d) model structure estimation effectiveness and assessment of the need for model order overdetermination, (e) effectiveness of the two-phase model validation and the model transformation/modal parameter extraction procedures. 3.1 Test Case I: Two degree-of-freedom system with closely-spaced modes In this case a two degree-of-freedom structural system of the form used in the first part of the paper [1], but with parameters selected such that its vibrational modes are closely spaced (Table 1), is identified from two-excitation two-response vibration displacement data records. The emphasis in this case is on the achievable accuracy and resolution, and hence the ability of the proposed approach to distinguish between the two modes when operating on data records corrupted by autocorrelated random noise at three different levels (N/S = 1 %, 10 %, and 32 %). In each N/S case an ARX model of order p=10 was used in the estimation [see Eq.(34)], and a two-dimensional ARMAX(2,1,4) model with unit time-delay was estimated as statistically adequate based on the procedure of subsection 2.3. The estimated adequate models were subsequently validated by using the procedure of subsection 2.4, and part of those results corresponding to the 10 % N/S ratio case are depicted in Figure 1, where the estimated model-based one-step-ahead predictions (within the validation subset consisting of 100 data samples) are observed to be, for both response signals, extremely close to their actual values. The estimated Q matrix is shown in Table 2 for each N/S ratio case, and is quite close to the identity matrix for the 1 % and 10 % N/S cases, while deviating from it in the 32 % case; a fact solely due to the relatively high level of noise, and not indicative of model inadequacy (see the discussion in 2.4). The ARMAX(2,1,4) models are thus accepted as accurate discrete-time representations of the structural and noise dynamics. Table 1 presents the results of a Monte Carlo analysis of the proposed approach for the 10 % N/S ratio case. The modal parameter estimates are evidently very good, with their sample mean values being very close to their theoretical counterparts (an indication of consistency) and their standard deviations being quite small. Modal analysis results for all the N/S ratio cases considered

are shown in Table 3, and are accompanied by comparisons with the deterministic ERA. The results of the proposed approach are obviously very good, as excellent accuracy and resolution are achieved in all three cases. The ERA results, obtained based on a 2 x 2 impulse response matrix and a cut-off singular value of 10-, are excellent in the low (1%) N/S ratio case, but the situation deteriorates significantly as the noise level increases (N/S = 10%), and, in the N/S = 32% case the method fails to identify one of the two modes (specifically the first one that is the weakest contributor to the vibration energy according to the Dispersion Analysis results) while providing 15 "extraneous" (in this case "computational") modes with Modal Amplitude Coherence (MAC) values greater than 0.50 (see Figure 2). The excellent performance characteristics of the proposed approach were finally verified with even shorter (N = 200, 500) data records, and such results indicating the achievable accuracy and resolution for the 10% N/S ratio case are presented in Table 4. 3.2 Test Case II: Three degree-of-freedom system In this case the three degree-of-freedom structural system of Figure 3, with physical and corresponding actual modal parameters indicated in Table 5, is identified from three-excitation threeresponse vibration displacement data records corrupted by uncorrelated noise at three different levels (N/S ratio of 1 %, 10 %, and 32 %). As the dispersion analysis results for the actual system indicate, mode 1 is dominant and accounts for more than 60% of the vibration signal energy (see Table 6 where the results corresponding to one of the system transfer functions are presented). In all cases an ARX model order of p = 10 was used in the estimation, and the modeling strategy indicated three-dimensional ARMAX(2,1,2) models as statistically adequate. The estimated adequate models were successfully validated, and the results of a Monte Carlo analysis of the proposed approach for the 10 % N/S ratio case are presented in Table 5. Once again, the sample mean values of the estimated modal parameters are very close to their theoretical counterparts, and their standard deviations quite small. Estimation results obtained from one data record at all three noise levels are presented in Table 6, where comparisons with results obtained by the stochastic Single-Input Single-Output (SISO) approach of Lee and Fassois [14] are also made. As it may be observed the results of the proposed MIMO approach are very good at all N/S ratios, and further improve those of the corresponding SISO method, especially for the modes with smaller dispersions and cases of higher (10% and 32%) N/S ratios. In Figure 4 the frequency response magnitude curves for the transfer function relating the v2 vibration measurement to the fi force excitation (Figure 3), as estimated by both the MIMO

and SISO approaches, are compared to the theoretical curve. In Figure 5 the frequency response magnitude curve of the transfer function v2/f2, as estimated by the proposed approach in the 32% N/S ratio case, is compared to its theoretical counterpart, and excellent accuracy is observed. In order to gain insight into the improvement that the proposed approach can achieve over the classical FDM method, the corresponding results obtained by the latter based on data records consisting of N = 1024 (that is slightly more than those used by the proposed approach) and N = 10240 samples each, are shown in Figure 6 (it is noted that the first and second FDM results were obtained by using a Hanning window and 7 and 20 averages, respectively). The first FDM result apparently is quite ambiguous and poor, whereas the second, although improved due to the very long data records and significantly more averaging used, it is still "wiggly" and ambiguous, especially in the low db regions, that is in the neighborhood of the valley, where several candidate frequencies seem to appear, and also before the first and after the second peaks, where a number of false candidate frequencies are also shown. Evidently none of these FDM results can match the smoothness and accuracy of that of Figure 5 that was obtained by the proposed approach, despite the fact that the latter was based on data records that were only N = 900 samples long. The performance characteristics of the proposed approach are finally examined with shorter (N = 200,500) data records corrupted at 10% N/S ratio, and excellent accuracy is observed from both the tabulated results (Table 7) and the estimated frequency response magnitude curves of the transfer function vl/f) (Figure 7). 3.3 Test Case III: Very lightly damped three degree-of-freedom system In this final case the three degree-of-freedom structural system of Figure 3, with physical parameters summarized in Table 8, is identified from three-excitation three-response vibration displacement data records corrupted by uncorrelated noise at three (1%, 10%, and 32%) N/S ratios. This case is especially interesting as the system is very lightly damped (all damping ratios are smaller than 0.2 %), and has two modes that are relatively closely spaced (at approximately 2.07 and 2.25 Hz). In all cases an ARX model order of p=10 was used in the estimation, and three-dimensional ARMAX(2,1,2) models were indicated as statistically adequate. The estimated adequate models were successfully validated, and the modal parameter estimation results corresponding to all three N/S ratios are summarized in Table 8. Evidently, the proposed approach very accurately estimates all three natural frequencies and provides very low damping ratio estimates that are very close

to their theoretical values for the lower N/S ratio case. As the N/S ratio, however, increases, the errors associated with the damping factors are increased, although the actual estimates are still very small (all smaller than 0.9%) and would be judged as very reasonable by any current standard. This behavior is to be, however, expected, as it is well-known that the damping factors of very lightly damped systems are very hard to accurately estimate in the presence of noise; a phenomenon that is, at least in- part, due to sensitivity problems associated with the discrete-tocontinuous transformation used in evaluating the modal parameters from the estimated discrete model [15]. 3.4 Discussion Based on the test case results of the previous subsections, as well as a number of additional simulations, the following general remarks may be made: (a) The proposed modal analysis approach offers high estimation accuracy and resolution from noise-corrupted and also short data records (data records consisting of as few as 200 samples were successfully used). The approach not only overcomes the well-known difficulties and limitations of the classical frequency domain and time-domain deterministic methods (such as the ERA), but also offers significant practical advantages and further improvements over earlier stochastic SISO methods such as [14]. (b) The model structure estimation scheme of the proposed approach appears to be effective in selecting the discrete model orders required for accurate modal parameter estimation. Excessive model order overdetermination, that is necessary with deterministic approaches, is avoided, and the dispersion-analysis-based methodology is very effective in assessing modal importance and distinguishing "extraneous" from "structural" modes. (c) The two-phase model validation procedure has been found to be effective in controlling model acceptance. (d) The mode shape estimation is quite accurate as long as the force excitation signals can be reasonably well approximated by trains of "impulsive" functions, so that the model form and the mode shape extraction procedure (implicit in the discrete-to-continuous model transformation) are appropriate. 4. EXPERIMENTAL RESULTS: MODAL ANALYSIS OF A THREE-SPAN BEAM In this section the proposed experimental modal analysis approach is used for the analysis of a three-span beam from laboratory data. A schematic diagram of the experimental set-up is shown in Figure 8. The beam was excited by band-limited white and zero-mean stochastic force signals

applied through an electromagnetic exciter, and measured by a force transducer. The vibration response signals were measured through accelerometers at two points on the beam, as shown in Figure 8. The frequency range of interest was selected as 0-150 Hz, and the force excitation and corresponding response signals were amplified, low-pass filtered with a cut-off frequency of 170 Hz, and sampled at a sampling frequency of 500 Hz by using a data acquisition system consisting of an IBM AT computer equipped with the DT2821 data acquisition board. Twenty data sets, each one consisting of 1000 samples with 900 used for estimation and the remaining 100 for model validation, were collected. In accordance with the discussion in subsection 2.1.3, since transfer inertance data were measured ARMAX(na, na - 1, nc) models with no time delay were considered. The model orders were then selected by using the procedure of subsection 2.3, and the corresponding results are shown in Figure 9, where an ARMAX(6,5,2) model is indicated as statistically adequate. In all estimation results the auxiliary ARX model order was selected as p = 20. The estimated model was subsequently validated by using the two-phase procedure of subsection 2.4; part of those results are, for one data set, presented in Figure 10, where the one-step-ahead predictions of both vibration signals are shown to be extremely close to their actual values. The ARMAX(6,5,2) model is thus accepted as an accurate discrete-time representation of the structural and noise dynamics. Since a bivariate (two-output) ARMAX(6,5,2) model has been selected as a discrete-time system representation, the order of the characteristic equation is 12, and up to 6 modes can be (and in fact are) estimated (Table 9). By considering the Dispersion Analysis results of Table 9, it is however evident that the dispersion percentages associated with the first, fifth, and sixth modes are very small, and almost all of the vibration signal energy is concentrated in the three remaining modes with frequencies approximately equal to 72.76, 97.46, and 139.92 Hz. The former modes are therefore treated as extraneous; a conclusion further supported by the fact that all of them are characterized by either very low or high (higher than the cut-off) frequencies and very high damping. The first mode is probably due to the accelerometer dynamics, whereas the fifth and sixth modes should be due to the dynamics of the cut-off filter and the other instruments used. Regarding the three actual (structural) modes, the Dispersion Analysis results indicate that the one at 97.46 Hz is dominant and accounts for about 70-80% of the vibration energy. In Figure 11 the frequency response curves of the transfer functions estimated by the proposed approach from one data set are compared to those obtained by the FDM based on the Fast Fourier Transform and a Hanning window. In accordance to observations made earlier (see subsection 3.2), the FDM results (especially those of the transfer function v2/f) appear quite inaccurate in the

low db region corresponding to the frequency range of 0-50 Hz. Otherwise, the two results are in good overall agreement; a fact primarily due to the rather low N/S ratio that characterized this experiment. It should be, nonetheless, noted, that in achieving this agreement, much longer data records (consisting of 10240 samples each, as opposed to 1024 with the proposed approach) still had to be used in conjunction with the FDM (where 18 data records, each consisting of 1024 samples, were in fact averaged with a 47.168% overlap factor). Finally, and in order to get a better feeling about the expected statistical variability of the results, a Monte Carlo analysis is performed based on all twenty measured data records. These results are presented in Table 10, and are evidently very good as the standard deviations of the modal parameter estimates are quite small for all three structural modes; a result that provides a quantitative answer to the question of repeatability of the test and estimation procedures. As a final remark it is mentioned that the relatively larger standard deviations associated with the estimated damping factors were expected, and are, to a certain extent, due to the lightly-damped nature of the beam (see the discussion in subsection 3.3). 5. CONCLUSION In this paper the problem of stochastic experimental modal analysis based on multiple-excitation multiple-response data was considered. The exact relationship between the actual structural and noise dynamics and their corresponding special-form multivariate ARMAX-type representation was studied for each one of the vibration displacement, velocity, and acceleration measurement cases, and an effective modal analysis approach developed. The proposed approach constitutes the first comprehensive scheme for the solution of this problem, as, unlike alternative methods, it is capable of operating on either displacement, or velocity, or acceleration vibration data records, and, in addition to overcoming the other main difficulties and limitations of current schemes, it also takes into account all major issues pertaining to the problem and attempts to minimize the various types of errors associated with the modal parameter estimates. Indeed, the proposed approach consists of a properly selected type of force excitation signals combined with mutually-compatible special-form multivariate ARMAX models, the Linear MultiStage parameter estimation algorithm, effective model structure estimation and validation procedures, as well as modal parameter extraction based on a compatible discrete-to-continuous model transformation and a physically meaningful Dispersion-Analysis-based procedure for the identification of "extraneous" modes, model reduction, and analysis. As a consequence it is characterized by the following main features and advantages:

1. Ability to operate on forced vibration data and handle arbitrary numbers of force excitation and vibration response (displacement, velocity, or acceleration) signals for effective and accurate modal analysis. 2. Unlike alternative frequency-domain and deterministic time-domain methods, the achievement of high accuracy and resolution with noise-corrupted data records. 3. Elimination of errors due to incompatibilities among the force excitation signal type, the discrete ARMAX system representation, and the modal parameter extraction procedure employed. 4. Accurate and effective estimated structural model analysis, including modal importance assessment, determination of the actual structural degrees of freedom, and modal parameter extraction. 5. Effective, and yet simplified, estimation of the appropriate discrete ARMAX system representation, so that problems requiring substantial operator judgement and training are largely eliminated. Indeed, the proposed approach requires no initial guess parameter values, is characterized by mathematically guaranteed algorithmic stability and the elimination of local extrema problems, and is equipped with a very systematic procedure for model structure estimation and model validation. 6. Modest computational complexity. The very good performance characteristics of the proposed approach were verified with a number of numerically simulated structural systems by using data records of various lengths and N/S ratios. Comparisons with the classical FDM and the deterministic ERA methods were also made, and the approach was finally used for the experimental modal analysis of a three-span beam from laboratory data.

REFERENCES 1. S.D. FASSOIS and J.E. LEE 1992 Journal of Sound and Vibration this issue. On the problem of stochastic experimental modal analysis based on multiple-excitation multiple-response data - part I: dispersion analysis of continuous-time structural systems. 2. S. HU, Y.B. CHEN and S.M. WU 1989 Proceedings of the 12th Biennial Conference on Mechanical Vibration and Noise, Montreal, Quebeck, Canada, 259-265. Multi-output modal parameter identification by vector time series modeling. 3. D. BONNECASE, M. PREVOSTO and A. BENVENISTE 1990 Proceedings of the 8th International Modal Analysis Conference, Kissimmee, Florida, 1, 382-388. Application of a multidimensional arma model to modal analysis under natural excitation. 4. J.N. JUANG and R.S. PAPPA 1985 AIAA Journal of Guidance, Control, and Dynamics 8, 620-627. An eigensystem realization algorithm for modal parameter identification and model reduction. 5. K.J. AiSTROM 1970 Introduction to Stochastic Control Theory. Academic Press. 6. R. BEN MRAD and S.D. FASSOIS 1991 ASME Journal of Vibration and Acoustics 113, 354-361. Recursive identification of vibrating structures from noise-corrupted observations, part I: identification approaches. 7. S.D. FASSOIS and J.E. LEE 1990 ASME Paper 90-WA/DSC-11. A linear multi-stage method for multivariate armax process identification, part I: the basic forms of the estimator. 8. S.D. FASSOIS and J.E. LEE 1990 ASME Paper 90-WA/DSC-12. A linear multi-stage method for multivariate armax process identification, part II: effective structure/parameter estimation and performance evaluation. 9. M. GEVERS and V. WERTZ 1987 Control and Dynamic Systems: Advances in Theory and Applications, edited by C.T. Leondes, Academic Press, 26, 35-86. Techniques for the selection of identifiable parametrizations for multivariable linear systems. 10. E.J. HANNAN 1976 Econometrica 44, 713-723. The identification and parametrization of armax and state space forms. 11. J.W. BREWER 1978 IEEE Transactions on Circuits and Systems 25, 772-781. Kronecker products and matrix calculus in system theory. 12. H. AKAIKE 1977 Applications of Statistics, edited by P.R. Krishnaiah, North-Holland Publishing Company, Amsterdam, 27-41. On entropy maximization principle. 13. S.M. PANDIT and S.M. WU 1983 Time Series and System Analysis with Applications.

John Wiley and Sons. 14. J.E. LEE and S.D. FASSOIS 1991 ASME Journal of Vibration and Acoustics, in press. Suboptimum maximum likelihood estimation of structural parameters from multiple-excitation vibration data. 15. S.D. FASSOIS, K.F. EMAN and S.M. WU 1990 ASME Journal of Dynamic Systems, Measurement, and Control 112, 1-9. Sensitivity analysis of the discrete-to-continuous dynamic system transformation.

APPENDIX: NOMENCLATURE A(B) autoregessive polynomial matrix A(i) i-th autoregressive matrix B backshift operator B(B) exogeneous polynomial matrix B(i) i-th exogeneous matrix C viscous damping matrix C(B) moving-average polynomial matrix C(i) i-th moving-average matrix d the total number of discrete estimated parameters e[t/8] (vector) model-based one-step-ahead prediction error f(t) (vector) force excitation G(s) receptance transfer matrix Gij(s) ij-th element of the receptance transfer matrix G(s) G(B) discrete-time receptance transfer matrix Ga(s) inertance transfer matrix Gaij(s) ij-th element of the inertance transfer matrix Ga(s) Ga(B) discrete-time inertance transfer matrix G,(s) mobility transfer matrix G,v(B) discrete-time mobility transfer matrix I, s x s identity matrix j imaginary unit (if not an index) K stiffness matrix M mass matrix m number of measured excitations n number of degrees of freedom N length of each sampled data record used in estimation n[t] discrete-time (vector) noise process na autoregressive order nb exogeneous order nc moving-average order p order of the ARX model used in LMS Rk residue matrix of the receptance transfer matrix R'k residue matrix of the mobility transfer matrix R"k residue matrix of the inertance transfer matrix s number of vibration measurement locations T sampling period v(t) (vector) vibration displacement 29

Va[t] discrete-time (vector) vibration acceleration vv[t] discrete-time (vector) vibration velocity y[t] discrete-time (vector) noise-corrupted vibration displacement Ya[t] discrete-time (vector) noise-corrupted vibration acceleration yV[t] discrete-time (vector) noise-corrupted vibration velocity wi[t] discrete-time (vector) white noise Ak k-th eigenvalue of the discrete-time transfer matrix (k k-th damping factor Ok k-th mode shape Wnk k-th natural frequency O the vector of the discrete model parameters E innovations covariance matrix.EXy cross-covariance between the signals {x} and {y} Ak k-th mode dispersion percentage matrix Conventions xi the i-th scalar component of the vector x x(t) the value at time t of the analog signal x x[t] the value at the discrete time t of the sampled signal x (x(t)} the signal x X(s) Laplace Transform of {x(t)} (superscript) denotes complex conjugate o (superscript) denotes quantities associated with the true system (not a model) denotes estimator/estimate F (subscript) indicates filtered quantity capital bold-face denotes matrix quantity lower-case bold-face denotes vector quantity Abbreviations AR Autoregressive ARMAX Autoregressive Moving-Average with exogeneous input(s) ARMAX(na, nb, nc) ARMAX model of AR, X, and MA orders na, nb, and nc, respectively ARX Autoregressive with exogeneous input(s) BIC Bayesian Information Criterion ERA Eigensystem Realization Algorithm FDM Frequency Domain Method LMS Linear Multi-Stage

MA Moving-Average MAC Modal Amplitude Coherence MIMO Multiple-Input Multiple-Output N/S noise-to-signal X exogeneous PE Prediction Error SISO Single-Input Single-Output TF Transfer Function 31

Theoretical Estimated Theoretical Estimated Parameters Parameters Parameters Parameters Natural 9.9793 9.9745 ~0.0118 Damping 0.1826 0.1818-0.0018 Frequencies 9.9274 9.9210 -0.0084 Factors 0.0480 0.0475~0.0006 (Hz) Theoretical Parameters Estimated Parameters (1.0, jO.O) (1.0, jO.0) Mode (-0.6576, jO.0117) (-0.6574~0.0233, jO.0145~0.0122) Shapes (1.0, jO.0) (1.0, jO.0) (1.5164, jO.0274) (1.5154~0.0165, jO.0301~0.0178) Sampling Interval: 0.025 sec m1 = 4.5 kg, m2 = 4.5 kg cl = 45 N-s/m, c2 = 35 N.s/m, c3 = 15 N-s/m k1 = 17500 N/m, k2 = 100 N/m, k3= 17500 N/m Table 1. Monte Carlo results obtained by the proposed approach. (Test case I; N/S = 10%). N/S = N/S= 10 % N/SS = 32 % Q 0.99992 0.00010 0.98637 0.00885 0.89017 0.05726 0.00001 0.99986 0.00499 0.98596 0.03718 0.88610 Table 2. The estimated Q matrix (Test Case I). Theoretical Estimated Parameters N/S = 1% N/S = 10 % N/S = 32'% Parameters Proposed ERA Proposed ERA Proposed ERA Natural 9.9793 9.9739 9.9697 9.9697 10.3747 9.9630 Frequencies 9.9274 9.9206 9.9201 9.9225 9.8538 9.9256 9.9608 (Hz)..... Damping 0.1826 0.1823 0.1819 0.1822 0.1876 0.1822 Factors 0.0480 0.0479 0.0475 0.0482 0.0486 0.0487 0.0490 Dispersion -35.77 -36.25 -36.51 -37.02 Percentages (%) 135.77 136.25 - 136.51 - 137.02 - (TF v2/fl) III. Table 3. Modal parameters of the system with closely-spaced modes estimated by both the proposed and the ERA methods at different N/S ratios. (Test case I). 32

Natural Damping Dispersion Percentages Frequencies Factors (%) (Hz) _ Mode 1 Mode 2 Theoretical 9.9793 0.1826 54.44 -35.77 45.56 135.77 Parameters 9.9274 0.0480 -35.77 16.10 135.77 83.90 No. of Data: 9.9574 0.1822 52.20 -35.57 1 F 47.80 135.57 200 9.9256 0.0487 [ -35.48 16.89 135.48 83.11 No. of Data: 9.9783 0.1865 [ 52.65 -37.10 1 [ 47.35 137.10 1 500 9.9233 0.0483 -36.99 16.93 J 136.99 83.07 No. of Data: 9.9687 0.1816 54.20 -37.13 1 F 45.80 137.13 900 9.9236 0.0486 -37.09 15.90 J 137.09 84.10 J'able 4. Modal parameters estimated by the proposed approach from data records of various lengths. (Test case I; N/S = 10 %). Theoretical Estimated Theoretical Estimated Parameters Parameters Parameters Parameters Natural 1.2926 1.2935~0.0015 Damping 0.0464 0.0464~0.0015 Frequencies 2.0622 2.0627~0.0036 Factors 0.0683 0.0687~0.0023 (Hz) 2.8295 2.8283~0.0036 0.0612 0.0611~0.0013 Theoretical Parameters Estimated Parameters (1.0, jO.O) (1.0, jO.O) Mode (1.3409, -j0.0181) (1.3416~0.0073, -jO.0205~0.0097) (0.7979, jO.0101) (0.7962~0.0055, jO.0087~0.0048) (1.0o, jo.) (1.0o, jo.o) Shapes (0.3218, jO.0154) (0.3202+0.0164, j0.0191~0.0127) (-0.8963, j0.0111) (-0.8925~0.0251, jO.0191~0.0247) (1.0, jo.0) (1.0, jO.0) (-1.1628, -j0.0303) (-1.1835~0.0346, -jO.0208~0.0335) (0.3492, -j0.0131) (0.3418~0.0201, -jO.0034~0.0109) Sampling Interval: 0.0884 sec ml = 1 kg, m2 = 1 kg, m3 = 2 kg cl = 0.6 N.s/m, c2 = 0.5 N-s/m, C3 = 0.6 N.s/m, C4 = 1.5 N.s/m, c5 = 0.7 N.s/m, c6 = 0.5 N-s/m k- = 100 N/m, k2 = 100 N/m, k3 = 100 N/m, k4 = 200N/m, ks = 0 N/m, k6 = 0 N/m Table 5. Monte Carlo results obtained by the proposed approach. (Test case II; N/S = 10%). 33

Theoretical Estimated Parameters N/S= 1% N/S 10 % N/S =32 % Parameters MIMO SISO MIMO SISO MIMO SISO Natural 1.2926 1.2928 1.2928 1.2929 1.2930 1.2934 1.2891 Frequencies 2.0622 2.0613 2.0612 2.0615 2.0846 2.0726 2.2104 (Hz) 2.8295 2.8275 2.8274 2.8263 2.8345 2.8171 2.9121 Damping 0.0464 0.0464 0.0468 0.0466 0.0461 0.0475 0.0513 Factors 0.0683 0.0686 0.0689 0.0684 0.0894 0.0727 0.0856 0.0612 0.0611 0.0618 0.0608 0.0630 0.0639 0.0739 Dispersion 87.78 87.75 87.77 87.51 87.60 85.22 88.33 Percentages (%) 1.72 1.69 1.70 1.98 2.53 2.89 4.35 (TF v2/fl) 10.50 10.55 10.53 10.50 9.88 11.89 7.32 Table 6. Modal parameter estimates obtained by the proposed MIMO and the stochastic SISO approach of [14] at different N/S ratios. (Test case II). Natural Damping Dispersion Frequencies Factors Percentages (%) (Hz) (TF V2/fl) Theoretical 1.2926 0.0464 87.78 Parameters 2.0622 0.0683 1.70 2.8295 0.0612 10.52 No. of Data: 1.2920 0.0468 87.97 200 2.0651 0.0706 1.71 2.8179 0.0658 10.32 No. of Data: 1.2934 0.0464 87.49 500 2.0611 0.0675 1.93 2.8216 0.0630 10.58 No. of Data: 1.2929 0.0466 87.52 900 2.0615 0.0684 1.98 2.8263 0.0608 10.50 Table 7. Modal parameters estimated by the proposed approach from data records of various lengths. (Test case II; N/S = 10 %). 34

Theoretical Estimated Parameters Parameters N/S = 1 % N/S = 10 % N/S = 32 % Natural -2.07358 2.07303 2.07559 2.07158 Frequencies 2.25055 2.25041 2.25141 2.25027 (Hz) 3.66498 3.66272 3.66560 3.62087 Damping 0.00061 0,00052 0.00085 0.00255 Factors 0.00124 0.00136 0.00132 0.00353 0.00035 0.00055 0.00663 0.00817 Dispersion 84.31 84.58 96.43 92.37 Percentages(%) -0.01 -0.01 0.12 -0.07 (TF of y2/fl) 15.70 15.43 3.45 7.70 Sampling Interval: 0.07 sec System: Figure 3 m l = 1 kg, m2 = 5 kg, m3 = 1 kg cl = 0.02 N-s/m, c2 = 0.02 N-s/m, cs = 0.02 N.s/m C4 = 0.0 N.s/m, c5 = 0.0 N.s/m, c6 = 0.0 N.s/m k1 = 100 N/m, k2 = 100 N/m, k3 = 200 N/m k4 = 0 N/m, k5 = 2200 N/m, k6 = 0 N/m Table 8. Modal parameters estimated by the proposed approach at different N/S ratios. (Test case III). Natural Damping Dispersion Percentages (%) Q Frequencies(Hz) Factors (TF v/lf) (TF v2/f) 5.5284 0.73290 -0.0410 0.0002 72.7574 0.00769 10.3254 17.0466 97.4616 0.00275 80.2673 71.4329 0.99815 -0.00030 139.9186 0.00213 9.6944 11.5192 -0.00146 0.99760 177.9703 0.47895 -0.2546 0.0009 233.5024 0.22033 0.0085 0.0002 Table 9. Modal parameters of the beam estimated by the proposed approach by using a single data record. (Beam experiment). Natural Damping Dispersion Percentages (%) Frequencies (Hz) Factors (TF vll/f) (TF v2/f) 72.9979~0.0934 0.00652~-0.00080 15.7311-~8.4643 19.4647~6.3263 97.4986~0.0434 0.00412~0.00054 72.1615~-7.1378 60.1651~9.1354 140.0033~0.0466 0.00220~0.00020 11.8184~-5.0296 19.5889~12.119 Residues (TF vl/f) (TF v2/f) Magnitude Phase Magnitude Phase 0.0850740.00521 36.28201~2.08339 0.08419+0.00589 38.20421+-2.05633 0.13867~0.00665 17.67108~0.76915 0.11896~::0.00785 19.69146~1.43482 0.06999~0.00432 -28.10695+1.14071 0.05206~0.00446 -23.87358~2.16576 Table 10. Monte Carlo results obtained by the proposed approach based on experimental beam data. (Beam experiment). 35

LIST OF FIGURES Figure 1. Validation of the estimated model: One-step-ahead predictions of the measured vibration displacement signals.: measured vibration, +: predictions. (Test case I; N/S = 10%). Figure 2. The modal amplitude coherence (MAC) values of the frequencies estimated by the ERA. (a) Theoretical; (b) N/S = 1%; (c) N/S = 10%; (d) N/S = 32 %. (Test case I). Figure 3. Three degree-of-freedom structural system. Figure 4. The frequency response magnitude curve as estimated by the proposed MIMO and the stochastic SISO approach of [14]. (a) N/S = 1 %; (b) N/S = 32 %. (Test case II; transfer function v2/fl). Figure 5. The frequency response magnitude curve of the transfer function v2/f2 as estimated by the proposed approach.: Theoretical, - - - - -: estimated. (Test case II; N/S = 32 %). Figure 6. The frequency response magnitude curve of the transfer function v2/f2. (a) Theoretical; estimated by the FDM based on N = 1,024 data points; (c) estimated by the FDM based on 10,240 data points. (Test case II; N/S = 32 %). Figure 7. The frequency response magnitude curve of vl/fl as estimated by the proposed approach from data sets of various lengths. -: Theoretical, -- - --: estimated from 900 samples, -- --: estimated from 500 samples, - -: estimated from 200 samples. (Test case II; N/S = 10 %). Figure 8. Schematic diagram of the experimental set-up. Figure 9. Model order estimation results. (Beam experiment). Figure 10. Validation of the estimated model: one-step-ahead predictions of the measured vibration signals.: measured vibration, +: predictions. (Beam experiment). Figure 11. The frequency response functions as estimated by the proposed and FDM approaches. (a) transfer function vl/fi; (b) transfer function v2/f. --: FDM,: proposed approach. (Beam experiment).

3.0 3.0 2.0 2.0 1.0 (N 1 +0 0. 0 - 0.0 )r~ I +I 4. -2. 0 -2. 0 -3.0 -3.0 900 920 940 960 980 1000 900 920 940 960 980 1000 DATA SAMPLE NUMBER DATA SAMPLE NUMBER (a) Output 1 (b) Output 2 Figure 1. Validation of the estimated model: One-step-ahead predictions of the measured vibration displacement signals.: measured vibration, +: predictions. (Test case I; N/S = 10%).

1.0 1.0 U U 0.0 0.0 0.0 10.0 20.0 0.0 10.0 FREQUENCY (Hz) FREQUENCY (Hz) (a) Theoretical (b) N/S= 1% 1.0 1.0 U I IIIU 0.0 10.0 20.0 0.0 10.0 20.0 FREQUENCY (Hz) FREQUENCY (Hz) (c) N/S = 10% (d) N/S = 32% Figure 2. The modal amplitude coherence (MAC) values of the frequencies estimated by the ERA. (a) Theoretical; (b) N/S = 1%; (c) N/S = 10%; (d) N/S = 32 %. (Test case I).

1a4 I,, A,,,<'A I "'""-f f 3 V k 6|v ~ ~~~~~~~~~~~~~~~~ 3 I, 3 1 WJs C 6 k 2 k3 k 4 m, I I m,-"n, ~ j c5 C2 -2 C3 C4 k s Figure 3. Three degree-of-freedom structural system.

FRce L4 0.0 0.0'"-20. 0' -20. 0 -D 40.0 -40. 0 -60.0 -60.0 Theoretical Theoretical - - Estimated by SISO -—: Estimated by MIMO -80.0 80.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 FREQUENCY (Hz) FREQUENCY (Hz) (a) N/S = 1% 0.0 0.0 -20.0 -20. 0 -40.0 40.0 -60.0 -60.0': Theoretical ~: Theoretical -- -: Estimated by SISO --: Estimated by MIMO -80.0 o', -80.0' 0.0 1.0 2.0 3.0 4.0 5.0 6.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 FREQUENCY (Hz) FREQUENCY (Hz) (b) N/S = 32% Figure 4. The frequency response magnitude curve as estimated by the proposed MIMO and the stochastic SISO approach of [14]. (a) N/S = 1 %; (b) N/S = 32 %. (Test case II; transfer function v2/fi).

I O. O -20.0 C, -30.0 40o. OA 0.0 1.0 2.0 3.0 4.0 o.0 *.0 FREQUENCY (Hz) Figure 5. The frequency response magnitude curve of the transfer function v2/f2 as estimated by the proposed approach. Theoretical,- - - estimated. (Test case II; N/S = 32 %).

0.0 10.0 20. m -20.0 -30.0 -40.0.'., I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 FREQUENCY (Hz) (a) 0.0 0.0 -10.0 — %-I 0. 0 -~ -10 0 I A O.|O -20.0:-20.0 -40.00.0 -50.0 -40.0', 0.0 1.0 2.0 3.0 4.0 5.0 6.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 FREQUENCY (Hz) FREQUENCY (Hz) (b) (c) Figure 6. The frequency response magnitude curve of the transfer function v2/f2. (a) Theoretical; estimated by the FDM based on N = 1,024 data points; (c) estimated by the FDM based on 10,240 data points. (Test case II; N/S = 32 %).

F. 0.0: Theoretical -- --: Estimated with 900 samples..- -~: Estimated with 500 samples - ---- Estimated with 200 samples -10. 0 -20.0 z <-30.0 -40.0 *,,, 0.0 1.0 2.0 3.0 4.0 5.0 6.0 FREQUENCY (Hz) Figure 7. The frequency response magnitude curve of vll/f as estimated by the proposed approach from data sets of various lengths. ---: Theoretical, — - —: estimated from 900 samples, -- —: estimated from 500 samples, - - - - -: estimated from 200 samples. (Test case II; N/S = 10 %).

Monitoring Fccelerometers Exciter Power Filter Signal - Amplifier Generator x exForce Transducer Beam Figure 8. Schematic diagram of the experimental set-up.

-4 ARMAX(l, k - l, k -6 _ -8 -10 -12 (3,2,3) (4,3,4) (5,4,5) (6,5,6) (7,6,7) MODEL ORDER -10.,~, [\#RMXAX(* k*-1l I) i-12 - -12 -13 -14 (6,5,6) (6,5,5) (6,5,4) (6,5,3) (6,5,2) (6,5,1) MODEL ORDER Figure 9. Model order estimation results. (Beam experiment).

.0.0 3.0 a_ 3.0. -2. 0 -2.0 -3.0 -1 - 3.0 900 920 940 960 980 1000 900 920 940 980 980 1000 DATA SAMPLE NUMBER DATA SAMPLE NUMBER (a) Output 1 (b) Output 2 Figure 10. Validation of the estimated model: one-step-ahead predictions of the measured vibration signals.: measured vibration, +: predictions. (Beam experiment).

(uamu.jadxa ureaa) q:peoiddpasodoid' — —' (Ij f/zn uo!l:)unj ijsusi (q) iTf/la uo!l:)unj ajsu'el (e) -saqnidd~e yjyjj pue pasodoid all~ q paemetui se suo~punj asuodsai X:uanbaij a'ql -IT Wud (//yf) g UO!3nUnf JTs=JtT (q) (f//t) T uoTun1 j,Jsl3T. (e) (zH) AON3nO3l J (zH) A3N3nf3l J 9,'OOG'00'00'o 0'O'0'OO'000'O'0 o,08-'09~~~I I~~09-' O~~~~~~~~~~~~~~~~~~'~r II ~~~~~~~C) C ~t,- z' ~ I,! z,0~ i - ~~~~~~~~~~~~Orm r 0 crs" t~~~ r? -'O'O0 ]'Or (ZH) AON3"lO3bJ (ZH) ADN3no3bJ'0o1'OOr'ogI'00'0og'C Oot'00C'ogI' oot'og 0 1111 \.....00~-..........~ -'OL'... _...,.....,...-,- Oor-..: 0OLQ g.,' 9C ___.:___._____-___ _!:9 t o tset r,~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 09 0909'0~ —'":' — I C J0 t'o cr,~~~~~~~~~~~~~~~~~~ ~.~~~~~~~~~u':'9 - "- o;'0 ) s -'~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~11~i~.~,

UNIVERSITY OF MICHIGAN 39015 02229 1473