CN101315649A - 含大量输入端口的微机电系统降阶建模方法 - Google Patents

Abstract

本发明公开了含大量输入端口的微机电系统降阶建模方法，根据结构设计的参数，建立结构的三维实体模型，并使用有限元方法进行模态分析，从分析结果中提取结构的质量、刚度矩阵，建立结构的二阶动力学行为方程后转换为一阶状态方程；将输入矩阵B按列平均拆分为k块，得到拆分后结构的状态方程；对于第i个子系统应用分块Arnoldi算法进行降阶，得到第i个子系统的宏模型状态方程；将所有子系统的宏模型按输出叠加以后，就可得到原始系统的整体宏模型。本发明使基于宏模型的MEMS系统级建模与仿真速度加快，提高MEMS的设计效率。

Description

含大量输入端口的微机电系统降阶建模方法

技术领域

本发明涉及一种微机电系统(MEMS)的建模方法，属微机电系统设计与模型降阶领域。

背景技术

随着MEMS技术的迅速发展，MEMS器件结构也变得越来越复杂，采用其原始模型直接进行系统级建模与仿真十分耗时，效率低下。为了提高仿真效率，一般对原始模型进行处理，采用一定的降阶方法，在维持输入输出特性基本不变的情况下使其自由度大幅缩减，建立能反应原始模型特性的低阶模型，也叫宏模型或缩聚模型，自由度缩减过程称为降阶建模过程。

一阶线性时不变MEMS结构的状态方程可表示如下：

$\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=\mathrm{Bu}\left(t\right)\\ y\left(t\right)=L^{T}x\left(t\right)\end{array}\right.---\left(1\right)$

其中：t是时间变量，x(t)∈R^N是状态变量，初始状态变量x(0)＝x0，u(t)∈R^m是输入向量，y(t)∈R^p是输出向量，C，G∈R^N×N是系统矩阵，B∈R^N×m是输入矩阵，L∈R^N×p是输出矩阵。一般系统矩阵的维数N非常大，复杂系统有时甚至高达数十万，对如此大规模的状态方程，采用普通有限元方法求解是相当费时且效率低下的。

降阶建模即将系统的原始状态方程采用一定降阶算法进行处理，在保证足够精度的条件下，使状态方程的维数大幅缩减，得到规模远小于原始系统的低阶模型，也称宏模型；而降阶建模过程通常也称为宏建模。宏模型可用如下状态方程表示：

