一种基于增强响应灵敏度的多体系统参数识别方法

1.本发明涉及参数识别领域，尤其涉及一种基于增强响应灵敏度的多体系统参数识别方法。


背景技术：

2.典型的多体系统包括运载器、航天器、机器人、各种机械装置等。多体系统广泛存在于工业生产与日常生活中，对于多体系统的模型建立、优化设计、测试诊断过程等，识别其系统参数都显得尤为重要。多体系统的参数识别主要是根据系统的响应识别出系统的几何参数和惯性参数。传统的参数识别方法主要缺陷在于形成的是有偏估计，识别准确度不高。


技术实现要素：

3.为了解决上述技术问题，本发明的目的是提供一种基于增强响应灵敏度的多体系统参数识别方法，准确识别多体系统的参数，并且具有良好的抗噪性。
4.本发明所采用的技术方案是：一种基于增强响应灵敏度的多体系统参数识别方法，包括以下步骤：
5.s1、基于绝对坐标方法构建多体系统的动力学方程，并建立待识别参数最小二乘优化方程；
6.s2、载入测量数据并基于增广法求解多体系统的动力学方程，计算得到系统响应和响应残差；
7.s3、基于中心差分求取多体系统的响应灵敏度；
8.s4、基于待识别参数最小二乘优化方程建立迭代参数的最小二乘方程；
9.s5、引入tikhonov正则化对迭代参数的最小二乘方程进行处理，结合响应残差和相应灵敏度计算得到迭代参数；
10.s6、根据迭代参数对待识别参数进行更新，得到更新后的待识别参数；
11.s7、基于更新后的待识别参数重新求解多体系统的动力学方程，计算得到新系统响应和新响应残差。
12.s8、引入置信域计算一致性指标并与预设置信因子比较，调整迭代参数；
13.s9、判断到迭代参数与更新后的待识别参数的比值超出误差容限内，输出更新后的待识别参数，否则返回步骤s2。
14.优选地，所述多体系统的动力学方程公式如下：
[0015][0016]
上式中，b∈r
k
表示一个系统参数组成的列向量，分别是n
×
1广义坐标系下的位移、速度和加速度，φ∈r
m
表示完整约束方程，是jacobian矩阵，与lagrange乘子λ∈r
m
一同表示由约束方程引起的系统的反作用力，m∈r
n
×
n
表示系统的广义
质量矩阵，而q∈r
n
表示系统受到的广义外力矢量。
[0017]
对于位置φ∈r
n
‑
m
和速度φ
′
∈r
n
‑
m
的系统的初始条件可以表示为如下形式：
[0018][0019]
利用传统增广法求解方程。由于式是指标3的微分代数方程，先对中的约束方程关于时间t求二阶导，将其化为指标1的微分代数方程，并写成如下的等价形式：
[0020][0021]
其中，
[0022]
由式中后两式直接增广计算，得到
[0023][0024][0025]
由的第一式和即构成常规的一阶常微分方程，
[0026][0027]
上式可以利用matlab的ode进行求解。
[0028]
然而，由于上述方法只考虑了加速度形式的约束方程，所得到的响应未能严格符合位移或者速度的约束方程，会引起约束违约的问题。为了进行违约修正，通过引入人工阻尼将不稳定的约束方程转化为
[0029][0030]
从而，的第三式变为
[0031][0032]
再利用matlab求解即可。
[0033]
对于阻尼参数的选取，提出了参数α和β的一种自动选取方法，能得到很好的数值稳定性，即，
[0034][0035]
其中，h为求解微分方程的离散差分格式的步长。
[0036]
优选地，在典型参数识别流程中，对于一组实际测量得到的数据待识别的系统参数为b∈r
k
应该是如下非线性最小二乘优化问题的解，所述待识别参数最小二乘优化方程的公式如下：
[0037][0038]
上式中，表示待识别参数b的可行域，表示测量数据，r
i
(b)表示通过求解动力
学方程得到的理论解，为长度是pl的列向量的形式，符号表示关于权重矩阵w的加权范数。
[0039]
方程中的r(b)通常是关于b的非线性隐函数，于是需要利用迭代法求解如上的非线性最小二乘问题。接下来的关键就在于如何从已知的参数确定一个合适的迭代更新δb。通常将在点处线性化，保留一阶项，即
[0040][0041]
优选地，所述响应灵敏度的形式如下：
[0042][0043]
上式中，表示响应灵敏度矩阵，表示对应的响应灵敏度。
[0044]
具体地，基于梯度法解决反问题，灵敏度分析是不可或缺的。由于直接利用传统增广法求解多体系统的灵敏度复杂耗时，故采用有限差分法求取。本文中心差分求取响应灵敏度，即，
[0045][0046]
其中，ε是一个充分小的正数；而q(t,b
i
+ε)和q(t,b
i
‑
ε)可以通过求解正问题得到。本文ε取0.001。
[0047]
优选地，至此，即可得到关于δb的一个线性最小二乘问题，所述迭代参数的最小二乘方程公式如下：
[0048][0049]
上式中，δb表示迭代参数，表示响应残差。
[0050]
然而，问题通常是不适定的，因此需要引入tikhonov正则化进行处理，所述引入tikhonov正则化对迭代参数的最小二乘方程进行处理，公式表示如下：
[0051][0052]
上式中，μ≥0为正则化参数，||
·
||2表示l2范数，从而可得，
[0053][0054]
上式中，i为单位矩阵。正则化参数μ的值通常通过l曲线法选取，将其记为
[0055]
优选地，通过线性化得到近似的线性最小二乘问题。要使线性化后得到的能较好地近似原非线性函数，需要在参数的邻域内，亦即置信域[]内。那么，对于式加入的置信域约束可以表示为，存在足够小的δ＞0，使得||δb||＜δ。
由此可见，解得的δb应该是δ的一个函数。为了选取合适的δ值，需要引入如下的一致性指标，所述一致性指标的公式如下：
[0056][0057]
良好的一致性指标应该使得比某一特定的正数大，亦即，
[0058][0059]
从而需要选取合适的δ使得一致性条件满足。
[0060]
满足置信域的条件下，响应r(b)通常是关于b的光滑函数，且至少满足因此在满足的条件下，应有
[0061][0062]
以及
[0063][0064]
从上可知，通过适当的增加正则化参数μ，可以使迭代更新δb
μ
变小，并且满足置信域条件。
[0065]
最后，更新待识别的参数b
k+1
＝b
k
+δb，判断到迭代参数与更新后的待识别参数的比值在误差容限内，输出更新后的待识别参数，即||δb||/||b
k+1
||≤tol，输出更新后的待识别参数。
[0066]
本发明方法及系统的有益效果是：本发明基于增广法准确对多体系统的动力学方程进行求解，准确得到系统响应，通过引入置信域构成增强响应灵敏度方法，结合系统响应对多体系统进行参数识别，具有良好的抗噪性、稳定性和准确性。
附图说明
[0067]
图1是本发明一种基于增强响应灵敏度的多体系统参数识别方法的步骤流程图；
[0068]
图2是本发明具体实施例平面曲柄连杆机构的结构示意图。
具体实施方式
[0069]
下面结合附图和具体实施例对本发明做进一步的详细说明。对于以下实施例中的步骤编号，其仅为了便于阐述说明而设置，对步骤之间的顺序不做任何限定，实施例中的各步骤的执行顺序均可根据本领域技术人员的理解来进行适应性调整。
[0070]
参照图1和图2，第一具体实施例以平面曲柄连杆机构为例说明识别过程，曲柄连杆机构广泛存在于机械领域，比如常见的有活塞发动机、往复泵、往复式活塞压缩机等。如图2所示，单缸内燃机的曲柄连杆机构部分的模型。l1和l2分别是曲柄和连杆的长度，m1,m2和m3分别是曲柄、连杆和滑块的质量。取系统的状态变量为q＝[x1,y1,θ1,x2,y2,θ2,x3]
t
。除了考虑该机构受重力作用外，为了充分激励参数灵敏度，在曲柄假设受到如下的转动副摩擦模型，
[0071][0072]
从而，系统的广义质量矩阵为
[0073]
m＝diag[m1,m1,j1,m2,m2,j2,m3]
[0074]
其中，j1和j2分别是曲柄和连杆的转动惯量。
[0075]
待识别的系统参数取为曲柄和连杆的转动惯量b＝[j1,j2]
t
。利用绝对坐标方法可以建立起多体系统动力学方程：
[0076][0077]
其中，分别是广义坐标系下的位移、速度和加速度。φ是完整约束方程，是jacobian矩阵，与lagrange乘子λ一同表示由约束方程引起的系统的反作用力。m是系统的广义质量矩阵，而q是系统受到的广义外力矢量。
[0078]
取系统的广义坐标为q＝[x1,y1,θ1,x2,y2,θ2,x3]
t
。从而，系统的广义质量矩阵为：
[0079]
m＝diag[m1,m1,j1,m2,m2,j2,m3]
[0080]
系统所受的外力矢量为：
[0081]
q＝[0，
‑
m1g,τ(t),0,
‑
m2g,0,0]
t
[0082]
系统的约束方程为：
[0083][0084][0085][0086][0087]
φ5＝x3‑
l
1 cosθ1‑
l
2 cosθ2＝0
[0088]
φ6＝l
1 sinθ1‑
l
2 sinθ2＝0
[0089]
φ＝[φ1；φ2；φ3；φ4；φ5；φ6]
[0090]
相应的jacobian矩阵为：
[0091][0092]
结合动力学方程，便得到平面曲柄连杆机构的动力学方程。
[0093]
设真实系统的参数如下，l1＝1m,以及m1＝m2＝1kg,m3＝2kg，k
f
＝
‑
1待识别
的转动惯量的参数为j1＝0.1,j2＝0.2，初始条件为测量数据的采样率为100hz，测量时间为10s。
[0094]
利用传统增广法求解其动力学方程。由于多体系统动力学方程是指标3的微分代数方程，先对中的约束方程关于时间t求二阶导，将其化为指标1的微分代数方程，并写成如下的等价形式：
[0095][0096]
其中，
[0097]
由式中后两式直接增广计算，得到：
[0098][0099][0100]
由的第一式和即构成常规的一阶常微分方程，
[0101][0102]
上式可以利用matlab的ode进行求解。
[0103]
然而，由于上述方法只考虑了加速度形式的约束方程，所得到的响应未能严格符合位移或者速度的约束方程，会引起约束违约的问题。为了进行违约修正，通过引入人工阻尼将不稳定的约束方程转化为：
[0104][0105]
从而，的第三式变为：
[0106][0107]
再利用matlab求解即可。
[0108]
对于阻尼参数的选取，即，
[0109][0110]
其中，h为求解微分方程的离散差分格式的步长。
[0111]
在典型参数识别流程中，设定初始的系统待识别参数b，定义权重矩阵w，定义误差容限tol，设定最大迭代步数nmax，设定置信域流程的最大迭代步数ntr。对于一组实际测量得到的数据待识别的系统参数为b∈r
k
应该是如下非线性最小二乘优化问题的解，
[0112][0113]
其中，为参数b的可行域；而r(b)＝[y(b,t1),y(b,t2),
…
,y(b,t
l
)]是通过求解动力学方程得到的理论解，它们都表示为长度为pl的列向量的
形式。符号表示关于权重矩阵w的加权范数。最适权重矩阵w
optimal
应该与测量数据误差的方差成反比。
[0114]
方程中的r(b)通常是关于b的非线性隐函数，于是需要利用迭代法求解如上的非线性最小二乘问题。接下来的关键就在于如何从已知的参数确定一个合适的迭代更新δb。通常将在点处线性化，保留一阶项，即：
[0115][0116]
其中，其中，是响应灵敏度矩阵，其形式如下，
[0117][0118]
具体地，基于梯度法解决反问题，灵敏度分析是不可或缺的。由于直接利用传统增广法求解多体系统的灵敏度复杂耗时，故采用有限差分法求取。本实施例中心差分求取响应灵敏度，即，
[0119][0120]
其中，ε是一个充分小的正数；而q(t,b
i
+ε)和q(t,b
i
‑
ε)可以通过求解正问题得到。本实施例ε取0.001。
[0121]
至此，即可得到关于δb的一个线性最小二乘问题，如下所示，
[0122][0123]
然而，问题通常是不适定的，因此需要引入tikhonov正则化进行处理，可得，
[0124][0125]
其中，μ≥0为正则化参数，||
·
||2代表l2范数。从而可得，
[0126][0127]
上式中，i为单位矩阵。正则化参数μ的值通常通过l曲线法选取，将其记为
[0128]
部分的分析是通过线性化得到近似的线性最小二乘问题。要使线性化后得到的能较好地近似原非线性函数，需要在参数的邻域内，亦即置信域[]内。那么，对于式加入的置信域约束可以表示为，存在足够小的δ＞0，使得||δb||＜δ。由此可见，解得的δb应该是δ的一个函数。为了选取合适的δ值，需要引入如下的一致性指标(agreement indicator)
[0129]
[0130]
良好的一致性指标应该使得比某一特定的正数大，亦即，
[0131][0132]
从而需要选取合适的δ使得一致性条件满足。
[0133]
满足置信域的条件下，可知，响应r(b)通常是关于b的光滑函数，且至少满足因此在满足的条件下，应有
[0134][0135]
以及
[0136][0137]
从上可知，通过适当的增加正则化参数μ，可以使迭代更新δb
μ
变小，并且满足置信域条件。
[0138]
更新待识别的参数b
k+1
＝b
k
+δb，判断到迭代参数与更新后的待识别参数的比值在误差容限内，输出更新后的待识别参数，即||δb||/||b
k+1
||≤tol，输出更新后的待识别参数。
[0139]
利用增强响应灵敏度方法识别系统参数，结果如表2所示。
[0140]
表2
[0141][0142]
以上是对本发明的较佳实施进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可做作出种种的等同变形或替换，这些等同的变形或替换均包含在本技术权利要求所限定的范围内。