# Modal Parameter Identification from Output Data Only: Equivalent Approaches

Authors: Joseph Lardies
Journal: Shock and Vibration
 (2015)

---

## Abstract

The problem of modal parameter identification from output data only is presented. To identify the modal parameters different algorithms are presented: the block Hankel matrix and its shifted version and the block observability and block controllability matrices and their shifted version. These algorithms are derived from properties of the subspace approach. It is shown in the paper that these algorithms give the same results even in the noisy data case. Numerical and experimental results are presented showing the effectiveness of the procedure. In particular a microsystem constituted of a perforated microplate is analysed.


---

## Body

## 1. Introduction

The traditional modal parameter identification techniques from input and output data have been well developed and are widely used in engineering. However, in operational modal analysis (OMA) the input data, or the excitation force, is not available as a measured signal and the need to develop techniques, which are capable of accurately extracting the modal parameters from output-only measurements, has arisen. In early 1990s the NExT (natural excitation technique) algorithm [1] was proposed using correlation functions of the random response of the analyzed structure. The correlation functions are expressed as a summation of decaying sinusoids and each decaying sinusoid has a damped natural frequency, a damping ratio, and a mode shape coefficient that is identical to one of the corresponding structural modes. A general autoregressive moving average (ARMA) model [2–4] has also been employed for operational modal parameter identification. The AR parameters describe the system dynamics while the MA part is related to external disturbances. Corresponding to the multioutput case, the ARMA model has been extended to a multidimensional ARMA model or ARMAV model, in which all the modal information is contained in the AR part. The AR coefficients can be obtained by the extended instrumental variable method [5] and the modal parameters are identified by the eigenvalue decomposition of the companion matrix of the AR polynomial. The efficiency of ARMAV algorithms has been proved with various excitation models, consisting of white or coloured noise, mixed with harmonics and nonstationarity [6, 7]. In particular, Spiridonakos and Fassois [8] have used a functional series vector time-dependent autoregressive moving average (FS-VTARMA) method and an experimental laboratory test of a bridge with a heavy passing vehicle has shown that this method is robust to nonstationarity of the excitation.The dynamic behaviour of a discrete mechanical system can be described by a matrix differential equation which can be converted in a discrete time state space model and numerous papers have been presented on system identification, acquiring the estimation of parameters from measured data. For a linear time-invariant system Ho and Kalman [9] introduced the minimal state space realization problem in which a Hankel matrix is constructed by a sequence of impulse response functions called the Markov parameters. Kung [10] proposed a concept combining singular value decomposition and minimal realization algorithm for the problem of retrieving sinusoidal processes from noisy measurements. Juang and Pappa [11] introduced an eigensystem realization algorithm (ERA) for modal parameter identification and model reduction for dynamical systems from test data. This algorithm is an extension of the Ho-Kalman procedure where two indicators have been developed to quantitatively identify the system and noise modes. A stochastic subspace identification (SSI) method has been developed by van Overschee and de Moor [12] and by Peeters and de Roeck [13]. The subspace method presented by these authors identifies the state space matrices by using robust numerical techniques such as QR-factorization, singular value decomposition (SVD), and least squares. The state space matrices are related to the modal parameters and the key concept of SSI is the projection of the row space of the future outputs into the row space of the past outputs.The fundamental problem in modal parameter identification by subspace methods is the determination of the state space matrix (or transition matrix) which characterizes the dynamics of the system. The purpose of this paper is to perform a comparison of four algorithms, derived from four properties of the subspace approach, to estimate the transition matrix. The aim is to prove that these algorithms are equivalent even in the noisy data case. The first algorithm uses properties of the shifted block controllability matrix; the second algorithm uses properties of shifted columns of the block Hankel matrix, and it is proved that these two decomposition techniques give the same modal parameters. The third algorithm uses properties of the shifted block observability matrix; the fourth algorithm uses properties of shifted rows of the block Hankel matrix, and it is also proved that these two decomposition techniques give the same modal parameters. The procedures presented in the paper can identify closely eigenfrequencies that cannot be identified by the traditional Fourier transform. In addition, we show that our procedure can be used to the modal parameter identification of a microsystem constituted of a perforated microplate. To our knowledge, it is the first time that subspace algorithms are used in the microsystem field.The paper is organized as follows. The second section is devoted to the presentation of four identification algorithms based on shifting properties of block Hankel, block controllability, and block observability matrices. Validity of the modal parameter identification procedures is presented in the third section with simulated and experimental tests in laboratory. The paper is briefly concluded in Section4.


## 2. Subspace Identification Methods

### 2.1. The Discrete State Space Representation

