Experiences with the Eigensystem Realization Algorithm

 

 

V.H. Balden^I; G.N. Nurick^II

^IDepartment of Mechanical Engineering, University of Cape Town, Private Bag, Rondebosch, 7700 South Africa
^IIDepartment of Mechanical Engineering, University of Cape Town

 

 

---

ABSTRACT

The Eigensystem Realization Algorithm (ERA) is a multi input / multi output time domain implementation of system realization theory. The input to the ERA is impulse response functions, typically from vibration analysis of a structure. The output of the ERA is realized state-space matrices that reconstruct the impulse response functions and contain the excited structure's modal characteristics. This article presents a brief mathematical background of the ERA and discusses the application to experimental vibration data. Practical issues of the Hankel matrix block dimensions and the number of singular values to retain are discussed. It is shown that the modal parameters converge only after massive over-specification of the Hankel matrix block dimensions and that the number of retained singular values do not determine the model order. Rather, the ERA user determines the model order by selecting successfully transformed and reliable modes that accurately reconstruct both the measured impulse and frequency response functions.

---

 

 

Nomenclature

A state transmission matrix, dimensions [2m × 2m]

B input influence matrix, dimensions [2m × 2m]

B_f input force influence matrix, dimensions [m × l]

C output influence matrix, dimensions [p × 2m]

D viscous damping matrix, dimensions [m × m]

D diagonal matrix of singular values

E Young's modulus

h (k) matrix of discrete time impulse reponse functions

H_rsHankel matrix (r block rows, s block columns)

I identity matrix

j complex operator

kdiscrete-time index

K stiffness matrix, dimensions [m × m]

Inumber of applied input forces

M mass matrix, dimensions [m × m]

m physical degrees of freedom of vibrating structure

N number of retained singular values

O zero matrix

P singular value decomposition matrix

pnumber of output variables measured

Q singular value decomposition matrix

q(t)displacement vector, dimensions [m × 1]

ttime

U_f(t) force vector, dimensions [l x 1]

Vobservable matrix

W controllable matrix

x(k) discrete-time state vector, dimensions [2m × 1]

x (t) continuous-time state vector, dimensions [2m × 1]

y(t) vector of measured variable, either displacement, velocity or acceleration, dimensions [p × 1]

Λ discrete-time diagonal eigenvalue matrix with elements λ

φ discrete-time eigenvector matrix with elements ψ

v Poisson ratio

τ dummy variable

ρmaterial density

α_i integer, such that α_i = i

ß_iinteger, such that ß_i = i

ωcircular frequency

∆t constant time interval

 

Subscripts

c continuous-time matrix

d discrete-time matrix

r number of block rows of Hankel matrix

s number of block columns of Hankel matrix

 

Superscripts

Λ ERA realised matrix

+ pseudo-inverse of matrix

T matrix transpose

 

Introduction

As both the precision required of structural models and the modal complexity of the structure themselves have increased, it has become more challenging to obtain accurate estimates of modal parameters. A central part of this challenge is improving the numerical procedure with which the modal parameters are identified from the vibration data. Current research efforts have focused primarily on time-domain, state-space system realization procedures. Juang¹ showed that most of the current time-domain modal identification algorithms may be reformulated in a unified way using the framework of system realization theory. The most notable state-space system realization procedure is the Eigensystem Realization Algorithm (ERA),² and its variants ERA/DC³ and OKID.⁴

To enhance the understanding of the practical application of the ERA, a brief overview of the algorithm is given. The algorithm commences by assembling the Hankel matrix and shifted Hankel matrix, from measured discrete-time impulse response functions. By minimising the error between the shifted and non-shifted Hankel matrices, the discrete-time state transition matrix may be defined.

At this point the state transition matrix is defined in terms of the shifted Hankel matrix and pseudo-inverses of both the observable and controllable matrices.

These pseudo-inverse matrices may be solved by the singular value decomposition of the Hankel matrix. An appropriate number of the most significant singular values are retained, which approximate the pseudo-inverse of both the observable and controllable matrices. Following substitution, the state transition matrix, output and input influence matrices may be solved.

Initially, both the continuous-time and discrete-time structural vibration models are formulated in state-space format. The purpose is to define the discrete-time impulse response function in terms of the discrete-time state-space format, which is used in the formulation of the ERA.

For the purpose of practically applying the ERA, experimental vibration measurements were obtained from the impact analysis of a simple beam structure. Brief descriptions of the data acquisition and post-processing techniques, used in acquiring the impulse response functions, are included.

Two practical issues concerning the influence of the Hankel matrix block dimensions and the number of retained singular values on the realized modal parameters are discussed. Finally, the ERA reconstructed impulse and frequency response functions are compared to the measured responses.

 

Continuous-Time State-Space Form

The vibration response of a structure is usually modelled as the linear matrix differential equation of the form

The appropriate output influence matrix, either C_displ, C_velor C_accel, each of dimension [p × m],is selected. For example, if the measured output variable, y(t), were velocity, then the appropriate output influence matrix is C_vel.

The continuous state-space representation of equation (1) is

