陈 敏， 刘 广， 汪 利， 吕中荣

(中山大学 航空航天学院，深圳 510275)

自20世纪60年代以来，人们日益关注机械、车辆、机器人的动力学特性。与此同时，电子计算机、数值计算方法的进步为这些复杂系统动力学的建模和分析提供了可能性[1]。在这样的背景下，人们开始研究由多个部件和运动副组成的动力学系统——多体系统。典型的多体系统包括运载器[2-3]、航天器[4-5]、机器人[6-7]、军事武器[8]、各种机械装置[9]等。正是由于多体系统广泛存在于工业生产与日常生活中[10]，对于多体系统的模型建立、优化设计、测试诊断过程等，识别其系统参数就显得尤为重要。多体系统的参数识别主要是根据系统的响应识别出系统的几何参数和惯性参数。从实际物理模型试验得到的数据可以用于识别只是已知其粗略大小的系统参数。

不同系统的组合、测量数据的要求、模型的误差以及优化的方法形成了不同的参数识别方法。对于多体系统，目前已经有很多识别其未知参数的优化方法。传统的非线性系统的识别方法有Gauss最小二乘法，Young提出的工具变量法，Fisher提出的最大似然估计法。其中，最小二乘法由于简易以及适用性广而被广泛利用。该方法的主要缺陷在于它所形成的是有偏估计，为了克服这一缺点，Schiehlen等[11]提出了在模型噪声是平稳随机过程的假设下，提出了基于协方差的参数识别方法。Serban等[12]利用最小二乘优化算法，并且通过直接差分法进行灵敏度分析，成功识别了简单的摆锤模型，最后利用矩阵的条件数提出了参数可识别性的数值测试手段。Bock等[13]发展了边值问题方法来估计微分代数方程的参数，主要是通过在有限差分或多重打靶法离散模型时引入一个可能的边界条件，并将其作为最小二乘目标函数的非线性约束。Uchida等[14]提出了一个新的办法来直接识别微分代数方程形式的多体系统方程，主要处理手法是将Lagrange乘子表示成关于时间的分段线性函数代入到动力学方程，通过标准的最小二乘算法就可以同时识别出系统的惯性参数和Lagrange乘子。Oberpeilsteiner等[15]利用伴随灵敏度方法结合二阶伴随系统识别多体系统的摩擦参数。他们与Lauss等之后将伴随灵敏度分析与Fourier分析相结合，提出了一个频域识别方法用于识别振荡多体系统，表现出很大的潜力和效率[16]。

本文利用Lü等[17]提出的增强响应灵敏度法用于识别多体系统问题。增强相应灵敏度方法首先通过时域测量数据与理论系统计算响应数据的误差的最小二乘估计，将参数识别问题转化为数学上典型的反问题，然后通过引入Tikhonov正则化解决该问题可能存在的不适定性，并通过置信域流程选取得到合适的正则化参数，再利用有限差分获得的一阶参数灵敏度数据，最终通过迭代法求解反问题得到系统的真实参数。由此可见，增强相应灵敏度方法只需利用到系统的动力学响应，在实际中只要试验测量时长足够，响应数据很容易通过传感器直接大量获取，这是其最明显的优势；再者，增强响应灵敏度法属于梯度法，收敛速度快，并且不像Newton法那样需要二阶灵敏度数据，再利用置信域流程，该方法计算效率更高，且对测量数据噪声不敏感；最后，该方法也已经成功地应用于振动系统[18]、混沌系统[19]、滞回系统[20]等的识别，均表现出了良好的识别精度、抗噪性和稳定性。

1 多体系统的动力学建模及求解

利用绝对坐标方法可以建立起多体系统动力学方程

(1)

对于位置φ∈Rn-m和速度φ′∈Rn-m的系统的初始条件可以表示为如下形式

(2)

利用传统增广法的Baumgarte方法求解方程式(1)。由于式(1)是指标3的微分代数方程，先对式(1)中的约束方程关于时间t求二阶导，将其化为指标1的微分代数方程，并写成如下的等价形式

(3)

其中

由式(3)中后两式直接增广计算，得到

(4)

(5)

由式(3)的第一式和式(4)即构成常规的一阶常微分方程

(6)

式(6)可以利用MATLAB软件的ode进行求解。