The subspace identification method assumes that the dynamic behaviour of a vibrating system can be described by a discrete time state space model [14, 15]: (1)zk+1=Azk+wkstate equation,yk=Czk+vkobservation equation,where zk is the unobserved state vector of dimension n, yk is the (m×1) vector of observations or measured output vector at discrete time instant k, wk contains the external nonmeasured force or the excitation which can be a random force, an impulse force, a step force and so on, and vk is a measurement noise term. A is the (n×n) transition matrix describing the dynamics of the system and C is the (m×n) output or observation matrix, translating the internal state of the system into observations. The subspace identification problem deals with the determination of the two state space matrices A and C using output-only data yk.The modal parameters of a vibrating system are obtained by applying the eigenvalue decomposition of the transition matrixA [14, 15]: (2)A=ΨΛΨ-1,where Λ=diag(λi), i=1,2,…,n, is the diagonal matrix containing the complex eigenvalues and Ψ contains the eigenvectors of A as columns. The eigenfrequencies fi and damping factors ξi are obtained from the eigenvalues which are complex conjugate pair:(3)fi=12πΔtln⁡λiλi∗24+cos-1⁡λi+λi∗2λiλi∗2,ξi=ln⁡λiλi∗2ln⁡λiλi∗2+4cos-1⁡λi+λi∗/2λiλi∗2with Δt being the sampling period of analyzed signals.The mode shapes evaluated at the sensor locations are the columns of the matrixC~ obtained by multiplying the output matrix C with the matrix of eigenvectors Ψ: (4)C~=CΨ.We propose four algorithms for the estimation of the transition matrix A, in order to identify the eigenfrequencies and damping factors of a vibrating system. However, we prove that the four algorithms are equivalent even in the noisy data case.


### 2.2. Determination of the Transition Matrix by Shifting Properties of the Block Controllability Matrix and by Shifting Properties of the Block Hankel Matrix

Define the(mp×1) future data vector as yk+=ykT,yk+1T,…,yk+p-1TT and the (mp×1) past data vector as yk-1-=yk-1T,…,yk-pTT, where the superscript T denotes the transpose operation. The (mp×mp) covariance matrix between the future and the past is given by(5)H=Eyk+yk-T=R1R2·RpR2R3·Rp+1····RpRp+1·R2p-1,where E denotes the expectation operator. H is the block Hankel matrix formed with the (m×m) individual theoretical autocovariance matrices Ri=Eyk+iykT=CAi-1G, with G=Exk+1ykT. The index p corresponds to the number of block lines or block columns needed to form the block Hankel matrix H.In practice, the autocovariance matrices are estimated fromN data points and are computed by Ri=N-1∑k=1Nyk+iykT; i=1,…,2p-1 and with these estimated autocovariance matrices we form the estimated block Hankel matrix.The block Hankel matrixH can be written as(6)H=CGCAG·CAp-1GCAGCA2G·CApG····CAp-1GCApG·CA2p-2G=CCA·CAp-1GAG⋯Ap-1G=OK.By identification we obtain the (mp×n) block observability matrix O and the (n×mp) block controllability matrix K:(7)O=CCA·CAp-1,K=GAG⋯Ap-1G.This last matrix can be written as(8)K=K1Ap-1G=GK2,where the block matrices K1 and K2 are (9)K1=GAG⋯Ap-2G,K2=AG⋯Ap-1G,where K1 and K2 are n×m(p-1) matrices obtained by deleting, respectively, the last and the first block column of the block controllability matrix K. It is easy to show that K2=AK1 and the transition matrix obtained by the deleted block column of the controllability matrix method can then be calculated via pseudoinverse:(10)AK=K2K1+.The eigenvalues of the transition matrix AK can be used to identify the modal parameters and we have (11)λAK=λK2K1+.A second method to determine the transition matrix is obtained by deleting block columns of the block Hankel matrix H. Let H1 and H2 be the matrices mp×m(p-1) obtained by deleting, respectively, the first and the last block column of the block Hankel matrix:(12)H1=R2R3·RpR3R4·Rp+1····Rp+1Rp+2·R2p-1=CAGCA2G·CAp-1GCA2GCA3G·CApG····CApGCAp+1G·CA2p-2G=CCA·CAp-1AGAG⋯Ap-2GH2=R1R2·Rp-1R2R3·Rp····RpRp+1·R2p-2=CGCAG·CAp-2GCAGCA2G·CAp-1G····CAp-1GCApG·CA2p-3G=CCA·CAp-1GAG⋯Ap-2G.From (12) we have (13)H1=OAK1=OK2,H2=OK1,and the transition matrix obtained by the deleted block column of the block Hankel matrix method can then be calculated via pseudoinverses of the observability and controllability matrices: (14)AH1=O+H1K1+.The eigenvalues of the transition matrix AH1 are(15)λAH1=λO+H1K1+=λH2+H1=λK1+O+OK2=λK1+K2=λK2K1+=λAK.Firstly the eigenvalues of the transition matrix obtained by shifting properties of the block controllability matrix are the same as those obtained by shifting properties of block columns of the block Hankel matrix and secondly the eigenvalues of the transition matrix are also the eigenvalues of the matrix (H2+H1). So, in practical conditions it is sufficient to form the block Hankel matrix H, to extract matrices H1 and H2 by deleting a block column from H and compute the eigenvalues (and eventually the eigenvectors) of the matrix (H2+H1).In the presence of noise, both algorithms use the singular value decomposition (SVD) of the block Hankel matrix to get the same performances. By retaining then dominant singular values and the corresponding singular vectors, we get(16)H=USVT=OKwith U, V(mp×n) being matrices of singular vectors and S(n×n) being diagonal matrix of singular values. By identification we choose K=VT and we perform the following matrix decomposition:(17)K=VT=V11TV2T=V1TV12T,where V11T and V12T are matrices formed, respectively, with the (p-1) first and (p-1) last block columns of the matrix VT. The eigenvalues of the transition matrix are (18)λAK=λK2K1+=λV12TV11T+and this relation can also be obtained by shifting properties of the block Hankel matrix. Indeed, we have (19)H1=USV12T,H2=USV11T,λAH1=λH2+H1=λUSV11T+USV12T=λV12TV11T+=λAK.We have proved that, in the state space approach, the modal parameters obtained by shifting properties of the block Hankel matrix are the same as those obtained by shifting properties of the block controllability matrix. In the next section we analyze the eigenvalues of the transition matrix estimated by shifting properties of the block observability matrix.


