具有解析式正解的Stewart衍生型并联机构的位移输入协调关系

叶鹏达，尤晶晶,2*，沈惠平，吴洪涛，茹 煜

(1. 南京林业大学 机械电子工程学院，江苏 南京 210037；2. 江苏省精密与微细制造技术重点实验室，江苏 南京 210016；3. 常州大学 现代机构学研究中心，江苏 常州 213016；4. 南京航空航天大学 机电学院，江苏 南京 210016)

1 引 言

1965年，STEWART首次提出含6条相同支链的机构，学者们称其为Stewart并联机构[1]。相比于串联机构，并联机构具有输出精度高、结构刚性好和承载能力强等优点，成为国内外机构学研究热点[2-3]。目前，大多数6自由度并联机构都是基于Stewart并联机构衍生而来，我们称之为Stewart衍生型并联机构，其在光学元器件、并联机床和六维加速度传感器等领域得到越来越多的应用[4]。

然而，由于拓扑构型较复杂，6自由度并联机构的精确实时控制一般很难实现。有研究发现[5-7]，基于并联机构的装置性能主要与其正向运动学问题有关。正向运动学的研究方法主要有解析法和数值法两种[8-9]。解析法主要包括消元法、Grobner基法和共形几何代数等方法。文献[10]基于计算机符号运算，利用矢量消元和Sylvester结式消元，得出一般6-4台体型Stewart并联机构位置正解的一元32次方程；文献[11]运用计算机代数系统中的分次字典序Grobner基算法，获得一般6-6平台型Stewart并联机构位置正解的一元20次方程且该机构最多有40组解的结论；文献[12]基于共形几何代数建立运动学方程，并通过构造9阶Sylvester结式，得到一般6-6台体型Stewart并联机构位置正解的一元40次方程。不难看出，这些方法最终都将并联机构位置正解表达为一个一元高次代数方程，其优点是可以得到全部解，但消元和推导过程繁琐，技巧性强，且计算量大，特别是针对台体型Stewart并联机构。数值法是另一种求解并联机构正向运动学的方法，主要运用Newton-Raphson法(N-R法)等数值逼近迭代思路求解非线性方程组。文献[13]针对6-6平台型Stewart并联机构提出了利用位置反解逐步迭代的方法求解位置正解的思路，保证了实时性，但每次只能迭代计算一组解；文献[14]针对6-3 Stewart平台结构提出了一种基于机构简化的位置正解数值方法，保证了实时性，但对机构结构具有一定的依赖性。不难看出，并联机构拓扑结构变得复杂时，数值法对于精确解的收敛性较差[15-17]。文献[18]引入传感器提高了模型算法的收敛性，但由于引入传感器而造成结构上的复杂性，最终对算法的收敛性并没有很大帮助。

考虑到并联机构位置正解的求解难度与机构耦合度指标有关[19]，并且冗余驱动具有减少奇异位形和增加有效工作空间等优点[20]，文献[21]提出了一种低耦合度冗余驱动Stewart衍生拓扑构型，由于尚未有与之相匹配的二重复合球铰链，因而限制了其结构模型的设计。鉴于此，本文设计了二重复合球铰链。除此之外，还设计了可转换主动、从动模式的移动副，该移动副可用于改变机构的冗余度，实现变冗余驱动。

从以上分析可以看出，目前针对冗余驱动并联机构的研究主要集中在机构内力协调，而对于移动副位移协调以及驱动模式优选方面的研究还较少，并且大多数研究的机构都是具有确定的冗余度，而对于变冗余度机构的研究相对较少。本文推导并化简了6个位移输入协调方程，为了便于对其进行验证，重新推导出机构的位置正解；然后，分别运用N-R法和Broyden法求解了协调方程；进一步研究不同冗余度对冗余协调算法数值结果的影响；最后，总结出冗余驱动模式优选的3点选取原则。

2 结构模型及正向运动学

2.1 结构模型

图1 Stewart衍生型并联机构结构模型

Stewart衍生型并联机构的结构模型和拓扑构型分别如图1和图2所示。该机构由1个边长为2(n+L)的静平台、1个边长为2n的动平台、6个二重复合球铰链、12个可转换主动、从动模式的移动副和12个一般球铰链组成。12个一般球铰链与静平台相连，6个二重复合球铰链与动平台相连。初始状态下，静平台与动平台的中心重合，并且姿态相同。当驱动N(6≤N≤11)个移动副时，动平台的位置和姿态发生变化，从而实现动平台的六维运动。

