一个$n$自由度的无阻尼振动系统可由下面二阶微分方程来刻画

$ M_a \ddot{q}(t)+K_a q(t)=0, $

(1.1)

其中$M_a, \ K_a \in {\bf R}^{n \times n}$分别是分析质量矩阵与刚度矩阵, $q(t)$是$n\times 1$位移向量.通常$M_a$是实对称正定矩阵, $K_a$是实对称半正定矩阵, 即$M_a=M_a^\top>0$, $K_a=K_a^\top\geq0$.方程(1.1)通常被称为有限元分析模型.设$q(t)=x {\rm e}^{{\rm i}\omega t}$是方程(1.1)的基础解, 则自然频率$\omega_i$和模态向量$x_i$满足下面的广义特征值方程

$ \lambda_i M_a x_i = K_a x_i, \ i=1, 2, \cdots, n, $

(1.2)

其中$\lambda_i=\omega_i^2$是第$i$个特征值, $x_i$是第$i$个特征向量.设

$ \Lambda=\left[ \begin{array}{cc} \Lambda_1~~&0 \\ 0 ~~& \Lambda_2 \\ \end{array} \right], \ \ X=[X_1, X_2], $

其中$\Lambda_1=\mbox{diag}\{\lambda_1, \cdots, \lambda_p\}, $ $\Lambda_2=\mbox{diag}\{\lambda_{p+1}, \cdots, \lambda_n\}, $ $ X_1=[x_1, \cdots, x_p], $ $X_2= [x_{p+1}, \cdots, x_n].$易见(1.2)式的$n$个特征值---特征向量关系式可以组装为如下的矩阵方程

$ M_aX\Lambda=K_aX. $

(1.3)

无阻尼振动模态的一个重要性质是它们关于质量矩阵的正交性, 即

$ X^\top M_a X=I_n. $

(1.4)

当对无阻尼振动系统施加一个主动控制力后, 方程(1.1)即为

$ M_a \ddot{q}(t)+K_a q(t)=Bu(t), $

(1.5)

其中$B\in{\bf R}^{n \times m}$是列满秩控制反馈矩阵, $u(t)\in{\bf R}^{m \times 1}$是控制向量.若取$u(t)$为

$ u(t)=-G\ddot{q}(t)-Fq(t), $

(1.6)

其中$G, F\in{\bf R}^{m \times n}$分别是加速度和位移反馈增益矩阵.把方程(1.6)代入方程(1.5)得到下面的闭环系统

$ (M_a+BG)\ddot{q}(t)+(K_a+BF)q(t)=0, $

(1.7)

其对应的特征结构方程为

$ (M_a+BG)Y\Sigma=(K_a+BF)Y, $

(1.8)

其中$\Sigma=\mbox{diag}\{\mu_1, \cdots, \mu_n\}$和$Y=[y_1, \cdots, y_n]$分别是闭环系统(1.7)对应的特征值和特征向量.修正一个有限元模型不仅要使得修正后的系数矩阵是对称的, 同时还要保持原模型中未知的或者不能测量的高阶模态数据是不变的, 即修正是不溢出的.用加速度和位移反馈控制同时修正质量和刚度矩阵的问题可以描述如下.

问题 1  设$M_a, K_a\in{\bf R}^{n \times n}$分别是分析质量矩阵和刚度矩阵. $B\in{\bf R}^{n \times m}$是列满秩控制矩阵, $\Sigma_1=\mbox{diag}\{\mu_1, \cdots, \mu_p\}\in {\bf R}^{p \times p}, \ Y_1=[y_1, \cdots, y_p]\in{\bf R}^{n \times p}$分别是测量的特征值和特征向量矩阵, 其中$p\ll n.$求矩阵$G$和$F$使得

$ M X_2\Lambda_2=KX_2, \ \ M Y_1 \Sigma_1=KY_1, \ \ \mbox{s.t.} \ \ M=M^\top, \ \ K=K^\top, $

(1.9)

其中$M=M_a+BG$, $K=K_a+BF$.通常, 问题1的解是不唯一的, 这样, 我们需进一步求解下面的最佳逼近问题.

问题 2  设${\mathcal{S}}_{E}$是问题1的解集.求$ (\hat{G}, \hat{F}) \in {\mathcal{S}}_{E}$使得

$ \| B \hat{G}\|^2+\| B \hat{F} \|^2 = \min\limits_{(G, F) \in {\mathcal{S}}_{E}} (\| BG \|^2+\|BF \|^2). $

(1.10)

一旦求得问题2的解, 则修正的质量矩阵和刚度矩阵可表示为

$ \hat{M}=M_a+ B\hat{G}, \ \ \hat{K}=K_a+B\hat{F}. $

(1.11)

有限元模型修正问题已经有相当多的讨论, 也建立了许多求解方法.有限元模型被广泛用于预测结构的动态特性.为了确保有限元模型的有效性, 工程师通常需要找到一种方法来修正现有的有限元模型以便修正后的模型能较好地反映出从物理结构测量的模态数据^[1].在过去的40多年里, 利用测量数据来修正无阻尼系统的质量矩阵与刚度矩阵已有许多方法. Baruch^[2]和Baruch和Bar-Itzhack^[3]利用拉格朗日乘子法导出了修正的刚度矩阵和柔度矩阵的解析解. Berman和Nagy^[4]给出了利用测量的正规模态和自然频率修正分析模型的质量矩阵和刚度矩阵的方法. Wei^[5-6]和Dai^[7]提出了基于约束极小化理论和最优逼近理论同时修正质量和刚度矩阵的方法. Carvalho等^[8]利用非完备的测量数据来修正刚度矩阵.对于阻尼结构系统修正问题, 在近些年也得到了较快的发展, Mao和Dai^[9]提出了一种用非完备的测量数据同时修正质量、阻尼和刚度矩阵的方法. Yuan和Dai^[10]给出了同时修正阻尼和刚度矩阵的一种数值解法. Bai等^[11]提出了一种牛顿型对偶优化方法来修正阻尼结构系统.在文献[12-13]中, 作者讨论了对称低秩情况下的质量、阻尼和刚度矩阵的修正问题.

