萤火虫算法求解多体系统动力学微分-代数方程*

张笑笑 丁洁玉 ，2†

（1.青岛大学数学与统计学院，青岛266071）（2.青岛大学计算力学与工程仿真研究中心，青岛266071）

引言

微分-代数方程是描述多体系统动力学的数学模型，它由微分方程和代数约束方程组成，其数值求解方法的研究是多体系统动力学研究的重要内容.保持约束稳定从而保证微分-代数方程求解的稳定是数值方法研究的重点［1］.近年来逐渐发展起来的一些方法，如：坐标缩并方法［2］、广义-α方法［3］、辛算法［4］、能量方法［5，6］、辛保持的变分数值积分方法［7］、Lie 群方法［8］等，均取得了较好的求解结果.

萤火虫算法（firefly algorithm，FA）是XS.Yang在2008年提出的一种新颖的智能优化算法［9］.该方法的思想是模拟萤火虫闪烁行为，通过萤火虫间的相互吸引来实现种群的迭代进化.近年来，研究人员对标准FA算法进行诸多改进，较好地解决算法收敛速度慢、求解精度低、易陷入局部最优等缺陷.例如：将变量从混沌空间变换到解空间然后再进行搜索的混沌萤火虫算法［10，11］；通过将FA与其他算法结合提出混合萤火虫算法［12-14］；基于离散空间的离散萤火虫算法［15］等.目前该类算法已经在诸多优化问题中得到了较好的效果，并且在机器学习、生产调度、机器人智能控制和图像处理等领域得到了广泛的应用.

本文主要通过Lagrange插值，结合Gauss数值积分方法将微分-代数方程求解问题转换成优化问题，然后对该优化问题进行萤火虫算法设计.最后对平面双连杆系统进行仿真实验.

1 多体系统动力学微分-代数方程

多体系统动力学方程通常为指标3的微分-代数方程组（DAEs）：

其中，qϵRn是广义坐标，q̇ϵRn是广义速度，λϵRd是Lagrange乘子，Φ为广义坐标q的约束方程，Φq为约束方程的Jacobi矩阵.

将约束方程Φ(q，t)=0求两阶导，方程（1）可以由指标3降为指标1，方程形式如下：

将仿真时间[0，T]平均划分为若干小区间[ti，ti+1]，i=0，1，…，N.在时间 [ti，ti+1]中对广义坐标q(t)进行Lagrange插值：

则

为保证位移约束、速度约束和加速度约束，插值函数q(t)、q̇(t)、q̈(t)需要满足如下 3d个约束方程：

通过方程组（7）将方程（2）中的变量进行缩并，3d个 变 量 (q1，...qd；，...；，...)可 以 由3(n-d)个 变 量 (qd+1，...qn；q̇d+1，...q̇n；q̈d+1，...q̈n)来表示 .然后将q(t)、q̇(t)、q̈(t)代入微分-代数方程（2）得到：

若对于∀t∈ [ti，ti+1]都有f(t)=0，则插值函数q(t)为式（1）在[ti，ti+1]上的精确解.显然非线性函数（8）一般不存在精确解，只能求解相对于插值函数q(t)尽可能满足（8）式接近于0的相对最优解.则可以转换成如下优化问题：

对上式进行Gauss数值积分，可以化成如下形式：

其中，As为Gauss积分系数，ts为Gauss积分节点.显然，式（10）是一个非凸优化问题［16］，已有的研究表明，传统的数学优化方法难以对此类问题进行有效求解，智能优化算法是求解此类问题的有效算法.萤火虫算法是模仿种群运动的一种智能算法，提供了求解DAE方程的新方法.

2 萤火虫算法设计

萤火虫算法包含一个有M个个体的萤火虫种群系统 X={X1，X2，X3，…，XM}，Xi∈RD.第i个个体在时间t处的位置：Xi(t)=(xi1(t)，xi2(t)，…，xiD(t)).萤火虫对彼此吸引的原因取决于两个因素，即自身亮度Li和吸引度β.其中，萤火虫的亮度取决于自身所在位置的目标值Li=φ(Xi).吸引度与亮度相关，愈亮的萤火虫吸引亮度比其弱的萤火虫往这个方向移动.吸引度与距离Ri，j成反比，距离越远的萤火虫吸引度越低.系统X的整体最优位置为Xg.迭代过程中，萤火虫的位置更新公式如下［17］：