图2 Stewart衍生型并联机构拓扑构型

如图3所示，为了降低机构的耦合度，设计了一种二重复合球铰链[25]，其主要由第一层铰链、第二层铰链和中心柱组成。第一层铰链由第一套筒、内叉和第一滑柱组成，第二层铰链由第二套筒、外叉和第二滑柱组成。中心柱与动平台相连，滑柱与移动副相连，两层铰链均具有3条互相垂直的转动轴线，并且转动轴线始终相交于点Bi(i=1,2,...,6)。

如图4所示，为了改变机构的冗余度，设计了一种可转换主动、从动模式的移动副[26]，其主要由导杆、内套筒、外套筒和转换套筒组成。导杆与内套筒之间滑动连接，内套筒与外套筒之间通过螺纹连接。当转换套筒和导杆通过螺栓螺母连接在一起时，随着内、外套筒之间发生相对转动，移动副的长度也发生改变，此时，该移动副处于主动模式；当拆掉螺栓螺母时，导杆与转换套筒分离，内套筒与导杆之间可相对滑动，此时，该移动副处于从动模式。

图3 二重复合球铰链结构模型

图4 移动副的三维图和剖面图

2.2 正向运动学

并联机构的正向运动学是指：已知驱动副的输入量，求解动平台的输出运动参数。惯性坐标系{O-XYZ}与静平台相连，动坐标系{M-UVW}与动平台相连，如图2所示，惯性坐标系的坐标原点位于初始状态下的动平台中心处。12个一般球铰链中心点在惯性坐标系中的笛卡尔坐标可表示成矩阵形式：

[b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12]=

(1)

式中：b1～b12为同名点的几何中心在惯性系内的笛卡尔坐标；n为动平台半边长；L为移动副初始长度。

12个移动副长度可表示为：

(2)

式中：Bi为同名点的几何中心在惯性系内的笛卡尔坐标；|·|为矢量的模。

点M为动平台的几何中心，根据机构拓扑构型并且运用立体几何知识可得出M点与Bi点之间满足如下的关系：

B1+B4=B2+B5=B3+B6=2M，

(3)

B1+M=B2+B3，

(4)

(5)

为了便于计算，选取M，B1，B2和B34个点为待求量：

(6)

根据12个移动副的长度约束关系并结合式(5)，建立动平台输入、输出量之间的位姿映射方程组。进一步分析后发现，同构方程两两相减，能够消去二次项，得到4组线性封闭方程：

(7)

(8)

(9)

(10)

其中：

再结合式(4)，得到剩余3个未知量的解析式：

(11)

至此，重新推导了新型并联机构的正向运动学，具体算法流程如图5所示。

图5 位姿正解算法流程图

为了验证位姿正解模型的正确性，在Mathematica中构建了Stewart衍生型并联机构的虚拟样机，如图6所示。其中，动平台半边长n和移动副初始长度L分别设置为15 mm和25 mm。虚拟仿真中，动平台为刚体，并且不考虑球面副、移动副的摩擦和间隙。不失一般性，取机构的初始状态和任意状态进行仿真，计算结果如表1所示。结果显示，位姿正解的计算值与准确值完全一致。

图6 Stewart衍生型并联机构的虚拟样机

表1 位姿正解模型验证

3 协调方程的推导及求解

3.1 协调方程的推导

将动平台视作刚性结构，运动过程中位置点之间的相对距离保持不变。由并联机构结构模型可知，移动副长度与位置点之间的尺度约束关系满足如下关系式：

|B1-b1|2=x12+

[y1-(n+L)]2+(z1+n)2=l12，

(12)

|B2-b3|2=(x2+n)2+

[y2-(n+L)]2+z22=l32，

(13)

|B3-b5|2=(x3-n)2+

y32+[z3+(n+L)]2=l52，

(14)

|B1-M|2=(x1-x0)2+

(y1-y0)2+(z1-z0)2=2n2，

(15)

|B2-M|2=(x2-x0)2+

(y2-y0)2+(z2-z0)2=2n2，

(16)

|B3-M|2=(x3-x0)2+

(y3-y0)2+(z3-z0)2=2n2，

(17)

(18)

(19)

(20)

