CN103530278A - 一种基于最小结构奇异子集的分块快速延展方法 - Google Patents

Abstract

本发明公开了一种基于最小结构奇异子集的分块快速延展方法，属于多物理场建模技术领域。本方法包括以下步骤：步骤一：将多物理场模型产生的稀疏高阶DAE系统转化为Σ矩阵；步骤二：对Σ矩阵进行块状上三角矩阵分解；步骤三：对每一独立的块状上三角矩阵采用分块快速延展方法获取它的最优偏移向量；步骤四：由每个矩阵局部最优偏移向量构成Σ矩阵的整体最优偏移向量，将DAE系统转化为ODE系统，并找出一致初始点值的约束方程。本方法通过将DAE系统转化为ODE系统，并找出一致初始点值的约束方程,使得快速准确地实现稀疏高指标DAE系统的指标约简，实现了直接获取DAE系统的最优偏移向量等结构信息。

Description

一种基于最小结构奇异子集的分块快速延展方法

技术领域

本发明属于多物理场建模技术领域，涉及一种基于最小结构奇异子集的分块快速延展方法。

背景技术

目前，针对一般结构的DAE系统主要的指标约简方法有：Gear方法，Pantelides方法和Pryce方法。

Gear方法是一种符号操作的DAE指标约简方法。该方法通过一系列符号操作找出DAE系统中需要微分的代数方程，反复执行这一过程，直至将DAE系统转化为微分方程(OrdinaryDifferential Equations,ODE)系统。但对于一般的非线性DAE系统，因微分了一些不必要的代数方程，Gear方法的时间复杂度很高，有时甚至无法实现。

Pantelides方法是一种基于MSS子集的指标约简方法。该方法首先将一阶DAE系统转化为对应的二部图，通过二部图的增广路算法来查找DAE系统的MSS子集；然后对MSS子集微分一次，并对DAE系统对应的二部图进行相应的删除和增加操作；反复执行直至DAE系统中找不到MSS子集。因为其相对比较低的时间复杂度，目前该方法被广泛应用于Dymola，MapleSim和MWorks等专业软件。但是针对高阶DAE系统，Pantelides方法不能直接处理，首先通过变量替换将其转化为一阶DAE系统，同时在一阶DAE系统指标约简中，又会不断引进新的变量来完成指标约简，从而不能直接获取DAE系统的结构信息。针对多物理场仿真模型容易产生的稀疏高阶DAE系统，虽然可能存在块状上或下三角结构，但该方法不易实现块状化处理。

Pryce方法是基于DAE系统Σ矩阵(Signature Matrix)的指标约简方法，从理论上等价于Pantelides方法，也被称为结构分析方法。该方法先将求解Σ矩阵最大横截面值的问题转化为指派问题，并通过求解这个指派问题来找到Σ矩阵值最大的横截面集合；然后通过求解指派问题的对偶问题(一般的线性规划问题)，得到方程的最优偏移值(Canonical Offset)和变量的最优偏移值；最后根据这些最优偏移值将DAE系统转化为ODE系统，并找出一致初始点值的约束方程。该方法能够更加直接地获取DAE系统结构信息(最优偏移值等)。针对存在块状上或下三角结构的稀疏高阶DAE系统，Pryce等人最近提出了分块处理的技术，并将他们的软件包DAESA植入到Matlab的工具箱中，虽然可以降低部分计算量，但他们有时只能获取系统的有效偏移值，不能获取最优偏移值。

发明内容

有鉴于此，本发明的目的在于提供一种基于最小结构奇异子集的分块快速延展方法，该方法用于解决在数学软件中DAE系统指标约简过程中的效率问题，在现有符号化简技术的基础上，快速准确地实现稀疏高指标DAE系统的指标约简，直接获取DAE系统的最优偏移向量等结构信息。

为达到上述目的，本发明提供如下技术方案：

一种基于最小结构奇异子集的分块快速延展方法，包括以下步骤：步骤一：将多物理场模型产生的稀疏高阶DAE系统转化为Σ矩阵；步骤二：对Σ矩阵进行块状上三角矩阵分解为：