其中，γ是光吸引系数；β0是最大吸引度；α∈ (0，1)是步长因子；α为扰动项；υ为动态视觉权重，越大寻优视野越大，越容易获得远距离信息；Ri，j为萤火虫i和j之间的距离.

迭代初期，希望对萤火虫种群个体进行较大扰动，以增强全局探索能力有利于跳出局部极值点；而在搜索后期，要求对个体减小扰动，以提高算法的搜索精度.因此，设计扰动项如下：

运用萤火虫算法求解区间[ti，ti+1]的插值函数q(t)，已知初始点q0=q(ti) ∈Rn，插值节点qj，j=1，2…m是待优化的节点，并且qj满足约束Φ(q，t)=0。.一般有e个约束方程，则维度d=(n-e)m.定义目标函数如下：

如果φ(qj)=0，则qj，j=1，2…m即为所求的最优值.由于目标函数φ为非线性函数，一般不存在解析解，所以φ(qj)只能尽可能地接近精确解，设置收敛精度φmin和最大迭代步数gmax，若达到收敛精度或达到最大迭代步数，则输出寻优结果.

萤火虫算法的求解步骤如下：

1）设置萤火虫算法的初始参数：种群规模M，萤火虫位置维度d，萤火虫位置的变化范围，最大迭代次数gmax，收敛精度φmin，光吸引系数γ；最大吸引度β0；步长因子α，初始扰动项=0，并随机生成初始种群，定义目标函数φ；

2）根据设定的目标函数，分别计算每个个体的适应度值，比较种群的适应度，确定种群的亮度最强个体Xg；

4）比较移动前后的萤火虫适应度值，若优于之前位置则完成位置更新，否则保持原位置不移动；

5）如果φ(Xg)<φmin或达到最大迭代步数gmax，则输出最优解gbest，否则重复步骤2）-步骤5），直到满足预设的终止条件.

3 数值算例

图1 平面双连杆机械臂Fig.1 Two-link planar manipulator

两点Gauss数值积分方法求解目标函数.

图2-图5为时间步长0.005时系统仿真结果.从连杆末端运动轨迹和系统能量变化可以看出，连杆运动过程稳定，没有出现明显偏移.位移约束、速度级约束和加速度级约束没有出现违约.

图2 连杆末端运动轨迹Fig.2 The trajectory of the end of rod

图3 系统总能量Fig.3 Total System Energy

图4 系统动能、势能Fig.4 System Kinetic Energy and Potential Energy

图5 双连杆机械臂约束Fig.5 Constraints of the two-link manipulator

采用相同时间步长h=0.01对平面双连杆系统进行仿真实验，结果比较见表2，其中FA-DAE表示本文的求解算法，RK4表示四阶Runge-Kutta方法，仿真时间为20s.

从表1可以看出，系统运动过程中精确保持位移约束、速度级约束和加速度级约束，系统总能量最大相对误差较小.随着时间步长减小，系统总能量最大相对误差和平均相对误差在减小，仿真结果也更加的准确.从表2可以看出，本文算法运行时间长，但在约束和能量方面保持较好.

表1 不同时间步长的误差结果比较Table 1 Comparison of errors for different time steps

表2 相同时间步长时各方法结果比较(h=0.01，t=20s)Table 2 Comparison of different methods with the same time step(h=0.01，t=20s)

4 小结

本文运用萤火虫优化算法求解多体系统动力学微分-代数方程，并对平面双连杆系统进行仿真实验.实验结果表明，在满足指标1、指标2和指标3约束的情况下，系统总能量的误差较小.萤火虫算法的优化精度高、算法设计简单且不需要目标函数的导数信息，对于求解多体系统微分-代数方程取得了较好效果，说明智能优化算法应用到多体动力学仿真中的可行性.然而，该方法存在计算量大和算法运行时间长的缺点.今后在有效地降低运行时间和提高算法精度方面是一个主要研究方向.