杨流松 姚文莉 薛世峰,†

1. 中国石油大学(华东)储运与建筑工程学院, 青岛 266580; 2. 青岛理工大学理学院, 青岛 266520;

多体系统是由多个刚体或柔体通过某种形式联结的复杂机械系统。考虑到系统中存在的运动约束, 一般会将多体系统动力学的数学模型表示为微分‒代数混合方程组(Differential-Algebraic Equations, DAEs)。由于不方便得到 DAEs 的解析解, 对多体系统进行动力学仿真时, 需通过一些数值分析方法得到其数值解, 其中应用最广的是拉格朗日乘子法[1]。假如一个多体系统受m个完整约束C作用, 广义坐标用q表示, 微分‒代数混合方程组描述的动力学模型[2]可简写为

其中,M(q,t)∊Rn×n是系统的广义质量矩阵,q∊Rn×1代表系统的广义坐标列阵,C q(q,t) ∊Rm×n为系统的雅可比矩阵,λ∊Rm×1代表拉格朗日乘子,Q(q,q.,t)∊Rn×1为系统的广义力列阵(包含外力和速度相关项)。对式(2)求导, 可以得到系统的速度和加速度的约束方程:

其中,

式(1)～(4)构成封闭的微分–代数混合方程组。

一些多体系统常常由于自身的几何构型的原因, 在运动过程中出现奇异位置, 系统的自由度发生突变, 导致运动出现分叉点[3]。利用微分‒代数混合方程组对其进行动力学仿真时, 在奇异位置处,约束方程(式(2))的雅可比矩阵的秩rank(Cq)﹤m, 给计算造成很大的困难, 导致仿真结果与实际情况不符, 甚至发生严重的失真。引起这些问题的直接原因是系统经过奇异位置时约束反力的突变, 根本原因在于奇异位置处约束方程允许运动的子空间扩展[4]。当系统离开奇异位置时, 扩展的子空间的运动不再满足约束方程, 致使系统的运动偏离实际运动。为了克服约束雅克比矩阵缺秩引起的约束反力突变, 一些数值方法会引入较大的冲击载荷来消除不在约束空间内的速度分量[5–6]。如果数值方法不够稳定, 过大的冲击载荷也可能导致仿真失败。Kurdila 等[7]给出一种修正的 Lagrange 方程, 适用于具有奇异位置的多体系统, 但由于存在较强的刚性,不具有普适性。Bayo 等[8]提出一种增广拉格朗日方法, 通过选择较大的惩罚因子, 对约束方程的违约进行限制, 最后利用迭代算法来求解动力学方程组。还有一些学者将增广拉格朗日方法应用到非完整多体系统中, 用于提高计算机模拟的实时仿真效率[9–10]。增广拉格朗日方程已成功地应用于机械系统的研究和仿真, 如重型机械模拟器[11]、生物力学[12]和车辆动力学的联合模拟设置[13]等。Aghili 等[11]和Phong 等[12]利用雅可比矩阵的零空间, 将微分‒代数方程组转换为常微分方程组。不论系统是否存在奇异位置或冗余约束, 利用零空间方法得到的微分方程组形式紧凑, 其系数矩阵始终保持正定状态,方程组具有唯一解。上述方法都可以处理冗余约束和具有奇异位置的动力学问题, 但都有局限性, 比如, 基于增广拉格朗日方程方法依赖于惩罚参数的选取, 零空间方法中广义逆的数值计算比较耗时,影响仿真效率。

近年来, 关于高斯原理的研究呈增长趋势。姚文莉等[13–17]将极值形式的高斯原理拓展到多刚体系统的单边接触动力学问题。Udwadia 等[18]对受约束的多体系统运动提出新的解释, 并基于高斯原理和广义逆矩阵推导出显示形式的动力学方程。刘延柱[19–20]将高斯原理应用到柔体和刚柔耦合系统, 推导高斯拘束函数的表示方法, 并探讨利用最小值模型解决多体系统动力学的数值方法。杨流松等[17,21]利用高斯原理, 对受冗余约束作用的多体系统进行建模, 将力学中的动力学问题转换为一个求函数极值的数学问题, 并采用智能优化方法进行求解, 可以大大地提高仿真的计算效率。

上述关于高斯原理的研究没有涉及多体系统的奇异位置问题。本文利用高斯原理, 建立具有奇异位置的多体系统的动力学优化模型, 直接采用数学优化方法求解, 将具有奇异位置的多体系统动力学问题转化为求函数极值的优化问题。目前, 优化问题的求解手段主要分为两类: 第一类是传统算法,对于凸问题求解效率高, 但不适用于多极值问题,所求解多为局部最优解; 另一类是智能算法, 此类算法具有一定的随机性, 适用于多峰优化问题, 鲁棒性强, 不易陷入局部最优。粒子群优化(particle swarm optimization, PSO)是一种基于种群的随机优化算法, 由于易于实施, 并且可以快速地收敛到最优解, 近年来受到广泛关注, 并应用于人工神经网络训练、功能优化、模糊控制和模式识别等领域。本文对 PSO 进行改进, 并将智能优化与传统优化方法相结合, 充分利用传统优化的快速收敛和智能优化优化的全局搜索特性, 求解高斯方法建立的动力学优化模型, 探讨其应用于求解动力学奇异位置问题的可行性。