### 2.3. Determination of the Transition Matrix by Shifting Properties of the Block Observability Matrix and by Shifting Properties of the Block Hankel Matrix

The(mp×n) observability matrix O has the form(20)O=O1CAp-1=CO2,where O1 and O2 are m(p-1)×n matrices obtained by deleting, respectively, the last and the first block row of the block observability matrix:(21)O1=CCA·CAp-2,O2=CACA2·CAp-1.It is easy to show that O2=O1A and the transition matrix obtained by the deleted block row of the observability matrix method is as follows:(22)AO=O1+O2.The eigenvalues of the transition matrix AO can be used to identify the modal parameters and we have (23)λAO=λO1+O2.Another method to determine the transition matrix is obtained by deleting block rows of the block Hankel matrix. Let H3 and H4 be the matrices m(p-1)×mp obtained by deleting, respectively, the first and the last block row of the block Hankel matrix:(24)H3=R2R3·Rp+1R3R4·Rp+2····RpRp+1·R2p-1=CAGCA2G·CApGCA2GCA3G·CAp+1G····CAp-1GCApG·CA2p-2G=CCA·CAp-2AGAG⋯Ap-1G,H4=R1R2·RpR2R3·Rp+1····Rp-1Rp·R2p-2=CGCAG·CAp-1GCAGCA2G·CApG····CAp-2GCAp-1G·CA2p-3G=CCA·CAp-2GAG⋯Ap-1G.From (24) we have (25)H3=O1AK=O2K,H4=O1Kand the transition matrix obtained by the deleted block row of the block Hankel matrix method is (26)AH3=O1+H3K+.The eigenvalues of the transition matrix AH3 are (27)λAH3=λO1+H3K+=λH4+H3=λK+O1+O2K=λO1+O2=λAO.The eigenvalues of the transition matrix obtained by shifting properties of the block observability matrix are the same as those obtained by shifting properties of block rows of the block Hankel matrix.In the presence of noise, we use the singular value decomposition of the block Hankel matrix to show that these two methods give the same performances. From (16) we choose O=U and we consider the following matrix decomposition:(28)O=U=U11U2=U1U12,where U11 and U12 are matrices formed, respectively, with the (p-1) first and (p-1) last block rows of the matrix of singular vectors U. The eigenvalues of the transition matrix are(29)λAO=λO1+O2=λU11+U12and this relation can also be obtained by shifting properties of the block Hankel matrix. Indeed, we have (30)H3=O2K=U12K,H4=O1K=U11K,λAH3=λH4+H3=λK+U11+U12K=λU11+U12=λAO.We have shown that, in the state space approach, the modal parameters obtained by shifting properties of the block Hankel matrix are the same as those obtained by shifting properties of the block observability matrix. Finally we note that the output or observation matrix C can be obtained by the first block row of the matrix O1 or U11.With estimates of the transition matrixA and the observation matrix C in hand, we compute the eigenvalues and eigenvectors of the transition matrix and we can identify eigenfrequencies, damping factors, and mode shapes of the vibrating system following (3) and (4). However, all the subspace modal identification algorithms have a serious problem of model order determination. When extracting physical or structural modes, subspace algorithms always generate spurious or computational modes to account for unwanted effects such as noise, leakage, residuals, nonlinearity. For these reasons, the assumed number of modes, or model order, is incremented over a wide range of values and we plot the stability diagram. The stability diagram involves tracking the estimates of eigenfrequencies and damping factors as a function of model order. As the model order is increased, more and more modal frequencies and damping ratios are estimated; hopefully, the estimates of the physical modal parameters are stabilized using a criterion based on the modal coherence of measured modes and identified modes [9]. Using this criterion we detect and remove the spurious modes. A numerical example and two experimental tests in laboratory are now presented to identify eigenfrequencies and damping factors of vibrating systems.


## 2.1. The Discrete State Space Representation