其中

步骤三：对每一独立的块状上三角矩阵Σⁱ，采用分块快速延展方法获取它的最优偏移向量；步骤四：由每个Σⁱ矩阵局部最优偏移向量构成Σ矩阵的整体最优偏移向量，将DAE系统转化为ODE系统，并找出一致初始点值的约束方程。

进一步，步骤三中的分块快速延展方法具体包括以下步骤：

1）存储块状上三角结构的Σ矩阵(符号矩阵)和对角子方阵阶数：输入n阶p块上三角结构的Σ矩阵

存储其每一对角子方阵阶数到对应的p维向量S:＝[s₁,s₂,...,s_p]_1×p中，且满足

2）初始化设定，令i＝1，n维p块参数向量md＝[md₁,md₂,..md_p]_1×p，

方程偏移值向量c＝[c₁,c₂,..c_p]_1×p，

和变量偏移值向量d＝[d₁,d₂,..d_p]_1×p，

3）存储中间矩阵MS：取MS＝M_ii，即Σ矩阵M的第i的子方阵；

4）运用含参数Σ矩阵约简方法，获取偏移向量c_i和d_i，[c_i,d_i]＝PSMRM(M_ii,md_i)，再取md_i＝d_i；

5）判断主算法计算是否结束：如果i＝p，输出c和d，结束；否则，i＝i+1；

6）非对角矩阵[M_i-1,i,…,M_i-1,p]更新：只要第i-1个子方阵M_i-1,i-1的方程偏移值c_i-1中存在元素c_i-1[j]＞0，则[M_i-1,i,…,M_i-1,p]矩阵的第j行中所有的非负元素都增加c_i-1[j]；

7）参数向量[md_i,md_i+1,..,md_p]更新：取更新后[M_i-1,i,…M_i-1,p]矩阵中各列的最大值，转步骤3）。

进一步，所述含参数Σ矩阵约简方法具体包括以下步骤：

1）输入n(＞1)阶Σ矩阵M和需要满足限制的n维参数向量md；

2）初始化设定：k＝0，方程n维偏移向量c＝[0,0,...,0]_1×n，变量偏移向量d＝[0,0,...,0]_1×n，方程集合

变量集合

和变量匹配方程n维数向量Matching＝[0,0,...,0]_1×n；

3）变量偏移向量d的初始更新：先取d为M中各列的最大值构成的向量，再取d为d和参数md对应各元素中的最大值来更新d中对应的元素值，即d＝max(d,md)；

4）方程主导数集合L选取：调用函数L＝FindLeadingDerivativeSet(M,d)，其中函数表示为L＝{(i,j)|M_i,j＝d_j,i＝1...n,j＝1...n}；

5）判断子算法计算是否结束：如果i＝n，输出c和d，结束；否则，i＝i+1，取PATHFOUND=FALSE；

6）标记变量向量初始化设定：Label＝[0,0,...,0]_1×n；

7）匹配查找：运用FindMatching(i,PATHFOUND,Label,Matching,L)；

8）判断是否找到匹配：如果找到匹配，PATHFOUND=TRUE，转步骤（2.5）；否则，找到集合X＝{j|Label[j]＝1}和关于X的MSS子集F＝{i}∪{Matching[j]|j∈X}；

9）对MSS子集进行微分，并更新方程偏移向量c和调整M矩阵：如果F仅有一个元素i，取md与M矩阵的第i行之间向量差的最小值为ac，c[i]＝c[i]+ac,对M矩阵的第i行中所有非负元素都增加ac；否则依次取i′∈F，c[i′]＝c[i′]+1，并对M矩阵的第i′行中所有非负元素增加1；

10）变量偏移向量d和方程主导数集合L更新：md取M矩阵各列的最大值，再更新d为max(d,md)，并取L＝FindLeadingDerivativeSet(M,d)，转步骤6）。

本发明的有益效果在于：本发明采用分块快速延展方法,获取它的最优偏移向量,根据由每个矩阵局部最优偏移向量构成矩阵的整体最优偏移向量，将DAE系统转化为ODE系统，并找出一致初始点值的约束方程,使得快速准确地实现稀疏高指标DAE系统的指标约简，实现了直接获取DAE系统的最优偏移向量等结构信息。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明所述方法流程图；

