Application of Particle Swarm Optimization on the Multi-body System Dynamics with Singular Positions

YANG Liusong1, YAO Wenli2,†, XUE Shifeng1,†

1. College of Pipeline and Civil Engineering, China University of Petroleum (East China), Qingdao 266580; 2. College of Science, Qingdao University of Technology, Qingdao 266520; † Corresponding authors, E-mail: ywenli1969@sina.com (YAO Wenli); sfeng@upc.edu.cn (XUE Shifeng)

Abstract Different from the traditional method, the mathematical optimization model is established with Gauss principle to handle the singular problems to deal with the singular problems. Traditional optimization method and intelligent optimization method (particle swarm optimization algorithm, PSO) are combined to solve the above optimization problem, which can fully utilize the fast convergence of the traditional optimization method and the characteristic of global searching of the intelligent algorithm. The numerical example is simulated by Lagrangian formulation, null space method and Gauss optimization method respectively. The simulation results show that Gauss optimization method has higher computational accuracy, keeps the stability of the numerical calculation and would not lead to simulation failure due to the sudden changes of the system degree of freedom, which validates the effectiveness and universality of the proposed method. Key words multi-body system; singular problem; Gauss principle; particle swarm optimization (PSO)

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

式(1)~(4)构成封闭的微分–代数混合方程组。一些多体系统常常由于自身的几何构型的原因, 在运动过程中出现奇异位置, 系统的自由度发生突变, 导致运动出现分叉点[3]。利用微分‒代数混合方程组对其进行动力学仿真时, 在奇异位置处,约束方程(式(2))的雅可比矩阵的秩rank(c )  m , 给

q计算造成很大的困难, 导致仿真结果与实际情况不符, 甚至发生严重的失真。引起这些问题的直接原因是系统经过奇异位置时约束反力的突变, 根本原因在于奇异位置处约束方程允许运动的子空间扩展[4]。当系统离开奇异位置时, 扩展的子空间的运动不再满足约束方程, 致使系统的运动偏离实际运动。为了克服约束雅克比矩阵缺秩引起的约束反力突变, 一些数值方法会引入较大的冲击载荷来消除不在约束空间内的速度分量[5–6]。如果数值方法不够稳定, 过大的冲击载荷也可能导致仿真失败。Kurdila等[7]给出一种修正的Lagrange 方程, 适用于具有奇异位置的多体系统, 但由于存在较强的刚性,不具有普适性。Bayo等[8]提出一种增广拉格朗日方法, 通过选择较大的惩罚因子, 对约束方程的违约进行限制, 最后利用迭代算法来求解动力学方程组。还有一些学者将增广拉格朗日方法应用到非完整多体系统中, 用于提高计算机模拟的实时仿真效率[9–10]。增广拉格朗日方程已成功地应用于机械系统的研究和仿真, 如重型机械模拟器[11]、生物力学[12]和车辆动力学的联合模拟设置[13]等。Aghili等[11]和

[12] Phong 等 利用雅可比矩阵的零空间, 将微分‒代数方程组转换为常微分方程组。不论系统是否存在奇异位置或冗余约束, 利用零空间方法得到的微分方程组形式紧凑, 其系数矩阵始终保持正定状态,方程组具有唯一解。上述方法都可以处理冗余约束和具有奇异位置的动力学问题, 但都有局限性, 比如, 基于增广拉格朗日方程方法依赖于惩罚参数的选取, 零空间方法中广义逆的数值计算比较耗时,影响仿真效率。近年来, 关于高斯原理的研究呈增长趋势。姚

[13–17]文莉等 将极值形式的高斯原理拓展到多刚体系统的单边接触动力学问题。Udwadia等[18]对受约束的多体系统运动提出新的解释, 并基于高斯原理和广义逆矩阵推导出显示形式的动力学方程。刘延柱[19–20]将高斯原理应用到柔体和刚柔耦合系统, 推导高斯拘束函数的表示方法, 并探讨利用最小值模型解决多体系统动力学的数值方法。杨流松等[17, 21]利用高斯原理, 对受冗余约束作用的多体系统进行建模, 将力学中的动力学问题转换为一个求函数极值的数学问题, 并采用智能优化方法进行求解, 可以大大地提高仿真的计算效率。上述关于高斯原理的研究没有涉及多体系统的奇异位置问题。本文利用高斯原理, 建立具有奇异位置的多体系统的动力学优化模型, 直接采用数学优化方法求解, 将具有奇异位置的多体系统动力学问题转化为求函数极值的优化问题。目前, 优化问题的求解手段主要分为两类: 第一类是传统算法,对于凸问题求解效率高, 但不适用于多极值问题,所求解多为局部最优解; 另一类是智能算法, 此类算法具有一定的随机性, 适用于多峰优化问题, 鲁棒性强, 不易陷入局部最优。粒子群优化(particle swarm optimization, PSO)是一种基于种群的随机优化算法, 由于易于实施, 并且可以快速地收敛到最优解, 近年来受到广泛关注, 并应用于人工神经网络训练、功能优化、模糊控制和模式识别等领域。本文对PSO进行改进, 并将智能优化与传统优化方法相结合, 充分利用传统优化的快速收敛和智能优化优化的全局搜索特性, 求解高斯方法建立的动力学优化模型, 探讨其应用于求解动力学奇异位置问题的可行性。

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