将式(7)～式(11)代入式(12)～式(20)，发现方程的最高次数达到8次，并且方程中所含高次项的项数也较多。为了降低方程的次数以及减少方程中高次项的项数，通过式(12)，式(15)相减、式(13)，式(16)相减、式(14)，式(17)相减、式(12)，式(13)，式(18)相减、式(13)，式(14)，式(19)相减、式(12)，式(14)，式(20)相减，得到如下6个位移输入协调方程：

f(X)=2x1x0+2y1y0+2z1z0-x02-

y02-z02+3n2+(n+L)2-2y1(n+L)+

2z1n-l12=0,

(21)

g(X)=2x2x0+2y2y0+2z2z0-x02-

y02-z02+3n2+(n+L)2-2y2(n+L)+

2x2n-l32=0,

(22)

v(X)=2x3x0+2y3y0+2z3z0-x02-y02-

z02+3n2+(n+L)2+2z3(n+L)-

2x3n-l52=0,

(23)

u(X)=-2x1x2-2y1y2-2z1z2-2x2n+

2y2(n+L)-2n2-2(n+L)2+

2y1(n+L)-2z1n-(2n2-l12-l32)=0,

(24)

w(X)=-2x2x3-2y2y3-2z2z3-2x2n+

2y2(n+L)-2n2-2(n+L)2-

2z3(n+L)+2x3n-(6n2-l32-l52)=0,

(25)

h(X)=-2x1x3-2y1y3-2z1z3-2z3(n+L)+

2x3n-2n2-2(n+L)2+2y1(n+L)-

2z1n-(2n2-l12-l52)=0.

(26)

将式(7)～式(11)代入式(21)～式(26)可知，消除了式(21)、式(22)、式(23)中的8次方项以及式(24)、式(25)、式(26)中的6次方项。化简前后的协调方程对比如表2所示。

表2 化简前后的协调方程对比

Tab.2 Comparison of the compatibility equations before and after simplification

方程特征项化简前化简后数目96最高次数88最高次项数459

为了验证协调方程的正确性，并联机构结构参数与2.2节中的参数一致，在Mathematica中仿真对比了6个协调方程在两种不同状态下的计算结果。表3列出了两种状态下协调方程的计算误差，将6个协调方程误差的平均值作为综合误差。误差δ的表达式为：

δ=|F(X)*-F(X)*|,

(27)

其中：

F(X)=

式中：上标“*”代表仿真结果；下标“*”代表准确结果；X为从动副长度。

由表3可知，当移动副位于初始状态时，计算结果没有产生误差；当移动副位于任意状态时，计算结果产生较小误差，该误差可能是由计算机浮点运算中的舍入误差和泰勒展开时的截断误差等因素造成，可以忽略不计。计算结果表明协调方程推导正确。

表3 协调方程的误差统计

3.2 协调方程的求解

并联机构12个移动副之间满足一定的位移输入协调关系，由于协调方程的次数均高于4次，给定N个驱动副长度，可通过数值方法计算出其它(12-N)个从动副长度。求解非线性方程组的数值方法主要包括不动点迭代法、牛顿法和拟Newton法等方法，本文选择其中的N-R法和Broyden法进行求解，其它方法将另文研究。

3.2.1 N-R法

N-R法是求解非线性方程组的经典方法，因其收敛速度快和自校正等优点，得到了广泛的应用。目前很多新的算法都是在此基础上改进而来。

协调方程可以表示为：

F(X)=0.

(28)

泰勒公式展开：

F(Xk)+

F′(Xk)(Xk+1-Xk)+O(|Xk+1-Xk|2)=0,

(29)

其中：

式中：k为迭代步数；F′(Xk)为Jacobian矩阵。

忽略二阶无穷小量可得到：

Xk+1=Xk-[F′(Xk)]-1F(Xk).

(30)

式(30)即为N-R法的迭代公式。N-R法是局部收敛的，只要初值选取得当，计算结果总能收敛到合理值。N-R法求解协调方程的算法流程如图7中实线所示。

3.2.2 Broyden法

Broyden法是基于改进的N-R法发展而来，是求解非线性方程组的重要方法之一。它克服了N-R法需要求导、求逆等缺点，将Jacobian矩阵简化为矩阵递推关系式，这样不仅简化了计算过程，同时还能保证方法的超线性收敛速度[27]。当然，在改进的N-R法成立的前提下，Broyden法是计算修正矩阵ΔAk的一种方法，其它方法还包括Broyden第二方法、DFP秩2算法和BFGS秩2算法等。