图2为分块快速延展方法流程图；

图3为分块快速延展方法BFPM与其他方法的耗时比较图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

图1为本发明所述方法流程图，如图所示，该基于最小结构奇异子集的分块快速延展方法，包括以下步骤：步骤一：将多物理场模型产生的稀疏高阶DAE系统转化为Σ矩阵；步骤二：对Σ矩阵进行块状上三角矩阵分解为：其中

步骤三：对每一独立的块状上三角矩阵Σⁱ，采用分块快速延展方法获取它的最优偏移向量；步骤四：由每个Σⁱ矩阵局部最优偏移向量构成Σ矩阵的整体最优偏移向量，将DAE系统转化为ODE系统，并找出一致初始点值的约束方程。

分块快速延展方法BFPM的原理是：对于块状上三角结构的Σ矩阵，按对角块的顺序由上往下依次运用含参数Σ矩阵约简方法（Signature Matrix Reduce Method with Parameter,PSMRM），从而逐渐从局部最优偏移向量扩展到Σ矩阵的整体最优偏移向量。其特征是：在任意给定参数向量的条件下，通过不断地查找Σ矩阵中的MSS子集，并更新偏移向量；直到找不到MSS子集时，变量偏移向量中各元素的值均不小于给定参数向量中与之对应元素的值。这一重要特性保证了分块快速延展方法的成功实现。

图2为分块快速延展方法流程图，下面通过实施例来对本方法进行说明。

以国产的盛通数控自动数控冲床STSK-Q2508为实施对象，利用Maplesim6.1软件建立仿真模型，通过它的冗余处理得到加工中心由163个微分或代数方程构成的DAE系统，并将DAE系统转化为Σ矩阵；进一步对Σ矩阵进行块状上三角矩阵分解，并可以表示为：

$M=\left[\begin{array}{cc}{M}^1& \\ & M^2\end{array}\right],$ 其中 $M^1=\left[\begin{array}{cc}0& M_{12}^1\\ & M_{22}^1\end{array}\right]$ 和 $M^2=\left[\begin{array}{cccccccc}0& 0& 0& 2& 2& 2& X& X\\ 0& 0& 0& 2& 2& 2& X& X\\ & & 0& 2& 2& 2& X& X\\ & & & 0& 0& 0& X& X\\ & & & & 0& 0& X& X\\ & & & & 0& 0& X& X\\ & & & & & & M_{55}^2& X\\ & & & & & & & M_{66}^2\end{array}\right];$

上述矩阵的空白处代表其元素值为-∞；且

和

的矩阵阶数都为52；X表示其中存在非负元素。

对M¹和M²分别运用分块快速延展方法，其中含参数Σ矩阵约简方法会运用到FindMatching(i,PATHFOUND,Label,Matching,L)方法。它的关键是查找在方程主导数集合L中是否存在以方程i为起始点的增广路径；如果存在，则PATHFOUND为True，并通过对增广路径取反为其中每一方程匹配每一变量；否则可以找到MSS子集F＝{i}∪{Matching[j]|j∈X}(它是关于X的MSS子集)。当分块快速延展方法结束时，可以得到Σ矩阵M的方程整体最优偏移向量c和变量整体最优偏移向量d，计算该方法的时间复杂度，并与其他方法比较，见图3,M1和M2分别表示M¹和M²，分块快速延展方法和Pantelides方法的效率均高于Pryce方法的效率，对于存在两个高阶数子方阵的Σ矩阵M²，分块快速延展方法的效率大约是Pantelides方法的一倍；最后根据方程整体最优偏移向量c，将DAE系统快速准确地转化为ODE系统，并找出一致初始点值的约束方程。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其做出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (3)

1.一种基于最小结构奇异子集的分块快速延展方法，其特征在于：包括以下步骤：

步骤一：将多物理场模型产生的稀疏高阶DAE系统转化为Σ矩阵；