1 建模方法

1.1 增广拉格朗日方法(ALF)

Bayo 等[5]根据多体系统动力学的哈密顿描述,提出增广拉格朗日方法。该方法具有普遍性, 并且可以应用到完整系统和非完整系统中。与传统的动力学模型不同, 增广拉格朗日方法通过增加一部分惩罚项, 使动力学模型转化为一组常微分方程, 表达形式如下:

其中,λ是拉格朗日乘子,κ是由惩罚因子组成一个对角方阵,Ω和μ是常系数对角方阵。这里的系数Ω和μ与约束稳定理论中的常系数有相同的作用, 但在增广拉格朗日方法中, 这两个参数有特别的物理意义, 分别代表与约束方程相关的自然频率和阻尼系数。

从式(5)可以看出, 与传统的拉格朗日方法不同, 增广拉格朗日方法不需要附加约束方程组, 它允许约束直接插入动力学方程中, 将最初的 DAEs转化成常微分方程组。实际上, 增广拉格朗日方法是将惩罚的约束方程合并到动力学方程中, 从而使约束条件得到满足。为了求解式(5), 需要先计算出乘子λ。

将式(3)和(4)代入式(5), 可以得到

比较式(1)和(5), 可以得到

根据式(7), 可以建立一个迭代程序, 计算未知乘子λ,公式如下:

其中,i代表迭代次数,λ0=0。

将式(8)左乘AT带入式(6), 可以得到关于广义加速度..q的迭代公式:

1.2 零空间方法(NS)

为了克服动力学系统中的奇异位置问题, Phong等[12]利用约束方程的雅克比矩阵的零空间, 推导出一种形式紧凑的动力学方程——零空间方法。假设H满足以下条件:

则H就是雅克比矩阵Cq的零空间。根据秩–零化度定理[22], 下列关系式始终成立:

对式(10)两边同时进行转置运算, 可得

用HT左乘式(1), 并利用式(11), 可得

合并式(4)和(12), 则微分代数混合方程可以用矩阵表示为如下零空间方程:

其中,

事实上, 不论系统受到冗余约束的作用或经过某些奇异位置引起自由度的突变, 由于 rank(HTM)+ rank(Cq)=n始终成立, 所以式(13)的方程数目始终与系统的广义坐标的数目保持相等。所以, 不论雅克比矩阵Cq是否缺秩, 式(13)始终是一个正定方程组, 系统加速度具有唯一解, 并且可以直接求出。此外, 虽然零空间可以简化系统的动力学方程,但同时也增加计算的负担。利用约束惯性矩阵的零空间, 系统的加速度可以显式地表达为

1.3 高斯优化方法

1.3.1 基于高斯原理建立的动力学优化模型

高斯原理是基本的变分原理。与其他积分变分原理相比, 作为一种以加速度为变量的微分变分原理, 高斯原理可以更方便地将多体系统的动力学问题转化为求有约束条件函数的最小值问题, 从而借助于数学的优化算法来处理[23]。

高斯原理可以描述为最小值形式: 在任意时刻,系统的真实运动与位置和速度相同, 但加速度不同的可能运动相比较, 其拘束函数G取极小值, 因此高斯原理也称为高斯最小拘束原理。其中, 拘束函数是系统真实运动偏离自由运动的度量。

假设系统由N个质点组成,mi与ir..分别是第i个质心的质量和加速度, 则拘束函数表达如下:

针对上述质点系, 选取n个坐标q1,q2,...,qn作为系统的广义坐标, 所选取的广义坐标彼此之间可以不独立。质点的矢径可用广义坐标表达为

对式(16)求导, 可得

将式(15)表示的拘束展开, 得到

将式(18)代入式(19)可得

令

引入下列矩阵

将式(20)用上述矩阵表示, 可以得到广义坐标形式的拘束函数:

则系统的广义加速度可以通过求解以下约束优化问题而得到:

1.3.2 奇异位置问题的求解难点

求解优化模型(式(22))常用的约束优化方法有很多, 包括传统的优化和智能优化算法。传统的优化中常用的有罚函数法、广义拉格朗日乘子法以及可行方向法。如果qC满秩, 则利用传统的拉格朗日乘子法得到式(22)的解是具有唯一性的。但是, 当系统存在奇异位置时, 在奇异位置附近约束方程(式(2))的雅克比矩阵是奇异的, 此时引用拉格朗日函数的约束最优性条件(K-T 条件)可以写为如下形式:

其中λ是引入的拉格朗日乘数。式(22)的全局最优解包含于式(23)的解空间。由于qC行不满秩, 求解式(23)可以得到无数个解。所以, 用传统的优化算法求解式(22)时, 容易陷入局部最优, 导致所求系统动力学的加速度出现分叉, 仿真与实际出现偏差,甚至失败。

因此, 如果采用传统的优化算法求解奇异约束条件下的优化问题, 难点在于如何在保证收敛效率的前提下跳出局部最优, 得到全局最优解。近年来,受人类智能、生物群体社会性和自热现象规律启发, 人们发明了很多智能算法, 用于解决复杂的高维数和多极值优化问题。智能算法不依赖于初值,算法独立于求解域, 具有全局优化的性能, 通用性强, 适用于并行处理。但是, 智能算法具有一定的随机性, 所以也存在一些缺陷, 比如收敛速度慢、局部搜索能力差以及控制的参数较多等。

本文充分利用传统优化的快速局部寻优能力和智能优化算法的全局搜索性能, 将二者相结合来求解多体系统的动力学优化模型, 克服因动力学问题的奇异性导致的计算精度和效率问题。本文将这种采用高斯原理进行建模, 并采用优化算法的动力学求解方法称为高斯优化方法(Gauss optimization method, GS)。

2 计算方法

2.1 粒子群算法简介

为了解决约束优化问题, 有学者开发了不同的确定算法和随机算法[24]。由于可行方向和广义梯度下降等确定性方法对目标函数的连续性和可微性给出较强的假设, 而在实际问题中这些条件很难全部满足, 所以它们的适用性有限。随着优化理论的发展, 一些智能算法得到迅速发展和广泛应用, 为非线性、多极值的复杂函数及组合优化问题提供了切实可行的解决方案, 如遗传算法、进化策略、进化规划和粒子群优化等。

粒子群优化算法作为一种启发式全局优化技术[25–26], 可以模拟鸟群随机搜寻食物的捕食行为。鸟群通过自身经验和种群之间的交流, 调整自己的搜寻路径, 从而找到食物。粒子群优化每次搜寻时都会根据自身经验(自身历史搜寻的最优位置)和种群交流(种群历史搜寻的最优位置), 调整自身搜寻方向和速度, 从而找到最优解。

假设在一个D维的目标搜索空间中, 有τ个颗粒组成一个群落, 其中第i个粒子的位置表示为一个D维的向量Xi＝(X i1,Xi2,...,XiD); 第i个颗粒的历史最优位置为Pi＝(Pi1,Pi2,...,PiD); 整个颗粒群迄今为止搜索到的最好位置记为Pg＝(Pg1,Pg2,...,PgD); 第i个颗粒的运动速度也是一个D维的向量Vi＝(Vi1,Vi2,...,ViD)。对于颗粒i, 第k+1 次迭代时的颗粒速度和位置可以表示为

其中,w是惯性权重,c1和c2是两个随机数, rand()产生一个(0, 1)之间的随机数,i表示第i个颗粒,k代表迭代的次数,d代表维度。本文中, 这些参数分别取为

2.2 算法流程

本文对经典的 PSO 进行改进, 在文献[25]的基础上, 增加一个传统无约束优化的可行解和上一时间步的最优解作为额外的颗粒, 添加到初始颗粒群中, 以便达到快速收敛和全局寻优的目的。针对奇异位置引起的求解空间的扩张, 本文将变异算子引入粒子群优化算法中, 进一步增加种群的多样性,使算法尽快跳出局部最优解, 防止算法过早收敛。新的粒子群算法如下:

其中ρ为变异因子, 取值为2。

对具有奇异位置的多体系统进行仿真的操作流程如下。

1) 预先给定仿真时间tend和时间步长Δt, 给出满足约束方程的初始条件: 速度和位移q0, 令t＝t0。

2) 计算系统的质量矩阵M, 包括外力和速度相关项的广义力矩阵Q, 约束方程的雅克比矩阵Cq以及约束方程的右端项η。

3) 用式(22)建立多体系统的动力学模型, 并按照下列步骤, 求解得到系统此时的加速度。

① 令k=0, 初始化满足约束条件的d维颗粒Xk,令所有颗粒的初始速度Vk＝0。

② 求得无约束目标函数的最优解x0, 将最后一个颗粒更新为Xτ＝x0; 如果k=0, 则Xτ-1＝x0, 否则,Xτ-1＝qt-1。

③ 将当前颗粒的最优位置赋值给Pk:

④ 评价每个颗粒的目标函数值f k(Xk), 找到最优值, 并赋值给:

⑤ 令k=k+1, 利用式(24)～(27)更新颗粒速度和位置。

