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: firstname.lastname@example.org 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. , Masjedian and Keshmiri , and Brincker and Ventura . 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 . 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  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 , 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  and Brincker  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 . 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  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 , 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  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 . 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 . 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 , 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 , 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 , 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  or Papoulis and Pillai , 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. , but has recently been generalized to the case of general damping and general Gaussian white noise input, Brincker  and Brincker and Ventura . 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 , 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 pﬃﬃﬃﬃﬃﬃﬃ 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 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ^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 ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ1 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  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.                        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.