步骤二：对Σ矩阵进行块状上三角矩阵分解为：

其中

步骤三：对每一独立的块状上三角矩阵Σⁱ，采用分块快速延展方法获取它的最优偏移向量；

步骤四：由每个Σⁱ矩阵局部最优偏移向量构成Σ矩阵的整体最优偏移向量，将DAE系统转化为ODE系统，并找出一致初始点值的约束方程。

2.根据权利要求1所述的一种基于最小结构奇异子集的分块快速延展方法，其特征在于：步骤三中的分块快速延展方法具体包括以下步骤：

1）存储块状上三角结构的Σ矩阵(符号矩阵)和对角子方阵阶数：输入n阶p块上三角结构的Σ矩阵

存储其每一对角子方阵阶数到对应的p维向量S:＝[s₁,s₂,...,s_p]_1×p中，且满足

2）初始化设定，令i＝1，n维p块参数向量md＝[md₁,md₂,..md_p]_1×p，

方程偏移值向量c＝[c₁,c₂,..c_p]_1×p，

和变量偏移值向量d＝[d₁,d₂,..d_p]_1×p，

3）存储中间矩阵MS：取MS＝M_ii，即Σ矩阵M的第i的子方阵；

4）运用含参数Σ矩阵约简方法，获取偏移向量c_i和d_i，[c_i,d_i]＝PSMRM(M_ii,md_i)，再取md_i＝d_i；

5）判断主算法计算是否结束：如果i＝p，输出c和d，结束；否则，i＝i+1；

6）非对角矩阵[M_i-1,i,…,M_i-1,p]更新：只要第i-1个子方阵M_i-1,i-1的方程偏移值c_i-1中存在元素c_i-1[j]＞0，则[M_i-1,i,…,M_i-1,p]矩阵的第j行中所有的非负元素都增加c_i-1[j]；

7）参数向量[md_i,md_i+1,..,md_p]更新：取更新后[M_i-1,i,…M_i-1,p]矩阵中各列的最大值，转步骤3）。

3.根据权利要求2所述的一种基于最小结构奇异子集的分块快速延展方法，其特征在于：所述含参数Σ矩阵约简方法具体包括以下步骤：

1）输入n(＞1)阶Σ矩阵M和需要满足限制的n维参数向量md；

2）初始化设定：k＝0，方程n维偏移向量c＝[0,0,...,0]_1×n，变量偏移向量d＝[0,0,...,0]_1×n，方程集合变量集合

和变量匹配方程n维数向量Matching＝[0,0,...,0]_1×n；

3）变量偏移向量d的初始更新：先取d为M中各列的最大值构成的向量，再取d为d和参数md对应各元素中的最大值来更新d中对应的元素值，即d＝max(d,md)；

4）方程主导数集合L选取：调用函数L＝FindLeadingDerivativeSet(M,d)，其中函数表示为L＝{(i,j)|M_i,j＝d_j,i＝1...n,j＝1...n}；

5）判断子算法计算是否结束：如果i＝n，输出c和d，结束；否则，i＝i+1，取PATHFOUND=FALSE；

6）标记变量向量初始化设定：Label＝[0,0,...,0]_1×n；

7）匹配查找：运用FindMatching(i,PATHFOUND,Label,Matching,L)；

8）判断是否找到匹配：如果找到匹配，PATHFOUND=TRUE，转步骤（2.5）；否则，找到集合X＝{j|Label[j]＝1}和关于X的MSS子集F＝{i}∪{Matching[j]|j∈X}；

9）对MSS子集进行微分，并更新方程偏移向量c和调整M矩阵：如果F仅有一个元素i，取md与M矩阵的第i行之间向量差的最小值为ac，c[i]＝c[i]+ac,对M矩阵的第i行中所有非负元素都增加ac；否则依次取i′∈F，c[i′]＝c[i′]+1，并对M矩阵的第i′行中所有非负元素增加1；

10）变量偏移向量d和方程主导数集合L更新：md取M矩阵各列的最大值，再更新d为max(d,md)，并取L＝FindLeadingDerivativeSet(M,d)，转步骤6）。

Images