A New Model Updating Method for Quadratic
Eigenvalue Problems
Yuen-Cheng Kuo
∗Wen-Wei Lin
†Shu-Fang Xu
‡Abstract
Finite element model updating of quadratic eigenvalue problems (QEPs) is proposed by Friswell, Inman and Pilkey 1998, to incorporate the measured model data into the finite element model to produce an adjusted finite element model on the damping and stiffness that closely match the experimental modal data. In this paper, we mainly develop an efficient numerical method for the finite element model updating of QEPs which needs only O(nk2) flops and is stored in
a sparse technique, where n is the size of coefficient matrices of the QEP and k is the number of the measured eigenpairs.
∗National Center for Theoretical Sciences Mathematics Division, Hsinchu, 300, Taiwan
†Department of Mathematics, National Tsinghua University, Hsinchu, 300, Taiwan
‡LMAM, School of Mathematical Sciences, Peking University, Beijing, 100871, China
1
Introduction
Given n × n real matrices M, C and K, the task of finding scalars λ ∈ C and nonzero vectors x ∈ Cn satisfying
Q(λ)x ≡ (λ2M + λC + K)x = 0 (1.1)
is known as the quadratic eigenvalue problem (QEP), which corresponds to solving the homogeneous second-order differential equation (see e.g., [10])
Mq(t) + C ˙q(t) + Kq(t) = 0.¨ (1.2)
The scalar λ and the associated vector x in (1.1) are called, respectively, eigenvalues and eigenvectors of the quadratic pencil Q(λ). It is known that the QEP has 2n finite eigenvalues, provided the leading matrix M is nonsingular. Recently the QEP has received much attention because its information has repeatedly arisen in many different discipline, including applied mechanics, electrical oscillation, vibro-acoustics, fluid mechanics and signal processing. A nice survey paper for the QEP can be found in [14] by Tisseur and Meerbergen. Vibrating systems, such as automotives, bridges, highways, buildings are described by distributed parameters. However, due to lack of viable computational methods to handle distributed parameter systems, a finite element method is generally used to discretize such systems to an analytical model (finite element model), namely,
Qa(λ) = λ2Ma+ λCa+ Ka, (1.3)
where Ma, Ca and Ka represent the mass, damping and stiffness, respectively, that all
are real n × n symmetric matrices with Ma being symmetric positive definite (Denoted
by Ma>0). See the book [9] by Friswell and Mottershead for details.
In the finite element model (1.3) for structural dynamics, the mass and stiffness are, in general, clearly defined by physical parameters. However, the damping for precise
dissipative effects is not well understood because it is a purely dynamics property that can not be measured statically. A common simplification is to assume proportional or modal damping, but it seems to be sufficient where damping levels are lower than 10% of critical [8].
Finite element model updating has emerged in the 1990’s as a significant subject to the design, construction, and maintenance of mechanical systems [9, 12]. Model updating, at its ambitious, is used to correct inaccurate analytical models by measured data, such as natural frequencies, damping ratios, mode shapes and frequency re-sponse function, which can usually be obtained by vibration test. In the past decade, Baruch/Bar-Itzak [1, 2], Bermann/Nagy [3, 4] and Wei [15, 17, 16] considered vari-ant aspects of finite element model updating by using measured data for undamped
structured systems (i.e. C = Ca = 0). In the works by Datta/Elhay/Ram/Sarkissian
[5, 6, 7], studies are undertaken toward a feedback design problem for second-order con-trol system. That consideration eventually leads to a partial eigenstructure assignment problem for the QEP. Recently, Friswell, Inman and Pilkey [8] proposed to incorporate the measured model data into the finite element model to produce an adjusted finite element model on the damping and stiffness with modal properties that closely match
the experimental modal data. That is, with M = Ma the penalty function
J = 1 2kN −1(K − K a)N−1k2F + 1 2µkN −1(C − C a)N−1k2F, (1.4a) is minimized, subject to MaΦΛ2+ CΦΛ + KΦ = 0, (1.4b) C> = C, K> = K, (1.4c) where N = M12
a, µis a weighting parameter, C and K are the updated damping and
and eigenvalue matrices, respectively. In practice, k is much less than the matrix size n. The solutions K and C of (1.4) are given by [8, 13] with,
K = Ka− 2MaRe(ΓΛΦ>+ ΦΓ>Λ)Ma (1.5) and C = Ca− 2 µMaRe(ΓΛΛΦ >+ ΦΛΓ> Λ)Ma, (1.6)
where ΓΛ ∈ Cn×k solves the 2nk linear equations
2MaRe(ΓΛΦ>+ ΦΓ>Λ)MaΦ + 2 µMaRe(ΓΛΛΦ > + ΦΛΓ>Λ)MaΦΛ =MaΦΛ2+ CaΦΛ + KaΦ. (1.7)
Generally, the size n in a finite element model (1.3) is quite large. It is impractical to solve a complex matrix ΓΛ for a large and dense 2nk × 2nk linear system in (1.7).
The purpose of this paper is to develop an efficient algorithm for the computation of the solutions C and K for (1.4) which is required only O(nk2) flops and is stored is a
sparse technique.
2
Solving a PD-IQEP.
To match the partial measured data of the spectrum information of a QEP, we consider to solve the partially described inverse quadratic eigenvalue problem (PD-IQEP): Let (Λ, Φ) ∈ Rk×k× Rn×k (k ≤ n) be a given pair of matrices, where
Λ = diag(λ[2]1 , . . . , λ [2] ` , λ2`+1, . . . , λk) (2.1a) with λ[2]j = αj βj −βj αj , βj 6= 0, for j = 1, . . . , `, and Φ = [ϕ1R, ϕ1I, . . . , ϕ`R, ϕ`I, ϕ2`+1, . . . , ϕk]. (2.1b)
Suppose Λ has only simple eigenvalues and Φ is of full column rank. Find a general form of symmetric matrices M , C and K with M being symmetric positive definite that satisfy the equation
MΦΛ2+ CΦΛ + KΦ = 0. (2.2)
Let Φ have the QR-decomposition Φ = Q R 0 ≡ [Q1, Q2] R 0 , (2.3)
where Q ∈ Rn×n is orthogonal with Q
1 ∈ Rn×k and R ∈ Rk×k is nonsingular. Partition
Q>M Q, Q>CQand Q>KQ by Q>M Q= M11 M12 M21 M22 , Q>CQ = C11 C12 C21 C22 , Q>KQ = K11 K12 K21 K22 , (2.4) where M11, C11 and K11 ∈ Rk×k, and M21, C21 and K21∈ R(n−k)×k.
A general solution of symmetric M , C, K with M being symmetric positive definite is given by the theorem in [11].
Theorem 2.1. [11] For a given matrix pair (Λ, Φ) as in (2.1), the general solution of PD-IQEP is given by
M11 : arbitrary fixed symmetric positive definite matrix, (2.5a)
C11= − M11S+ S>M11+ R−>DR−1 , (2.5b) K11 = S>M11S+ R−>DΛR−1, (2.5c) K21 = K12> = − M21S2+ C21S , (2.5d) where S = RΛR−1 and D= diag ξ1 η1 η1 −ξ1 , . . . , ξ` η` η` −ξ` , ξ2`+1, . . . , ξk (2.6)
in which ξi and ηi are arbitrary real numbers. Furthermore, M21 = M12>, C21 = C12>,
C22 = C22> and K22 = K22> are chosen arbitrary, and M22 = M22> is chosen so that
M22− M21M11−1M12 is symmetric positive definite.
3
Solving the optimization problem
(1.4).
In this section we shall develop an efficient algorithm for solving the optimization problem described in (1.4). We first solve two optimization problems. Let (Λ, Φ) ∈ Rk×k×Rn×k be given in (2.1), D and R be given in (2.6) and (2.3), respectively. Denote
R−1 = [r1, . . . , rk] = r11 . . . r1k . .. ... 0 rkk . (3.1)
Optimization Problem I. Given A = [a1, . . . , ak], B = [b1, . . . , bk] ∈ Rk×k and let
x= (ξ1, η1, . . . , ξ`, η`, ξ2`+1, . . . , ξk) >
(3.2) correspond to the matrix D in (2.6). Minimize
f(x) = µkA + R−>DR−1k2 F + kB − R −>Λ>DR−1k2 F = k X j=1 fj(x) (3.3a) for x, where fj(x) = µkaj + R−>Drjk22+ kbj − R−>Λ>Drjk22, j = 1, . . . , k. (3.3b) Note that ξ η η −ξ u v = u v −v u ξ η
. The vector Drj in (3.3b) can be rewritten
by
where (i) j = 2s, s 6 `, Γj = diag r1j r2j −r2j r1j , . . . , r2s−1,j r2s,j −r2s,j r2s−1,j , 0, . . . , 0 ∈ Rk×k, (ii) j = 2s + 1, s < `, Γj = diag r1j r2j −r2j r1j , . . . , r2s−1,j r2s,j −r2s,j r2s−1,j , r2s+1,j, r2s+1,j,0, . . . , 0 ∈ Rk×k, (iii) j > 2`, Γj = diag r1j r2j −r2j r1j , . . . , r2`−1,j r2`,j −r2`,j r2`−1,j , r2`+1,j, . . . , rj,j,0, . . . , 0 ∈ Rk×k .
Substituting (3.4) into (3.3b) we compute ∇fj(x) = ∂fj ∂x1 , . . . , ∂fj ∂xk > = 2µ R−>Γ j > aj+ R−>Γjx − 2 R−>Λ>Γ j > bj− R−>Λ>Γjx . (3.5) Consequently, ∇f (x) = k X j=1 ∇fj(x) = 2 k X j=1 h µ R−>Γj > aj + µΓ>j R >R−1Γ jx− Γ>j ΛR −1b j+ Γ>jΛ R >R−1Λ>Γ jx i . (3.6) Setting ∇f (x) = 0 we derive the linear equation for x
where G= k X j=1 h µΓ> j R >R−1Γ j + Γ>jΛ R >R−1Λ>Γ j i (3.8a) and b= k X j=1 Γ> j ΛR −1b j− µΓ>jR −1a j . (3.8b)
Since the function f (x) in (3.3a) must have an optimum, the linear system of (3.7) is consistent, and therefore, x is solvable.
Optimization Problem II. Given E, F ∈ R(n−k)×k and S = RΛR−1. Minimize
g(X) = µkE − Xk2F + kF + XSk2F (3.9)
for X ∈ R(n−k)×k.
Let
x= vec (X) , e = vec (E) , f = vec (F ) . (3.10)
Here “ vec ” denotes the vectorization of the given matrix columnwisely. By (3.10), the function g(X) in (3.9) becomes
g(x) = g(X) = µke − xk22+ kf + S>⊗ I n−k xk22. (3.11) Then ∇g(x) = 2−µ(e − x) + (S ⊗ In−k) f + S>⊗ In−k x. (3.12)
Here ⊗ denotes the Kronecker product of two matrices. The matrix form of (3.12) can be written by ∂ ∂Xg(X) = 2 −µE + µX + F S>+ XSS>. (3.13) By setting ∂ ∂Xg(X) = 0, we then solve X = (µE − F S>)(µI + SS>)−1. (3.14)
We now solves the optimization problem (1.4). Redefine Ca := M −1 2 a CaM −1 2 a , Ka:= M −1 2 a KaM −1 2 a , (3.15a) C := M−12 a CM −1 2 a , K := M −1 2 a KM −1 2 a , (3.15b) Φ := M12 aΦ. (3.15c)
Let Q = [Q1, Q2] be the orthogonal matrix such that Φ = Q[R>,0]> with R
nonsingu-lar. Then the problem (1.4) becomes
min = 1 2µkCa− Ck 2 F + 1 2kKa− Kk 2 F = 1 2µ Q >C aQ− C11 C12 C21 C22 2 F + 1 2 Q >K aQ− K11 K12 K21 K22 2 F , (3.16)
where Cij = Q>i CQj and Kij = Q>i KQj, i, j = 1, 2, are undetermined. Note that
M = In. Set
C22= Q>2CaQ2, K22 = Q>2KaQ2. (3.17)
Then from (2.5) the optimization problem (3.16) is equivalent to the following two optimization problems: min =µ Q> 1CaQ1− C11 2 F + Q> 1KaQ1− K11 2 F (3.18) =µ Q> 1CaQ1+ S + S>+ R−>DR−1 2 F + Q> 1KaQ1− S>S− R−>DΛR−1 2 F and min =µ Q> 2CaQ1− C21 2 F + Q> 2KaQ1− K21 2 F (3.19) =µ Q> 2CaQ1− C21 2 F + Q> 2KaQ1+ C21S 2 F.
Hence, the optimal solutions C11 and K11 of (3.18) are solved by the Optimization
Problem I by setting
A= Q>1CaQ1+ S + S>, B = Q>1KaQ1− S>S. (3.20)
The optimal solutions C21and K21of (3.19) are solved by the Optimization Problem
II by setting E = Q> 2CaQ1 and F = Q>2KaQ1. (3.21) Reset M = Ma, C = M12 aQ C11 C > 21 C21 C22 Q>M12 a, K = M 1 2 aQ K11 K > 21 K21 K22 Q>M12 a (3.22)
which solve the optimization problem (1.4). The steps of computation for solving (1.4) are summarized in the following algorithm.
Algorithm 3.1. Given Qa(λ) = λ2Ma+ λCa+ Ka and (Λ, Φ) ∈ Rk×k × Rn×k as in
(2.1). The optimal solutions C and K of (1.4) are computed by
step 1. Set Ca:= M −1 2 a CaM −1 2 a , Ka := M −1 2 a KaM −1 2 a , Φ := M 1 2 aΦ ;
step 2. Compute the QR-factorization of Φ :
Φ = [Q1, Q2] R 0 , and S = RΛR−1; step 3. Compute C22 = Q>2CaQ2, K22= Q>2KaQ2;
where G= k X j=1 Γ> j h µ R>R−1+ Λ R>R−1Λ>iΓ j, b= k X j=1 Γ> j ΛR−1vj− µR−1uj , Γj = diag r1,j r2,j −r2,j r1,j , · · · , r2`−1,j r2`,j −r2`,j r2`−1,j , r2`+1,j,· · · , rk,j , [u1,· · · , uk] = Q>1CaQ1+ S + S>, [v1,· · · , vk] = Q>1KaQ1− S>S, (r1,j,· · · , rk,j)>= R−1ej; step 5. Compute C11 = − S + S>+ R−>DR−1 , K11= S>S+ R−>DΛR−1, where D = diag ξ1 η1 η1 −ξ1 , · · · , ξ` η` η` −ξ` , ξ2`+1,· · · , ξk , and compute C21= Q>2 µCaQ1− KaQ1S> (µI + SS>)−1, K21 = −C21S; step 6. C = M12 aQ C11 C12 C21 C22 Q>M12 a, K = M 1 2 aQ K11 K12 K21 K22 Q>M12 a, where Q = [Q1, Q2].
Remark 3.1. (i) In a finite element model, the size of the analytical matrices Ma, Ca
and therefore, it is easily invertible. In practice, the number of the measured data of eigenpairs is much less than the size of the finite element model, i.e., k n. The orthogonal matrix Q = [Q1, Q2] in step 2 of Algorithm 3.1 can be computed and stored
in the form of a diagonal matrix plus a low rank updating by Householder transforma-tions. Suppose the sparse matrix Ca or Ka times a vector needs O(n) flops. Then, the
computational cost of Algorithm 3.1 can be easily estimated by O(nk2) flops.
(ii) Using Algorithm 3.1 to solve the optimization problem (1.4) is different from using (1.7) to solve (1.4). The latter needs to solve a large, but possibly dense nk × nk linear system as in (1.7) which is impractical in a finite element model updating process when n is sufficiently large.
4
Numerical results.
A set of pseudo simulation data was provided by The Boeing Company for testing purpose. The dimension of matrices Ma, Ca and Ka are 42 × 42.
Test I. We first test that the Algorithm 3.1 computes the optimal solution of the optimization problem (1.4). Choosing an eigenmatrix pair (Λa,Φa) ∈ R14×14× R42×14
of the analytical model Qa= λ2Ma+ λCa+ Ka, i.e., MaΦaΛ2a+ CaΦaΛa+ KaΦa = 0,
the Algorithm 3.1 should theoretically give the optimal solution C = Ca and K = Ka.
The numerical result of the relative errors computed by Algorithm 3.1 are estimated by kC − CakFa kCakFa ' 10−7, kK − KakFa kKakFa ' 10−10,
where k · kFa = kM
−12 a (·)M
−12 a kF.
Test II. Now we are given the measured eigenvalues
{λmj}14j=1 = { − 0.60939 ± 37.365ι, −0.73496 ± 36.707ι, −2.8832 ± 31.992ι,
− 1.8907 ± 61.437ι, −1.9112 ± 54.181ι, −2.2785 ± 39.639ι, (4.1)
− 5.0387, −4.3416}
and the measured mode shapes vj ∈ Rs, j = 1, . . . , 14. The measured eigenvectors ϕj
is estimated by
ϕj = DjDej†vj, j = 1, . . . 14, (4.2)
where Dj is defined by Dj = [λ2mjMa+ λmjCa+ Ka]−1Ba with Ba ∈ Rn×t (t < s), and
the matrix eDj consistent of the first s rows of Dj and the superscript “ † ” denotes the
pseudo inverse. We construct the eigenmatrix pair (Λ, Φ) ∈ R14×14× R42×14 as in (2.1)
associated with (4.1) and (4.2). The Algorithm 3.1 computes the new updated matrices M = Ma, C and K with µ = 0.1, 1.0 and 10, which minimizes the optimization problem
(1.4). We define the relative residual by
res = kM ΦΛ 2+ CΦΛ + KΦk 2 kM ΦΛ2k 2+ kCΦΛk2+ kKΦk2 . (4.3)
The numerical results are shown in Table 4.1.
Table 4.1 relative residuals
µ 0.1 1.0 10
res 1.4725e-014 1.4826e-014 1.4859e-014
5
Conclusions.
We have developed an efficient numerical algorithm for finite element model updating of quadratic eigenvalue problems. This method can serve as a fast and reliable manner
for updating the analytical model. It was shown to be insightful in a simple pseudo test suit provided by The Boring Company.
References
[1] M. Baruch, Optimization procedure to correct stiffness and flexibility matrices using vibration data. AIAA Journal, 16(11):1208–1210, 1978.
[2] M. Baruch and I. Y. Bar-Itzack, Optimal weighted othogonalization of measured modes. AIAA Journal, 16(4):346–351, 1978.
[3] A. Berman, Comment on ‘optimal weighted othogonalization of measured modes’. AIAA Journal, 17(8):927–928, 1979.
[4] A. Berman and E. J. Nagy, Improvement of a large analytical model using test data. AIAA Journal, 21(8):1168–1173, 1983.
[5] B. N. Datta, Finite element model updating, eigenstructure assignment and eigen-value embedding techniques for vibrating systems, mechanical systems and signal processing. Special Volume on Vibration Control, 16:83–96, 2002.
[6] B. N. Datta, S. Elhay, Y. M. Ram and D. R. Sarkissian, Partial eigenstructure assignment for the quadratic pencil. Journal of Sound and Vibration, 230:101–110, 2000.
[7] B. N. Datta and D. R. Sarkissian, Theory and computations of some inverse eigen-value problems for the quadratic pencil. in Structured Matrices in Mathematics, Computer Science, and Engineering. I, Contemp. Math. 280, AMS, Providence, RI, pages 221–240, 2001.
[8] M. I. Friswell, D. J. Inman and D. F. Pilkey, The direct updating of damping and stiffness matrices. AIAA Journal, 36(3):491–493, 1998.
[9] M. I. Friswell and J. E. Mottershead, Finite element model updating in structural dynamics. Kluwer Academic Publishers, 1995.
[10] I. Gohberg, P. Lancaster and L. Rodman, Matrix Polynomials. Academic Press, New York, 1982.
[11] Y. C. Kuo, W. W. Lin and S. F. Xu, On a general solution of
partially described inverse quadratic eigenvalue problems and its
applica-tions. NCTS Tech-Report 2005-1-2, National Tsing-Hua Uni., Taiwan.
http://math.cts.nthu.edu.tw/Mathematics/preprints/prep2005-1-001.pdf.
[12] J. E. Mottershead and M. I. Friswell, Model updating in structural dynamics: A survey. Journal of Sound and Vibration, 167(2):347–375, 1993.
[13] D. F. Pilkey, Computation of a Damping Matrix for finite Element Model Updat-ing. Dissertation of Virginia Polytech. Inst. and State Uni., 1998.
[14] F. Tisseur and K. Meerbergen, The quadratic eigenvalue problem. SIAM Review, 43:253–286, 2001.
[15] F.-S. Wei, Structural dynamic model identification using vibration test data. 7th
IMAC, Las Vegas, Nevada, pages 562–567, 1989.
[16] F.-S. Wei, Mass and stiffness interaction effects in analytical model modification. AIAA Journal, 28(9):1686–1688, 1990.
[17] F.-S. Wei, Structural dynamic model improvement using vibration test data. AIAA Journal, 28(1):175–177, 1990.
