1. 引言
多体系统动力学为武器、航空、航天、车辆、机器人、精密机械等多个工程机械系统的动态性能评价和优化设计提供了强有力的工具 [1] 。在对多体运动系统进行建模时，由于选择的广义坐标的相关性，系统的广义质量矩阵往往容易发生奇异，虽然在经典的无约束动力学系统中，奇异质量矩阵并不常见，但奇异质量矩阵常常出现在一些受到约束的复杂多体系统之中。当系统的广义坐标数目超过了系统所需要的最小值，系统的质量矩阵最容易出现奇异。一般来说，系统模型的灵活性越大，奇异质量矩阵出现的可能性就越大 [2] 。
另外，在对多体系统进行编程自动化的过程中引入冗余约束几乎是不可避免的, 而且如果系统存有奇异构型，系统的自由度可能会随系统的运动而发生变化，从而导致冗余约束的产生。总之，在独立约束的前提下来解决多体系统动力学问题很难实现。对于存在冗余约束的多体系统，需要剔除约束方程中的多余约束保证约束方程的独立性 [3] 。判断约束独立性的方法有很多，其中，增广拉格朗日公式声称可以解决多体系统中的冗余等奇异问题 [4] 。
但是包括增广拉格朗日公式在内的一些数值方法都不能解决具有奇异构型的多体动力学问题。这些方法之所以失效，就是因为系统运动到接近奇异构型的位置时，系统的广义质量矩阵就会发生奇异，约束方程违约增大，从而导致了加速度的突变。
因此，质量奇异矩阵问题在复杂多体系统动力学的研究中是不可忽略的问题，但对于该问题的研究并不多见，其中F.E. Udwadia及R.E. Kalaba提出的显式方程可用于质量奇异矩阵的情形 [1] ，但该方程在应用到非理想多体系统时会存在问题，所以其有效性还有待商榷 [5] 。
基于高斯最小拘束原理来研究刚体系统的动力学问题的优化方法，可以将动力学问题通过寻找函数极值的变分方法直接得出运动规律，而无需建立运动微分方程 [6] 。本文采用广义坐标形式的高斯最小拘束原理来研究经典情况下的难点奇异性问题，通过引入广义逆，建立了广义质量矩阵奇异情形下的高斯最小拘束原理，研究了针对奇异性问题的优化方法的数值求解策略。并将数值模拟结果同经典的拉格朗日方程方法进行对比分析，旨在从建模方式出发，研究解决经典动力学框架下的难点问题的方法。
2. 广义坐标形式的高斯最小拘束原理
建立多体系统动力学模型时通常需要采用笛卡尔广义坐标或铰链坐标的形式，若采用质点笛卡尔形式的高斯拘束 [6] 来建立多体系统动力学的优化模型会比较繁琐，下面将引入广义坐标形式的高斯最小拘束原理。
研究受约束的系统，其广义坐标可表达为
，此组广义坐标可以是不独立的。
系统动能的表达式为：
其中
，
，
分别为广义速度的二次项,一次项和零次项；为系数矩阵，M为系统的广义质量矩阵。
可据此定义下列矩阵：
，
2.1. 广义质量矩阵矩阵A为非奇异的情形 [7]
此时高斯最小拘束可以取为：
(1)
即在时刻
，
，
固定的情况下，在所有满足约束的可能运动中，实际运动所对应的高斯最小拘束最小。
在上述高斯拘束G中，Q为广义坐标所对应的广义力矩阵，g为含广义坐标及广义速度的矩阵。即
，
，
，
其中
， 而
为系数矩阵
的第一类Christoffel符号。
2.2. 广义质量矩阵为奇异矩阵的情形
在对多刚体系统进行动力学建模时，由于选择的广义坐标的相关性，系统的广义质量矩阵M往往容易发生奇异，若M是奇异方阵，我们引入M的广义逆
来代替
[8] ，则此种情形下的高斯拘束可以写作：
(2)
即在时刻
,
,
固定的情况下，在所有满足约束的可能运动中，实际运动所对应的高斯最小拘束最小。
2.3. 讨论
1) 在上述定理(1)及(2)中，广义加速度并未要求是独立的，它的适用性很广，对完整系统、非完整系统及单、双边约束都是成立的；
2) 这样多体系统动力学问题可以基于上述高斯最小拘束原理转化为有约束或无约束的求最小值的优化问题；
3) 方程(1)及方程(2)可分别用于质量非奇异矩阵的优化形式的动力学建模；
3. 基于高斯最小拘束原理的数值求解策略
对于含约束的上述形式的高斯最小拘束原理，根据函数求极值的必要条件，可由上式推导出微分代数方程组的形式，该形式与第一类拉格朗日方程得到的微分代数方程组等价，本文不采用微分代数方程组的求解，而是直接采用有约束或无约束的优化的方法来解决。其目标函数为关于广义加速度的二次非线性函数，本文选用智能优化算法来求解。具体求解步骤如框图1所示。首先对系统用高斯最小拘束原理建模，之后对系统的广义质量矩阵进行判断，看其是否奇异，分别用(1)式和(2)式建模。在对系统进行数值模拟时，采用智能优化算法——遗传算法借助MATLAB对系统的运动规律进行仿真。
4. 算例
算例的几何描述如下图2所示，多体系统由1、2、3杆组成，可在铅锤平面内运动，1杆的受到
的初始力矩，初始位置
，
。系统参数如表1。
该平面四连杆系统是单自由度系统，多体系统自动建模时往往采用铰链法，如取
为系统的广义坐标，因此会出现多余的广义坐标。
由图示中所给参数可以看出，系统的几何构型为平行四边形，得到如下约束条件：
(3)
则系统的动能表达式可得到系统的广义质量矩阵M。
Figure 1. The block diagram of the numerical solution of the singular mass matrix
图1. 质量奇异矩阵的优化数值求解框图
Figure 2. The plane four-bar linkage
图2. 平面四连杆的机构参数
Table 1. Mechanical parameters of planar four link
表1. 平面四连杆的机构参数
，此质量矩阵为奇异方阵，需要采用广义逆来替代通常意义的逆矩阵，
，为系统的广义力矩阵。
4.1. 采用第一类拉格朗日方程建模
(4)
将(3)式微分两次后与(4)联立，且也需要引入相应的广义逆矩阵，可以求得相应的广义加速度及拉格朗日乘子，通过积分得到下个时间积分步的广义坐标及广义速度。用MATLAB进行数值仿真可以得到1杆y方向的位移曲线。
4.2. 优化方法处理
确定高斯拘束(2)，采用图1框图所示的计算步骤进行计算，用MATLAB处理后的结果如下图所示。
4.3. 讨论
为了讨论两种方法在解决广义质量矩阵问题方面的优劣，对两种结果做如下处理，
，
其中，分别为两种不同建模方法所得到的1杆质心沿y轴的运动曲线函数。图5为两种方法中对于杆1质心在y方向的差值图。图3和图4是针对系统质量矩阵奇异情况下分别采用第一类拉格朗日方程及高斯最小拘束原理对1杆质心在y方向的运动轨迹分析图，从图5可以看出两种方法的误差曲线的数量级为，可以认为两种方法是等效的，但因采用优化方法建立的模型不需要引入拉格朗日乘子，因此计算效率更高。
Figure 3. The motion curve of the 1 rod centroid in y direction
图3. 1杆y方向的运动曲线
Figure 4. The motion curve of the 1 rod centroid in y direction
图4. 1杆质心y方向运动曲线
Figure 5. The difference of two methods in y direction
图5. 两种方法在y方向的差值图
5. 结论
质量矩阵奇异是在多体系统动力学在自动建模中经常遇到的难点问题，本文从建模方式出发，通过引入广义逆矩阵，将该类问题从通常的求解微分代数方程转变为采用求函数极值的优化问题，并建立了求解该问题的数值求解方案，算例表明文中给出的方法可以解决质量矩阵的奇异性问题，同时因为运算过程中无需引入拉格朗日乘子，从而使其具有更高的计算效率。
基金项目
本文作者感谢国家自然科学基金课题11272167的支持。