The subspace identification method assumes that the dynamic behaviour of a vibrating system can be described by a discrete time state space model [14, 15]: (1)zk+1=Azk+wkstate equation,yk=Czk+vkobservation equation,where zk is the unobserved state vector of dimension n, yk is the (m×1) vector of observations or measured output vector at discrete time instant k, wk contains the external nonmeasured force or the excitation which can be a random force, an impulse force, a step force and so on, and vk is a measurement noise term. A is the (n×n) transition matrix describing the dynamics of the system and C is the (m×n) output or observation matrix, translating the internal state of the system into observations. The subspace identification problem deals with the determination of the two state space matrices A and C using output-only data yk.The modal parameters of a vibrating system are obtained by applying the eigenvalue decomposition of the transition matrixA [14, 15]: (2)A=ΨΛΨ-1,where Λ=diag(λi), i=1,2,…,n, is the diagonal matrix containing the complex eigenvalues and Ψ contains the eigenvectors of A as columns. The eigenfrequencies fi and damping factors ξi are obtained from the eigenvalues which are complex conjugate pair:(3)fi=12πΔtln⁡λiλi∗24+cos-1⁡λi+λi∗2λiλi∗2,ξi=ln⁡λiλi∗2ln⁡λiλi∗2+4cos-1⁡λi+λi∗/2λiλi∗2with Δt being the sampling period of analyzed signals.The mode shapes evaluated at the sensor locations are the columns of the matrixC~ obtained by multiplying the output matrix C with the matrix of eigenvectors Ψ: (4)C~=CΨ.We propose four algorithms for the estimation of the transition matrix A, in order to identify the eigenfrequencies and damping factors of a vibrating system. However, we prove that the four algorithms are equivalent even in the noisy data case.


## 2.2. Determination of the Transition Matrix by Shifting Properties of the Block Controllability Matrix and by Shifting Properties of the Block Hankel Matrix

Define the(mp×1) future data vector as yk+=ykT,yk+1T,…,yk+p-1TT and the (mp×1) past data vector as yk-1-=yk-1T,…,yk-pTT, where the superscript T denotes the transpose operation. The (mp×mp) covariance matrix between the future and the past is given by(5)H=Eyk+yk-T=R1R2·RpR2R3·Rp+1····RpRp+1·R2p-1,where E denotes the expectation operator. H is the block Hankel matrix formed with the (m×m) individual theoretical autocovariance matrices Ri=Eyk+iykT=CAi-1G, with G=Exk+1ykT. The index p corresponds to the number of block lines or block columns needed to form the block Hankel matrix H.In practice, the autocovariance matrices are estimated fromN data points and are computed by Ri=N-1∑k=1Nyk+iykT; i=1,…,2p-1 and with these estimated autocovariance matrices we form the estimated block Hankel matrix.The block Hankel matrixH can be written as(6)H=CGCAG·CAp-1GCAGCA2G·CApG····CAp-1GCApG·CA2p-2G=CCA·CAp-1GAG⋯Ap-1G=OK.By identification we obtain the (mp×n) block observability matrix O and the (n×mp) block controllability matrix K:(7)O=CCA·CAp-1,K=GAG⋯Ap-1G.This last matrix can be written as(8)K=K1Ap-1G=GK2,where the block matrices K1 and K2 are (9)K1=GAG⋯Ap-2G,K2=AG⋯Ap-1G,where K1 and K2 are n×m(p-1) matrices obtained by deleting, respectively, the last and the first block column of the block controllability matrix K. It is easy to show that K2=AK1 and the transition matrix obtained by the deleted block column of the controllability matrix method can then be calculated via pseudoinverse:(10)AK=K2K1+.The eigenvalues of the transition matrix AK can be used to identify the modal parameters and we have (11)λAK=λK2K1+.A second method to determine the transition matrix is obtained by deleting block columns of the block Hankel matrix H. Let H1 and H2 be the matrices mp×m(p-1) obtained by deleting, respectively, the first and the last block column of the block Hankel matrix:(12)H1=R2R3·RpR3R4·Rp+1····Rp+1Rp+2·R2p-1=CAGCA2G·CAp-1GCA2GCA3G·CApG····CApGCAp+1G·CA2p-2G=CCA·CAp-1AGAG⋯Ap-2GH2=R1R2·Rp-1R2R3·Rp····RpRp+1·R2p-2=CGCAG·CAp-2GCAGCA2G·CAp-1G····CAp-1GCApG·CA2p-3G=CCA·CAp-1GAG⋯Ap-2G.From (12) we have (13)H1=OAK1=OK2,H2=OK1,and the transition matrix obtained by the deleted block column of the block Hankel matrix method can then be calculated via pseudoinverses of the observability and controllability matrices: (14)AH1=O+H1K1+.The eigenvalues of the transition matrix AH1 are(15)λAH1=λO+H1K1+=λH2+H1=λK1+O+OK2=λK1+K2=λK2K1+=λAK.Firstly the eigenvalues of the transition matrix obtained by shifting properties of the block controllability matrix are the same as those obtained by shifting properties of block columns of the block Hankel matrix and secondly the eigenvalues of the transition matrix are also the eigenvalues of the matrix (H2+H1). So, in practical conditions it is sufficient to form the block Hankel matrix H, to extract matrices H1 and H2 by deleting a block column from H and compute the eigenvalues (and eventually the eigenvectors) of the matrix (H2+H1).In the presence of noise, both algorithms use the singular value decomposition (SVD) of the block Hankel matrix to get the same performances. By retaining then dominant singular values and the corresponding singular vectors, we get(16)H=USVT=OKwith U, V(mp×n) being matrices of singular vectors and S(n×n) being diagonal matrix of singular values. By identification we choose K=VT and we perform the following matrix decomposition:(17)K=VT=V11TV2T=V1TV12T,where V11T and V12T are matrices formed, respectively, with the (p-1) first and (p-1) last block columns of the matrix VT. The eigenvalues of the transition matrix are (18)λAK=λK2K1+=λV12TV11T+and this relation can also be obtained by shifting properties of the block Hankel matrix. Indeed, we have (19)H1=USV12T,H2=USV11T,λAH1=λH2+H1=λUSV11T+USV12T=λV12TV11T+=λAK.We have proved that, in the state space approach, the modal parameters obtained by shifting properties of the block Hankel matrix are the same as those obtained by shifting properties of the block controllability matrix. In the next section we analyze the eigenvalues of the transition matrix estimated by shifting properties of the block observability matrix.