[5] Bayo 等 根据多体系统动力学的哈密顿描述,提出增广拉格朗日方法。该方法具有普遍性, 并且可以应用到完整系统和非完整系统中。与传统的动力学模型不同, 增广拉格朗日方法通过增加一部分

惩罚项, 使动力学模型转化为一组常微分方程, 表达形式如下:

Mq  +  T  ( C   2 C  2 C )  qt  Q, (5) q其中,  是拉格朗日乘子,  是由惩罚因子组成一个对角方阵, 和 是常系数对角方阵。这里的系数和 与约束稳定理论中的常系数有相同的作用, 但在增广拉格朗日方法中, 这两个参数有特别的物理意义, 分别代表与约束方程相关的自然频率和阻尼系数。

从式(5)可以看出, 与传统的拉格朗日方法不同, 增广拉格朗日方法不需要附加约束方程组, 它允许约束直接插入动力学方程中, 将最初的DAES转化成常微分方程组。实际上, 增广拉格朗日方法是将惩罚的约束方程合并到动力学方程中, 从而使约束条件得到满足。为了求解式(5), 需要先计算出乘子 。将式(3)和(4)代入式(5), 可以得到其中, 第一次迭代时, q 满足 Mq  Q。整个迭代

0过程一直运行到满足条件 q   q  为止,  为用

i 1 i户自定义的误差容限, 本文取 =10  6 。在整个动力学问题的求解过程中, 即便系统经过奇异位置, 从而导致雅克比矩阵Cq缺秩, 但只要选择合适的惩罚因子 , 式(9)的系数矩阵 M + CT  Cq 就可以始终保

q持正定和对称, 从而使得系统的加速度具有唯一解。所以, 可以将增广拉格朗日方法作为一种建模方法, 用来求解具有奇异位置的多体系统动力学问题。

通常情况下, 惩罚因子 取值越大, 约束方程的违约越小, 但不能过大, 否则 M + CT  Cq 会变得

q病态。但是, 从式(8)可以看出, 在迭代过程中, 乘子可以对约束方程不满足的地方进行补偿, 所以在用增广拉格朗日方程求解动力学问题时, 惩罚因子通常取106~107之间的常数, 本文取106。

1.2 零空间方法(NS)

[12]为了克服动力学系统中的奇异位置问题, Phong等 利用约束方程的雅克比矩阵的零空间, 推导出一种形式紧凑的动力学方程——零空间方法。假设H 满足以下条件:事实上, 不论系统受到冗余约束的作用或经过某些奇异位置引起自由度的突变, 由于 rank ( H T M)  rank (C )=n 始终成立, 所以式(13)的方程数目始q终与系统的广义坐标的数目保持相等。所以, 不论雅克比矩阵Cq是否缺秩, 式(13)始终是一个正定方程组, 系统加速度具有唯一解, 并且可以直接求出。此外, 虽然零空间可以简化系统的动力学方程,但同时也增加计算的负担。利用约束惯性矩阵的零空间, 系统的加速度可以显式地表达为

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

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

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

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

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

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

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

因此, 如果采用传统的优化算法求解奇异约束条件下的优化问题, 难点在于如何在保证收敛效率的前提下跳出局部最优, 得到全局最优解。近年来,受人类智能、生物群体社会性和自热现象规律启发, 人们发明了很多智能算法, 用于解决复杂的高维数和多极值优化问题。智能算法不依赖于初值,算法独立于求解域, 具有全局优化的性能, 通用性强, 适用于并行处理。但是, 智能算法具有一定的随机性, 所以也存在一些缺陷, 比如收敛速度慢、局部搜索能力差以及控制的参数较多等。本文充分利用传统优化的快速局部寻优能力和智能优化算法的全局搜索性能, 将二者相结合来求解多体系统的动力学优化模型, 克服因动力学问题的奇异性导致的计算精度和效率问题。本文将这种采用高斯原理进行建模, 并采用优化算法的动力学求解方法称为高斯优化方法(Gauss optimization method, GS)。

