0 引　言

少自由度并联机构由于其高刚度、高速度、高负载能力和无误差累积等特点，近年来在工业、航海、航空等领域得到广泛的应用^{[1 − 2]}。特别是在舰载平台上，船舶在海上作业时，由于波浪运动而带来的船体升沉、横摇和纵摇，会给过驳吊装作业带来风险和隐患^[3]。因此，必须通过并联机构对船舶在这3个自由度的扰动进行补偿，保证海上过驳作业的高效与稳定。

荷兰Barger Master公司基于3UPS-3SS并联机构研制代号为BM-T700的海上波浪补偿平台，依靠动平台的运动补偿横摇、纵摇和升沉3个自由度扰动，但该补偿平台存在体积大、力传递性能差等问题^[4]。YONG 等^[5]开发了一种8PSS-UP并联式舰船波浪补偿机构，对其进行正、逆运动学分析，并利用速度雅克比矩阵验证其工作空间内无奇异点，但其平均最小奇异值较小，机构整体力传递性能较差。强红宾等^[6]提出一种具有解析解的1T2R变异3UPS-PUU-2SS并联稳定机构，利用空间闭环矢量方程构建逆运动学模型与约束方程，求解参数化表达式及速度雅克比矩阵，并根据其局部条件数评估工作空间灵巧度，但该机构在工作空间内的全局灵巧度较差。张伟中等^[7]对1T2R并联机构提出一种尺度综合多目标优化设计方法，并采用多目标粒子群优化算法得到目标函数的Pareto最优解集。彭斯洋等^[8]提出一种非对称1T2R型UPS-RPU-PU并联机构，基于第二代非支配排序遗传算法（NSGA-II），以姿态能力、全域传递指标、传递波动指标为目标进行优化设计，最终得到该并联机构参数的Pareto最优解集。上述算法对双目标优化问题求解效果接近，但对三目标及以上则存在解集收敛性差、种群多样性低、耗时较长等问题。毕晓君等^[9]针对带约束的高维多目标优化问题设计一种基于参考点的约束支配的第三代非支配排序遗传算法（NSGA-III），相对于其他算法的解集具有更好的收敛性和分布性。

针对舰船载3UPS-3SS并联机构存在体积大、力传递性能差等问题，本文以3UPS-3SS并联机构为对象，通过空间闭环矢量方程法，建立机构的逆运动学模型和约束方程组，进而求解逆运动学的参数化表达式及速度雅克比矩阵；根据速度雅克比矩阵的局部条件数，求解该并联机构的力传递性能。在此基础上，以可达工作空间（下文简称工作空间）体积、机构体积和力传递性能为优化目标，建立带约束多目标优化模型，采用NSGA-III算法得到了在约束条件下较优的机构尺度参数，为该机构的控制和实际应用奠定理论基础。

1 结构描述

图1和图2所示分别为3UPS-3SS并联机构初始位置三维模型和机构简图，该机构由动平台、定平台、3条驱动支链和3条约束支链组成。定平台固定在船体上；3条驱动支链A_iB_i（i = 1,2,3）连接动平台与定平台，为UPS构型。驱动支链驱动元件为液压驱动单元（P副），其在定平台端连接件为虎克铰（U副），在动平台端连接件为球铰（S副），驱动杆A_iB_i长度为l_i，且A_iB_i初始长度为动、定平台初始高度差h。3个等长约束杆F_iD_i两端分别与动平台、定平台通过球铰（S副）连接，为SS构型的约束支链。铰点F_i位于动平台下方h₁处，即|C_iG_i| = h₁。

如图2所示，在定平台上建立定坐标系O₁x_Ly_Lz_L（记为{L}），x_L轴和y_L轴正向分别沿A₁A₂和O₁A₃方向，z_L轴垂直于定平台。在动平台上建立运动坐标系Ox_Py_Pz_P（记为{P}），x_P轴和y_P轴正向沿B₁B₂和OB₃方向，z_P轴垂直于动平台。

驱动杆与定平台连接铰点为A_i，与动平台连接铰点为B_i。如图3所示，铰点A_i相连组成等边三角形，外接于半径为R_a的圆内，圆心为O₁，三角形各边对应的圆心角为θ_a；铰点B_i相连组成等边三角形，外接于半径为R_b的圆内，圆心为O，三角形各边对应的圆心角为θ_b，其中θ_a = θ_b。该机构的特征在于，当其处于初始位姿时，动平台平行于定平台，定坐标系原点O₁是动平台中心O在定平台上的投影，驱动杆轴线垂直于定平台，3个约束杆在一平行于定平台的平面内。如图3所示，B、F、D在一条直线上，且|B_iF_i| = l_n，|F_iD_i| = l_m。

