基于动力学的拟人肩关节动载协调分配优化研究*

王 林，李研彪*，孙 鹏，罗怡沁，徐梦茹，郑 航

(1.浙江工业大学 机械工程学院，浙江 杭州 310032；2.浙江工业大学 特种装备制造与先进加工技术教育部重点实验室，浙江 杭州 310032)

0 引 言

并联机构具有结构紧凑、承载能力强、运动惯性小等优点，故广泛应用于各种拟人关节[1-6]。由于并联机构的动力学方程是非线性多输入多输出系统，存在多种分配组合来实现运动[7-10]。在稳定外载荷作用下，合理地优化分配各驱动可以有效降低能耗、减少驱动力。

目前，主要从驱动力矩最优和能耗最优两个角度来进行动载协调分配，通过构建力矩分配模型和能量分配模型来实现优化[11-15]。其中，分配方式有加权最小二乘法、规划载荷分配系数等。但上述优化方法均属于单目标优化，只考虑驱动力矩或者能耗，且均从瞬时状态考虑，未衡量机构的整个运动过程。因此，优化结果可能存在驱动力矩、速度等方面的波动或突变，这种波动、突变现象不利于机构的稳定运行。同时，只考虑力矩和能耗最小，也可能出现机构的运动时间较长，不利于提高机构的运动效率。

针对上述动载协调分配优化的不足，本文将考虑机构性能、时间、能耗和力矩波动4个因素，提出一种动载协调分配优化方法，并基于动力学模型，定义综合性能指标，采用Dijkstra算法优化求解性能最优轨迹，最后用遗传算法求解得到最优广义时间。

1 肩关节机构的运动学分析

本文研究的肩关节机构以球面5R并联机构为原型，肩关节的模型图如图1所示。

图1 肩关节的模型图

该机构由定平台、动平台和连接两者的二条支链组成。运动副均为回转副，且各轴线汇交于O点。其中，轴线OB1与轴线OC1、轴线OA2与轴线OC2、轴线OC1与轴线OC2相互垂直。

肩关节的机构简图如图2所示。

图2 肩关节的机构简图α1—平面A1OZ和平面A2OZ所在平面的夹角，α1=90°；α2—Z轴与轴线OA1的夹角，α2=60°；α3—轴线OA1与轴线OB1的夹角，α3=70°；d1—回转副C1、C2与机构中心O点的距离，d1=70 mm；d2—回转副B1与机构中心O点的距离，d2=90 mm；d2—回转副A1、A2与机构中心O点的距离，d2=185 mm

本研究建立定坐标系{O-XYZ}原点与机构中心O重合，Z轴沿OC1轴线方向，Y轴沿OA2轴线方向，X轴满足右手螺旋定则。建立动坐标系{O-X1Y1Z1}原点与机构中心O重合，Z1轴沿OC1轴线方向，X1轴沿OC2轴线方向，Y1轴满足右手螺旋定则。当定坐标系和动坐标系重合时，该机构处于初始位姿。

本研究采用Z-Y-X型的欧拉角描述动平台的姿态，动平台绕Z1轴旋转角度为α，绕Y1轴旋转角度为β，绕X1轴旋转角度为γ，其中α=0 rad。关节驱动器1输入角度为θ1，关节驱动器2输入角度为θ2。

本研究根据肩关节机构的几何关系，建立矢量约束方程，化简可得位置反解：

(1)

式中：

其中，cαi=cosai，sai=sinai，(i=1,2,3)。

将式(1)两边对于时间t求导，速度反解可得：

(2)

2 肩关节机构的动力学分析

为了便于建立动力学模型，本研究将肩关节机构的各连杆和动平台均视为刚体，同时忽略各运动副之间的摩擦力和轴类零件的回转运动。考虑惯性力、外力作用，建立肩关节机构的动力学模型。

2.1 惯性力分析

本研究采用拉格朗日方程计算肩关节机构的惯性力，将系统惯性力转换到广义欧拉坐标q=[γβ]T上。其计算过程如下：

肩关节机构的动能E包括E1和E2(其中：E1—平台的动能；E2—各连杆的动能)。

(3)

Ip表示过质心坐标系的转动惯量矩阵，即：

Ip=RIHRT

(4)

式中：R—动平台欧拉角所对应的旋转矩阵。

综合式(3,4)，可得动能E1:

(5)

动能E2的求解如下：

(6)

式中：I1—连杆A1B1的转动惯量；I2—连杆A2C2的转动惯量；I3—连杆B1C1的转动惯量。

则肩关节机构的动能为：

(7)

设定坐标系中OXY面为重力零势能面，由于动平台的质心与机构中心O点重合，且动平台始终绕质心旋转，动平台势能变化为零。则机构的总势能V为：

(8)