(7)

从而，式(3)的第三式变为

(8)

再利用MATLAB软件求解即可。

对于阻尼参数的选取，参考文献[22]提出了参数α和β的一种自动选取方法，能得到很好的数值稳定性，即

(9)

式中，h为求解微分方程式(3)所用的离散差分格式的步长。

2 多体系统的参数识别

2.1 反问题的建立

(10)

(11)

(12)

其中各项如下所示

至此，即可得到关于Δb的一个线性最小二乘问题，如下所示

(13)

然而，问题式(13)通常是不适定的，因此需要引入Tikhonov正则化进行处理，可得

(14)

(15)

2.2 时域灵敏度分析

基于梯度法解决反问题式(13)，灵敏度分析是不可或缺的。由于直接利用传统增广法求解多体系统的灵敏度复杂且耗时巨大，故采用有限差分法求取。本文中心差分求取响应灵敏度，即

(16)

式中：ε为一个充分小的正数；而q(t,bi+ε)和q(t,bi-ε)可以通过求解正问题式(1)得到。本文ε取0.000 1。

2.3 引入置信域

(17)

(18)

从而需要选取合适的δ使得一致性条件式(18)满足。

(19)

以及

(20)

从上可知，通过适当的增加正则化参数μ，可以使迭代更新Δbμ变小，并且满足置信域条件式(18)。

2.4 增强响应灵敏度法

通过以上几部分的分析，至此，建立起了用于识别多体系统参数的增强响应灵敏度法的实现过程。具体的算法流程，如表1所示。

表1 增强响应灵敏度的算法流程Tab.1 Algorithmic procedure for enhanced response sensitivity approach

3 数值算例

这个部分以平面连杆机构和空间连杆机构为例进行参数识别。系统的动力学方程利用MapleSim软件的Multibody Library建模提取，通过第1章所述的方法利用MATLAB软件求解。

表1中的各项参数的取值分别为tol=1×10-6，γ=1.414,ρcr=0.5,Nmax=1 000,Ntr=20。测量数据均由MapleSim仿真得到，若考虑噪声，均以服从标准正态分布的随机噪声进行模拟，表示为如下形式

(21)

3.1 平面曲柄连杆机构

考虑如图1所示的平面曲柄连杆机构，l1和l2分别是曲柄和连杆的长度，m1,m2和m3分别是曲柄、连杆和滑块的质量。取系统的状态变量为q=[x1,y1,θ1,x2,y2,θ2,x3]T。除了考虑该机构受重力作用外，为了充分激励参数灵敏度，在曲柄假设受到如下的转动副摩擦模型

图1 平面曲柄连杆机构Fig.1 A planar slider-crank system

(22)

从而，系统的广义质量矩阵为

M=diag[m1,m1,J1,m2,m2,J2,m3]

(23)

式中，J1和J2分别为曲柄和连杆的转动惯量。

待识别的系统参数取为b=[J1,J2]T。系统所受的外力矢量为

Q=[0,-m1g,τ(t),0,-m2g,0,0]T

(24)

系统的约束方程为

(25)

相应的Jacobian矩阵为

(26)

利用传统增广法，结合MATLAB软件的ode113求解结果与MapleSim软件仿真结果比较，如图2所示；其误差如图3所示。可以看出，引入Baumgarte稳定性理论的传统增广法得到的结果与仿真结果相比的误差相当小，误差也没有发散的趋势，说明正问题的求解流程是正确。接下来采用MapleSim软件的仿真数据加上不同水平的噪声来模拟实际的测量数据，选取连杆的平动位移数据x2和y2进行识别。如图4所示是噪声水平为20%的情况下的测量数据。取参数初始值b0=[1,1]T，利用增强响应灵敏度法识别系统参数，结果如表2所示。20%噪声水平下识别过程中目标函数的变化情况,如图5所示。可以看出，增强响应灵敏度法识别系统参数均在二十几个迭代步内收敛，且即使在20%噪声污染的情况下识别结果的误差在可以接受的范围内，该方法表现出了良好的收敛速度和抗噪性。

(a) 连杆水平方向的位移响应