The solution for the continuous state vector, x(t),with initial condition x(t₀), has been given by Franklin et al.⁷ as

 

Discrete-Time State-Space Form

The discrete-time state vector, x(k), may be obtained by evaluating equation (3) at equally spaced discrete intervals of time, i.e. t = (k + l)∆t and t₀ = k∆t. Further, assuming that u_f(τ) is constant over the interval k∆t < τ < (k + 1)∆t and has the value u_f(k∆t), then equation (3) may be rewritten as⁷

 

Modal Parameters

It is assumed that has a complete set of linearly independent eigenvectors φ = (ψ₁,ψ₂, .... ψ_N) with corresponding eigenvalues (λ₁, λ₂,...., λ_N), which are not necessarily distinct.

Then the realized discrete-time state transition matrix, may be decomposed as = Λ = diag (λ₁, λ₂,.... λ_N).

Then the discrete-time state-space model, as defined in equation (5), may be transformed via x = φ_z. Similarly, the realized triple state-space matrices would transform to

In this form the discrete-time modal damping rates and damped natural frequencies are, respectively, the real and imaginary parts of the diagonal matrix Λ. The matrix , after transformation from the discrete to continuous time domain defines the mode shapes.

 

Implementation of the ERA

The Hanke1 block matrix, as defined in equation (9), is assembled from the measured impulse response matrix. Following the singular value decomposition of the Hankel matrix, it is necessary that the ERA user select an appropriate number of the largest singular values to retain.

Considerable savings in both computational expense and computer memory storage may be effected by noting that, the shifted Hankel block matrix may be formed by the appropriate augmenting and truncating of the assembled Hankel matrix.

The realized state-space matrices are calculated using equation (14) and following transformation the modal parameters may be extracted. The algorithm was implemented on a personal computer using the C programming language.

 

Experimental Technique

The experimental apparatus used to perform the impact analysis, for application of the ERA, is shown in Figure 1.

 

 

The apparatus consisted of a vertically suspended uniform rectangular steel beam of dimensions 48.1mm × 19.2mm × 500mm and typical material properties E = 204GPa, v = 0.3, and ρ= 7860kg/m³. Two soft rubber slings supported the beam to simulate free-free boundary conditions.

An impact hammer fitted with a force transducer (PCB Piezotronics Model 208-B03) was used to excite the beam, perpendicular to the suspension direction. Three accelerometers (Wilcoxon Research Model 736) captured the vibration response of the beam, including a co-located driving point measurement for mass normalisation of the mode shapes.

The respective output signals were passed through a low-pass filter with a band limit of 8 kHz and an attenuation of 30 dB per octave. Following sampling at a frequency of 40 KHz, 4 096 samples were stored for postprocessing.

Twenty-five captured sets of vibration data were used to calculate the ensemble average frequency response functions and coherence vectors. The frequency response functions, for each input/output pair, were calculated by computing the ratio of the cross-spectrum between the input and output to the power-spectrum of the input, as described by Halvorsen and Brown.¹⁰

The measured discrete-time impulse response matrix was calculated from the inverse Fourier transform of the frequency response functions.

 

Discussion

Experimental Results

The measured discrete-time impulse response function for accelerometer #3, as labelled in Figure 1, is given in Figure 2.

 

 

The amplitude of the impulse response function decays in an exponential envelope. After ≈1000 discrete samples or ≈25ms the amplitude reverts to a virtual steady-state response, a phenomenon known as 'ringing'. Of interest are the measured impulse response samples containing high frequency content, which occur during the exponential decay envelope. The extent of the exponential decay envelope was estimated to be the initial ≈25ms of the impulse response function as shown in Figure 3. The corresponding frequency response function and phase vector, for accelerometer #3, are shown in Figure 4. Notice the presence of weakly excited and highly attenuated modes at frequencies greater than the filter cut-off frequency of 8KHz.

 

 

 

 

Hankel Matrix Block Dimensions

The influence of the Hankel matrix block dimensions, r and 5, on the realised poles was investigated. Figure 5 shows the result of plotting, in descending order, the normalised singular values for various Hankel block dimensions.

 

 

From Figure 5, it could be suggested that the Hankel matrix of smallest block dimensions yields an 'efficient' method of retaining the most significant singular values for the least block dimensions.

However, Peterson¹¹ states that most modal parameters, specifically poles, converge only after massive over-specification of the Hankel matrix block dimensions, while some poles converge at faster rates than others.

Further, poles that have converged can occasionally split into two or more closely repeated poles as the Hankel matrix block dimensions are increased to ensure convergence of other poles.

Convergence of the high frequency poles was especially noticeable when comparing the realized poles for the r = 100, s = 100 and the r = 300, s = 1 200 cases.

This may be expected as only the initial 200 discrete samples are used in the realisation for the r = 100, s = 100 case, while 1 500 samples are used for the r = 300, s = 1 200 case. From inspection of the measured impulse response fimction, Figure 3, it would be expected that there is valuable high frequency content after 5ms or 200 samples, which would not be present in the r = 100, s = 100 realization.