⑥ 更新每个颗粒的历史最优位置与群体最优颗粒位置:

4) 对加速度进行积分, 得到下一时刻的速度+Δt和位移qt+Δt, 并令t＝t+Δt。

5) 若t﹥tend, 则程序终止, 否则, 转到步骤2。

2.3 3 种方法对比

尽管本文提到的3 种方法所建立的动力学方程的形式及求解过程存在差异, 但它们在求解非奇异问题方面都具有可行性。虽然零空间方法和增广拉格朗日方法在求解多体系统中的奇异问题方面可行, 但其求解速度不能满足实际仿真的需要[8,12]。比如, 虽然零空间方法得到的动力学方程是一个正定方程组, 但用奇异值分解方法计算约束雅克比矩阵的零空间非常耗时。增广拉格朗日方法将一个动力学问题转化成一个非线性问题, 并通过循环迭代来求解, 当系统存有奇异位置时, 迭代不一定收敛,并且十分依赖参数的选择。

与增广拉格朗日方法和零空间方法相比, 在求解具有奇异位置的动力学问题时, 高斯优化方法具有以下独特的优势。

1) 不需要进行复杂的奇异值分解计算, 可以大大地提高求解效率。

2) 不需要选择合适的惩罚因子, 更具有普适性。

3) 将动力学问题转变为优化问题, 并通过数学优化理论进行求解, 扩展了求解方式。

4) 引入粒子群算法进行数值求解, 并借鉴遗传算法的思想, 引入变异因子来扩展求解空间, 在初始种群中添加与时间相关的颗粒, 丰富了种群的多样性。

3 数值算例

如图 1 所示, 曲柄滑块机构由两个质量和长度相同的均质杆件和一个滑块构成, 在重力作用下,在铅锤平面内运动, 重力加速度向下。假设x轴水平,y轴垂直向上, 滑块只能沿x轴滑动, 笛卡尔坐标为(x, 0)。两个杆的质量均为m1=6 kg, 长度均为l=1 m, 滑块的质量为m2=2 kg。1θ和2θ分别是两根杆绕x轴旋转的角度。系统有一个自由度, 取q=[x,θ1,θ2]T作为系统的广义坐标。作用于系统的两个完整约束表示如下:

图1 曲柄滑块机构Fig. 1 Slider-crank mechanism

对约束方程分别求一次和二次导数, 可以得到速度和加速度约束方程:

约束方程的雅可比矩阵及右端项可以写成如下形式:

约束方程的零空间为

同时, 可以得到式(21)中的广义质量矩阵和广义力列阵, 分别为

从图 2 可以看出, 零空间方法在仿真进行到 5 s左右时, 计算结果发生突变, 滑块的位移开始发散,不再做周期运动, 与曲柄滑块机构的实际运动不符。高斯优化方法在整个仿真过程中都表现得比较稳定, 而且约束方程只发生极小的违约(10–20)。由于增广拉格朗日方程中含 3 个不确定系数, 不同的取值对计算结果会有不同的影响, 本文取两组不同的参数对曲柄滑快机构进行动力学仿真, 分别为 µ=10, Ω=1 (case 1)和 µ=5, Ω=3 (case 2)。从图3 可以看出, 当增广拉格朗日方程取 case 1 中参数时, 系统的仿真比较精确稳定; 当取 case 2 中参数时, 仿真在 8 s 左右开始失效。因此, 增广拉格朗日方程是否保持平稳的仿真结果依赖于选取的参数,这会增加动力学计算的复杂度。高斯优化方法可以平稳地经过奇异位置, 并且只发生极小的违约, 说明与零空间方法和增广拉格朗日方法相比, 在不进行违约修正的前提下, 该方法更适合应用于含有奇异位置的多体系统的动力学求解。

图2 高斯优化与零空间方法的对比Fig. 2 Time history of displacement of point B and constraint violation obtained by GS and NS

图3 高斯优化与增广拉格朗日方程的对比Fig. 3 Time history of displacement of point B and constraint violation obtained by GS and ALF

4 结论

本文以广义坐标形式的高斯原理作为建模方法, 从数学优化的角度对动力学奇异问题进行求解。在求解极值问题的最优解时, 对传统的颗粒流算法进行改进, 将用传统优化方法得到的动力学方程无约束最优解以及前一时刻的最优解作为额外的颗粒添加到初始颗粒群, 并充分利用传统优化方法的快速收敛和智能优化的全局搜索特性, 达到提高颗粒寻优的计算效率的目的。分别采用高斯优化方法、零空间方法和增广拉格朗日方程, 对一个具有奇异位置的多体系统进行数值仿真, 结果表明高斯优化方法不仅具有较高的计算精度, 而且可以长时间保持数值计算的稳定, 不会因为系统自由度的突变而导致仿真失败。