## 2.3. Determination of the Transition Matrix by Shifting Properties of the Block Observability Matrix and by Shifting Properties of the Block Hankel Matrix

The(mp×n) observability matrix O has the form(20)O=O1CAp-1=CO2,where O1 and O2 are m(p-1)×n matrices obtained by deleting, respectively, the last and the first block row of the block observability matrix:(21)O1=CCA·CAp-2,O2=CACA2·CAp-1.It is easy to show that O2=O1A and the transition matrix obtained by the deleted block row of the observability matrix method is as follows:(22)AO=O1+O2.The eigenvalues of the transition matrix AO can be used to identify the modal parameters and we have (23)λAO=λO1+O2.Another method to determine the transition matrix is obtained by deleting block rows of the block Hankel matrix. Let H3 and H4 be the matrices m(p-1)×mp obtained by deleting, respectively, the first and the last block row of the block Hankel matrix:(24)H3=R2R3·Rp+1R3R4·Rp+2····RpRp+1·R2p-1=CAGCA2G·CApGCA2GCA3G·CAp+1G····CAp-1GCApG·CA2p-2G=CCA·CAp-2AGAG⋯Ap-1G,H4=R1R2·RpR2R3·Rp+1····Rp-1Rp·R2p-2=CGCAG·CAp-1GCAGCA2G·CApG····CAp-2GCAp-1G·CA2p-3G=CCA·CAp-2GAG⋯Ap-1G.From (24) we have (25)H3=O1AK=O2K,H4=O1Kand the transition matrix obtained by the deleted block row of the block Hankel matrix method is (26)AH3=O1+H3K+.The eigenvalues of the transition matrix AH3 are (27)λAH3=λO1+H3K+=λH4+H3=λK+O1+O2K=λO1+O2=λAO.The eigenvalues of the transition matrix obtained by shifting properties of the block observability matrix are the same as those obtained by shifting properties of block rows of the block Hankel matrix.In the presence of noise, we use the singular value decomposition of the block Hankel matrix to show that these two methods give the same performances. From (16) we choose O=U and we consider the following matrix decomposition:(28)O=U=U11U2=U1U12,where U11 and U12 are matrices formed, respectively, with the (p-1) first and (p-1) last block rows of the matrix of singular vectors U. The eigenvalues of the transition matrix are(29)λAO=λO1+O2=λU11+U12and this relation can also be obtained by shifting properties of the block Hankel matrix. Indeed, we have (30)H3=O2K=U12K,H4=O1K=U11K,λAH3=λH4+H3=λK+U11+U12K=λU11+U12=λAO.We have shown that, in the state space approach, the modal parameters obtained by shifting properties of the block Hankel matrix are the same as those obtained by shifting properties of the block observability matrix. Finally we note that the output or observation matrix C can be obtained by the first block row of the matrix O1 or U11.With estimates of the transition matrixA and the observation matrix C in hand, we compute the eigenvalues and eigenvectors of the transition matrix and we can identify eigenfrequencies, damping factors, and mode shapes of the vibrating system following (3) and (4). However, all the subspace modal identification algorithms have a serious problem of model order determination. When extracting physical or structural modes, subspace algorithms always generate spurious or computational modes to account for unwanted effects such as noise, leakage, residuals, nonlinearity. For these reasons, the assumed number of modes, or model order, is incremented over a wide range of values and we plot the stability diagram. The stability diagram involves tracking the estimates of eigenfrequencies and damping factors as a function of model order. As the model order is increased, more and more modal frequencies and damping ratios are estimated; hopefully, the estimates of the physical modal parameters are stabilized using a criterion based on the modal coherence of measured modes and identified modes [9]. Using this criterion we detect and remove the spurious modes. A numerical example and two experimental tests in laboratory are now presented to identify eigenfrequencies and damping factors of vibrating systems.


## 3. Applications

### 3.1. Simulated Data