式中：zi—各连杆质心坐标的Z轴值；mi—各连杆的质量。

建立拉格朗日方程可得：

(9)

式中：L=E-V，FI—惯性力。

将式(7，8)代入上式，化简可得：

(10)

2.2 外力分析

设作用在动平台上的外力为F，均可简化为过旋转中心O点的力矩Ms：

Ms=F×r

(11)

式中：r—外力F作用点到旋转中心O点的矢量。

2.3 动力学模型建立

综合式(10，11)，根据虚功原理，将惯性力FI和外力Ms映射到相应关节上的驱动力矩：

(12)

式中：JT—力雅克比矩阵。

根据式(12)建立的动力学模型可知，肩关节机构驱动力与惯性力和外力两者有关。其中，惯性力大小受速度、加速度和姿态影响，而外力仅与姿态有关。

2.4 动力学仿真验证

给出肩关节机构的结构参数如表1所示。

表1 肩关节机构结构参数

外力F=[1 1]N·m，并给出一组动平台的运动方程：

(13)

基于上述动力学模型，将肩关节机构的结构参数和运动方程代入式(12)中，并利用Matlab软件计算得到肩关节机构的关节驱动力矩。

为了验证动力学模型的正确性，本研究利用ADAMS动力学仿真软件对肩关节机构进行动力学仿真，且添加各类约束条件和外力，使得仿真与理论计算的环境保持一致。笔者将仿真得到的驱动力矩与理论计算值进行比较，如图3所示。

图3 关节驱动力矩理论值与仿真值

通过对比可得：驱动力矩的理论值与仿真值两者基本相等，故验证了动力学模型的正确性。

3 肩关节机构的性能分析

3.1 动力学性能评价指标

肩关节机构不同的运动状态会影响驱动力的大小，其中，速度、加速度会直接影响机构惯性力的大小。而肩关节机构一般作低速运动，速度较小，故忽略速度对惯性力的影响，只考虑加速度对惯性力的影响，因此式(10)可简化为：

(14)

(15)

将式(15)求导，化简可得：

(16)

式中：λ—矩阵DTD的特征值。

由式(16)可知：在工作空间W内，λ随机构姿态变化而变化，且λ值越小表明由加速度引起的惯性力越小，因此将km作为动力学传递性能评价指标：

(17)

km数值越大，表示肩关节机构的动力学传递性能越好，则km在工作空间W的全域值为：

(18)

且肩关节机构存在多个输入力矩，其差值越大则动力学传递性能越差，因此将kmc作为动力学传递均衡性能评价指标：

(19)

kmc数值越接近于1，表示肩关节机构的输入力矩偏差越小，其动力学传递均衡性能越好，则kmc在工作空间W的全域值为：

(20)

3.2 力映射性能评价指标

由动力学模型可知，惯性力和外力通过虚功原理，将其映射为关节的驱动力矩。其中，力雅克比矩阵JT代表该映射关系，其值受机构姿态影响，即：

τ=JTτF

(21)

由于在工作空间W内，rank(J)=2，力雅克比矩阵JT可以奇异值分解，存在正交阵U∈R2×2和V∈R2×2，使：

JT=UΛV

(22)

设欧拉坐标系上的广义力τF为单位向量可得：

τTU(ΛΛT)-1UTτ=1

(23)

当广义力τF为单位矩阵时，关节驱动力矩分布在椭圆上；σi值越大，表明广义力转换到关节驱动力矩的效率越高，因此将kJ作为力传递性能评价指标[17-20]：

kJ=σ2

(24)

kJ数值越大，表示肩关节机构的力传递性能越好。则kJ在工作空间W的全域值为：

(25)

当σ1=σ2时，关节驱动力矩分布在圆上，关节驱动力矩之间的差值最小，因此将kJc作为力传递均衡性能评价指标：

(26)

kJc数值越接近于1，表示肩关节机构的力传递均衡性能越好，关节驱动力矩之间的偏差越小。则kJc在工作空间W的全域值为：

(27)

4 肩关节机构的动载协调分配优化

在稳定外载荷作用下，要求肩关节机构从起点运动到终点，但由于其运动轨迹和运动状况的不确定性，存在多组运动方式。本文先依据性能指标得到性能最优的轨迹，再通过遗传算法优化求解最优广义时间，使得肩关节机构在能耗最低、时间最快、性能最好的情况下完成运动。

4.1 性能最优路径规划

为综合考虑动力学性能和力映射性能因素，笔者采用加权求和法，将上述多性能指标转换成单性能指标。

通过式(18，20，25，27)，分别求解出各性能指标在全域范围内最大值kimax和最小值kimin，则不同姿态下的各性能指标可表示为：

(28)

式中：Ki—性能指标在全域变换范围内的比值，Ki值越大，机构性能越好。

