基于稀疏矩阵技术的拉格朗日乘子快速求解.docxVIP

基于稀疏矩阵技术的拉格朗日乘子快速求解 1 基于稀疏矩阵的算例 多体系统动力学广泛应用于航空航天、机械、机器人等领域。近年来, 为了满足多体系统动力学仿真的需要, 发展了多种数值仿真算法 以上算法都是对传统增广法进行改进, 提高了算法的速度和精度, 使之适合快速或实时仿真的需要。但同时, 这些改进算法也增加了算法的复杂性, 给编程实现带来了难度, 且一些并行算法只能在并行计算机上运行, 这限制了算法的应用范围。针对这些问题, 本文以传统增广法为基础, 从提高算法计算效率入手, 给出了一种基于稀疏矩阵计算技术的改进算法。算法修改了拉格朗日乘子的求解方程, 并利用系数矩阵的稀疏性进行排序和符号分解, 获得一个最佳的计算数据结构来求解拉格朗日乘子, 将微分代数方程转化为易求解的常微分方程, 提高了计算效率。针对违约现象, 建立了违约修正方程, 再次利用获得的数据结构求解修正值, 直接修正每一步计算值, 使计算结果具有较高的精度。 1 超定方程的产生 完整约束的多体系统动力学基本方程可写成如下形式: 式中:M (q, t) ∈? 系统除了满足约束方程 (1b) 以外, 还应满足以下派生约束方程: 方程组 (1) 为一组超定方程。为求解该组方程, 直接增广法常有如下计算过程: 由方程 (1a) 可得: 将式 (2) 代入方程 (1d) , 整理得: 其中:A=φ 求解方程 (3) 可获得拉格朗日乘子, 于是将DAE方程转化为易求解的ODE方程: 对ODE方程求解方法的研究已经非常成熟, 有许多优良的算法 这一计算过程常常因为计算量过大而限制了应用。 2 采用稀疏矩阵技术的改进和扩大算法 2.1 模型求解过程 在采用传统增广法作为仿真算法时, 重复的拉格朗日乘子求解将占用大量的计算时间。如果能提高其计算效率, 则可以大量减少仿真过程的计算量。传统算法是利用线性方程组 (3) 来求解拉格朗日乘子, 该方程系数矩阵A的稀疏程度较低, 甚至对于一些机械系统根本不存在稀疏性 考察线性方程组: 此方程组可由方程 (1a) 和 (1d) 联立并去掉广义力矢量得到。求解该方程组也可以获得λ。与方程 (3) 相比, 方程 (5) 的优点在于系数矩阵的稀疏程度高, 且稀疏程度不受机械系统拓扑结构的影响。因此, 我们可以用稀疏矩阵计算技术来修改求解过程。为了详细说明系数矩阵结构, 将广义质量矩阵M中质量矩阵M′和惯量矩阵J′分开表达, 并用B表示方程 (5) 的系数矩阵, 于是: 利用方程 (5) 求解拉格朗日乘子的过程可以分为四个步骤: (1) 对系数矩阵B进行排序, 确定最佳矩阵结构; (2) 对排序后的矩阵进行LDL (3) 数值LDL (4) 回代求解。 对于约束不变的机械系统来讲, 尽管在整个仿真过程中, 雅克比矩阵 对稀疏矩阵排序等同于对其相应的无向图进行排序, 交换矩阵的行和列等价于重新标号无向图的顶点 符号分解是依据排序后得到的顺序进行的一次虚拟分解。这次分解并不进行具体元素的数值计算, 只是记录分解时填充元的位置, 同时也将分解时消去元素的位置依据修改顺序记录下来, 建立了一组针对LDL 具体的仿真预处理过程可以概括如下: (1) 首先给出系数矩阵B对应的无向图G= (V, E) ; (2) 对无向图G采用改进RCM排序算法进行排序, 并记录排序后的矩阵; (3) 根据获得的最佳分解顺序, 取出当前要进行操作的主行 (主列) ; (4) 找出将要被消去的行 (列) , 利用主行 (主列) 的相应元素对这些行 (列) 进行消去, 确定修改元素的位置并依次记录下来; (5) 重复步骤 (3) 和 (4) , 直至分解结束。 下面举例说明仿真预处理过程。如图2所示为一太阳展开帆板的结构简图, 帆板简化为一个四铰链和五杆机构, 建立机构的约束雅克比矩阵如式 (6) 所示。 根据约束雅克比矩阵, 建立方程的系数矩阵B如下。 图3为矩阵B对应的无向图。利用上述的改进RCM排序算法对无向图进行排序, 即对矩阵B进行排序, 获得排序后的矩阵B′。 然后, 将排序后的矩阵进行一次虚拟分解, 得到填充元和分解时消去元素的位置并依次存储, 存储的数据结构将在整个仿真过程中都可以利用。 在上述求解拉格朗日乘子的两个方程中, 尽管方程组 (5) 的数量大于方程组 (3) (约为两倍) , 但方程组 (5) 的系数矩阵稀疏程度不受系统拓扑结构的影响, 且稀疏程度高, 可采用稀疏矩阵技术提高计算效率。在计算量上, 方程 (5) 经过仿真预处理后, 只有数值求解在仿真循环内, 求解拉格朗日乘子约需O (2m) 的计算量, 而方程 (3) 的整个求解过程均在仿真循环内, 求解拉格朗日乘子约需O (m 2.2 仿真效果仿真 违约现象是直接增广法固有的缺陷,