Broyden首先对改进的N-R法提出了一种修正Ak的计算方案，它每步迭代只需计算n个分量函数值及O(n2)次算术运算，大大降低了计算量，其迭代公式为：

(31)

其中：

Sk=Xk+1-Xk

yk=F(Xk+1)-F(Xk).

(32)

式中：SkTBkyk≠0。逆Broyden秩1公式每步总计算量为n个分量函数值和O(n2)次算术运算。

式(32)即为Broyden法的迭代公式。Broyden法求解协调方程的算法流程如图7中虚线所示。

图7 两种迭代算法流程图

观察图7可以发现，与N-R法相比，Broyden法在迭代过程中不需要计算Jacobian矩阵，而是用当前函数值代替导数，避免了每次迭代过程中的求导，降低了计算量。

为了对比数值结果的精度和效率，将两种算法编写成Mathematica程序，并在软件中通过Timing指令获取算法的计算时间。不失一般性，任意给定动平台一组位姿：x0=0.5 mm，y0=0 mm，z0=0 mm；λ1=-0.05，λ2=0.05，λ3=0.05，迭代精度控制为1.0×10-6，仿真结果如表4所示。表5列出了许可初值偏差从20%变化到25%时，计算结果收敛到给定位姿所运行的迭代步数和计算时间。

表4 协调方程数值结果精度对比

表5 协调方程数值结果效率对比

Tab.5 Efficiency comparison of numerical results of the compatibility equations

许可初值偏差/%N-R迭代步数时间/msBroyden迭代步数时间/ms205187.510156.3215203.110156.3225203.111156.3235203.111156.3246218.811171.9256250.012187.5平均值5.3210.910.8164.1

观察并分析表4，5可知， N-R法数值结果的精度优于Broyden法至少3倍，但效率较低。从迭代步数和计算时间可以看出，在相同许可初值偏差情况下，虽然Broyden法收敛到满足精度要求的给定位姿所运行的迭代步数较多，但由于每一步迭代不需要构造Jacobian矩阵，因此计算量较少，计算时间仅为N-R法的78%左右。

将两种算法的许可初值偏差均取为30%，其迭代收敛性如图8所示，其中综合杆长表示计算杆长的平均值。为了更清晰地反映两种算法的收敛性，选取图中计算值与准确值较接近的部分曲线进行局部放大。从图8中可以看出，使用Broyden法计算得到的杆长收敛曲线围绕准确值的扰动较小，并且在迭代步数达到4步时，计算值与准确值基本吻合。因此，Broyden法的收敛速度较快，即收敛性较优。

图8 两种迭代算法的收敛性

4 协调方程数值结果的影响因素

通过切换移动副的工作模式，从而改变机构的冗余度，实现6种驱动模式，如图9所示，其中R代表机构的冗余度。协调方程数值结果的影响因素主要包括许可初值偏差和并联机构结构参数，针对各种驱动模式，分别对比分析这两者对冗余协调算法的影响。

图9 6种驱动模式

4.1 许可初值偏差

N-R法和Broyden法在计算协调方程时都需要给定算法迭代初值，如果初值选取不当，就会造成计算结果发散，失去有效性。并联机构结构参数与2.2节中的参数一致，分别给定算法迭代精度为1.0×10-6与1.0×10-9，仿真结果如图10和图11所示。本文以一组结构参数为例，给出了确定最佳冗余度的方法和思路，为类似并联机构驱动模式的选取提供参考。

图10 许可初值偏差与迭代步数的关系(1.0×10-6)

图11 许可初值偏差与迭代步数的关系(1.0×10-9)

观察并分析图10和图11可知：

(1)随着许可初值偏差的增大和迭代精度的提高，迭代步数呈上升趋势。R=3时，冗余协调算法对初值的依赖程度较低；R=1时，Broyden法对初值的依赖程度较高；R=5时，迭代步数都较低，即实时性较优。

(2)冗余度为0，1，3时，迭代步数对迭代精度的敏感程度都较高；特别地，对于N-R法，当R=0时，迭代步数对迭代精度较敏感，并且对许可初值偏差的要求较高，究其原因，这可能是由于Jacobian矩阵维数的增加导致计算量变大，以及更易产生奇异等原因造成的。