此外, 在结构动力学问题中, 特征值配置和特征结构配置问题也有很多潜在的应用.在工程中, 振动通常是需控制的, 以防止机器损坏或过早疲劳.对于振动控制已有很多方法, 其中改变固有频率或增加阻尼是常用的方法. Saad^[14]给出了一阶系统部分极点配置的投影算法.对二阶控制系统, Xu和Qian^[15]研究了基于状态反馈控制的鲁棒部分特征值配置问题. Cai等^[16]考虑了高阶控制系统的部分极点配置问题. Cai和Xu^[17]讨论了对称特征值嵌入问题, 其中不需要的特征值被新的特征值所取代, 而其余的特征值保持不变.文献[18-22]也对部分鲁棒极点配置问题进行了讨论.在文献[23]中, Zhang等通过加速度和位移反馈控制技术讨论了无阻尼结构系统的特征结构配置问题.然而, 他们的方法并不能确保修正后的质量和刚度矩阵是对称的. Yuan等^[24]考虑了基于加速度和位移反馈的无阻尼结构系统的特征结构配置问题, 据此给出了修正质量和刚度矩阵的一种直接方法, 但他们的方法并不能保证修正是不溢出的.本文提出了一种有效的迭代方法来求解问题1和2.通过选择合适的初始矩阵对, 可以得到问题2的极小Frobenius范数解$(\hat{G}, \hat{F})$.最后给出的两个数值例子显示提出的方法是可靠有效的.

本文记${\bf R}^{m \times n}$为$m \times n$矩阵空间, ${\bf R}^{n \times n}$中的对称矩阵全体记为${\bf SR}^{n \times n}$. $A^{\top}, $ $\mbox{tr}(A)$和$\mathcal{R}(A)$分别表示矩阵$A$的转置, 迹及列向量所展成的空间. $I_n$表示$n$阶单位矩阵.矩阵空间$ {\bf R}^{m \times n}$的内积定义为$(A, B) =\mbox{tr}(B^\top A), \ \ \forall A, B \in {\bf R}^{m \times n}, $则${\bf R}^{m \times n}$是一个Hilbert空间, 此内积诱导的矩阵范数${\| \cdot \|}$即为Frobenius范数.给定矩阵$A=[a_{ij}] \in {\bf R}^{m \times n}$, $ B=[b_{ij}] \in {\bf R}^{p \times q}$, $A$和$B$的Kronecker积定义为$A \otimes B=[a_{ij}B]\in {\bf R}^{mp \times nq}.$此外, 对$m\times n$矩阵$A=[a_1, a_2, \cdots, a_n]$, 其中$a_j, j=1, \cdots, n, $是$A$的第$j$个列向量, 拉直函数$\mbox{vec}(A)$定义为$\mbox{vec}(A)=[a_1^\top, a_2^\top, \cdots, a_n^\top]^\top. $设$A, B$和$X$是三个适当维数的矩阵, 则可得如下的等式^[25]

$ {\rm{vec}}(AXB)=(B^\top\otimes A){\rm{vec}}(X). $

由(1.3)和(1.4)式可得

$ M_a^{-1}=X_1X_1^\top+X_2X_2^\top, \ \ M_a^{-1}K_a M_a^{-1}=X_1\Lambda_1 X_1^\top+X_2 \Lambda_2 X_2^\top. $

(2.1)

由(1.3)、(2.1)式和$X_2$是列满秩阵, 则方程(1.9)的第一个等式即为

$ BGP=BFL, $

(2.2)

其中$P=M_a^{-1}K_a M_a^{-1}-X_1\Lambda_1 X_1^\top, \ L=M_a^{-1}-X_1X_1^\top.$在工程实践中, 方程(1.9)中的矩阵$X_2$和$\Lambda_2$通常是未知的, 注意到在方程(2.2)中没有出现矩阵$X_2$和$\Lambda_2$.为了使$BF$和$BG$是对称的, 令$F=\tilde{S}B^\top$, $G=\tilde{H}B^\top$, 其中$\tilde{S}, \ \tilde{H}$是对称矩阵.则

$ BG=B\tilde{H}B^\top, $

(2.3)

$ BF=B\tilde{S}B^\top. $

(2.4)

设$B$的QR -分解如下

$ B=Q\left[ \begin{array}{c} R \\ 0 \\ \end{array} \right], $

(2.5)

其中$Q=[Q_1, Q_2]$是$n\times n$正交矩阵$(Q_1\in {\bf R}^{n\times m})$, $R$是$m\times m$非奇异矩阵.把方程(2.5)代入方程(2.3)可得

$ B\tilde{H}B^\top=Q_1R\tilde{H}R^\top Q_1^\top. $

设$R\tilde{H}R^\top=H$, 则

$ BG=B\tilde{H}B^\top =Q_1HQ_1^\top, \ \ G=R^{-1}HQ_1^\top, $

(2.6)

其中$ H \in {\bf SR}^{m\times m}$是对称矩阵.

同理把方程(2.5)代入方程(2.4)并记$R\tilde{S}R^\top=S$, 则得

$ BF=B\tilde{S}B^\top =Q_1SQ_1^\top, \ \ F=R^{-1}SQ_1^\top, $