2 计算方法2.1 粒子群算法简介

为了解决约束优化问题, 有学者开发了不同的

[24]确定算法和随机算法 。由于可行方向和广义梯度下降等确定性方法对目标函数的连续性和可微性给出较强的假设, 而在实际问题中这些条件很难全部满足, 所以它们的适用性有限。随着优化理论的发展, 一些智能算法得到迅速发展和广泛应用, 为非线性、多极值的复杂函数及组合优化问题提供了切实可行的解决方案, 如遗传算法、进化策略、进化规划和粒子群优化等。粒子群优化算法作为一种启发式全局优化技术[25–26], 可以模拟鸟群随机搜寻食物的捕食行为。鸟群通过自身经验和种群之间的交流, 调整自己的搜寻路径, 从而找到食物。粒子群优化每次搜寻时都会根据自身经验(自身历史搜寻的最优位置)和种群交流(种群历史搜寻的最优位置), 调整自身搜寻方向和速度, 从而找到最优解。

假设在一个D维的目标搜索空间中, 有 τ个颗粒组成一个群落, 其中第i个粒子的位置表示为一个D维的向量 X  ( Xi 1, Xi , ..., X ) ; 第i个颗粒的

i 2 id历史最优位置为 P  ( P , P , ..., P ) ; 整个颗粒群

i i 1 i 2 id迄今为止搜索到的最好位置记为 P  ( P , P , ..., g g1 g2 P ); 第i个颗粒的运动速度也是一个D维的向量gd V  ( Vi 1, Vi , ..., VID )。对于颗粒i, 第k+1次迭代时的

i 2

颗粒速度和位置可以表示为

2.2 算法流程