To prove the effectiveness of the identification procedure based on the subspace analysis we consider a two-DOF system with very closely spaced modes. The parameters of the system aref1=31 Hz, f2=31.5 Hz, ξ1=0.01, and ξ2=0.02. Figure 1 shows the free response of the system where a Gaussian white noise has been added: the generated data were corrupted by a random noise. The sampling frequency is 100 Hz and 300 time samples are used in the simulation. We convert to the frequency domain this time response by taking the discrete Fourier transform of the noisy signal. It is impossible to identify the two frequencies components by using the FFT as shown in Figure 2, where the power spectral density has been plotted. In our identification procedure, we plot the stabilization diagram on eigenfrequencies and damping factors. Figure 3 shows stabilization diagrams using shifting properties of the block controllability matrix with the modal coherence indicator: spurious modes have been eliminated and only physical modes are present. Figure 4 shows stabilization diagrams using shifting properties of the block observability matrix with the modal coherence indicator. Our procedure can separate closely spaced modes and the mean values on identified eigenfrequencies and damping factors obtained by an average over the orders of the stabilization diagram are f^1=31 Hz, f^2=31.5, ξ^1=0.01, and ξ^2=0.0199. A very satisfactory estimation on eigenfrequencies and damping factors has been obtained using simulated data. These values have been estimated using the observability identification procedure; however, very similar results are obtained if we use the controllability identification procedure. Two experimental tests in laboratory are presented in the following sections.Figure 1
Time response of the two-DOF system.Figure 2
Spectrum for the two-DOF system.Figure 3
Stabilization diagram on eigenfrequencies and damping factors for the two-DOF system using shifting properties of the block controllability matrix.Figure 4
Stabilization diagram on eigenfrequencies and damping factors for the two-DOF system using shifting properties of the block observability matrix.


### 3.2. Modal Parameter Identification of a Clamped Beam

Figure5 shows the first experimental system tested in laboratory. It is a simple clamped horizontal rectangular cantilever beam with five measurement locations equally spaced. A Gaussian excitation is applied transversally at the free end of the beam and Figure 6 shows a typical response of an accelerometer. Only the time responses of all accelerometers are used for modal parameter identification of the clamped beam which is excited with an unmeasured random force. The sampling frequency of signals is 1280 Hz and 8192 data points are collected for each channel. Figure 7 shows the stabilization diagram on eigenfrequencies and damping factors and Table 1 gives the mean values of these identified modal parameters obtained by an average over the orders of the stabilization diagram. For comparison purposes the theoretical values on eigenfrequencies are also computed. These values are obtained from the mechanical characteristics of the beam: Young’s modulus E=2.1×105 MPa, mass density ρ=9000 kg·m−3, length of the beam L=0,56 m, and thickness of the beam e=97×10-4 m. According to the hypothesis of clamped-free beam, the flexural eigenfrequencies result from the expression derived from the Euler-Bernoulli theory:(31)fi=αi2e43πL2Eρ.Using boundary conditions of the beam we obtain the equation cos(αi)cosh(αi)=-1 and by resolution of this equation we determine the values of coefficients αi. We note that the eigenfrequencies are very well identified; the maximum value of the relative error is only −0.61%. Another experimentation test in laboratory is presented in the next section.Table 1
Natural eigenfrequencies and damping factors of the experimental beam.
Modes
Theoretical eigenfrequenciesfi  (Hz)
Identified eigenfrequenciesf^i  (Hz)
Relative error on eigenfrequencies (%)
Identified damping factor (%)
1
24.13
24.17
0.16
1.06
2
150.81
150.58
−0.15
0.13
3
422.29
419.72
−0.61
0.04Figure 5
Schematic of the clamped beam in laboratory.Figure 6
Time response of an accelerometer.Figure 7
Stabilization diagram on eigenfrequencies and damping factors for the clamped beam using shifting properties of the block observability matrix.


### 3.3. Modal Parameter Identification of a Perforated Microplate