2 运动学分析

首先，结合关节坐标与旋转矩阵，构建逆运动学模型；然后，通过约束支链关系建立方程组，求得位置逆解；最后，基于动平台与驱动杆速度关系，解出速度雅克比矩阵表达式。

2.1 位置逆解

通过Roll-Pitch-Yaw(RPY)角定义变换矩阵R，描述动坐标系{P}对定坐标系{L}的转动：先绕{L}的x_L轴转γ，再绕y_L轴转β，最后绕z_L轴转α，则旋转矩阵为：

$ \begin{split} &{{\boldsymbol{ R}} = {\boldsymbol{R}}(z,\alpha ){\boldsymbol{R}}(y,\beta ){\boldsymbol{R}}(x,\gamma ) = }\\ &{\left( {\begin{array}{*{20}{c}} {\cos \alpha \cos \beta }&{\cos \alpha \sin \beta \sin \gamma - \sin \alpha \cos \gamma }&{\cos \alpha \sin \beta \cos \gamma + \sin \alpha \sin \gamma } \\ {\sin \alpha \sin \beta }&{\sin \alpha \sin \beta \sin \gamma + \cos \alpha \cos \gamma }&{\sin \alpha \sin \beta \cos \gamma - \cos \alpha \sin \gamma } \\ { - \sin \beta }&{\cos \beta \sin \gamma }&{\cos \beta \cos \gamma } \end{array}} \right)} \end{split} 。$

(1)

定平台U副中心点A_i在{L}下的位置矢量为：

$ {}^L{{A}_i} = {\left( {\begin{array}{*{20}{c}} {{A_{ix}}}&{{A_{iy}}}&{{A_{iz}}} \end{array}} \right)^{\text{T}}} 。$

(2)

动平台S副球铰中心位置B_i在{P}下的位置矢量为：

$ ^PB_i=\left(\begin{array}{*{20}{c}}B_{ix} & B_{iy} & B_{iz}\end{array}\right)^{\text{T}}。$

(3)

动坐标系原点O在{P}下的位置矢量在{L}下为：

$ ^LO_P=\left(\begin{array}{*{20}{c}}^LO_{Px} & ^LO_{Py} & ^LO_{Pz}\end{array}\right)^{\text{T}}。$

(4)

则铰点B_i在{L}下的位置矢量为：

$ {}^L{{B}_i} = {}^L{{O}_P} + {R}{}^P{{B}_i}。$

(5)

约束支链在动平台上的球铰铰点F_i在{P}下的位置矢量为：

$ ^PF_i=\left(\begin{array}{*{20}{c}}F_{ix} & F_{iy} & F_{iz}\end{array}\right)^{\text{T}}。$

(6)

则F_i在{L}下的位置矢量为：

$ ^LF_i=^LO_P+R^PF_i。$

(7)

约束支链在定平台上的球铰铰点D_i在{L}下的位置矢量为：

$ ^LD_i=\left(\begin{array}{*{20}{c}}D_{ix} & D_{iy} & D_{iz}\end{array}\right)^{\text{T}}。$

(8)

则3个约束杆的约束方程为：