综合考虑上述4个性能指标的影响，构造综合性能指标函数：

(29)

式中：δi—目标比重系数。

δi值越大表示该衡量指标越重要。本文δi均等于1，将各性能指标视为相同比重，且Kmin值越小，机构的综合性能越好。

给定机构起点q0=[0 0]T，终点qt=[-1.5 1]T，在起点至终点区域，均匀采样n个控制点(每个控制点代表一种运动姿态，即动平台运动经过的点)，且每个控制点对应一个综合性能指标值。从起点运动到终点，建立路径规则。从起点出发，可以到达相邻3个控制点(控制点1、控制点2和控制点4)；之后，再从某一个控制点出发，可以达到下一个相邻的3个控制点；以此类推，最终达到终点。

路径规则如图4所示。

图4 路径规则

采用上述方式，可以得到数个由控制点组成的路径。根据式(29)，笔者将一条路径上所有控制点的综合性能指标值求和，代表该条路径的综合性能。采用Dijkstra算法优化求解出从起点到终点的最优路径，满足路径最短、综合性能最好。

性能最优路径如图5所示。

图5 性能最优路径

由于规划出的路径为折线，无法实现光滑轨迹运动，基于上述性能最优路径的结果，本研究采用B样条曲线拟合方法构造动平台运动轨迹，并得到动平台运动方程：

(30)

式中：u—路径中第u个控制点，u∈[0，15]。

4.2 广义时间优化

根据上述建立的性能最优轨迹可知，动平台的运动轨迹与参数u有关。通过建立参数u与时间t的函数，可以进一步确定机构的速度特性。因此，建立参数u关于时间t的函数：

u(t)=a0+a1t+a2t2+a3t3+a4t4

(31)

式中：ai—时间函数的系数。

由于机构需要满足起点和终点的位置约束条件，且在起点和终点处速度为零，故建立时间约束方程：

(32)

式中：u0=0，u1=15；T—机构完成运动的时间。

同时，考虑肩关节机构的关节驱动器输出速度和力矩存在最大值，故建立关节输出约束方程：

(33)

为实现肩关节机构在最短时间内完成运动，建立时间优化目标函数：

(34)

机构在运动过程中，可能存在时间很短而能耗很大或者关节力矩变化很大的情况，属于不理想的工作状况。因此，需要考虑能耗和关节力矩因素，故建立能耗优化目标和力矩波动优化目标：

(35)

(36)

式中：Fτ—力矩波动优化目标，为机构关节运动的力矩变化率，衡量关节力矩的波动情况。

综合考虑时间、能耗和力矩波动，建立广义时间优化目标函数：

(37)

则优化过程如下：首先根据式(30，31)计算出动平台的角度和角速度；其次通过式(1，2)得到关节驱动角度和角速度；然后根据式(13)动力学模型计算出关节驱动力矩；通过式(37)得到广义时间优化目标函数；采用遗传算法对时间函数式(31)进行优化求解，得到最优的广义时间。设定种群个数为100，基因数为80，迭代代数为110，变异概率为0.1，交叉概率为0.5。

广义时间优化结果如图6所示。

图6 广义时间优化结果

由图6可知：经过数代运算后收敛，得到优化结果为：时间T=2.03 s；时间函数系数a0=0，a1=0，a2=8.85，a3=-1.58，a4=-0.48。将优化结果代入式(31)得到时间函数，并根据式(1，2，12，30)，通过Matlab软件计算，得到优化后关节驱动力矩、关节角速度。

关节驱动力矩如图7所示。

图7 关节驱动力矩

关节角速度如图8所示。

图8 关节角速度

由图7和图8可知：关节驱动力矩和关节角速度均为光滑曲线，不存在突变现象，说明优化后，关节驱动器能够实现平稳运行。在起点和终点时刻，关节输出角速度均接近于零，符合实际运动状况。两个关节的驱动力矩和角速度值均处于同一数量级，数值大小接近，说明该优化方法能够把机构惯性力和外力合理地分配到两个关节上，说明采用上述动载协调分配优化方法，肩关节机构能够在时间最短、能耗最低、性能最好的条件下完成运动。

5 结束语

本研究结合拉格朗日方程和虚功原理，建立了肩关节机构的动力学模型，并用动力学仿真软件验证了上述动力学模型的正确性；然后基于动力学模型，构建了动力学性能指标和力映射性能指标，采用加权求和法将各性能指标转换为综合性能指标，并利用Dijkstra算法求解得到综合性能最优的轨迹；最后基于性能最优轨迹，考虑时间、能耗和力矩波动3个因素，建立广义时间优化目标函数，采用遗传算法优化时间函数，得到最优的广义时间。

通过上述动载协调分配优化，使得肩关节机构能够在时间最短、能耗最小、性能最好的条件下完成运动。