A microelectromechanical system (MEMS) can be constituted by a perforated rectangular plate [16] as shown in Figure 8. The main purpose of perforations is to reduce the damping and spring forces acting in the MEMS due to the fluid flow inside and around the microstructure. Generally, the modeling problem is quite complicated since the damping force acting in the moving MEMS depends on the 3D fluid flow in the perforations and also around the structure. In this paper a single degree of freedom of a spring-mass-damper model is used to study the microplate behavior. The model is constituted by the following parameters: the plate mass m concentrated in the central plate, the damping coefficient c, and the stiffness coefficient k which are constituted of fluidic and nonfluidic components. The second order differential equation describing the dynamic behavior of the microplate is(32)mz¨t+cz˙t+kzt=gt,where g(t) is the external force acting in the microplate. In our experimental test this excitation force is step function and Figure 9 shows the time response of the microplate center. The dynamic measurements are conducted in the time domain by means of a laser vibrometer. Only this time response of the structure is used in the identification procedure where the sampling frequency of signals is 2 MHz and 3000 time samples are considered. Table 2 shows the dimensions and material properties of the microplate. The perforated microplate area is given by A=a2-s0Nh=3.127×10-8 m2 and its mass is m=ρAhc=3.814×10-9 kg. The microplate stiffness is given by k=m2πf2 and the damping coefficient is c=4πmFξ, where F is the resonance frequency of the perforated microplate and ξ is the damping factor. These two modal parameters are obtained by an average over the orders of the stabilization diagram presented in Figure 10. To our knowledge, it is the first time that the subspace algorithm is applied to a microsystem and in particular to a perforated microplate.Table 2
Dimensions and properties of the perforated microplate.
Plate sidea
Gap thicknessh
Plate thicknesshc
Hole sizes0
Plate densityρ
Number of holesNh
185.96μm
3μm
6.312μm
7.19μm
19320 kg·m−3
64Figure 8
Schematic of the perforated microplate.Figure 9
Time response of the microplate to a step force actuation.Figure 10
Stabilization diagram on resonance frequency and damping factor for the experimental perforated microplate using shifting properties of the block observability matrix.Table3 shows the identified microplate parameters and Figure 11 shows the comparison between the measured time response of the structure and the reconstructed response obtained from the identified modal parameters. An important part of the identification process is to assess the quality of the estimated model obtained from the identified parameters and as a measure of model quality we use(33)Q=1-∑k=1Nyk-y^k2∑k=1Nyk2which is a normalized estimation error where increasing values of Q indicate better fit and 1 indicates a perfect fit of the model to the validation data. In this expression yk are the measured data and y^k are the estimated data. For the perforated microplate we obtain Q=0.817 which is a satisfactory value. The decay function of the response is expressed by the relation x(t)=Aexp(-Dt), where A is obtained by initial conditions and D=2πfξ is computed from the identified modal parameters. This decay function is plotted in Figure 11. Note that the coefficients A and D=c/(2m) can also be evaluated directly by interpolating the measured response of the microplate and in this case we can obtain the value of the damping coefficient c from D and m.Table 3
Parameters identification of the perforated microplate.
Resonance frequencyf
Damping factorξ
Stiffness coefficientk
Damping coefficientc
22810 Hz
1.62%
78.33 N·m−1
17.71× 10−6 N·s·m−1Figure 11
Comparison between the measured time response of the microplate (in blue) and the reconstituted time response of the microplate (in red).


## 3.1. Simulated Data

To prove the effectiveness of the identification procedure based on the subspace analysis we consider a two-DOF system with very closely spaced modes. The parameters of the system aref1=31 Hz, f2=31.5 Hz, ξ1=0.01, and ξ2=0.02. Figure 1 shows the free response of the system where a Gaussian white noise has been added: the generated data were corrupted by a random noise. The sampling frequency is 100 Hz and 300 time samples are used in the simulation. We convert to the frequency domain this time response by taking the discrete Fourier transform of the noisy signal. It is impossible to identify the two frequencies components by using the FFT as shown in Figure 2, where the power spectral density has been plotted. In our identification procedure, we plot the stabilization diagram on eigenfrequencies and damping factors. Figure 3 shows stabilization diagrams using shifting properties of the block controllability matrix with the modal coherence indicator: spurious modes have been eliminated and only physical modes are present. Figure 4 shows stabilization diagrams using shifting properties of the block observability matrix with the modal coherence indicator. Our procedure can separate closely spaced modes and the mean values on identified eigenfrequencies and damping factors obtained by an average over the orders of the stabilization diagram are f^1=31 Hz, f^2=31.5, ξ^1=0.01, and ξ^2=0.0199. A very satisfactory estimation on eigenfrequencies and damping factors has been obtained using simulated data. These values have been estimated using the observability identification procedure; however, very similar results are obtained if we use the controllability identification procedure. Two experimental tests in laboratory are presented in the following sections.Figure 1
Time response of the two-DOF system.Figure 2
Spectrum for the two-DOF system.Figure 3
Stabilization diagram on eigenfrequencies and damping factors for the two-DOF system using shifting properties of the block controllability matrix.Figure 4
Stabilization diagram on eigenfrequencies and damping factors for the two-DOF system using shifting properties of the block observability matrix.


## 3.2. Modal Parameter Identification of a Clamped Beam

Figure5 shows the first experimental system tested in laboratory. It is a simple clamped horizontal rectangular cantilever beam with five measurement locations equally spaced. A Gaussian excitation is applied transversally at the free end of the beam and Figure 6 shows a typical response of an accelerometer. Only the time responses of all accelerometers are used for modal parameter identification of the clamped beam which is excited with an unmeasured random force. The sampling frequency of signals is 1280 Hz and 8192 data points are collected for each channel. Figure 7 shows the stabilization diagram on eigenfrequencies and damping factors and Table 1 gives the mean values of these identified modal parameters obtained by an average over the orders of the stabilization diagram. For comparison purposes the theoretical values on eigenfrequencies are also computed. These values are obtained from the mechanical characteristics of the beam: Young’s modulus E=2.1×105 MPa, mass density ρ=9000 kg·m−3, length of the beam L=0,56 m, and thickness of the beam e=97×10-4 m. According to the hypothesis of clamped-free beam, the flexural eigenfrequencies result from the expression derived from the Euler-Bernoulli theory:(31)fi=αi2e43πL2Eρ.Using boundary conditions of the beam we obtain the equation cos(αi)cosh(αi)=-1 and by resolution of this equation we determine the values of coefficients αi. We note that the eigenfrequencies are very well identified; the maximum value of the relative error is only −0.61%. Another experimentation test in laboratory is presented in the next section.Table 1
Natural eigenfrequencies and damping factors of the experimental beam.
Modes
Theoretical eigenfrequenciesfi  (Hz)
Identified eigenfrequenciesf^i  (Hz)
Relative error on eigenfrequencies (%)
Identified damping factor (%)
1
24.13
24.17
0.16
1.06
2
150.81
150.58
−0.15
0.13
3
422.29
419.72
−0.61
0.04Figure 5
Schematic of the clamped beam in laboratory.Figure 6
Time response of an accelerometer.Figure 7
Stabilization diagram on eigenfrequencies and damping factors for the clamped beam using shifting properties of the block observability matrix.