(2.7)

其中$ S \in {\bf SR}^{m\times m}$是对称矩阵.

现在求解问题1等价于找对称矩阵$S$和$H$使得

$ Q_1SQ_1^\top P-Q_1HQ_1^\top L=0, $

(2.8)

$ Q_1SQ_1^\top Y_1 \Sigma_1 -Q_1HQ_1^\top Y_1=K_aY_1-M_a Y_1 \Sigma_1. $

(2.9)

为了方便, 令$A_1=Q_1^\top P, \ B_1=Q_1^\top L, \ A_2=Q_1^\top Y_1 \Sigma_1, \ B_2=Q_1^\top Y_1$和$C_2=Q_1^\top (K_a Y_1-M_a Y_1 \Sigma_1), $则方程(2.8)和(2.9)等价于下面的方程

$ SA_1-HB_1=0, $

(2.10)

$ SA_2-HB_2=C_2, $

(2.11)

$ Q_2^\top (K_aY_1-M_aY_1 \Sigma_1)=0. $

(2.12)

假设方程(2.10)和(2.11)是相容的, 且条件(2.11)成立, 容易验证$S_E$是非空的闭凸集.由最优逼近定理^[26]可知问题2存在唯一解.由

$ \|BG\|^2+\|BF\|^2=\|Q_1SQ_1^\top\|^2+\|Q_1HQ_1^\top\|^2\\ =\|Q^\top Q_1SQ_1^\top Q\|^2+\|Q^\top Q_1HQ_1^\top Q\|^2\\ =\left\| \left[ \begin{array}{cc} S ~~ 0 \\ 0 ~~ 0 \\ \end{array} \right] \right\|^2+\left\| \left[ \begin{array}{cc} H ~~ 0 \\ 0 ~~ 0 \\ \end{array} \right] \right\|^2\\ =\|S\|^2+\|H\|^2. $

所以求解问题1和2等价于找方程(2.10)和(2.11)的极小Frobenius范数对称解.下面我们给出求解矩阵方程(2.10)和(2.11)的对称解的迭代算法.

算法 2.1   1) 输入矩阵$M_a, K_a, B, Y_1, \Sigma_1$, 任意选取初始对称矩阵$\tilde{S_1}\in{{\bf SR}}^{m \times m}, $ $\tilde{H_1}\in{{\bf SR}}^{m \times m}$.

2) 由(2.5)式计算矩阵$B$的QR -分解.

3) 计算矩阵

$ P=M_a^{-1}K_a M_a^{-1}-X_1\Lambda_1 X_1^\top, \ L=M_a^{-1}-X_1X_1^\top, \ \\A_1=Q_1^\top P, \ B_1=Q_1^\top L, $

$ A_2=Q_1^\top Y_1 \Sigma_1, \ B_2=Q_1^\top Y_1, \ C_2=Q_1^\top (K_a Y_1-M_a Y_1 \Sigma_1), $

$ S_1=R\tilde{S_1}R^\top, \ H_1=R\tilde{H_1}R^\top. $

4)

$ E_1=H_1B_1-S_1A_1, \ F_1=C_2-S_1A_2+H_1B_2, $

$ R_1=\left[ \begin{array}{cc} E_1~~ 0 \\ 0 ~~ F_1 \\ \end{array} \right], $

$ M_1=E_1A_1^\top+F_1A_2^\top, \ N_1=-(E_1B_1^\top+F_1B_2^\top), $

$ U_1=\frac{1}{2}(M_1+M_1^\top), \ V_1=\frac{1}{2}(N_1+N_1^\top). $

5) 如果$R_k=0$, 或者$R_k\neq0 $但$ U_k=0, \ V_k=0$, 则停止; 否则进行下一步.

6)

$ a_k=\frac{\|R_k\|^2}{\|U_k\|^2+\|V_k\|^2}, \ S_{k+1}=S_k+a_kU_k, \ H_{k+1}=H_k+a_kV_k, $

$ E_{k+1}=H_{k+1}B_1-S_{k+1}A_1, \ F_{k+1}=C_2-S_{k+1}A_2+H_{k+1}B_2, $

$ {R_{k + 1}} = \left[ {\begin{array}{*{20}{c}} {{E_{k + 1}}}&0\\ 0&{{F_{k + 1}}} \end{array}} \right], $

$ M_{k+1}=E_{k+1}A_1^\top+F_{k+1}A_2^\top, \ N_{k+1}=-(E_{k+1}B_1^\top+F_{k+1}B_2^\top), $

$ b_k=\frac{\mbox {tr}(M_{k+1}U_k)+\mbox {tr}(N_{k+1}V_k)}{\|U_k\|^2+\|V_k\|^2}, $

$ U_{k+1}=\frac{1}{2}(M_{k+1}+M_{k+1}^\top)-b_kU_k, \ V_{k+1}=\frac{1}{2}(N_{k+1}+N_{k+1}^\top)-b_kV_k. $

由算法2.1易知, $S_k\in{{\bf SR}}^{m \times m}$和$H_k\in{{\bf SR}}^{m \times m}, \ k=1, 2, \cdots.$

引理2.1  由算法$2.1$产生的数列$\{R_i\}, \{U_i\}, \{V_i\}, \{M_i\}$和$\{N_i\}, $满足关系式

$ {\rm{tr}}\left(R_{i+1}^{\top}R_j\right)={\rm{tr}}\left(R_i^{\top}R_j\right)-a_i\left[{\rm{tr}}(U_iM_j)+{\rm{tr}}(V_iN_j)\right], \ \ i, j=1, 2, \cdots . $

(2.13)

证 由算法2.1可得

