1081286517733034

Original Article
Modal participation in multiple input
Ibrahim time domain identification
Mathematics and Mechanics of Solids
1–13
Ó The Author(s) 2017
Reprints and permissions:
sagepub.co.uk/journalsPermissions.nav
DOI: 10.1177/1081286517733034
journals.sagepub.com/home/mms
Rune Brincker
Civil Engineering Department, Technical University of Denmark (DTU), Lyngby, Denmark
Peter Olsen
Department of Engineering, Aarhus Universitet, Aarhus, Denmark
Sandro Amador
Centre for Oil and Gas - DTU, Danmarks Tekniske Universitet, Lyngby, Denmark
Martin Juul
Department of Engineering, Aarhus Universitet, Aarhus, Denmark
Abdollah Malekjafarian
School of Civil Engineering, University College Dublin, Dublin, Ireland
Mohammad Ashory
Mechanical Engineering Department, Semnan University, Semnan, Iran
Received 9 March 2017; accepted 31 August 2017
Abstract
The Ibrahim time domain (ITD) identification technique was one of the first techniques formulated for multiple output
modal analysis based on impulse response functions or general free decays. However, the technique has not been used
much in recent decades due to the fact that the technique was originally formulated for single input systems that suffer
from well-known problems in case of closely spaced modes. In this paper, a known, but more modern formulation of the
ITD technique is discussed. In this formulation the technique becomes multiple input by adding some Toeplitz matrices
over a set of free decays. It is shown that a special participation matrix can be defined that cancels out whenever the system matrix is estimated. The participation matrix becomes rank deficient if a mode is missing in the responses, but if any
mode is present in one of the considered free decays, the participation matrix has full rank. This secures that all modes
will be contained in the estimated system matrix. Finally, it is discussed how correlation functions estimated from the
operational responses of structures can be used as free decays for the multiple-input ITD formulation, and the estimation errors of the identification technique are investigated in a simulation study with closely spaced modes. The simulation study shows that the multiple-input formulation provides estimates with significantly smaller errors on both mode
shape and natural frequency estimates.
Corresponding author:
Rune Brincker, Civil Engineering Department, Technical University of Denmark (DTU), Building 118, 2800 Lyngby, Denmark
Email: runeb@byg.dtu.dk
2
Mathematics and Mechanics of Solids 00(0)
Keywords
Modal participation, Ibrahim time domain, multiple-input formulation, operational modal analysis, closely spaced modes
1. Introduction
In this paper, the Ibrahim time domain (ITD) is studied in relation to applications in operational modal
analysis (OMA), i.e. cases where the excitation forces are unknown and only the responses have been
measured and, thus, forms the basis for the identification. An overview of OMA techniques can be
found in papers by Zhang et al. [1], Masjedian and Keshmiri [2], and Brincker and Ventura [3].
The ITD technique was developed in the 1970s and is one of the first techniques developed in the
modal community for identification of multiple-output systems, i.e. where advantage is taken of information from several measurement channels at the same time. Introduced by Ibrahim [7–9], the technique was from the beginning aimed at applications in OMA where free decays were obtained from the
random responses by the random decrement technique [8].
Shortly after the introduction of the ITD technique, the polyreference technique was introduced by
Vold and co-workers [10,11] at the beginning of the 1980s and the ERA technique was introduced by
Juang, Pappa, and co-workers in the mid-1980s [12–14]. Because these techniques were formulated for
multiple-input applications that handle several free decays at the same time, they performed better in
cases of closely spaced modes. Owing to their better efficiency in these cases, they quickly became more
popular than the ITD technique.
In 1986 Fukuzono [15] proposed a modified version of the ITD technique that allows taking multiple
inputs into account. However, this multiple-input version of the ITD technique has never become widely
used. In the studies by Malekjafarian et al. [16,17] and Olsen and Brincker [18], the capability of the classical ITD to identify closely spaced modes has been investigated.
In these contributions, apart from studying the ITD capabilities, ways have been sought to find a
multiple-input ITD formulation that is closely related to the original one. These studies have indicated
the possibility to add some matrices that are central in the ITD technique (the Toeplitz matrices in the
ITD formulation) and, in this way, extend the capabilities of the ITD technique to simultaneously handle several free decays at the same time. Similar Toeplitz matrices are also introduced in Brincker and
Ventura [3] and Brincker [19] where a modern multiple-input formulation of the ITD technique is given
that will be further studied in this paper.
An important result of the work presented in this paper, is that one of the matrices in the matrix product defining the Toeplitz matrices – a participation matrix to the Toeplitz matrices – become rank deficient whenever the modal participation factor of a mode vanishes. It is shown that this participation
matrix can be added over a set of different inputs. This assures the full rank of the Toeplitz matrices and
that all modes are present in the eigenvalue decomposition of system matrix that is estimated with the
ITD technique.
It is then discussed how the ITD techniques can be applied to identify the modal parameters from the
correlation functions estimated from a set of random responses. The idea behind this approach is that
each column – or row – in the correlation function matrix represents a free decay, and thus constitutes a
basis for multiple-input identification in OMA. In order to distinguish physical modes from noise modes
when using the ITD technique in OMA applications, it in this study a common modal participation factor for the whole correlation function matrix is being used as defined in Brincker and Ventura [3].
Finally, the accuracy of the multiple-input formulation is investigated in a simulation study with
closely spaced modes. In this investigation the exact theoretical solutions are being used for the correlation functions so the obtained results for the accuracy is not influenced by errors from estimating the
correlation functions.
2. Background theory
The classical equation of motion for a single particle of mass m, attached to a spring of stiffness k, and
a viscous damper width damping c, is given by m€y + c_y + ky = f , where f is the external force on the
particle.
Brincker et al.
3
In Cook [20] it is shown how this equation can be extended from a particle to continuum dynamics of
a volume V and with a surface S using the variational principle of virtual work, that for any discretization into a set of finite elements turns into the following matrix equation for the discretized system
_ + Ky(t) = f(t)
M€y(t) + Cy(t)
ð1Þ
where f(t) is a vector of external forces corresponding to discretized set of FE degrees of freedom (the
DOFs) nd of the entire structure, t the continuous time variable, and y(t) 2 <nd the vector of displacements. The matrices M, C, and K contain the discretized mass, damping, and stiffness properties of the
solutions yn (t) to the homogenous verstructure. It is shown in Ewins [21], that there exist nd possible
sion of this equation of the form yn (t) = cn bn eln t + cn bn eln t , where bn 2 C nd are the mode shape vectors,
ln 2 C the poles and cn 2 C a participation coefficient corresponding to the nth vibration mode. It is
well known that mode shapes, poles and modal participation factors occurs in complex conjugate pairs
in case of general non-proportional damping. Considering the first N modes, it is also well known, that
the discrete time (sampled) version of a free decay solution can be written as
y(k) = c1 b1 q1 ðk Þ + c1 b1 q1 ðk Þ + c2 b2 q2 ðk Þ + c2 b2 q2 ðk Þ + + cN bN qN ðk Þ + cN bN qN ðk Þ
ð2Þ
where k 2 @ denotes any discreet time instant tk = kDt, y(k) 2 <nd is the displacement response vector
and qn ðk Þ = eln tk the modal coordinate of the nth vibration mode at the discreet time instant k. The operator (:) denotes the complex conjugate.
2.1. Classical ITD (SIMO formulation)
In order to define the framework for the definition of the participation matrix in the following section,
we will revisit the original ITD technique in a form very close to the original, but using a formalism
found in Brincker and Ventura [3] which makes it easier to generalize to the multiple-input case. For
more details about the original formulation of the ITD technique, refer, for instance, to Ibrahim [22].
Assuming N modes occurring in complex conjugate pairs in equation (2), we have the 2N terms consisting of the mode shape vectors bn , the modal participation factors cn , and the modal coordinates
qn ðk Þ = eln tk , and their corresponding complex conjugates bn , cn and qn ðk Þ = eln tk . These modal properties can be gathered in matrices and vectors with following structure
8
9
3
2
q1 ðk Þ >
0
0
c1 0 0
>
>
>
>
>
>
>
7
6 0 c1 0
0
ð
k
Þ
q
>
>
1
<
=
7
h
i 6
7
6
.
.
.
.
n
.
.
.
.
,
q
ð
k
Þ
=
B = ½ b1 b1 bN bN ,
cnn = 6 0 .
ð3Þ
7
. .
0 7
. >
>
6
>
>
>
>
4 0 0 cN 0 5
>
q ðk Þ >
>
>
: N
;
0 0 0
0 cN
qN ðk Þ
where B 2 C nd × 2N is the mode shape matrix containing the 2N mode shapes bn , n cnn 2 C 2N × 2N the
diagonal matrix containing the 2N modal participation factors cn and qðk Þ 2 C 2N a vector with the 2N
modal coordinates qn ðk Þ. From the definitions of equation (3), equation (2) can be rewritten in matrix
notation, as
h
i
y(k) = B n cnn q(k)
ð4Þ
Writing down equation (4) for all Ns sampled time intervals, k, and combining the resulting equations,
yields
h
i
Y1 = B n cnn L
ð5Þ
where L = ½ q(1) q(2) qðNs Þ 2 C2N × Ns and Y1 = ½ yð1Þ yð2Þ yðNs Þ 2 <nd × Ns are
matrices containing the modal coordinates and the free decay response vectors and measured at all Ns
discreet time instants k, respectively. The classical ITD formulation operates with three time constants
4
Mathematics and Mechanics of Solids 00(0)
Dt1 , Dt2 and Dt3 describing different time delays of the involved data. The last time constant Dt3 is introduced in order to obtain a delayed version of equation (5). Delaying all the measured responses in equation (5), we obtain the corresponding delayed response matrix Y2 expressed by
h
i
n cnn L
ð6Þ
Y2 = B
= el1 Dt3 b1 el1 Dt3 b elN Dt3 bN elN Dt3 b 2 Cnd × 2N is the modified mode shape matrix.
where B
1
N
Equations (5) and (6) can be combined in one single matrix equation, as
h
i
h
i
Y1
B
H1 =
= n cnn L = C n cnn L
ð7Þ
Y2
B
with H1 2 C2nd × Ns now gathering the free decays responses Y1 and their delayed counterparts Y2 . Again
delaying the responses, but now with the time delay Dt1 , we obtain the response matrix H2 that can be
expressed similarly to equation (7), as
h
ih
i
h
ih
i
B
H2 = n m n n n c n n L = C n m n n n c n n L
ð8Þ
B
where n mn n 2 C2N × 2N is diagonal matrix containing the discrete quantities mn = eln Dt1 . Assuming that
the number of modes
N is equal to the number of DOFs nd , C becomes a 2N × 2N matrix. By isolating
n
the matrix cnn L from equations (7) and (8), and combining the resulting equations, the identity
n
1
mnn C1 H2 = C1 H1 is obtained, leading to the equation that is central in the ITD technique:
H2 = SH1
where S 2 C2N × 2N is a system matrix whose eigenvalue decomposition is obtained by
h
i
S = C n mn n C1
ð9Þ
ð10Þ
Traditionally, equation (9) is solved by a pre-multiplication with HT1 . Assuming that the resulting square
matrix H1 HT1 has full rank, it can be inverted to give the least square (LS) estimate of S:
^ 1 = H2 HT (H1 HT )1
S
1
1
ð11Þ
The matrix HT1 (H1 HT1 )1 is also known as the LS version of the pseudo inverse H+
1 of H1 . Similarly, we
can also post-multiply equation (9) with HT2 to obtain the estimate
^ 2 = H2 HT (H1 HT )1
S
2
2
ð12Þ
^ 1 and S
^ 2 are found, the system matrix can be obtained as an average between the
Once the estimates S
two estimates
^ 2 )=2
^ = (S
^1 + S
ð13Þ
S
which is also known as the double LS estimate in the ITD technique [24]. In a modern formulation the
classical ITD can be derived based on a Hankel matrix. Defining the Hankel matrix with four block
rows and np columns, we have
3
2
y(1) y(2) y(np )
6 y(2) y(3) y(np + 1) 7
7
ð14Þ
H=6
4 y(3) y(4) y(np + 2) 5
y(4) y(5) y(np + 3)
Assuming equidistant sampling with the time step Dt, the time delay Dt3 = Dt and the time delay
Dt1 = 2Dt, we see that the above mentioned matrices H1 and H2 are the upper and the lower part of the
Brincker et al.
5
Hankel matrix H. This defines the quantities mn = eln 2Dt = m2n , where the discrete time poles are defined
as mn = eln Dt . Once the eigenvalues and eigenvectors are estimated, the continuous time poles are found
from the well-known relation
ln =
log (mn )
Dt
ð15Þ
where log (:) stands for the natural logarithm. The mode shapes are found as the upper or lower part of
the eigenvectors in C. It is worth noticing that a minimum size formulation can be obtained by defining
the Hankel matrix
2
3
y(np )
y(1) y(2) H = 4 y(2) y(3) y(np + 1) 5
ð16Þ
y(3) y(4) y(np + 2)
and then defining the two overlapping Hankel matrices H1 and H2 , where H1 consists of the two upper
block rows and H2 of the two lower block rows of H. In this specific case, the two Hankel matrices are
only delayed one time step, so that Dt1 = Dt. This can be generalized to a Hankel matrix H with an arbitrary number of block rows larger than 3 and to Hankel matrices H1 and H2 with arbitrary overlap.
For simplicity, in the following we restrict our analysis to Hankel matrices of the form given by equation (14) and to Hankel matrices H1 and H2 with zero overlap. This means that in the following we use
the time delays Dt1 = 2Dt, Dt2 = Dt and Dt3 = Dt, where Dt is the sampling time step.
The maximum number of modes in the model is limited to the number of degrees of freedom nd in
the experiment. The number of rows in the matrices H1 and H2 is 2nd , thus the size of the system matrix
is 2nd × 2nd . Therefore the number of eigenvalues is 2nd and the number of modes is N = nd . In order to
increase the number of modes, Ibrahim has recommended adding some of the already measured degrees
of freedom to the response vector – also called pseudo measurements. In classical ITD the pseudo measurements are delayed by the time delay Dt2 . However, in a modern formulation the pseudo measurement can just be seen as the above mentioned generalization where the matrices H1 and H2 can have an
arbitrary number of block rows.
The important characteristics of the ITD technique is the least square solutions in equations (11) and
(12) for estimation of the system matrix, and the possibility to use the double least square estimate given
by equation (13). This approach is known to give a smaller bias which in many cases might be important for a good estimation of the damping. In the following we will try to keep these characteristics of
the ITD, while generalizing the technique to multiple-input cases.
2.2. MIMO participation matrix
One of the problems with the classical formulation of the ITD, as presented above, is that it is single
input, i.e. only a single free decay is allowed at the same time. However, we will use the fact, see for
instance Brincker and Ventura [3], that this problem can be removed by adding the following matrices
that all appear in equations (11) and (12):
T11 = H1 HT1 ;
T12 = H1 HT2
T21 = H2 HT1 ;
T22 = H2 HT2
ð17Þ
over a set of free decays representing different inputs . These commonly known matrices are square consisting of four blocks defined by the two block rows in H1 and H2 . Since they are constant along the
block diagonal they are denoted block Toeplitz matrices. In this case we have only two blocks in the
diagonal, but the Toeplitz property become clearer in case the Hankel matrices are generalized to more
block rows. In this case the resulting Toeplitz matrices become constant along all the block diagonals.
We will start out by proving that the Toeplitz matrices can be added over a set of free decays without
changing the physics of the system. To this end, it is useful to realize that each column in the Hankel
matrices H1 and H2 is a free decay response of the discrete time state space response vector
6
Mathematics and Mechanics of Solids 00(0)
u(k) =
y(k)
y(k + 1)
ð18Þ
From equations (14) and (17) we see that we can obtain the Toeplitz matrix T11 by the summation of the
outer products of the state space response vectors
T11 =
np
X
u(k)uT (k)
ð19Þ
k =1
We then see from equation (4) using that the time delay Dt3 is equal to the sampling time, thus
el1 Dt3 = mn , and using a similar approach as when we arrived at equation (7), that
h
i
h
i
B n
u(k) = cnn q(k) = C n cnn q(k)
ð20Þ
B
h
i
h
i
is now B
= B n mnn , with n mnn denoting a diagonal matrix containwhere the earlier defined matrix B
ing the discrete time poles mn . By inserting equation (20) into equation (19), the Toeplitz matrix T11 is
becomes
T11 =
np
X
np
h
i
h
i
h
i X
i
h
C n cnn q(k)qT (k) n cnn CT = C n cnn
q(k)qT (k) n cnn CT
k=1
k=1
h
i
h
i
= C n cnn M n cnn CT = CPCT
ð21Þ
where the matrices M and P are given, respectively, by
np
X
M=
q(k)qT (k)
P=
h
k =1
n
i
h
n
cnn M cnn
i
ð22Þ
Similarly, we have for the Toeplitz matrix T21
T21 = H2 HT1
=
np
X
h
i
u(k + 2)uT (k) = C n m2nn PCT
ð23Þ
k =1
^ 1 using equation
The matrix P is a matrix that is canceled out whenever we estimate the system matrix S
(11) because from equations (21) and (23) we have (assuming full rank of the involved matrices)
^ 1 = H2 HT (H1 HT )1 = T21 T1
S
1
1
11
h
i
h
i
= C n m2nn PCT (CPCT )1 = C n m2nn C1
ð24Þ
This provides a more clear background for the modal decomposition shown before in equation (11). As
it appears from equation (22), the matrix P is determined by system poles and the modal participation
coefficients cn . It plays a central role because a mode can only be estimated if its modal participation
contributes significantly to this matrix. The other matrices in the matrix products of the expressions for
T11 and T21 , given by equations (21) and (23), always have full rank. Therefore the rank of the resulting
Toeplitz matrices is only limited by a possible rank deficiency of the matrix P.
It is easy to realize that, in case of non-repeated poles, the matrix M will have full rank if enough data
points are used in the time averaging process, i.e. if the data count variable np is large enough. Thus, if
enough data is being
used in the averaging process, the rank of P is only limited by the modal participa
tion matrix n cnn .
Brincker et al.
7
If an input is being used where a mode is not excited, i.e. cn = 0, it follows, from equation (22), that
the corresponding row and column of P is zero, yielding a rank deficient matrix. As a result, the rank of
the Toeplitz matrices reduces accordingly. In practice due to the noise that is always present in real data,
the modal participation does not need to be exactly zero in order to cause problems. If the modal participating is small and comparable to the modal participation of the noise modes, it might not be possible
to identify the mode, because it acts like a noise mode.
Because the matrix P play this role and might become rank deficient due to a lack of participation of
one or more modes, it is fair to say that this matrix describes the modal participation to the Toeplitz
matrices. Using the so defined participation matrix P, it is possible to realize that Toeplitz matrices from
two or more free decays can be added without changing the physics of the system, and that the addition
of the Toeplitz matrices might help in the identification of modes that are not well represented in some
of the free decays .
For simplicity, let us consider two free decays y1 (k) and y2 (k) where both free decays might have missing modes, but where all modes are well represented in the two free decays considered together. If we
add the Toeplitz matrices from the individual free decays then we have
T11 = C P1 CT + C P2 CT = C(P1 + P2 )CT = CPC
h
i
h
i
h
i
h
i
T21 = C n m2nn P1 CT + C n m2nn P2 CT = C n m2nn (P1 + P2 )CT = C n m2nn PCT
ð25Þ
where P1 and P2 are the modal participation matrices for the Toeplitz matrices of the individual free
decays and where P = P1 + P2 is the participation matrix for the added Toeplitz matrices. We see, using
equation (24), that we can do nearly whatever we want with the participation matrix, because this matrix
does not disturb the physics of the system, as long as it has full rank. We also realize, that if the matrices
P1 and P2 are rank deficient due to missing modes in the individual free decays, but the resulting matrix
P has full rank, then the system matrix will have full rank and, as a consequence, all modes can be estimated. This is also the case for cases with closely spaced modes and repeated poles.
Similarly the Toeplitz matrices T12 and T22 are given by
T12 = H1 HT2
=
T22 = H2 HT2 =
np
X
k =1
np
X
h
i
u(k)uT (k + 2) = CP n m2n n CT
h
i h
i
u(k + 2)uT (k + 2) = C n m2nn P n m2nn CT
ð26Þ
k =1
Again we see that contributions from several free decays corresponding to different inputs can be added
^ 2 gets full rank, and
to assure that the Toeplitz matrices and the corresponding system matrix estimate A
that all modes are present in the modal decomposition of the system matrix.
In Fukuzono [15], it was suggested to combine several free decays y1 (k), y2 (k) , , ynd (k), in the
matrix
ð27Þ
Y(k) = y1 ðk Þ y2 ðk Þ ynd ðk Þ
and then assemble the Hankel matrix
2
Y(1)
6 Y(2)
H=6
4 Y(3)
Y(4)
Y(2)
Y(3)
Y(4)
Y(5)
3
Y(np )
Y(np + 1) 7
7
Y(np + 2) 5
Y(np + 3)
ð28Þ
Here the number of DOFs nd might be different from the earlier mentioned number of DOFs in the
finite element (FE) model. It follows from the rules of matrix multiplication that this corresponds exactly
to adding the contribution from the different free decays when forming the Toeplitz matrix, as given by
equation (25).
8
Mathematics and Mechanics of Solids 00(0)
2.3. Correlation function matrices as free decays
The idea of using correlation function (CF) matrices for the OMA identification is that the columns – or
the rows – of the CF matrix can be interpreted as free decays of the system loaded by white noise. First,
let us consider the general form of a free day as given by equation (2) but now we write it as a summation
over the N modes
y(k) =
N X
cn bn eln tk + cn bn eln tk
ð29Þ
n=1
If we have several free decays y1 (k), y2 (k) , the modal participation of mode n of the first decay is
c1n , the modal participation of mode n of the second decay is c2n , etc. so that the modal participation of
mode n for the different free decays is described by the column vector cTn = fc1n , c2n , , cnd n g. The
response matrix given by equation (27) is then
Y(k) =
N X
ln tk
bn cTn eln tk + bn cH
ne
ð30Þ
n=1
It is common practice, for instance as described in Bendat and Piersol [4], to estimate the CF matrix
based on the definition
Ry (t) = E y(t)yT (t + t)
ð31Þ
It is well known that because of the stationarity of the response the correlation function matrix has the
symmetry properties
RTy (t) = Ry ( t)
ð32Þ
However, some authors, e.g. Brandt [5] or Papoulis and Pillai [6], use the alternative definition for the
CF matrix
~ y (t) = E y(t + t)yT (t)
R
ð33Þ
Using the stationary properties of the response and equations (31)–(33) we have
~ y (t) = E y(t)yT (t t) = RT (t)
R
y
ð34Þ
and we see that the alternative correlation function matrix is just the transpose of the classical definition.
The analytical solution for the CF matrix of the random responses of a structure was originally due to
James et al. [24], but has recently been generalized to the case of general damping and general Gaussian
white noise input, Brincker [19] and Brincker and Ventura [3]. The general solution for time lags greater
than zero is
Ry + (t) = 2p
N X
ln t
gn bTn eln t + gn bH
ne
ð35Þ
n=1
where gn are the modal participation vectors given by
gn = Bn Gx bn =an
in which Bn is a weighted sum of the residual matrices An = bn bTn =an
N X
As
As
Bn =
+
ln ls ln ls
s=1
ð36Þ
ð37Þ
Brincker et al.
9
The matrix Gx is the white noise input spectral density that is flat in the Nyquist band and zero elsewhere
and an is the generalized modal mass. We see that if the classical formulation of the correlation function
given by equation (31) is being used, then the CF matrix for positive time lags given by equation (35) is
not a free decay of the system, because the modal participation vectors are replacing the mode shapes in
equation (35). However, if we use the alternative formulation of the CF matrix given by equation (33),
that can also be obtained taking the transpose of the solution in equation (35), we get
~ y + (t) = 2p
R
N X
ln t
bn gTn eln t + bn gH
ne
ð38Þ
n=1
We see that this is indeed a free decay of the system with the modal participation vectors
gn =
1
cn
2p
ð39Þ
As a result, we can use the transposed form of the classical CF matrix as a basis for multiple-input
OMA defining the input matrix given in equation (29) as
~ y + (k) = RT (k)
Y(k) = R
y+
ð40Þ
Once the modal identification has been performed as described in the preceding section, we now have
some estimates of mode shapes and poles. The modal participation of each mode can then be found as
follows. Assuming the estimated mode shape matrix A = ½an contains the complex conjugate pairs of the
mode shapes in equation (35) (now the theoretical mode shapes bn are replaced by the estimated mode
shapes an ) we extend the summation to the set of 2N modal parameters and we sample the CF matrix at
the discrete time lags t = kDt and obtain the expression
Ry (k) = 2p
2N
X
gn aTn eln kDt
ð41Þ
n=1
This is a series of outer products that can be expressed by the matrix equation
h
ik
Ry (k) = 2p G n mn n AT
h
n
ð42Þ
i
where the diagonal matrix mnn as before contains the discrete time poles and the matrix G contains
the modal participation vectors G = ½gn ; both matrices are holding a set of complex conjugate pairs for
each mode. In Brincker and Ventura [3], it is shown that arranging the CF matrix in the single block row
Hankel matrix
H = ½ Ry (0)
Ry (1) Ry (2)
and the modal information in the matrix
h
i0
h
i1
n
M = n mnn AT
mnn AT
h
n
mnn
i2
A
T
Ry (np ) h
n
ð43Þ
mnn
i np
A
T
ð44Þ
the participation matrix can be obtained as
^ = 1 HM+
G
2p
where M+ is the pseudo inverse of M. The resulting fit of the CF matrix is
h
ik
T
^ nm
^ y (k) = 2p G
R
nn A
ð45Þ
ð46Þ
10
Mathematics and Mechanics of Solids 00(0)
The participation vector can be used to determine a participation factor for each mode. This is important in practical identification because noise modes normally have low modal participation factors, thus
a high participation factor
indicates that the corresponding mode might have physical importance.
lt
l t
Just
as
ce
+
c
e
can
also be thought of as harmonic with an amplitude that is proportional to
pffiffiffiffiffiffiffi
jcj = cc or as having an energy proportional to cc . Since energy is additive, we can take one of the
^ and compute its absolute scalar measure pn by means of the follow^ n in G
modal participation vectors g
ing expression
qffiffiffiffiffiffiffiffiffiffiffi
^H
^n
pn = g
ð47Þ
ng
A relative modal participation factor can then be calculated as
pn =
p2n
pT p
ð48Þ
where the column vector
p contains the absolute modal participation factors as defined by equation (47)
pT = p1 , p2 , pnp . It follows from the definition that the factors pn adds up to 100% over all the
considered modes.
3. Simulation study
The following simulation case illustrates the performance of the multiple-input, multiple-output
(MIMO) formulation of ITD compared to the original formulation. The case is based on the three
degree of freedom (3 DOF) system given in Table 1. The system is constructed so that the natural frequencies of the first mode f1 = 1Hz and of the third mode f3 = 2Hz are fixed. The natural frequency f2
of the second mode varies from f1 to f3 according to f2 = f1 + Df where Df 2 ½0; 1Hz. The damping ratio
is 1% for all three modes. The variation of mode 2 results in systems with closely spaced modes when
mode 2 is close either to mode 1 or mode 3, and with well separated modes when mode 2 is in the mid
region between mode 1 and mode 3.
The 3 DOF system is loaded by white noise in all three DOFs. As input for the identification the exact
theoretical solution for the CF matrix given in equation (39) is being used. That is to be sure that the
simulation study is not influenced by errors from estimation of the correlation functions.
The exact CF matrix is computed for 40 time lags with a time step of 0.1 s corresponding to four
cycles of mode 1 and eight cycles of mode 2 reducing the responses to approximately 78% and 60%,
respectively, of the initial value. For the SIMO formulation of ITD the first column of the CF matrix is
being used and for the MIMO formulation the full matrix is being used.
For each variation of mode 2 a new CF matrix is computed using equation (38) where the participation vectors are obtained from equation (36) and where the constant input spectral density matrix Gx is
taken as a random matrix that is forced to be positive definite by taking the absolute value of its
eigenvalues.
Table 1. System parameters used in the simulation study.
Type:
Mode 1
fixed
Mode 2
variable
Mode 3
fixed
Natural frequency (Hz)
f1 = 1
f2 = 3
Damping ratio
Mode shape vector
z1 = 0:01
2
3
1
b1 = 4 0 5
1
f2 = f1 + Df
Df 2 ½0; 1
z2 = 0:01
2
3
1
b2 = 4 1 5
1
z3 = 0:01
2 3
1
b3 = 4 1 5
1
Brincker et al.
11
Figure 1. The relative error of the estimated frequencies.
Table 2. The relative error e of the estimated natural frequencies.
ITD SIMO
ITD MIMO
Mean (%)
Standard deviation (%)
Mean (%)
Standard deviation (%)
Mode 1
Mode 2
Mode 3
0.014
0.021
0.003
0.003
0.066
0.293
0.006
0.009
0.133
0.206
0.037
0.039
Using the solution for the CF matrix given by equation (38) leads to estimated modal parameters that
will be almost exact. Therefore, a small amount of noise is added to the exact CF matrix. In this study is
used a white noise with a standard deviation of 0.1% of the maximal value of each individual correlation
function. This assures a reasonable variability of the estimated modal parameters.
In order to compare the SIMO and MIMO formulations the relative error of the estimated natural
frequencies is calculated as
abs fi fj
eij =
ð49Þ
fi
where fi are the exact natural frequencies of the system and fj are the estimated natural frequencies.
The results of the simulations study for the errors on the estimated natural frequencies are shown in
Figure 1. As it appears from Figure 1, the relative errors of the estimated frequencies for the SIMO case
are significantly larger than the errors for the MIMO case. This is especially true in case of closely spaced
modes, i.e. when mode 2 gets close to mode 1 or mode 3. The mean and standard deviation of the relative estimation error are shown in Table 2. We see significantly larger estimation errors for the SIMO
formulation compared to the MIMO one.
In order to evaluate the correlation between the mode shape vectors the generalized angles aij between
the exact mode shape vectors bj and the estimated mode shape vectors ai are obtained as
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1
0v
u u aH bj 2
i
1 @t
A
aij = cos
ð50Þ
H
aH
i ai bj bj
12
Mathematics and Mechanics of Solids 00(0)
Figure 2. The generalized angle a between the estimated and the system mode shape vectors.
The results of the simulations study for the generalized angles between exact and estimated mode
shapes are shown in Figure 2. As it appears from Figure 2, the generalized angles for the SIMO case are
significantly larger than the angles for the MIMO case. This is especially true in case of closely spaced
modes, i.e. when mode 2 gets close to mode 1 or mode 3.
4. Conclusions
A multiple-input formulation of the ITD technique that is based on addition of Toeplitz matrices over
a set free decays has been investigated, and it has been shown how a special modal participation matrix
can be defined so that the contribution of each free decay is added to form the resulting participation
matrix. The contribution to this matrix is rank deficient if a mode is not present in the free decay, but it
is shown that adding the Toeplitz matrices over all free decays will assure the full rank of the resulting
system matrix if all modes are present in one of the considered free decays.
In case of OMA applications, it is possible to use the CF matrix in the form where it represents free
decays of the system (normally the transposed of the CF matrix), and to obtain a modal participation
factor for the whole CF matrix that makes it easier to distinguish between physical modes and noise
modes. Finally, the performance of the proposed multiple-input ITD technique has been illustrated on
a case with three degrees of freedom, where the natural frequency of one of the modes are being varied,
so that the performance of the technique can be studied in cases of closely and well separated modes.
The results clearly show that both natural frequencies and mode shapes are estimated more accurately
when the multiple-input formulation is being used. Specifically – as it normally observed in practice –
using the multiple-input formulation it is easier to identify closely spaced modes.
Funding
This work was supported by the Centre for Oil and Gas – DTU/Danish Hydrocarbon Research and Technology Centre
(DHRTC), DTU, Denmark and by Aarhus University, Denmark.
References
[1]
Zhang, L, Brincker, R, and Andersen, P. An overview of operational modal analysis: major developments and issues. In:
Proceedings of the 1st International Operational Modal Analysis Conference (IOMAC), Copenhagen, Denmark, 26–27
April 2005.
Brincker et al.
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
13
Masjedian, MH, and Keshmiri, M. A review on operational modal analysis researches: classification of methods and
applications. In: Proceedings of the 3rd International Operational Modal Analysis Conference (IOMAC), Portonova, Spain,
4–6 May 2009.
Brincker, R, and Ventura, CE. Introduction to operational modal analysis. Wiley, 2015.
Bendat, JS, and Piersol, AG. random data: analysis and measurement procedures. 4th ed. Wiley, 2010.
Brandt, A. Noise and vibration analysis: signal analysis and experimental procedures. Wiley, 2011.
Papoulis, A, and Pillai, SU. Probability, random variables and stochastic processes. 4th ed. McGraw-Hill Europe, 2002.
Ibrahim, SR, and Milkulcik, EC. A method for direct identification of vibration parameters from the free response. Shock
Vib Bull 1977; 47: 183-196.
Ibrahim, SR. Random decrement technique for modal identification of structures. J Spacecraft Rockets 1977; 14: 696–700.
Ibrahim, SR. Modal confidence factor in vibration testing. J Spacecraft Rockets 1978; 15: 313–316.
Vold, H, Kundrat, J, Rocklin, GT, et al. A multi-input modal estimation algorithm for mini-computers. SAE paper 1982820194, 1982.
Vold, H, and Rocklin, GT. The numerical implementation of a multi-input modal estimation method for mini-computers.
In: Proceedings of the International Modal Analysis Conference (IMAC), Location?, Date? 1982, pp.542–548.
Juang, JN, and Pappa, RS. An eigen system realization algorithm for modal parameter identification and modal
reduction. J Guidance Control Dyn 1985; 8: 620–627.
Pappa, RS, Elliott, KB, and Schenk, A. Consistent-mode indicator for the eigensystem realization algorithm. J Guidance
Control Dyn 1993; 16: 852–858.
Pappa, RS. Eigensystem realization algorithm, user’s guide for VAX/VMS computers. Technical memorandum 109066,
NASA, 1994.
Fukuzono, K. Investigation of multiple-reference Ibrahim time domain modal parameter estimation technique. MSc Thesis,
University of Cincinnati, USA, 1986.
Malekjafarian, A, Brincker, R, Ashory, MR, et al. Identification of closely spaced modes using Ibrahim time domain
method. In: Proceedings of the International Operational Modal Analysis Conference (IOMAC), Istanbul, 2011.
Malekjafarian, A, Brincker, R, Ashory, MR, et al. Identification of closely spaced modes using Ibrahim time domain
method: experimental results. In: Proceedings of the International Modal Analysis Conference (IMAC), Garden Grove,
February 11–14, 2013.
Olsen, P, and Brincker, R. Using random response input in Ibrahim time domain. In: Proceedings of the 31st International
Modal Analysis Conference (IMAC), Garden Grove, February 11–14, 2013.
Brincker, R. On the application of correlation function matrices in OMA. J Mech Syst Signal Process 2017; 87: 17–22.
Cook, RD. Concepts and applications of finite element analysis. John Wiley & Sons, 2007.
Ewins, DJ. Modal testing: theory and practice. Vol. 15. Letchworth, UK: Research Studies Press, 1984.
Ibrahim, SR. Fundamentals of time domain identification. In: Silva, JMM, and Maia, NMM (eds.) Modal analysis and
testing. NATO Science Series, Series E: Applied Sciences – Vol 363, Kluwer Academic Publishers, 1999.
James, GH, Carne, TG, and Lauffer, JP. The natural excitation technique (NExT) for modal parameter extraction from
operating structures. Int J Anal Exp Modal Anal 1995; 10: 260–277.
Ibrahim, SR. Double least squares approach for use in structural modal identification. AIAA J 1986; 24: 3.