(3)随着许可初值偏差的增大，N-R法的迭代步数变化较平缓，即稳定性较优，而对于Broyden法，当许可初值偏差接近最大许可初值偏差时，迭代步数产生较大波动，从波动的整个趋势来看，呈上升态势。

4.2 结构参数

对于不同的驱动模式，冗余协调算法的最大许可初值偏差与并联机构结构参数有关。不失一般性，取5～30 mm分别作为动平台半边长和移动副初始长度的变化范围，均以5 mm作为变化步长，依次获取每种结构参数所对应的最大许可初值偏差，仿真结果如图12所示。其中黄色代表N-R法，蓝色代表Broyden法。

图12 结构参数与最大许可初值偏差的关系

Fig.12 Relation ship between structural parameters and maximum admissible initial deviation

观察并分析图12可知：

(1)从整体趋势来看，随着动平台边长的增大和移动副初始长度的减小，最大许可初值偏差呈上升趋势。

(2)R=5时，最大许可初值偏差的区间宽度都较小，分别为19%和0%，此时，最大许可初值偏差对结构参数变化的敏感程度较低；R=2时，最大许可初值偏差的区间宽度都较大，分别为1 218%和1 174%，此时，最大许可初值偏差对结构参数变化的敏感程度较高，即在该驱动模式下，通过合理选取结构参数，能够增大冗余协调算法的最大许可初值偏差。

(3)R=0时，最大许可初值偏差在移动副初始长度为15 mm左右出现凸峰，究其原因，这可能是由于算法在计算Jacobian矩阵时较易产生奇异所导致。

为了进一步定量地对比驱动模式对冗余协调算法产生的影响，运用区间分析理论，引入并定义扰动适应性能评价指标，将6种驱动模式下的扰动适应性能的平均值作为对应算法的综合扰动适应性能。扰动适应性能的表达式为：

(33)

显然，扰动适应性能P值越小，表明冗余协调算法的抗扰动性能越强。表6列出了不同驱动模式下冗余协调算法的扰动适应性能。

表6 冗余协调算法的扰动适应性能对比

Tab.6 Comparison of perturbation adaptive performance of redundant coordination algorithm

驱动模式(R)扰动适应性能P/1N-RBroydenⅠ(5)0.360Ⅱ(4)2.051.84Ⅲ(3)0.600.43Ⅳ(2)4.394.54Ⅴ(1)0.911.00Ⅵ(0)5.082.76综合扰动适应性能2.231.76

观察并分析表6可知：

(1)冗余度为1，3，5时，冗余协调算法的抗扰动性能较强，即最大许可初值偏差对结构参数变化的敏感程度较低，从移动副的布局方式来看，冗余度为奇数时，至少存在1个二重复合球铰链，其两端连接的移动副由驱动副和从动副混连构成，该构成形式能够实现冗余协调算法的扰动补偿。

(2)特别地，当R=5时，Broyden法的扰动适应性能为0，在该驱动模式下，结构参数的变化对最大许可初值偏差不产生影响。从综合扰动适应性能来看，Broyden法优于N-R法1.27倍，即Broyden法的抗扰动性能较优。

5 结 论

以具有解析式正解的Stewart衍生拓扑构型为研究对象，为了降低机构的耦合度，发明了一种二重复合球铰链，还给出了可转换主动、从动模式的移动副的一种设计方案，该移动副可以改变机构的冗余度，从而实现驱动模式的切换。构造了位置正解全解析算法，并验证了位置正解模型的正确性，该方法同样适用于动平台上含3个以上二重复合球面副、球铰中心不局限于动平台棱边中点且耦合度小于2的台体型并联机构位置正解的求解。

根据位置点之间的几何关系推导并化简了6个位移输入协调方程，分别对比分析了N-R法和Broyden法的数值性态，研究发现，N-R法的精度优于Broyden法至少3倍，Broyden法的计算时间约为N-R法的78%，Broyden法的综合扰动适应性能优于N-R法的1.27倍。同时，研究了许可初值偏差和结构参数对冗余协调算法的影响。

Stewart衍生型并联机构在具体应用场合需要合理选取驱动模式，其选取原则为：要允许有较大的最大许可初值偏差时，最佳冗余度为3；要使最大许可初值偏差对结构参数变化较敏感时，最佳冗余度为2；对实时性和抗扰动性能要求较高时，最佳冗余度为5。

本文仅讨论了N-R法和Broyden法的数值性态，对于其他算法是否会得到不同结论，该问题将另文研究。