$ E_{i+1} = H_{i+1}B_1-S_{i+1}A_1 =(H_i+a_iV_i)B_1-(S_i+a_iU_i)A_1 \\=E_i+a_i(V_iB_1-U_iA_1). $

同理可得

$ F_{i+1}=F_i+a_i(V_iB_2-U_iA_2). $

所以

$ {\rm{tr}}\left(R_{i+1}^{\top}R_j\right) ={\rm{tr}}\left(R_i^{\top}R_j\right)\\+a_i{\rm{tr}}\left( \left[ \begin{array}{cc} B_1^\top V_iE_j-A_1^\top U_iE_j&0 \\ 0&B_2^\top V_iF_j-A_2^\top U_iF_j\\ \end{array} \right] \right) \\ ={\rm{tr}}\left(R_i^{\top}R_j\right)+a_i\left[{\rm{tr}}\left(B_1^\top V_iE_j-A_1^\top U_iE_j\right) + {\rm{tr}}\left(B_2^\top V_iF_j-A_2^\top U_iF_j \right)\right] \\ ={\rm{tr}}\left(R_i^{\top}R_j\right)+a_i\left[{\rm{tr}}\left(V_i(E_jB_i^\top+F_jB_2^\top)\right)-{\rm{tr}}\left(U_i(E_jA_i^\top+F_jA_2^\top)\right)\right] \\ ={\rm{tr}}\left(R_i^{\top}R_j\right)-a_i\left[{\rm{tr}}(U_iM_j)+{\rm{tr}}(V_iN_j)\right]. $

证毕.

引理 2.2  由算法$2.1$产生的数列$\{R_i\}, \{U_i\}, \{V_i\}, $满足关系式

$ {\rm{tr}}(R_i^\top R_j)=0, \ {\rm{tr}}(U_i U_j)+{\rm{tr}}(V_i V_j)=0, \ \ i, j=1, 2, \cdots, \ i\neq j. $

(2.14)

证   因为$\mbox {tr}(R_{i}^{\top}R_j)=\mbox {tr}(R_j^{\top}R_{i})$, 所以我们只需证明

$ {\rm{tr}}(R_i^\top R_j)=0, \ {\rm{tr}}(U_i U_j)+{\rm{tr}}(V_i V_j)=0, \ i>j. $

下面用数学归纳法证明引理2.2.当$j=1 $, $i=2$时, 由引理2.1和算法2.1可得

$ {\rm{tr}}(R_2^{\top}R_1) =\|R_1\|^2-a_1{\rm{tr}} (U_1M_1+V_1N_1)\\ =\|R_1\|^2-a_1 {\rm{tr}}\left(U_1\frac{M_1+M_1^\top}{2}+V_1\frac{N_1+N_1^\top}{2}\right)\\ =\|R_1\|^2-a_1{\rm{tr}} (U_1U_1+V_1V_1)\\ =\|R_1\|^2-a_1\left(\|U_1\|^2+\|V_1\|^2\right)\\ =0, $

$ {\rm{tr}}(U_2U_1)+\mbox {tr}(V_2V_1) ={\rm{tr}}\left[\left(\frac{1}{2}(M_2+M_2^\top)-b_1U_1\right)U_1\right]\\+\mbox {tr}\left[\left(\frac{1}{2}(N_2+N_2^\top)-b_1V_1\right)V_1\right] \\ ={\rm{tr}}(M_2U_1-b_1U_1U_1)+\mbox {tr}(N_2V_1-b_1V_1V_1)\\ ={\rm{tr}}(M_2U_1)+\mbox {tr}(N_2V_1)-b_1(\|U_1\|^2+\|V_1\|^2)\\ =0. $

假设(2.14)式在$j=i-1>1$时成立, 即

$ {\rm{tr}}(R_i^\top R_{i-1})=0, \ {\rm{tr}}(U_i U_{i-1})+{\rm{tr}}(V_i V_{i-1})=0, $

则