The comments of Peterson and the observed convergence of the poles imply that the largest practical Hankel matrix block dimensions should be chosen. However, due to computation restrictions this may not be practically possible. It is therefore suggested that the Hankel block dimensions be selected to ensure that the sum of r and s is at least sufficient to capture those discrete impulse response samples contained in the initial exponential decay envelope.

 

Number of Retained Singular Values

Using the Hankel matrix block dimensions s = 1 200 and r = 300 the effect of varying the number of retained singular values, N in equation (14), was considered. Each retained singular value represents a prospective mode.

Typically, ≈75% of the retained singular values are eliminated as they fail to successfully transform from the discrete to the continuous-time domain or the mode shape may not be scaled to achieve a mass-normalised mode shape.

For the remaining ≈25% of singular values, a modal accuracy indicator is required to access the reliability of the prospective modes. In this study the 'modal-amplitude coherence'¹ indicator was used, although the poles were easily identified. The realized poles are listed in Table 1, for varying number of retained singular values.

 

 

The tabulated poles demonstrate the phenomenon of split-modes, as described by Peterson,¹¹ for both the fifth and sixth modes. The results also suggest that the reliable modal parameters do converge as the number of retained singular values is increased. For this reason the retained modes were taken from the N = 400 case.

A minimum order realization was obtained by retaining only the appropriate rows and columns of the triple state-space matrices [Λ, φ^-1 B_d C_d^φ], corresponding to the four poles and four split poles.

 

Reconstruction of the Measured Vibration Data

The minimum order state-space matrices were used to reconstruct both the impulse and frequency response functions for accelerometer #3, as labelled in Figure 1.

The initial ≈1000 discrete samples of the measured and reconstructed impulse response functions are shown in Figure 3. The difference between corresponding discrete samples is also plotted in Figure 3, as an indication of the error in the reconstructed response.

There is a difference for the initial ≈200 discrete samples. It is assumed that this indicates high frequency content in the measured response, which is not present in the ERA reconstruction.

The difference following the initial «200 discrete samples is small and considered negligible. Figure 4 shows good correlation between the measured and ERA reconstructed phase and frequency response function. In the filter bandwidth, the location and magnitude of the both the poles and zeros are in close agreement. The presence of weakly excited high-frequency modes, above the filter band limit may be observed, which the ERA reconstruction does not simulate. These high frequency modes are thought to be responsible for the difference in the initial «200 discrete samples of the impulse response.

 

Conclusions

The modal parameters, especially the higher frequency modes, converge only after massive over-specification of the Hankel matrix block dimensions. The Hankel block dimensions should be selected to ensure that the sum of r and s is at least sufficient to capture those discrete impulse response samples contained in the initial exponential decay envelope.

The ERA user should select the greatest possible number of singular values to retain, as the modal parameters converge as the number of retained singular values is increased. A minimum order realization is possible if the ERA user selects successfully transformed and reliable modes that accurately reconstruct both the impulse and frequency response functions. Split modes should be included in the minimum order realization.

Using the minimum order realization, the theoretical reconstruction is capable of simultaneously satisfying both the measured frequency and impulse response functions to a high degree of accuracy.

 

References

1. Juang JJ 1987. Mathematical Correlation of Modal Parameters Identification methods via System Realization Theory. Journal of Modal Analysis, 1986, pp.1-18.

2. Juang JJ & Pappa RS 1983. An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction. Journal of Guidance, Control, and Dynamics, 8, pp.620-627.         [ Links ]

3. Juang JN, Cooper JE & Wright JR 1988. An Eigensystem Realization Algorithm using Data Correlations (ERA/DC) for Modal Parameter Identification. Control Theory and Advanced Technology, 4, pp.5-14.         [ Links ]

4. Juang JN, Phan M., Hort a LG & Longman RW 1991. Identification of Observer/Kalman Filter Markov Parameters: Theory and Experiments. Proceedings of the AIAA Guidance, Navigation, and Control Conference, pp.1195-1207.

5. Ho BL & Kalman RE 1965. Effective construction of linear state variable models from Input/Output data. Proceedings of the 3rd Annual Allerton Conference on Circuit and System Theory, pp.449-459.

6. Zeiger HP & McEwen AJ 1974. Approximate linear realizations of given dimension via Ho's algorithm. IEEE Trans. Automatic Control, AC19, 1974, p.153.         [ Links ]

7. Franklin GF, Powell JD & Workman ML 1990. Digital Control of Dynamic Systems. 2nd edn. Addison-Wesley Publishing Company.

8. Gilbert EG 1963. Controllability and Observability in Multivariable Control Systems. SIAM Journal of Control, 1, pp.128-151.         [ Links ]

9. Kalman RE 1963. Mathematical Description of Linear Dynamic Systems. SIAM Journal of Control, 1, pp.152-192.         [ Links ]

10. Halvorsen WG & Brown DL 1977. Impulse Technique for Structural Frequency Response Testing. Sound and Vibration, 11, pp.8-21.         [ Links ]

11. Peterson LD 1995. Efficient Computation of the Eigensystem Realization Algorithm. Journal of Guidance, Control, and Dynamics, 18, pp.395-403.         [ Links ]

 

 

Received November 1997
Final version June 1998