$ \left\{\begin{aligned} & \left(D_{1x}-F_{1x}\right)^2+\left(D_{1y}-F_{1y}\right)^2+\left(D_{1z}-F_{1z}\right)^2=l_m^2，\\ & \left(D_{2x}-F_{2x}\right)^2+\left(D_{2y}-F_{2y}\right)^2+\left(D_{2z}-F_{2z}\right)^2=l_m^2，\\ & \left(D_{3x}-F_{3x}\right)^2+\left(D_{3y}-F_{3y}\right)^2+\left(D_{3z}-F_{3z}\right)^2=l_m^2。\\ \end{aligned}\right. $

(9)

机构第i条驱动杆长度矢量为l_i，其在定坐标系{L}下可由第i条驱动支链的球铰中心位置^LB_i和虎克铰中心位置^LA_i计算得到，如下式：

$ l_i=^L {B_i}-^L {A_i}。$

(10)

联立式(5)与式(10)得：

$ l_i=^L {O_P}+R^P {B_i}-^L {A_i}。$

(11)

则机构驱动杆长度为：

$ \left|l_i\right|=\sqrt{l_1^2+l_2^2+l_3^2}。$

(12)

式中：l₁、l₂、l₃分别为驱动杆在x、y、z方向上的实际长度。

机构驱动杆理论长度L_i为：

$ L_i=\sqrt{(B_i-A_i)^{\text{T}}(B_i-A_i)}。$

(13)

使用最小二乘法，建立目标函数求解驱动杆长度：

$ F(x,y,z,\theta_x,\theta_y,\theta_z)=\sum\limits_{i=1}^3[l_i-L_i]^2。$

(14)

式中：(x, y, z, θ_x, θ_y, θ_z)为驱动杆在定坐标系{L}下的位姿；θ_x、θ_y、θ_z分别为驱动杆绕x_L、y_L、z_L轴的旋转角度。

2.2 速度逆解

根据动平台的位姿，求得动平台的速度v_p和角速度ω_p，进而求解3个驱动杆的伸缩速度。

B_i的速度矢量，一方面可以通过动平台的速度和角速度求得，如下式：

$ v_{B_i}=\omega_p\times\left(R^PB_i\right)+v_p。$

(15)

另一方面可以通过驱动杆的速度$ \dot{{{L}}_{{i}}} $和角速度ω_li求得，如下式：

$ {{v}_{{B_i}}} = {{\omega }_{{l_i}}} \times {{l}_i} + {\dot L_i}{{n}_i} 。$

(16)

在式(15)、式(16)两边同时点乘驱动杆轴向单位向量n_i，化简可得驱动杆的伸缩速度为：

$ {\dot L_i} = \left[ {\begin{array}{*{20}{c}} {{{n}_i}}&{\left( {{{R}^P}{{B}_i}} \right) \times {{n}_i}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{v}_p}} \\ {{{\omega }_p}} \end{array}} \right]。$

(17)

机构速度雅克比矩阵为：

$ {J = }\left[ {\begin{array}{*{20}{l}} {{{n}_1}}&{\left( {{{R}^P}{{B}_1}} \right) \times {{n}_1}} \\ {{{n}_2}}&{\left( {{{R}^P}{{B}_2}} \right) \times {{n}_2}} \\ {{{n}_3}}&{\left( {{{R}^P}{{B}_3}} \right) \times {{n}_3}} \end{array}} \right] 。$

(18)

将3条驱动杆的伸缩速度写为向量形式，可得各驱动杆伸缩速度、动平台速度与角速度的关系为：

$ \dot{L}_i={J}\cdot\dot{T} 。$

(19)

式中：$ \dot{{T}} $为动平台的速度和角速度，$ \dot{{T}}=\left[\begin{array}{c}{v}_{{p}} \\ {w}_{{p}}\end{array}\right] $。

3 性能分析 3.1 工作空间分析

工作空间是评价并联机构运动学性能的重要指标^[10]。设置工作空间的快速极坐标搜索步骤，从起始角开始，极径由0 mm增加，当伸缩杆长度l、转动关节的转动角度达到极限值时，极径的端点即为工作空间的边界点。国内近海以4级海况为主，在该海况下舰船升沉、横摇和纵摇运动的稳定补偿^[11]选取结构参数如表1所示。

设置工作空间的约束条件为：

$ \left\{\begin{aligned} & l_{\min}=l_d+l_s，\\ & l_{\max}=l_d+l_s+l_s，\\ & \theta_{ai}=\theta_{a\max}，\\ & \theta_{bj}=\theta_{b\max}。\\ \end{aligned}\right. $

(20)

式中：l_max、l_min分别为驱动杆最大、最小长度；θ_ai为虎克铰转角大小；θ_bj(j = 1, 2, …, 9)为球铰转角大小。

使用Matlab软件通过快速极坐标法绘制工作空间如图4所示。可知，该机构的工作空间较大，并非为规则形状，不可对其体积进行显式表达。本文采用包围盒法^[12]将工作空间包围成规则形状后对其体积V_W进行计算，计算式为：

$ V_W=\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)\cdot\left(y_{\mathrm{max}}-y_{\mathrm{min}}\right)\cdot\left(z_{\mathrm{max}}-z_{\mathrm{min}}\right)。$

(21)

式中：x_max、y_max、z_max分别为工作空间在x、y、z轴方向的最大坐标值；x_min、y_min、z_min分别为工作空间在x、y、z轴方向的最小坐标值。

3.2 机构物理尺寸

当并联机构能适应舰船运动范围时，缩小机构尺寸能降低成本，便于安装。一般选择最小机构体积作为优化参数^[13]。本文定义并联机构初始位姿时的体积为其物理尺寸。机构尺寸示意如图5所示。

初始位姿时，并联机构为一个三棱柱与长方体的组合体，其动平台、定平台圆外接圆半径分别为R_b、R_a，高度为h。机构体积表达式为：