$ {\rm{tr}}\left(R_{i+1}^{\top}R_i\right) =\|R_i\|^2-a_i{\rm{tr}}(U_iM_i+V_iN_i)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;=\|R_i\|^2-a_i{\rm{tr}}\left(U_i\frac{M_i+M_i^\top}{2}+V_i\frac{N_i+N_i^\top}{2}\right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;=\|R_i\|^2-a_i{\rm{tr}}\left[U_i(U_i+b_{i-1}U_{i-1})+V_i(V_i+b_{i-1}V_{i-1})\right]\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;=\|R_i\|^2-a_i(\|U_i\|^2+\|V_i\|^2)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;=0, $

${\rm{tr}}(U_{i+1}U_i)+{\rm{tr}}(V_{i+1}V_i)\\ ={\rm{tr}}\left[\left(\frac{1}{2}(M_{i+1}+M_{i+1}^\top)-b_iU_i\right)U_i\right]+ {\rm{tr}}\left[\left(\frac{1}{2}(N_{i+1}+N_{i+1}^\top)-b_iV_i\right)V_i\right]\\ = {\rm{tr}}(M_{i+1}U_i)+{\rm{tr}}(N_{i+1}V_i)-b_i(\|U_i\|^2+\|V_i\|^2) \\ =0. $

由数学归纳法可知, 对所有的正整数$j>0$都有

$ {\rm{tr}}(R_{j+1}^\top R_{j})=0, \ {\rm{tr}}(U_{j+1} U_{j})+{\rm{tr}}(V_{j+1} V_{j})=0. $

假设当$j=1, 2, \cdots, i-1$时, (2.14)式成立, 即

$ {\rm{tr}}(R_i^\top R_j)=0, \ {\rm{tr}}(U_i U_j)+{\rm{tr}}(V_i V_j)=0, $

则

$ {\rm{tr}}(R_{i+1}^\top R_j) ={\rm{tr}}(R_i^\top R_j)-a_i{\rm{tr}}(U_iM_j+V_iN_j) \\ =-a_i{\rm{tr}}\left(U_i\frac{M_j+M_j^\top}{2}+ V_i\frac{N_j+N_j^\top}{2}\right) \\ =-a_i{\rm{tr}}\left[U_i(U_j+b_{j-1}U_{j-1})+V_i(V_j+b_{j-1}V_{j-1})\right]\\ =-a_i{\rm{tr}}\left[U_iU_j+V_iV_j+b_{j-1}(U_iU_{j-1}+V_iV_{j-1})\right]\\ =0. $

$ {\rm{tr}}(U_{i+1}U_j) ={\rm{tr}}\left[\left( \frac{1}{2}(M_{i+1}+M_{i+1}^\top)-b_iU_i \right) U_j\right] \\ ={\rm{tr}} \left(M_{i+1}U_j\right)-b_i{\rm{tr}}(U_iU_j). $

同理可得

$ {\rm{tr}}(V_{i+1}V_j)={\rm{tr}} \left(N_{i+1}V_j\right)-b_i{\rm{tr}}(V_iV_j). $

则

$ {\rm{tr}}(U_{i+1} U_j)+{\rm{tr}}(V_{i+1} V_j) \\={\rm{tr}}(M_{i+1}U_j)+{\rm{tr}}(N_{i+1}V_j)-b_i[{\rm{tr}}(U_iU_j)+{\rm{tr}}(V_iV_j)]\\ ={\rm{tr}}(M_{i+1}U_j)+{\rm{tr}}(N_{i+1}V_j)\\ =\frac{1}{a_j} \left[{\rm{tr}}(R_j^\top R_{i+1})-{\rm{tr}}(R_{j+1}^\top R_{i+1})\right]\\ =0. $

由归纳原理即得引理2.2.

引理 2.3  若$(S^*, H^*)$是方程$(2.10)$和$(2.11)$的任意对称解.则对任意的对称矩阵对$(S_1, H_1)$, 其中$S_1, H_1 \in {\bf SR}^{m \times m}$, 可得

$ {\rm{tr}}\left[(S^*-S_i) U_i\right]+{\rm{tr}}\left[(H^*-H_i) V_i\right]=\|R_i\|^2, \ \ \ i=1, 2, \cdots, $

(2.15)

其中数列$ \{S_i\}, \{H_i\}, \{U_i\}, \{V_i\}$和$\{R_i\}$由算法$2.1$生成.

证 我们用数学归纳法证明引理2.3.当$i=1$时, 有

$ {\rm{tr}}\left[(S^*-S_1)U_1\right]+{\rm{tr}}\left[(H^*-H_1)V_1\right]\\ ={\rm{tr}}\left[(S^*-S_1)M_1\right]+{\rm{tr}}\left[(H^*-H_1)N_1\right]\\ ={\rm{tr}}\left[(S^*-S_1)(E_1A_1^\top+F_1A_2^\top)\right]-{\rm{tr}}\left[(H^*-H_1)(E_1B_1^\top+F_1B_2^\top)\right]\\ ={\rm{tr}}(E_1^\top E_1+F_1^\top F_1)\\ =\|R_1\|^2. $

假设当$i=t$时, (2.15)式成立, 即

$ {\rm{tr}}\left[(S^*-S_t) U_t\right]+{\rm{tr}}\left[(H^*-H_t) V_t\right]=\|R_t\|^2. $

因为

$ \;\;{\rm{tr}}\left[(S^*-S_{t+1})U_t\right]+{\rm{tr}}\left[(H^*-H_{t+1})V_t\right]\\ ={\rm{tr}}\left[(S^*-S_t-a_tU_t)U_t\right]+{\rm{tr}}\left[(H^*-H_t-a_tV_t)V_t\right]\\ ={\rm{tr}}\left[(S^*-S_t)U_t\right]+{\rm{tr}}\left[(H^*-H_t)V_t\right]-a_t\left(\|U_t\|^2+\|V_t\|^2\right)\\ =0. $

当$i=t+1$时, 有

$ {\rm{tr}}\left[(S^*-S_{t+1}) U_{t+1}\right] ={\rm{tr}}\left[(S^*-S_{t+1}) \left(\frac{1}{2}(M_{t+1}+M_{t+1}^\top)-b_tU_t\right)\right]\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;={\rm{tr}}\left[(S^*-S_{t+1})M_{t+1}\right]-b_t{\rm{tr}}\left[(S^*-S_{t+1})U_t\right]. $

同理可得

$ {\rm{tr}}\left[(H^*-H_{t+1}) V_{t+1}\right] ={\rm{tr}}\left[(H^*-H_{t+1})N_{t+1}\right]-b_t{\rm{tr}}\left[(H^*-H_{t+1})V_t\right]. $

则

$ \;\;\;{\rm{tr}}\left[\left((S^*-S_{t+1}) U_{t+1}\right)\right]+{\rm{tr}}\left[\left((H^*-H_{t+1}) V_{t+1}\right)\right]\\ ={\rm{tr}}\left[\left((S^*-S_{t+1}) M_{t+1}\right)\right]+{\rm{tr}}\left[\left((H^*-H_{t+1}) N_{t+1}\right) \right]\\ =\|R_{t+1}\|^2. $

由归纳原理即得引理2.3.

定理2.1  设方程$(2.10)$和$(2.11)$是相容的.对任意的初始矩阵对$(S_1, H_1)$, $S_1, H_1 \in {\bf SR}^{m \times m}$, 则在不计舍入误差的情况下, 经有限步迭代就可得到方程$(2.10)$和$(2.11)$的一个对称解.

证  若$R_l\neq 0, \ l=1, 2, \cdots, m(n+p)$.由引理2.2可知$U_l\neq 0$或$V_l\neq 0$.由算法2.1可以计算$R_{ m(n+p)+1}$和$(S_{ m(n+p)+1}, H_{ m(n+p)+1})$.由引理2.4可知

$ {\rm{tr}}(R_{ m(n+p)+1}^\top R_t)=0, \ t=1, 2, \cdots, m(n+p) $

和

$ {\rm{tr}}(R_j^\top R_i)=0, \ \ i, j=1, 2, \cdots, m(n+p), \ i\neq j. $

因此, $\{R_1, R_2, \cdots, R_{ m(n+p)}\}$构成矩阵空间$\Phi$的一组正交基, 其中

$ \Phi=\left\{T\left|T=\left[ \begin{array}{cc} X ~~& 0 \\ 0 ~~& Y \\ \end{array} \right], \ X\in {\bf R}^{m \times n}, \ Y\in {\bf R}^{m \times p}\right. \right\}, $

这意味着$R_{ m(n+p)+1}=0$, 即$(S_{ m(n+p)+1}, H_{ m(n+p)+1})$是方程$(2.10)$和$(2.11)$的对称解.

引理2.4^[27] 设相容线性方程组$Ax=b$有解$\tilde{x} \in \mathcal{R}(A^\top)$, 那么$\tilde{x}$是此线性方程组的唯一极小Frobenius范数解.

引理2.5^[28] 方程$(2.10)$和$(2.11)$有对称解当且仅当下面方程组

$ \begin{array}{l} SA_1-HB_1=0, \\ SA_2-HB_2=C_2, \\ A_1^\top S-B_1^\top H=0, \\ A_2^\top S-B_2^\top H=C_2^\top \end{array} $

(2.16)

是相容的.

应用Kronecker积和拉直函数, 可知矩阵方程组(2.16)等价于

$ \left[\begin{array}{cc} A_1^\top \otimes I_m, ~~&-B_1^\top \otimes I_m \\ A_2^\top \otimes I_m, ~~&-B_2^\top \otimes I_m \\ I_m \otimes A_1^\top, ~~&-I_m \otimes B_1^\top \\ I_m \otimes A_2^\top, ~~& -I_m \otimes B_2^\top \end{array} \right]\left[\begin{array}{c} {\rm{vec}}(S) \\ {\rm{vec}}(H) \end{array} \right]=\left[\begin{array}{c} 0 \\ {\rm{vec}}(C_2)\\ 0 \\ {\rm{vec}}(C_2^\top) \end{array} \right]. $

设$W_1\in {\bf R}^{m\times p}, Z_1 \in {\bf R}^{m\times n}$是任意矩阵, 则

$ \left[\begin{array}{cc} {\rm{vec}}(Z_1A_1^\top+W_1A_2^\top+A_1Z_1^\top+A_2W_1^\top) \\ {\rm{vec}}(-(Z_1B_1^\top+W_1B_2^\top+B_1Z_1^\top+B_2W_1^\top)) \end{array} \right]\\ = \left[\begin{array}{cccc} A_1 \otimes I_m, ~~ A_2 \otimes I_m, ~~I_m \otimes A_1, ~~ I_m \otimes A_2 \\ -B_1 \otimes I_m, ~~-B_2 \otimes I_m, ~~-I_m \otimes B_1, ~~-I_m \otimes B_2 \end{array} \right]\left[\begin{array}{c} {\rm{vec}}(Z_1)\\ {\rm{vec}}(W_1) \\ {\rm{vec}}(Z_1^\top)\\ {\rm{vec}}(W_1^\top) \end{array} \right] \\ = \left[\begin{array}{cc} A_1^\top \otimes I_m, ~~-B_1^\top \otimes I_m \\ A_2^\top \otimes I_m, ~~-B_2^\top \otimes I_m \\ I_m \otimes A_1^\top, ~~-I_m \otimes B_1^\top \\ I_m \otimes A_2^\top, ~~ -I_m \otimes B_2^\top \end{array} \right]^\top \left[\begin{array}{c} {\rm{vec}}(Z_1)\\ {\rm{vec}}(W_1) \\ {\rm{vec}}(Z_1^\top)\\ {\rm{vec}}(W_1^\top) \end{array} \right]\\ \in\mathcal{R}\left( \left[\begin{array}{cc} A_1^\top \otimes I_m, ~~-B_1^\top \otimes I_m \\ A_2^\top \otimes I_m, ~~-B_2^\top \otimes I_m \\ I_m \otimes A_1^\top, ~~-I_m \otimes B_1^\top \\ I_m \otimes A_2^\top, ~~ -I_m \otimes B_2^\top \end{array} \right]^\top \right). $

若选取

$ S_1=Z_1A_1^\top+W_1A_2^\top+A_1Z_1^\top+A_2W_1^\top, $

(2.17)

$ H_1=-(Z_1B_1^\top+W_1B_2^\top+B_1Z_1^\top+B_2W_1^\top), $

(2.18)

由算法2.1得到的$S_i, H_i, \ i=2, 3, \cdots, $满足

$ \left[\begin{array}{l} {\rm{vec}}(S_i) \\ {\rm{vec}}(H_i) \end{array} \right] \in \mathcal{R}\left( \left[\begin{array}{cc} A_1^\top \otimes I_m, ~~&-B_1^\top \otimes I_m \\ A_2^\top \otimes I_m, ~~&-B_2^\top \otimes I_m \\ I_m \otimes A_1^\top, ~~&-I_m \otimes B_1^\top \\ I_m \otimes A_2^\top, ~~& -I_m \otimes B_2^\top \end{array} \right]^\top \right). $

由$R\tilde{H}R^\top =H_1, $ $ R\tilde{S}R^\top =S_1$可得

$ \tilde S = {R^{ - 1}}({Z_1}A_1^ \top + {W_1}A_2^ \top + {A_1}Z_1^ \top + {A_2}W_1^ \top ){R^{ - \top }}, $

(2.19)

$ \tilde{H}=-R^{-1}(Z_1B_1^\top+W_1B_2^\top+B_1Z_1^\top+B_2W_1^\top)R^{-\top}. $

(2.20)

综上可得下面结果.

定理2.2  设$M_a, K_a\in{\bf R}^{n \times n}$分别是分析质量矩阵和刚度矩阵, $B\in {\bf R}^{n \times m}$是列满秩矩阵, $\Sigma_1=\mbox{diag}\{\mu_1, \cdots, \mu_p\}\in {\bf R}^{p \times p}, \ Y_1=[y_1, \cdots, y_p]\in{\bf R}^{n \times p}$分别是测量的特征值和特征向量矩阵.令$A_1=Q_1^\top (M_a^{-1}K_a M_a^{-1}-X_1\Lambda_1 X_1^\top), \ B_1=Q_1^\top (M_a^{-1}-X_1X_1^\top), $ $\ A_2=Q_1^\top Y_1 \Sigma_1, \ B_2=Q_1^\top Y_1 \ \mbox{和} \ C_2=Q_1^\top (K_aY_1-M_a Y_1 \Sigma_1)$ .假设方程是(2.10)和(2.11) 相容的, 且条件$(2.12)$成立.据$(2.19)$和$(2.20)$式选取初始矩阵对, 其中$W_1\in {\bf R}^{m\times p}, \ Z_1 \in {\bf R}^{m\times n}$是任意的矩阵, 特别地, 可以取$W_1=0, \ Z_1=0$, 那么在不计舍入误差的情况下, 经过有限迭代后可得方程$(2.10)$和$(2.11)$的唯一极小Frobenius范数对称解$(\hat{H}, \hat{S})$, 从而问题$2$的唯一解即为

$ \hat{G}=R^{-1}\hat{S}Q_1^\top, \ \hat{F}=R^{-1}\hat{H}Q_1^\top, $

(2.21)

修正的质量矩阵和刚度矩阵由$(1.11)$式给出.

算法 2.2   1) 输入矩阵$M_a, K_a, B, \Lambda_1, X_1, Y_1, \Sigma_1 $.