$\left\{\begin{array}{c}{C}_r{\stackrel{\·}{x}}_r\left(t\right)+G_r{x}_r\left(t\right)=B_{r}u\left(t\right)\\ y_r\left(t\right)={L_r}^T{x}_r\left(t\right)\end{array}\right.---\left(2\right)$

与方程(1)对应，t是时间变量，x_r(t)∈R^r是降阶后的状态变量，u(t)∈R^m是输入向量，y(t)∈R^p是输出向量，C，G∈R^r×r是降阶后的系统矩阵，B∈R^r×m是降阶后的输入矩阵，L∈R^r×p是降阶后的输出矩阵。一般情况下，由于宏模型系统矩阵的维数r＜＜N，使得其求解效率大幅提高，这对于有效地提高MEMS的设计效率具有重要意义。对于多输入多输出系统(MIMO)，即m＞1且p＞1，当前普遍采用的是分块Arnoldi降阶算法完成降阶建模，所得宏模型的规模与输入端口数量m及传递函数矩匹配数量r的乘积成正比，因此随着输入端口数量m的及对宏模型精度要求的提高，宏模型的规模将迅速增大，同时还将破坏原始模型系统矩阵的稀疏特性，使得宏模型的求解效率有时甚至比原始模型还低。

IBM公司Peter Feldmann针对含大量输入输出端口的电路系统提出了一种基于矩阵奇异值分解与传递函数矩匹配相结合的降阶方法——SVDMOR法，该方法在降阶过程中认为各电路输入端口之间并不是相互独立的，它们之间存在某种相关性，即输入矩阵B为低秩矩阵，这种假设在电路系统中符合现实，但是在MEMS结构中，输入矩阵B为列满秩矩阵，输入端口之间不存在相关性，采用SVDMOR方法显然无法得到理想结果。IMTEK的Peter Benner教授与Lihong Feng等针对含大量输入端口的系统提出了叠加Arnoldi算法，将输入矩阵B拆分为m个列向量，使原始MIMO系统变为m个单输入单输出(SISO)子系统，然后对这些SISO系统分别采用标准Arnoldi算法进行降阶，并将得到的m个宏模型进行叠加以提取原始系统的整体宏模型。这种降阶方法的缺点之一是降阶次数与输入矩阵B的列数相等，降阶频繁，效率低下；其次得到的系统整体宏模型是m个单输入单输出(SISO)子系统宏模型的叠加，其规模依旧庞大。

发明内容

为了克服现有技术求解效率低下的不足，本发明提供一种微机电系统降阶建模方法，能够实现含大量输入端口的MEMS结构的降阶建模。

本发明解决其技术问题所采用的技术方案包括以下步骤：

第一步：建立MEMS结构的状态方程。根据结构设计的参数，采用三维造型工具建立结构的三维实体模型，并使用有限元方法进行模态分析，从分析结果中提取结构的质量、刚度矩阵，建立结构的二阶动力学行为方程：

$\left\{\begin{array}{c}M\stackrel{\·\·}{q}\left(t\right)+D\stackrel{\·}{q}\left(t\right)+\mathrm{Kq}\left(t\right)=\mathrm{Pu}\left(t\right)\\ y\left(t\right)=E^{T}q\left(t\right)\end{array}\right.---\left(3\right)$

式中：t为时间，q(t)为节点位移向量，P为载荷矩阵，u(t)为输入向量，y(t)为输出向量，E为输出矩阵，M，D，K分别为结构的质量、阻尼和刚度矩阵，MEMS结构一般采用Rayleigh结构阻尼模型，即D＝αM+βK，其中α、β是和结构模态阻尼比相关的参数，具体取值可参考结构动力学文献。

二阶动力学行为方程经矩阵等价变换后，可转换为如式(1)所示的一阶状态方程，

$\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=\mathrm{Bu}\left(t\right)\\ y\left(t\right)=L^{T}x\left(t\right)\end{array}\right.$

其中： $x\left(t\right)=\left[\begin{array}{c}q\left(t\right)\\ \stackrel{\·}{q}\left(t\right)\end{array}\right],$ $C=\left[\begin{array}{cc}D& M\\ W& 0\end{array}\right],$ $G=\left[\begin{array}{cc}K& 0\\ 0& -W\end{array}\right],$ $B=\left[\begin{array}{c}P\\ 0\end{array}\right],$ $L=\left[\begin{array}{c}E\\ 0\end{array}\right].$

W是一个N×N的非奇异矩阵，实际应用中通常假设W为单位矩阵，即W＝I，在结构力学中，由于M、D、K均为对称正定矩阵，也可以取W＝M，这样上述状态方程的系统矩阵C、G保持了原二阶状态方程的系数矩阵的对称特性。由上式可以看出，二阶动力学行为方程转换为一阶状态方程后，方程的系数矩阵的规模翻倍。

第二步：将输入矩阵B按列分块。输入矩阵B的列通常对应着结构的输入端口，根据结构特点将B按列平均拆分为k块，如果m/k为整数，则每一块B的列数均为m/k，否则最后一块的列数＜m/k的整数部分。拆分后结构的状态方程可表示为：

$\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{B}_i{u}_i\left(t\right)\\ y\left(t\right)=L^{T}x\left(t\right)\end{array}\right.---\left(4\right)$

式(4)表明原始状态方程被拆分为k个子系统描述的状态方程。其中B_i为拆分后第i个子系统的输入矩阵，u_i(t)为相应第i个输入向量。

第三步：采用分块Arnoldi算法对所有子系统分别降阶。对于第i个子系统，以A＝-G^-1C、R＝G^-1B_i作为输入，应用分块Arnoldi算法进行降阶。则第i个子系统的宏模型状态方程可描述为：

$\left\{\begin{array}{c}{C}_{\mathrm{ir}}{\stackrel{\·}{x}}_{\mathrm{ir}}\left(t\right)+G_{\mathrm{ir}}{x}_{\mathrm{ir}}\left(t\right)=B_{\mathrm{ir}}{u}_i\left(t\right)\\ y_i\left(t\right)={L_{\mathrm{ir}}}^T{x}_{\mathrm{ir}}\left(t\right)\end{array}\right.,i\∈(1,k)---\left(5\right)$

其中： $C_{\mathrm{ir}}=V_i^T{\mathrm{CV}}_i,$ $C_{\mathrm{ir}}=V_i^T{\mathrm{GV}}_i,$ $B_{\mathrm{ir}}=V_i^{T}B,$ $L_{\mathrm{ir}}=V_i^{T}L,$ t是时间变量，x_ir(t)是降阶后的状态变量，u_i(t)是输入向量，y_i(t)是输出向量，C_ir、G_ir为系统矩阵，B_ir是降阶后的输入矩阵，L_ir是降阶后的输出矩阵，r为传递函数匹配矩的数量。

设子系统输入矩阵B_i的列数为c，则宏模型的系统矩阵维数w＝r×c。且一般远小于原始系统矩阵维数，因而所有子系统模型得到降阶。

第四步：将所有子系统的宏模型按输出进行叠加，提取系统的整体宏模型。

叠加属性是线性系统的固有属性，其原理描述如下：如果存在映射关系u₁(t)→y₁(t)，且u₂(t)→y₂(t)，则映射关系αu₁(t)+βu₂(t)→αy₁(t)+βy₂(t)成立，其中u_i(t)(i∈1，2...)是输入信号，y_i(t)(i∈1，2...)是相应的输出响应，α、β是任意两个常数。

对于式(5)所表示子系统的宏模型，根据叠加属性，可得原始系统整体宏模型的输出为：

$y\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{y}_i\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{L_{\mathrm{ir}}}^T{x}_{\mathrm{ir}}\left(t\right)---\left(6\right)$

其物理意义描述如下：对于一个含有大量输入端口的MEMS结构，对应于状态方程中某一块的任意某些端口的输入都会对所有端口的输出产生贡献，而每个端口的输出为所有端口输入贡献的代数和。

将所有子系统的宏模型按输出叠加以后，就可得到原始系统的整体宏模型，其状态方程为：

$\left\{\begin{array}{c}{C}_{\mathrm{ir}}{\stackrel{\·}{x}}_{\mathrm{ir}}\left(i\right)+G_{\mathrm{ir}}{x}_{\mathrm{ir}}\left(t\right)=B_{\mathrm{ir}}{u}_i\left(t\right)i\∈(1,k)\\ y\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{y}_i\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{L_{\mathrm{ir}}}^T{x}_{\mathrm{ir}}\left(t\right)\end{array}\right.---\left(7\right)$

将宏模型的状态方程采用硬件描述语言进行编码，即可将其用于MEMS的系统级快速建模与仿真。

本发明的有益效果是：由于采用分块叠加Arnoldi算法提取含有大量输入端口MEMS结构的整体宏模型，使原始模型的自由度规模得到大幅缩减，从而使基于宏模型的MEMS系统级建模与仿真速度加快，提高MEMS的设计效率。

下面结合附图和实施例对本发明进一步说明。

附图说明

图1是本发明提出的分块叠加Arnoldi降阶建模流程图。

图2是一种变截面折叠梁示意图；

图中，1-变截面折叠梁的粗梁，2-变截面折叠梁的细梁，3-变截面折叠梁的纵梁。

具体实施方式

下面结合一种变截面折叠梁的宏模型提取过程对本发明进行进一步说明。变截面折叠梁是MEMS设计中的一种典型结构，广泛应用于微加速度计、微陀螺、微型光栅等MEMS器件的设计中。参照附图2，其结构由五根粗梁1、二十二根细梁2和十二根纵梁3构成。变截面折叠梁的宏模型提取过程包括如下步骤：

第一步：使用三维造型工具建立变截面折叠梁的几何模型，采用有限元软件ANSYS，选择beam188材料和Solid92单元进行模态分析，其有限元模型含有两个约束面，7718个节点，其中包含两个输入输出节点，由于beam188的每个节点包含6个自由度，即分别为沿x、y、z坐标轴的平动和绕坐标轴的转动，因此变截面折叠梁共包含12个输入端口，属于本发明研究的包含大量输入端口的范畴之内，另外，模型的质量刚度矩阵的规模为7718×6＝46308阶。提取有限元模型的质量、刚度矩阵，并通过矩阵变换得到模型的一阶状态方程：

$\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=\mathrm{Bu}\left(t\right)\\ y\left(t\right)=L^{T}x\left(t\right)\end{array}\right.$

由于二阶动力学行为方程转换为一阶状态方程后，系数矩阵的规模加倍，因此，C，G∈R^{(46038×2)×(46038×2)}，B∈R^{(46038×2)×12}，L∈R^{(46038×2)×12}，u(t)∈R¹²，y(t)∈R¹²。

第二步：根据变截面折叠梁的输出特点，可将输入矩阵B按节点自由度划分，即将输入矩阵的划分为两个6列的子矩阵，得到结构的状态方程为：

$\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=B_1{u}_1\left(t\right)+B_2{u}_2\left(t\right)\\ y\left(t\right)=L^{T}x\left(t\right)\end{array}\right.$

其中：B₁，B₂∈R^{(46038×2)×6}，u₁(t)，u₂(t)∈R⁶。该状态方程结构的原始状态方程划分为两个子系统。

第三步：对步骤二所提取的变截面折叠梁状态方程，以A＝-G^-1C、R＝G^-1B₁及A＝-G^-1C、R＝G^-1B₂作为输入，分别应用分块Arnoldi算法对两个子系统进行降阶，得到子系统的宏模型状态方程为：

$\left\{\begin{array}{c}{C}_{\mathrm{ir}}{\stackrel{\·}{x}}_{\mathrm{ir}}\left(t\right)+G_{\mathrm{ir}}{x}_{\mathrm{ir}}\left(t\right)=B_{\mathrm{ir}}{u}_i\left(t\right)\\ y_i\left(t\right)={L_{\mathrm{ir}}}^T{x}_{\mathrm{ir}}\left(t\right)\end{array}\right.,i\∈(1,2)$

设矩匹配数量为r＝2，则子系统宏模型系数矩阵维数均为2×6＝12。即：C_ir，G_ir∈R^12×12，B_ir，L_ir∈R^12×6。

第四步：将步骤三所获得的子系统宏模型按输出进行叠加，得到原始系统的整体宏模型，其状态方程描述如下：

$\left\{\begin{array}{c}{C}_{\mathrm{ir}}{\stackrel{\·}{x}}_{\mathrm{ir}}\left(t\right)+G_{\mathrm{ir}}{x}_{\mathrm{ir}}\left(t\right)=B_{\mathrm{ir}}{u}_i\left(t\right)i\∈(1,2)\\ y\left(t\right)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{y}_i\left(t\right)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{L_{\mathrm{ir}}}^T{x}_{\mathrm{ir}}\left(t\right)\end{array}\right.$

将状态方程采用MAST硬件描述语言进行编码，即可得到含有两个输入输出节点共12个端口的变截面折叠梁宏模型。

Claims (1)

1、含大量输入端口的微机电系统降阶建模方法，其特征在于包括下述步骤：

(a)根据结构设计的参数，建立结构的三维实体模型，并使用有限元方法进行模态分析，从分析结果中提取结构的质量、刚度矩阵，建立结构的二阶动力学行为方程 $\left\{\begin{array}{c}M\stackrel{\·\·}{q}\left(t\right)+D\stackrel{\·}{q}\left(t\right)+\mathrm{Kq}\left(t\right)=\mathrm{Pu}\left(t\right)\\ y\left(t\right)=E^{T}q\left(t\right)\end{array}\right.,$ 式中t为时间，q(t)为节点位移向量，P为载荷矩阵，u(t)为输入向量，y(t)为输出向量，E为输出矩阵，M，D，K分别为结构的质量、阻尼和刚度矩阵；二阶动力学行为方程经矩阵等价变换后，转换为一阶状态方程 $\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=\mathrm{Bu}\left(t\right)\\ y\left(t\right)=L^{T}x\left(t\right)\end{array}\right.,$ 其中 $x\left(t\right)=\left[\begin{array}{c}q\left(t\right)\\ \stackrel{\·}{q}\left(t\right)\end{array}\right],$ $C=\left[\begin{array}{cc}D& M\\ W& 0\end{array}\right],$ $G=\left[\begin{array}{cc}K& 0\\ 0& -W\end{array}\right],$ $B=\left[\begin{array}{c}P\\ 0\end{array}\right],$ $L=\left[\begin{array}{c}E\\ 0\end{array}\right];$

(b)将输入矩阵B按列平均拆分为k块，如果m/k为整数，则每一块B_i的列数均为m/k，否则最后一块的列数＜m/k的整数部分；拆分后结构的状态方程表示为 $\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{B}_i{u}_i\left(t\right)\\ y\left(t\right)=L^{T}x\left(t\right)\end{array}\right.,$ 其中B_i为拆分后第i个子系统的输入矩阵，u_i(t)为相应第i个输入向量；

(c)对于第i个子系统，以A＝-G^-1C、R＝G^-1B_i作为输入，应用分块Arnoldi算法进行降阶，则第i个子系统的宏模型状态方程为 $\left\{\begin{array}{c}{C}_{\mathrm{ir}}{\stackrel{\·}{x}}_{\mathrm{ir}}\left(t\right)+G_{\mathrm{ir}}{x}_{\mathrm{ir}}\left(t\right)=B_{\mathrm{ir}}{u}_i\left(t\right)\\ y_i\left(t\right)={L_{\mathrm{ir}}}^T{x}_{\mathrm{ir}}\left(t\right)\end{array}\right.,$ 其中i∈(1，k)， $C_{\mathrm{ir}}=V_i^T{\mathrm{CV}}_i,$ $G_{\mathrm{ir}}=V_i^T{\mathrm{GV}}_i,$ $B_{\mathrm{ir}}=V_i^{T}B,$ $L_{\mathrm{ir}}=V_i^{T}L,$ r为传递函数匹配矩的数量；

(d)将所有子系统的宏模型按输出叠加以后，就可得到原始系统的整体宏模型，其状态方程为 $\left\{\begin{array}{c}{C}_{\mathrm{ir}}{\stackrel{\·}{x}}_{\mathrm{ir}}\left(t\right)+G_{\mathrm{ir}}{x}_{\mathrm{ir}}\left(t\right)=B_{\mathrm{ir}}{u}_i\left(t\right)i\∈(1,k)\\ y\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{y}_i\left(t\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{L_{\mathrm{ir}}}^T{x}_{\mathrm{ir}}\left(t\right)\end{array}\right..$

Images