$ V=\frac{R_a}{2}\left(1+\sin\frac{\theta_a}{4}\right)\left(3R_a\cos\frac{\theta_a}{4}+l_m+l_n\right)h。$

(22)

式中：θ_a = 120°。

3.3 力传递性分析

力雅克比矩阵描述的是驱动杆受力与动平台受力之间的关系，如下式：

$ {\boldsymbol{F}}_{F}={H}\cdot{\tau}_{j} 。$

(23)

式中：F_F为动平台所受到的力和力矩；τ_j为驱动杆的驱动力；H为力雅克比矩阵。

在机构学中，力雅克比矩阵H为速度雅克比矩阵J的转置J^T，则式(23)可变换为：

$ \boldsymbol{\tau}_j=\boldsymbol{H}^{-1}\cdot{F}_F=\left(\boldsymbol{J}^{\text{T}}\right)^{-1}\cdot F_F。$

(24)

对式(24)中的速度雅克比矩阵J进行奇异值分解，可得到驱动杆的驱动力τ_j的取值范围为：

$ \frac{ \left\| F_F \right\| }{\sigma_{\max}\left(J\right)}\leqslant \left\| \tau_j \right\| \leqslant\frac{ \left\| F_F \right\| }{\sigma_{\min}\left(J\right)}。$

(25)

驱动杆驱动力上限由速度雅克比矩阵J的最小奇异值σ_min(J)决定。σ_min(J)减少时，上限增大，机构接近奇异位姿，力传递效率下降；反之，上限减小，机构趋近各向同性，力传递性能提升。因此，σ_min(J)可作为局部力传递性能指标（LTI）。

LTI只反映机构在特定位姿下的力传递性能，无法对机构整体力传递性能进行评价。为弥补LTI指标的不足，引入全域力传递指标（GTI）。从传递角的定义出发，将工作空间局部力传递性能指标大于0.7的区域定为优质传递空间。优质传递空间与整体工作空间的比值为GTI，其表达式为：

$ {\psi _{{\text{GTI}}}} = \frac{{\displaystyle\int\limits_{{V_G}} {{\text{d}}W} }}{{\displaystyle\int\limits_{{V_W}} {{\text{d}}W} }}。$

(26)

式中：W为机构的整体工作空间；V_G为优质传递空间的体积；V_W为整体工作空间的体积。

ψ_GTI取值在0～1之间，越趋近于1，表示机构的力传递性能越好；越趋近于0，机构越接近奇异状态，力从驱动端至末端执行器的传递效率降低。

根据以上定义和3.1节参数，搜索得到雅可比矩阵最小奇异值的分布图，如图6所示。可知，机构在目标工作空间内绝大部分的区域的σ_min(J) > 0.3，平均最小奇异值$ {\bar{{\sigma}}}_{\text{min}}\left({J}\right) $为0.28，对应ψ_GTI为0.46，说明该机构整体力传递性能较差，尤其工作空间边界区域σ_min(J)在0～0.1之间，接近奇异位姿。

4 机构尺度参数优化 4.1 建立优化模型

机构整体性能受机构尺度参数影响。在该机构尺度优化问题中，尺度参数均为长度变量。该优化设计问题的决策变量为：

$ x=\left(x_1,x_2,\cdot\cdot\cdot,x_D\right)^{\text{T}}=\left(R_a,R_b,h_1,l_m,l_n,l_d,l_s\right)^{\text{T}}。$

(27)

式中：D为决策变量维数。

为提升机构在工作空间内的整体力传递性能，使得机构在尽可能多的区域内都具备较高的力传递效率，将平均最小奇异值作为优化目标。在机构设计优化过程中，若机构尺度参数足够大，即机构物理体积足够大，其工作空间就越大，同时可能降低系统的刚性和力传递效率。因此，最大化工作空间体积、最小化机构物理尺寸和最大化机构平均最小奇异值为一组互相矛盾的目标，可作为尺度参数优化时的3个目标。该优化问题的目标函数可设置为：