2) 由式(2.5)计算矩阵$B$的QR-分解.

3) 计算矩阵$A_1=Q_1^\top (M_a^{-1}K_a M_a^{-1}-X_1\Lambda_1 X_1^\top), $ $ B_1=Q_1^\top (M_a^{-1}-X_1X_1^\top), $ $ A_2=Q_1^\top Y_1 \Sigma_1, $ $ B_2=Q_1^\top Y_1, $ $ C_2=Q_1^\top (K_aY_1-M_a Y_1 \Sigma_1).$

4) 如果方程$(2.10)$和$(2.11)$是相容的且满足条件$(2.12)$, 则进行下一步; 否则停止计算且问题无解.

5) 根据式$(2.19)$和$(2.20)$选取初始矩阵对, 根据算法2.1计算方程$(2.10)$和$(2.11)$的极小Frobenius范数对称解.

6) 根据式(2.21)计算$\hat{G}, \hat{F}$.

7) 根据式(1.11)计算$\hat{M}, \hat{K}$.

在这一节, 我们给出两个数值例子来验证算法的有效性, 所有的测试均使用计算软件MATLAB R2011b.

例 3.1  考虑一个$6$ -自由度系统, 其质量矩阵, 刚度矩阵和控制矩阵分别为

$ M_a=\left[ \begin{array}{cccccc} 1.56&0.66&0.54&-0.39&0&0 \\ 0.66&0.36&0.39&-0.27 &0 &0 \\ 0.54&0.39&3.12&0&0.54 &-0.39 \\ -0.39&-0.27&0&0.72&0.39&-0.27 \\ 0&0&0.54&0.39&3.12&0 \\ 0&0&-0.39&-0.27&0&0.72 \\ \end{array} \right], $