(b) 连杆竖直方向的位移响应图2 平面曲柄连杆机构的MATLAB求解结果(ode113)与MapleSim仿真结果Fig.2 The calculation results of calculation(ode113) and the simulation results of MapleSim for planar crank-slider system

(a)

(b)图3 平面曲柄连杆机构的计算结果(ode113)与仿真结果的误差Fig.3 The error between the calculation results(ode113) and the simulation results of planar crank-slider system

(a) 受噪声污染的连杆水平方向的位移响应

(b) 受噪声污染的连杆竖直方向的位移响应图4 20%高斯噪声污染的测量数据Fig.4 Measured response polluted by Gaussian noise at level 20%

图5 20%高斯噪声水平下目标函数随迭代步数的变化曲线Fig.5 Objective function curves with respect to the iteration steps at 20% Gaussian noise level

表2 平面曲柄连杆机构的参数识别结果Tab.2 Parameter identification results of the planar crank-slider system

3.2 空间曲柄连杆机构

图6 空间曲柄连杆机构Fig.6 Spatial slider-crank system

表3 空间曲柄连杆机构的各项参数及其取值Tab.3 The parameters of the spatial crank-slider system and their values

取q=[s,θ,α,β]T为广义坐标，由MapleSim软件建模提取得到系统的质量矩阵为

(27)

系统受到的广义外力矩阵为

(28)

系统的约束矩阵为

(29)

其对应的Jacobian矩阵为

(30)

由此，结合方程式(1)便建立起了空间连杆机构的动力学方程。同样地，利用传统增广法，结合MATLAB软件的ode113求解结果与MapleSim软件仿真结果比较，如图7所示；其误差如图8所示，可以看出计算结果和仿真结果相当吻合，说明正问题的求解方法是正确的。利用MapleSim软件在0～10 s内的仿真数据q作为测量数据，采样率为200 Hz，利用式(21)加入不同水平的高斯噪声模拟，待识别的参数取为b=[Ix1,Iy1,Iy2,Iz2]T，参数初值取为[1,4.2,1,4.2]T，同样地，利用增强响应灵敏度法识别系统参数，结果如表4所示。可以看出，增强响应灵敏度法识别系统参数能够成功识别出系统的惯性参数，且即使在20%噪声污染的情况下识别结果的误差依旧在可以接受的范围内，同样地表现出了良好抗噪性。

(a) 滑块的位移响应

(b) 曲柄的位移响应

(c) 万向铰α角位移响应

(d) 万向铰β角位移响应图7 空间曲柄连杆机构的计算(ode113)与仿真结果对比Fig.7 The comparison between the calculation results(ode113) and the simulation results of spatial crank-slider system

(a)

(b)

(c)

(d)图8 空间曲柄连杆机构的计算结果(ode113)与仿真结果的误差Fig.8 The error between the calculation results(ode113) and the simulation results of spatial crank-slider system

表4 不同噪声水平下空间连杆机构的识别结果Tab.4 The parameter identification results of spatial crank-slider system on different noise levels

值得注意的是，由迭代步数可以看出，本算例在不同噪声水平下的收敛步数差异较大。对于无噪声的情况，收敛速度很慢，是因为系统的参数灵敏度并没有被外力充分激励，导致灵敏度过小，使得收敛速度大大降低。合适的系统激励可以通过参考文献[24-25]等的方法得到，但只适用于系统响应可以视作与参数线性相关的情况，并不存在适用于一般系统的情况[26]。

有噪声的情况下，收敛速度较快得益于置信域的引入。实际上，若不引入置信域，5%噪声的情况下，一般的灵敏度算法会在47步循环后得到结果[0.145 2,3.956 9,0.359 5,3.423 7]T，此时算法远远未满足循环终止的容限tol，与真实参数的差值将会越来越大，这是由于缺少了置信域流程的抗噪性，且系统的参数灵敏度过低，而噪声污染的影响相比之下很大，从而无法正确地收敛到真实参数。

4 结 论

本文利用了增强灵敏度方法识别多体系统的参数，通过两个算例表现出了良好的识别精度、收敛速度以及抗噪性。然而，实际实现过程其实存在着很多挑战，多体系统的特性会带来很多问题，诸如系统的最优激励问题，系统参数的独立性和可识别性等，未来的工作将着力解决以上问题，使得该方法能更好地应用于多体系统的参数识别。