本文对经典的PSO进行改进, 在文献[25]的基础上, 增加一个传统无约束优化的可行解和上一时间步的最优解作为额外的颗粒, 添加到初始颗粒群中, 以便达到快速收敛和全局寻优的目的。针对奇异位置引起的求解空间的扩张, 本文将变异算子引入粒子群优化算法中, 进一步增加种群的多样性,使算法尽快跳出局部最优解, 防止算法过早收敛。新的粒子群算法如下:其中  为变异因子, 取值为 2。对具有奇异位置的多体系统进行仿真的操作流程如下。1) 预先给定仿真时间tend 和时间步长t , 给出满足约束方程的初始条件: 速度 q0 和位移 q0 , 令t  t0 。2) 计算系统的质量矩阵M, 包括外力和速度相关项的广义力矩阵Q, 约束方程的雅克比矩阵Cq以及约束方程的右端项。3) 用式(22)建立多体系统的动力学模型, 并按照下列步骤, 求解得到系统此时的加速度q 。t

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作为系统的广义坐标。作用于系统的两个完整约束表示如下:  l cos  l cos2  x 1 C  。 l sin  l sin2  1对约束方程分别求一次和二次导数, 可以得到速度和加速度约束方程:  x  l sin( 1)  lsin( 2) 2  C  ,   cos( )  cos( )  1 1 2 2  x  l(sin( )   sin( )  )     l(cos( ) + cos( ) )。  cos( )   cos( )      sin(    sin(  约束方程的雅可比矩阵及右端项可以写成如下形式:  l sin  l sin2 1 1  ,  cos cos 0  1 2

, 系统处于奇异位置, 约

1 2束方程的雅可比矩阵的秩rank(cq )  1, 而在其他位置时 rank(cq )  2。由此可见, 在奇异位置处, 约束方程的雅可比矩阵奇异, 矩阵的秩发生突变。系统的初始条件如下:将式(29)~(36)分别带入式(5)、(14)和(22), 可以得到用增广罚函数法、零空间方法与高斯优化方法建立的动力学模型。取步长t =0.002 s, 以龙哥库塔算法作为积分器, 在不采用违约修正方法的前提下, 分别用3种方法对曲柄滑块机构进行时长为10 s的动力学仿真, 计算结果如图2和3所示, 其中

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

这会增加动力学计算的复杂度。高斯优化方法可以平稳地经过奇异位置, 并且只发生极小的违约, 说明与零空间方法和增广拉格朗日方法相比, 在不进行违约修正的前提下, 该方法更适合应用于含有奇异位置的多体系统的动力学求解。

4 结论

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

参考文献

[1] Fasano A, Marmi S. Analytical mechanics: an introduction. Oxford: Oxford University Press, 2006 [2] 潘振宽, 洪嘉振. 多体系统动力学微分/代数方程组数值方法. 青岛大学学报(自然科学版), 1996, 9(1): 83–96 [3] 许永生, 齐朝晖. 具有奇异位置的多体系统动力学方程的改进算法. 动力学与控制学报, 2006, 4(2): 109–113 [4] González F, Dopico D, Pastorino R, et al. Behaviour

of augmented Lagrangian and Hamiltonian methods for multibody dynamics in the proximity of singular configurations. Nonlinear Dynamics, 2016, 85(3): 1491–1508 [5] Bayo E, Ledesma R, Arbor A. Augmented Lagrangian and mass-orthogonal projection methods for constrained multibody dynamics. Nonlinear Dynamics, 1996, 9: 113–130 [6] García de Jalón J, Gutiérrez-lópez M D. Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces. Multibody System Dynamics, 2013, 30(3): 311–341 [7] Kurdila A J, Junkins J L, Hsu S. Lyapunov stable penalty methods for imposing holonomic constraints in multibody system dynamics. Nonlinear Dynamics, 1993, 4(1): 51–82 [8] Bayo E, Garcia De Jalon J, Serna M A. A modified Lagrangian formulation for the dynamic analysis of constrained mechanical systems. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 183–195 [9] Bayo E, Jalón J G D, Avello A, et al. An efficient computational method for real time multibody dynamic simulation in fully Cartesian coordinates. Computer Methods in Applied Mechanics and Engineering, 1991, 92(3): 377–395 [10] Potosakis N, Paraskevopoulos E, Natsiavas S. Application of an augmented Lagrangian approach to multibody systems with equality motion constraints. Nonlinear Dynamics, 2020, 99(1): 753–776 [11] Aghili F, Piedboeuf J C. Simulation of motion of

constrained multibody systems based on projection operator. Multibody System Dynamics, 2003, 10(1): 3–16 [12] Phong D V, Hoang N Q. Singularity-free simulation of closed loop multibody systems by using null space of Jacobian matrix. Multibody System Dynamics, 2012, 27(4): 487–503 [13] 姚文莉, 戴葆青. 广义坐标形式的高斯最小拘束原理及其推广. 力学与实践, 2012, 36(6): 779–785 [14] 姚文莉. 基于广义坐标形式的高斯最小拘束原理的多刚体系统动力学建模. 北京大学学报(自然科学版), 2016, 52(4): 708–712 [15] Yao Wenli, Yang Liusong, Song Kewei, et al. Optimization method for dynamics of non-holonomic system based on Gauss’ principle. Acta Mechanica Sinica, 2020, 36(5): 1133–1141 [16] 姚文莉, 刘彦平, 杨流松. 基于高斯原理的非理想系统动力学建模. 力学学报, 2020, 52(4): 945–953 [17] Yang Liusong, Xue Shifeng, Yao Wenli. Application of Gauss principle of least constraint in multibody systems with redundant constraints. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, 2020, 235(1): 150–163 [18] Udwadia F E, Kalaba R E. Response to Bucy’s comment

on a paper by Udwadia and Kalaba. Mathematical, Physical and Engineering Sciences, 1996, 452: 1055–1056 [19] 刘延柱. 杆网系统基于高斯原理的动力学建模. 动力学与控制学报, 2018, 16(4): 289–294 [20] 刘延柱. 关于高斯原理与动力学的优化法建模//中国力学大会会议论文集(CCTAM 2019). 北京: 中国学术期刊电子出版社, 2019: 4094–4096 [21] 杨流松, 姚文莉. 基于高斯最小拘束原理的广义质量矩阵奇异性问题研究. 力学研究, 2017, 6(1): 56– 62 [22] Meyer C D. Matrix analysis and applied linear algebra. Philadelphia: SIAM, 2001 [23] Gauß C F. Über ein neues allgemeines grundgesetz der mechanik. Journal für die Reine und Angewandte Mathematik, 1829, 4: 232–235 [24] Kochenderfer M J, Wheeler T A. Algorithms for optimization. Cambridge: MIT Press, 2019 [25] Kennedy J. Particle swarm optimization. Boston: Springer, 2017: 967–972 [26] Eberhart R, Kennedy J. A new optimizer using particle swarm theory // Proceeding of the 6th International Symposium on Micromachine and Human Science. Nagoya, 1995: 39–43

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

(a) B点位移随时间的变化; (b) 约束方程的时间历程图 2高斯优化与零空间方法的对比Fig. 2 Time history of displacement of point B and constraint violation obtained by GS and NS

(a) B 点位移随时间的变化; (b) 约束方程的时间历程图 3高斯优化与增广拉格朗日方程的对比Fig. 3 Time history of displacement of point B and constraint violation obtained by GS and ALF