$ K_a=\left[ \begin{array}{cccccc} 12&18&-12&18&0&0 \\ 18&36&-18&8&0&0 \\ -12&-18&24&0&-12&18 \\ 18&18&0&72&-18&18 \\ 0&0&-12&-18&24&0 \\ 0&0&18&18&0&72 \\ \end{array} \right], $

$ B=\left[ \begin{array}{ccc} -0.0037 & 0.0903 & 0.3687 \\ -0.0018 & 0.0319 & 0.0312 \\ -0.0047 &-0.1127 & -0.9881 \\ 0.0007 &-0.0624 & 0.1236 \\ -0.0015 &-0.1789 & 1.3131 \\ 0.0005 & 0.0322 & 0.2278 \\ \end{array} \right]. $

开环系统特征值为$\Lambda = {\mbox {diag}}\{0.0363, 1.4365, 11.4697, 58.1668, 206.023, 818.8383\}$, 其中$\Lambda_1 =\mbox {diag}\{0.0363, 1.4365, 11.4697\}$, $\Sigma_1 =\mbox {diag}\{0.0333, 1.3168, 10.5139 \}$,

$ X_1=\left[ \begin{array}{ccc} 0.5636&-0.5669&-0.5662 \\ -0.0862&0.3014&0.5000 \\ 0.3082&0.2401&0.3714 \\ -0.0820&0.1864&-0.0798 \\ 0.0933&0.3344 &-0.4212 \\ -0.0566&-0.1111&-0.1045 \\ \end{array} \right], \ \ \\Y_1=\left[ \begin{array}{ccc} 0.5144 &-0.5175&-0.5168\\ -0.0787 & 0.2751& 0.4565\\ 0.2814 & 0.2191& 0.3391\\ -0.0748 & 0.1702& -0.0728\\ 0.0852 & 0.3053& -0.3845\\ -0.0517 &-0.1014& -0.0954 \end{array} \right]. $