$ \min \left\{ \begin{aligned} &{f_1}\left( {x} \right) = 1/{V_W} ，\\ &{f_2}\left( {x} \right) = V ，\\ &{f_3}\left( {x} \right) = 1/{{\bar \sigma }_{\min }}\left( {J} \right) 。\\ \end{aligned} \right. $

(28)

为保证机构能适应舰船运动范围、避免干涉、避免机构奇异以及优化力传递性能，引入以下约束条件：

1）决策变量取值范围限制。为保证机构能适应舰船运动范围，将决策变量取值限制在一定范围内，如表2所示。

2）虎克铰、球铰转角限制。为避免干涉，将虎克铰、球铰转角限制在一定范围内，整理得：

$ \left\{ \begin{aligned} &{g_j}\left( x \right) = {\theta _{bj}} - {\theta _{b\max }} \leqslant 0 ，\\ &{g_{i + 9}}\left( x \right) = {\theta _{ai}} - {\theta _{a\max }} \leqslant 0 。\\ \end{aligned} \right. $

(29)

3）避免奇异位形，以保证机构具有较优的传递性能。因此，引入约束条件为：

$ {g_{13}}\left( x \right) = {\sigma _{{\text{min,min}}}} = \min {\sigma _{{\text{min}}}}\left( {J} \right) \geqslant 0.2。$

(30)

式中：σ_min,min为工作空间中最小奇异值中的最小值。

4.2 NSGA-III多目标优化方法

并联机构尺度参数优化是一个多目标优化问题，DEB等^[14]提出的NSGA-III算法是一种多目标优化算法，其原理如图7所示。NSGA-III算法摒弃了NSGA-II算法中基于拥挤距离的选择机制，引入参考点的概念，有效降低计算代价的同时，改善算法的收敛性和种群的多样性，加快搜索解决方案的速度，适合解决高维多目标（3个及以上）优化问题。NSGA-III算法流程如图8所示。

NSGA-III算法包括以下基本步骤：

步骤1　迭代次数、种群规模、变异概率、交叉概率等相关参数设置，参考点与种群初始化；

步骤2　对父代种群P_t进行选择、交叉和变异等操作得到子代Q_t并分别计算Q_t中所有个体的适应度；

步骤3　合并父代子代种群得到新种群R_t并通过非支配排序将其划分为等级不同的非支配层，如F₁，F₂，…，F_l，…，F_n；

步骤4　从F₁非支配层开始，依次往下选择个体保留至下代种群S_t中直到S_t中个体数量大于等于种群规模N，若此时被选择个体所在层为F_l，则舍弃F_l后剩余的非支配层级中的个体。若此时S_t中个体数量等于N，则此次迭代结束，并将S_t作为下一次迭代的父代种群P_t+1；若此时S_t中个体数量大于N，则将非支配层F_l之前的所有个体存放到下一次迭代的父代种群P_t+1中，然后基于参考点的种群选择策略，即根据最短距离分配。选取F_l中部分个体存放到P_t+1中直至P_t+1的种群规模大小为N，同时结束此次迭代。

步骤5　判断是否达到所设定的迭代次数，若是则终止迭代，输出最优种群，否则回到步骤2。

4.3 多目标优化及结果分析

NSGA-III算法参数设定为：种群规模200，交叉概率0.5，变异概率0.5，迭代次数80。搜索得到200个Pareto最优解，如图9所示。

本文为三目标优化，采用TOPSIS法^[15]通过评价对象与理想目标的接近程度进行排序，确定各目标对总体的影响。其评价指标计算方法为：

$ C=\frac{D^-}{D^-+D^+}。$

(31)

式中：D⁺为解到最优解的距离；D⁻为解到最劣解的距离。

根据以上评价指标计算，可知在综合考虑物理体积、工作空间体积和平均最小奇异值样本16是最优解，对应的机构尺度参数为{R_a，R_b，h₁，l_m，l_n，l_s，l_d}={1201，1388，259，1264，346，387，190}。

根据以上参数，搜索得到优化后雅可比矩阵最小奇异值的分布如图10所示。可知，机构在目标工作空间内绝大部分的区域的σ_min(J) > 0.8，平均最小奇异值$ {\bar{\sigma}}_{\min}\left({J}\right) $为0.86，对应ψ_GTI为1，说明该机构整体力传递性能较好。

机构优化前后性能参数对比如表3所示。

5 结　语

1）对舰船载3UPS-3SS并联机构进行运动学分析，推导出运动学反解和速度表达式，并求得速度雅克比矩阵。针对该机构体积大、力传递性能差等问题，使用Matlab软件采用快速极坐标搜索法绘制机构工作空间，并对机构物理体积、工作空间体积、力传递性能进行分析。

2）采用NSGA-III算法对机构物理体积、工作空间体积、平均最小奇异值进行优化，在工作空间体积损失率仅为14.01%的情况下，物理体积优化了17.58%，平均最小奇异值优化了207.14%，提高了经济性的同时提升了机构的力传递性能，有效保证了海上过驳作业的高效与稳定。