## 3.3. Modal Parameter Identification of a Perforated Microplate

A microelectromechanical system (MEMS) can be constituted by a perforated rectangular plate [16] as shown in Figure 8. The main purpose of perforations is to reduce the damping and spring forces acting in the MEMS due to the fluid flow inside and around the microstructure. Generally, the modeling problem is quite complicated since the damping force acting in the moving MEMS depends on the 3D fluid flow in the perforations and also around the structure. In this paper a single degree of freedom of a spring-mass-damper model is used to study the microplate behavior. The model is constituted by the following parameters: the plate mass m concentrated in the central plate, the damping coefficient c, and the stiffness coefficient k which are constituted of fluidic and nonfluidic components. The second order differential equation describing the dynamic behavior of the microplate is(32)mz¨t+cz˙t+kzt=gt,where g(t) is the external force acting in the microplate. In our experimental test this excitation force is step function and Figure 9 shows the time response of the microplate center. The dynamic measurements are conducted in the time domain by means of a laser vibrometer. Only this time response of the structure is used in the identification procedure where the sampling frequency of signals is 2 MHz and 3000 time samples are considered. Table 2 shows the dimensions and material properties of the microplate. The perforated microplate area is given by A=a2-s0Nh=3.127×10-8 m2 and its mass is m=ρAhc=3.814×10-9 kg. The microplate stiffness is given by k=m2πf2 and the damping coefficient is c=4πmFξ, where F is the resonance frequency of the perforated microplate and ξ is the damping factor. These two modal parameters are obtained by an average over the orders of the stabilization diagram presented in Figure 10. To our knowledge, it is the first time that the subspace algorithm is applied to a microsystem and in particular to a perforated microplate.Table 2
Dimensions and properties of the perforated microplate.
Plate sidea
Gap thicknessh
Plate thicknesshc
Hole sizes0
Plate densityρ
Number of holesNh
185.96μm
3μm
6.312μm
7.19μm
19320 kg·m−3
64Figure 8
Schematic of the perforated microplate.Figure 9
Time response of the microplate to a step force actuation.Figure 10
Stabilization diagram on resonance frequency and damping factor for the experimental perforated microplate using shifting properties of the block observability matrix.Table3 shows the identified microplate parameters and Figure 11 shows the comparison between the measured time response of the structure and the reconstructed response obtained from the identified modal parameters. An important part of the identification process is to assess the quality of the estimated model obtained from the identified parameters and as a measure of model quality we use(33)Q=1-∑k=1Nyk-y^k2∑k=1Nyk2which is a normalized estimation error where increasing values of Q indicate better fit and 1 indicates a perfect fit of the model to the validation data. In this expression yk are the measured data and y^k are the estimated data. For the perforated microplate we obtain Q=0.817 which is a satisfactory value. The decay function of the response is expressed by the relation x(t)=Aexp(-Dt), where A is obtained by initial conditions and D=2πfξ is computed from the identified modal parameters. This decay function is plotted in Figure 11. Note that the coefficients A and D=c/(2m) can also be evaluated directly by interpolating the measured response of the microplate and in this case we can obtain the value of the damping coefficient c from D and m.Table 3
Parameters identification of the perforated microplate.
Resonance frequencyf
Damping factorξ
Stiffness coefficientk
Damping coefficientc
22810 Hz
1.62%
78.33 N·m−1
17.71× 10−6 N·s·m−1Figure 11
Comparison between the measured time response of the microplate (in blue) and the reconstituted time response of the microplate (in red).


## 4. Conclusion

In modal parameter identification the key point is to determine the relationship between the system parameters and the measured dynamical data. Because in practice the system input data are often unavailable, attention has been paid to system parameter identification when only output measurements are available. Such parameters can then be used for fault detection, structural health monitoring, and model validation. In this paper, we have proposed four algorithms for modal parameter identification, and we have proved that the shifted controllability matrix algorithm gives the same modal parameters as the shifted block column of the block Hankel matrix algorithm and that the shifted observability matrix algorithm gives the same modal parameters as the shifted block row of the block Hankel matrix algorithm. These relationships have been established in the exact data case and in the noisy data case and are derived from properties of the subspace identification approach. Numerical and experimental results have shown the effectiveness of the procedure presented in the paper in modal parameter identification.



---

*Source: 124932-2015-11-02.xml*