当误差$\|R_k\|\leq 1e-010$时, 迭代停止.根据算法2.2并选取初始矩阵对$\tilde{S}_1=0$, $\tilde{H}_1=0$, 经迭代46步后, 可得问题2的唯一解

$ \hat{G}=\left[ \begin{array}{cccccc} -0.3192&-0.1511&-0.4080 & 0.0637& -0.1331 & 0.0434\\ 0.0085& 0.0030& -0.0106 &-0.0059& -0.0169 & 0.0030\\ -0.0334& -0.0028& 0.0896 &-0.0112& -0.1191 &-0.0207 \end{array} \right], $

$ \hat{F}=\left[ \begin{array}{cccccc} -1.1287 &-0.5344&-1.4430& 0.2254&-0.4706& 0.1536\\ 0.6997 & 0.2476& -0.8732& -0.4835& -1.3868& 0.2493\\ 0.0006 & 0.0000& -0.0015& 0.0002& 0.0020& 0.0004 \end{array} \right]. $

最优修正的质量和刚度矩阵为

$ \hat{M}=M_a+B\hat{G}= \left[ \begin{array}{rrrrrr} 1.5679 & 0.6602& 0.5144&-0.3863 & 0.0342 & 0.0057\\ 0.6602 & 0.3598& 0.3876&-0.2695 & 0.0031 & 0.0005\\ 0.5144 & 0.3876& 3.1850&-0.0087 & 0.4486 &-0.4051\\ -0.3863 &-0.2695& -0.0087& 0.7207 & 0.4005 &-0.2679\\ 0.0342 & 0.0031& 0.4486& 0.4005 & 3.2365 & 0.0211\\ 0.0057 & 0.0005& -0.4051& -0.2679 & 0.0211 & 0.7235 \end{array} \right], $

$ \hat{K}=K_a+B\hat{F}= \left[ \begin{array}{rrrrrr} 11.9486& 17.9815&-11.9437& 18.0338 & 0.0934 &-0.0168\\ 17.9815& 35.9933&-17.9807& 18.0120 & 0.0330 &-0.0059\\ -11.9437&-17.9807& 23.9188&-0.0405 &-12.1190& 18.0222\\ 18.0338& 18.0120&-0.0405& 71.9769 &-18.0658& 18.0117\\ 0.0934& 0.0330&-12.1190&-18.0658 &23.8086 &0.0338\\ -0.0168& -0.0059& 18.0222& 18.0117 & 0.0338 &71.9938 \end{array} \right]. $

此时误差为

$ \|\hat{M}X_2\Lambda_2-\hat{K}X_2\|= 6.8566e-012, \ \|\hat{M}Y_1\Sigma_1-\hat{K}\Sigma_1\|= 1.5724e-011. $

上式显示, 指定的特征值和特征向量已嵌入新模型$\hat{M} Y_1\Sigma_1=\hat{K}Y_1$, 而修正系统的剩余谱数据与原模型的剩余谱数据是一致的, 也即修正是不溢出的, 并且修正的矩阵$\hat{M}, \ \hat{K}$仍是对称正定的.

例 3.2  设分析质量矩阵和刚度矩阵如下

$ M_a=I_n, \ K_a=\left[ \begin{array}{cccccc} 1&~~-1~~&0 &~~ \cdots ~~& 0&~~0 \\ -1&2&-1&\cdots&0&~~0 \\ 0&-1&2&\cdots&0&~~0 \\ \vdots&\vdots&\vdots&\ddots&\vdots&~~\vdots \\ 0&0&\cdots&-1&2&~~-1 \\ 0&0&\cdots&0&-1&~~1 \\ \end{array} \right]\in {\bf R}^{n \times n}. $

实验仿真模型为$M=M_a$和$K=0.85K_a$, 并取仿真模型的前五阶模态数据作为测量实验数据.在试验中, 选取控制反馈矩阵$B=K_aY_1-M_aY_1\Sigma_1$, 初始迭代矩阵为$\tilde{S}_1=0, \tilde{H}_1=0$, 且当误差$\|R_t\|\leq 1.0e-010$时, 认为矩阵对$(S_t, H_t)$是方程$(2.10)$和$(2.11)$的极小Frobenius范数对称解.数值结果见下表 3.1.

其中"It."表示迭代次数, $ RR(\Sigma_1, Y_1)=\|\hat{M}Y_1\Sigma_1-\hat{K}Y_1\|, $ $ RR(\Lambda_2, X_2)=\|\hat{M}X_{2}\Lambda_{2}-\hat{K}X_{2}\|, $ $\hat{M}=M_a+B\hat{G}, $ $ \hat{K}=K_a+ B\hat{F}.$有趣的是, 从表 3.1看到, 迭代步数并不一定随着矩阵阶数的提高而增加.

该文给出了一种利用加速度和位移反馈控制技术来修正有限元模型的迭代方法.该方法通过选取合适的初始迭代矩阵, 经有限步迭代即可求得问题2的最优逼近解$(\hat{G}, \hat{F})$, 且可以保证修正是不溢出的.最后, 两个数值例子验证了算法的有效性.