随着机器人技术的迅速发展和社会需求的不断提高，机器人与人共融、协同作业，将成为下一代机器人的本质特征，也是当前国内外学者研究的热点问题。然而，在人机协同作业过程中，如何确保人与机器的安全是共融机器人技术研发中亟待解决的关键和难点问题之一。

索驱动并联机器人是一种特殊的并联机器人，其动平台是由绳索代替刚性杆进行驱动^[1-3]。索并联机器人因具有优良的运动性能、较高的负载能力以及较大的工作空间等优点，受到了学者的青睐^[4-6]。此外，索并联机器人可以通过调节绳索的索力来改变系统的刚度^[7]，这使研制刚度自适应调整的变刚度机器人成为可能，也为解决共融机器人在“人机交互”过程中的人机安全问题提供了强有力的技术支撑。

刚度问题是改善和提高机器人性能极为重要的方面。文[8-11]通过研究机器人末端刚度与姿态之间的关系，得出了其最优刚度姿态和增加刚度的方法。杜敬利等^[12-13]在充分考虑悬索因自重而产生垂度以及存在弹性变形等因素的条件下，对机器人的刚度进行了有效分析。Yeo等^[2]提出一种新型的变刚度索驱动机器人，通过在其驱动绳索上安装变刚度设备以大幅度改变机器人的刚度。Wen等^[14]提出一种基于梯度投影法的索力优化分配算法，在避免绳索虚牵的同时对机器人的刚度进行有效调整。Wang等^[15]针对八索驱动的6自由度悬吊系统，提出一种基于刚度优化的力位混合控制方法应用于飞机风洞实验测试系统。上述研究主要是为了保证索机器人的位置精度，从提高其刚度的角度进行相关研究，并没有研究其刚度的精确控制。而从人机交互的角度出发，需要对机器人的可控刚度特性进行充分研究，并对其刚度进行精确控制。

为此，本文根据八索并联机器人特点，对其运动学和静力学进行分析，推导其可控刚度和固有刚度矩阵，并建立静刚度模型，研究其可控刚度特性问题并给出索力多边形的计算方法，在索力可行域内对索力进行分配。该方法可以在系统满足静态平衡和绳索不出现虚牵的情况下，对其索力进行调整，实现对系统刚度的有效控制。最后，通过实例计算验证了本方法的正确性和有效性。

1 静态分析 1.1 机构描述

图 1a是拟研究的八索驱动并联机器人实验平台。对其进行机构建模，如图 1b所示，动平台在8根绳索的冗余驱动下可实现空间6自由度运动，其结构为4-4构型，即8根绳索以上四、下四的对称方式进行布局，出索点位于固定框架的8个顶点位置，连接点分别位于动平台的上、下表面4个定点位置。其中，A_i(i=1, 2, …, 8)表示绳索出索点，B_i(i=1, 2, …, 8)表示绳索与动平台的连接点，其尺寸参数如表 1所示。

A_i和B_i(i=1, 2, …, 8)形成的长方体尺寸分别为1.580 m×1.040 m×1.780 m和0.124$\sqrt 2 $ m×0.124$\sqrt 2 $ m×0.085 m。

根据绳索悬链线模型^[16]，绳索的实际长度s可表示为参数λ的函数，且当λ趋近于0时，绳索实际长度为直线长度，其中λ=q₀/f_h，q₀为绳索的单位质量，f_h为绳索拉力f的水平分量。由于本文所研究机器人的长宽高均小于2 m，索跨度小，单位重量q₀较小，λ趋近于0，故可忽略索因自重而产生的变形，用直线模型代替绳索悬链线模型^[17]。

1.2 位置与静力学分析

八索驱动并联机器人位置分析的目的是建立其动平台位姿和绳索长度之间的关系。如图 1b所示，建立八索并联机器人的固定坐标系O-xyz和运动坐标系O′-x′y′z′，并分别用O和O′表示，其中，O建立在框架底面中心处，其z轴垂直于底面，x//A₁A₂；O′建立在动平台上表面中心处，其z′轴垂直于动平台上表面，x′//B₁B₂。因此，动平台的位置可用O′的原点在O下的位置矢量p=[x_o′, y_o′, z_o′]^T表示，姿态可用O′相对于O的旋转Euler角[α, β, γ]^T表示，故动平台的位姿在O下可表示为

$ \mathit{\boldsymbol{x}} = {\left[ {{x_{o'}},{y_{o'}},{z_{o'}},\alpha ,\beta ,\gamma } \right]^{\rm{T}}}. $

(1)

根据矢量环方程，建立绳索长度向量与动平台位姿之间的关系如下：

$ {\mathit{\boldsymbol{l}}_i} = {\mathit{\boldsymbol{a}}_i} - \mathit{\boldsymbol{p}} - \mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}, $

(2)

其中：l_i表示从B_i到A_i绳索矢量；a_i表示A_i在O下的位置矢量；^o′b_i表示B_i在O′下的位置矢量；R表示O′相对于O的旋转矩阵，i=1, 2, …, 8。

其中：

$ \mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {{\rm{c}}\gamma \cdot {\rm{c}}\beta }&{{\rm{c}}\gamma \cdot {\rm{s}}\beta \cdot {\rm{s}}\alpha - {\rm{s}}\gamma \cdot {\rm{c}}\alpha }&{{\rm{c}}\gamma \cdot {\rm{s}}\beta \cdot {\rm{c}}\alpha + {\rm{s}}\gamma \cdot {\rm{s}}\alpha }\\ {{\rm{s}}\gamma \cdot {\rm{c}}\beta }&{{\rm{s}}\gamma \cdot {\rm{s}}\beta \cdot {\rm{s}}\alpha + {\rm{c}}\gamma \cdot {\rm{c}}\alpha }&{{\rm{s}}\gamma \cdot {\rm{s}}\beta \cdot {\rm{c}}\alpha - {\rm{c}}\gamma \cdot {\rm{s}}\alpha }\\ { - {\rm{s}}\beta }&{{\rm{c}}\beta \cdot {\rm{s}}\alpha }&{{\rm{c}}\beta \cdot {\rm{c}}\alpha } \end{array}} \right]. $

为了便于描述，本文采用cθ代表cosθ，sθ代表sinθ。

因此，当动平台位姿x确定，绳索长度‖l_i‖可通过式(2)计算获得。进一步可得绳索单位矢量：

$ {\mathit{\boldsymbol{u}}_i} = {\mathit{\boldsymbol{l}}_i}/\left\| {{\mathit{\boldsymbol{l}}_i}} \right\|. $

(3)

式(2)进一步可得机构运动学方程：

$ \mathit{\boldsymbol{x}} = f\left( \mathit{\boldsymbol{l}} \right). $

(4)

其中f是非线性矢量方程。式(4)对时间求导得：

$ \mathit{\boldsymbol{\dot I = }} - {\mathit{\boldsymbol{J}}^{\rm{T}}} \cdot \mathit{\boldsymbol{\dot x}}. $

(5)

其中：İ=[İ₁, İ₂, …, İ₈]^T表示绳索速度矢量，J∈R^6×8表示结构Jacobian矩阵，$\mathit{\boldsymbol{\dot x}} = {\left[{\mathit{\boldsymbol{\dot p}}, \mathit{\boldsymbol{\omega }}} \right]^{\rm{T}}} = $${\left[{{{\dot x}_{o'}}, {{\dot y}_{o'}}, {{\dot z}_{o'}}, \dot \alpha, \dot \beta, \dot \gamma } \right]^{\rm{T}}}$表示动平台的速度矢量。

其中：

$ \mathit{\boldsymbol{J}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{u}}_1}}& \cdots &{{\mathit{\boldsymbol{u}}_8}}\\ {\left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_1}} \right) \times {\mathit{\boldsymbol{u}}_1}}& \cdots &{\left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_8}} \right) \times {\mathit{\boldsymbol{u}}_8}} \end{array}} \right]. $

根据虚功原理可得八索机构静力学平衡方程：

$ \mathit{\boldsymbol{JT}} + \mathit{\boldsymbol{F}} = {\bf{0}}, $

(6)

其中：T=[t₁, t₂, …, t₈]^T，t_i表示第i根绳索张紧力大小，F=[f_e, m_e]^T，且f_e和m_e分别表示外界作用在动平台上的合力和合力矩(包括动平台自身重力)。

2 静态刚度模型

研究机器人的刚度，是研究变刚度机器人的基础，同时也是实现人机共融的关键。动平台在外扰力作用下的抗变形能力，是衡量系统刚度的标准之一。因此，建立动平台受到微小外扰力δF和动平台产生微小位姿δx之间的关系：

$ \delta \mathit{\boldsymbol{F}} = \mathit{\boldsymbol{K}} \cdot \delta \mathit{\boldsymbol{x}}, $

(7)

其中K表示系统刚度矩阵。

由式(6)和(7)可得：

$ \mathit{\boldsymbol{K}} = \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial \mathit{\boldsymbol{x}}}} = - \left( {\frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial \mathit{\boldsymbol{x}}}}\mathit{\boldsymbol{T}} + \mathit{\boldsymbol{J}}\frac{{\partial \mathit{\boldsymbol{T}}}}{{\partial \mathit{\boldsymbol{x}}}}} \right) = {\mathit{\boldsymbol{K}}_1} + {\mathit{\boldsymbol{K}}_2}. $

(8)

由式(8)可以看出，系统刚度矩阵由2部分组成：K₁与系统的结构矩阵变换和索力相关，在某位姿处可通过索力对其进行控制，称为可控刚度，也是本文主要研究内容；K₂与系统结构和动平台位姿相关，称为固有刚度。

2.1 可控刚度

可控刚度是实现机器人刚度自适应调整的关键，本文将通过引入线矢量和微分变换的方式对其矩阵进行推导，由式(8)进一步可得：

$ {\mathit{\boldsymbol{K}}_1} = - \mathit{\boldsymbol{HT}}. $

(9)

其中：H=[H¹, H², …, H⁶]∈R^6×6×8表示三维Hessian矩阵，H^j∈R^6×8是矩阵H的第j个子阵。

利用矩阵J对向量x的第j个元素求偏导，得到

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}^j} = \frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = }\\ {\left[ {\begin{array}{*{20}{c}} {\frac{{\partial {\mathit{\boldsymbol{u}}_1}}}{{\partial {\mathit{\boldsymbol{x}}_j}}}}& \cdots &{\frac{{\partial {\mathit{\boldsymbol{u}}_8}}}{{\partial {\mathit{\boldsymbol{x}}_j}}}}\\ {\frac{{\partial \left[ {\left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_1}} \right) \times {\mathit{\boldsymbol{u}}_1}} \right]}}{{\partial {\mathit{\boldsymbol{x}}_j}}}}& \cdots &{\frac{{\partial \left[ {\left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_8}} \right) \times {\mathit{\boldsymbol{u}}_8}} \right]}}{{\partial {\mathit{\boldsymbol{x}}_j}}}} \end{array}} \right].} \end{array} $

(10)

为简化矩阵H推导过程，将绳索单位向量u_i写成方向余弦形式

$ {\mathit{\boldsymbol{u}}_i} = {\left[ {c{\alpha _i},c{\beta _i},c{\gamma _i}} \right]^{\rm{T}}}. $

(11)

其中α_i、β_i、γ_i分别表示u_i与x轴、y轴和z轴夹角。

由式(10)和(11)可得：

$ \frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = \frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\alpha _i}}}\frac{{\partial {\alpha _i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} + \frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\beta _i}}}\frac{{\partial {\beta _i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} + \frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\gamma _i}}}\frac{{\partial {\gamma _i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}}, $

(12)

$ \begin{array}{*{20}{c}} {\frac{{\partial \left[ {\left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right) \times {\mathit{\boldsymbol{u}}_i}} \right]}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = \frac{{\partial \left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right)}}{{\partial {\mathit{\boldsymbol{x}}_j}}} \times }\\ {{\mathit{\boldsymbol{u}}_i} + \left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right) \times \frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}}.} \end{array} $

(13)

由于u_i是单位向量，式(12)进一步可得：

$ \begin{array}{*{20}{c}} {\frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = {{\left[ { - {\rm{s}}{\alpha _i},0,\frac{{{\rm{c}}{\alpha _i} \cdot {\rm{s}}{\alpha _i}}}{{{\rm{c}}{\gamma _i}}}} \right]}^{\rm{T}}} \cdot \frac{{\partial {\alpha _i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} + }\\ {{{\left[ {0, - {\rm{s}}{\beta _i},\frac{{{\rm{c}}{\beta _i} \cdot {\rm{s}}{\beta _i}}}{{{\rm{c}}{\gamma _i}}}} \right]}^{\rm{T}}} \cdot \frac{{\partial {\beta _i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}}.} \end{array} $

(14)

由式(2)和(3)得：

$ {l_i} \cdot {\mathit{\boldsymbol{u}}_i} = {\mathit{\boldsymbol{a}}_i} - \mathit{\boldsymbol{p}} - \mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}. $

(15)

由于R′=ω×R，式(15)两边对时间求导得：

$ \begin{array}{*{20}{c}} {{{\dot I}_i} \cdot {\mathit{\boldsymbol{u}}_i} + {l_i} \cdot {{\mathit{\boldsymbol{\dot u}}}_i} = - \mathit{\boldsymbol{\dot p}} - \omega \times \left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right) = }\\ {\left[ { - \mathit{\boldsymbol{I}}\left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right) \times } \right]\mathit{\boldsymbol{\dot x}} = {\mathit{\boldsymbol{G}}_i}\mathit{\boldsymbol{\dot x}}.} \end{array} $

(16)

其中G_i∈R^3×6表示代数矩阵。

由式(5)得：

$ {{\mathit{\boldsymbol{\dot I}}}_i} = - \mathit{\boldsymbol{J}}_i^{\rm{T}}\mathit{\boldsymbol{\dot x}}. $

(17)

其中J_i^T表示矩阵J^T第i行向量。

将式(17)写成向量分量形式，可得：

$ {{\dot I}_i}{\mathit{\boldsymbol{u}}_i} = - \left[ {\begin{array}{*{20}{c}} {{\rm{c}}{\alpha _i} \cdot \mathit{\boldsymbol{J}}_i^{\rm{T}}}\\ {{\rm{c}}{\beta _i} \cdot \mathit{\boldsymbol{J}}_i^{\rm{T}}}\\ {{\rm{c}}{\gamma _i} \cdot \mathit{\boldsymbol{J}}_i^{\rm{T}}} \end{array}} \right]\mathit{\boldsymbol{\dot x}} = {\mathit{\boldsymbol{Q}}_i}\mathit{\boldsymbol{\dot x}}. $

(18)

其中Q_i∈R^3×6表示代数矩阵。

将式(18)代入式(16)得：

$ \begin{array}{*{20}{c}} {{l_i} \cdot {{\left[ { - {{\dot \alpha }_i} \cdot {\rm{s}}{\alpha _i}, - {{\dot \beta }_i} \cdot {\rm{s}}{\beta _i}, - {{\dot \gamma }_i} \cdot {\rm{s}}{\gamma _i}} \right]}^{\rm{T}}} = }\\ {\left[ {{\mathit{\boldsymbol{G}}_i} - {\mathit{\boldsymbol{Q}}_i}} \right]\mathit{\boldsymbol{\dot x}}.} \end{array} $

(19)

其中[G_i-Q_i]∈R^3×6表示代数矩阵。

由式(19)可得：

$ \left\{ \begin{array}{l} \frac{{\partial {\alpha _i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = - \frac{1}{{{l_i} \cdot {\rm{s}}{\alpha _i}}}{\left[ {{\mathit{\boldsymbol{G}}_i} - {\mathit{\boldsymbol{Q}}_i}} \right]_{1,j}}\\ \frac{{\partial {\beta _i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = - \frac{1}{{{l_i} \cdot {\rm{s}}{\beta _i}}}{\left[ {{\mathit{\boldsymbol{G}}_i} - {\mathit{\boldsymbol{Q}}_i}} \right]_{2,j}} \end{array} \right.. $

(20)

其中[G_i-Q_i]_{1, j}和[G_i-Q_i]_{2, j}分别表示矩阵[G_i-Q_i]第1行第j列和第2行第j列元素。

R·^o′b_i对时间求导得：

$ {\left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right)^\prime } = \mathit{\boldsymbol{\omega }} \times \left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right) = {\mathit{\boldsymbol{M}}_i}\mathit{\boldsymbol{\dot x}}. $

(21)

其中M_i∈R^3×6表示代数矩阵。

由式(21)得：

$ \frac{{\partial \left( {\mathit{\boldsymbol{R}}{ \cdot ^{o'}}{\mathit{\boldsymbol{b}}_i}} \right)}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = {\mathit{\boldsymbol{M}}_{i,j}}. $

(22)

其中M_{i, j}表示矩阵M_i第j列向量。

将式(14)、(20)和(22)代入式(12)、(13)和(10)可得矩阵H，根据式(9)和T可计算刚度矩阵K₁。

2.2 固有刚度

固有刚度虽然不受索力控制，但对系统整体刚度依然有较大影响，本小节将给出其矩阵的推导过程，由式(8)进一步可得：

$ {\mathit{\boldsymbol{K}}_2} = - \mathit{\boldsymbol{J}}\frac{{\partial \mathit{\boldsymbol{T}}}}{{\partial \mathit{\boldsymbol{l}}}}\frac{{\partial \mathit{\boldsymbol{l}}}}{{\partial \mathit{\boldsymbol{x}}}}. $

(23)

绳索的索力与变形之间的关系：

$ \partial {t_i} = \frac{{{E_i} \cdot {A_i}}}{{{l_{oi}}}} \cdot \partial {l_i}, $

(24)

其中：E_i、A_i和l_oi分别表示绳索的弹性模量，横截面积和静态长度。

由式(24)得：

$ \frac{{\partial \mathit{\boldsymbol{T}}}}{{\partial \mathit{\boldsymbol{l}}}} = {\rm{diag}}\left( {\frac{{{E_1} \cdot {A_1}}}{{{l_{o1}}}},\frac{{{E_2} \cdot {A_2}}}{{{l_{o2}}}}, \cdots ,\frac{{{E_8} \cdot {A_8}}}{{{l_{o8}}}}} \right). $

(25)

由式(5)得：

$ \frac{{\partial \mathit{\boldsymbol{l}}}}{{\partial \mathit{\boldsymbol{x}}}} = - {\mathit{\boldsymbol{J}}^{\rm{T}}}. $

(26)

将式(25)和(26)代入(23)得：

$ {\mathit{\boldsymbol{K}}_2} = \mathit{\boldsymbol{J}} \cdot {\rm{diag}}\left( {\frac{{{E_1} \cdot {A_1}}}{{{l_{o1}}}},\frac{{{E_2} \cdot {A_2}}}{{{l_{o2}}}}, \cdots ,\frac{{{E_8} \cdot {A_8}}}{{{l_{o8}}}}} \right) \cdot {\mathit{\boldsymbol{J}}^{\rm{T}}}. $

(27)

利用上述公式可分别求出K₁和K₂，进而求出系统刚度矩阵K。

2.3 索力分配

已知动平台所受外力，如何在满足系统平衡条件下，有效调整各绳索索力是索力分配研究的重要内容，也是控制系统刚度的基础。由节2.1可知，通过改变绳索索力可以调整八索并联机器人系统刚度。因此，本文提出并采用以下方法实现对索力的分配。

由式(6)得：

$ \mathit{\boldsymbol{T}} = {\mathit{\boldsymbol{J}}^ + }\left( { - \mathit{\boldsymbol{F}}} \right) + \mathit{\boldsymbol{N}}\lambda = {\mathit{\boldsymbol{t}}_p} + {\mathit{\boldsymbol{t}}_n}. $

(28)

其中：J⁺=J^T(JJ^T)^-1表示J的伪逆矩阵；N=null(J)是列满秩的8×2矩阵，N的两个列向量是J的化零空间(null space)；λ=[λ₁  λ₂]^T是任意二维向量；J⁺(-F)是式(6)的一个最小最佳二乘解，Nλ是式(6)的齐次通解。

2.3.1 索力多边形

索力多边形所形成的区域称为索力可行域，是研究索力分配及控制的基础，本节将给出其定义及计算方法。

令Σ⊂ R⁸是式(6)的二维解向量空间，其表示如下：

$ \Sigma = \left\{ {\mathit{\boldsymbol{T}}\left| {\mathit{\boldsymbol{JT}} + \mathit{\boldsymbol{F}} = 0} \right.} \right\}. $

(29)

令Ω⊂ R⁸是由可行索力组成的八维超立方体，其表示如下：

$ \mathit{\Omega } = \left\{ {\mathit{\boldsymbol{T}}\left| {{t_i} \in \left[ {{t_{\min }},{t_{\max }}} \right],1 \le i \le 8} \right.} \right\}. $

(30)

其中: t_min和t_max分别表示绳索所允许的最小值和最大值，本文分别取10 N和1 000 N(机构采用8根相同绳索)。交集Λ=Σ∩Ω是一个凸多边形^[18]，称为索力多边形。Λ是由式(6)解中满足t_min≤t_i≤t_max部分组成，如图 2所示，八索并联机构的动平台在位姿x=[0.1, 0.1, 0.8, 5, 5, 5]处，承受外载荷F=[50, 50, 50, 10, 10, 10]时所形成的索力多边形。

由式(28)得，索力多边形可由如下2m个不等式定义：

$ {t_{\min }} - {\mathit{\boldsymbol{t}}_p} \le \mathit{\boldsymbol{N}}\lambda \le {t_{\max }} - {\mathit{\boldsymbol{t}}_p}. $

(31)

当索力取最小值t_min或最大值t_max时，随着λ的变化，每个不等式将通过一条边界直线定义一个半平面。式(31)所定义的2m个半平面相交可形成索力多边形Λ。

索力边界直线可用

$ L = \left\{ {{L_{ij}}\left| {i = 1,2, \cdots ,m;j = 1,2} \right.} \right\} $

表示，索力直线交点为C₂^2m-m=4C₂^m。机构驱动索两两之间不相互平行，若出现平行，则机构发生奇异。

2.3.2 索力多边形计算

根据索力多边形和索力边界直线特点，基于Graham‘s Scan扫描法设计索力多边形计算算法，计算原理为：通过循环计算索力边界直线交点所形成的凸包，逐步计算出索力多边形。计算流程如图 3所示。

步骤1 计算索力边界直线交点，并删除重复点，得点序列A：< p₀, p₁, …, p_n>。

步骤2 计算A中最低且最左点，作为基点p₀。

步骤3 对序列 < p₁, p₂, …, p_n>中各点相对于p₀的极角按从小到大进行排序，相应的，点序列更新为A′：< p′₁, p′₂, …, p′_n>(对于极角相同点，根据与p₀的距离排序)，如图 4所示。

步骤4 将p₀、p₁^′先后压入栈S。

步骤5 赋值i=1。

步骤6 赋值i=i+1。

步骤7 令当前点p=p_i^′。

步骤8 连接S顶端两个元素p_top-1和p_top，得到有向线段L。

步骤9 若p在L右边，则栈顶元素出栈，转步骤8，否则执行下一步。

步骤10 将p压入栈。

步骤11 若p=p_n^′，则执行下一步，否则转步骤6。

步骤12 若栈S内元素连线包含实边，则执行下一步，否则转步骤14。

步骤13 若A包含不满足实边的点，则将不满足点删除，并转步骤2，否则转步骤15。

步骤14 若A封闭，则令A=A-S，转步骤2，否则执行下一步。

步骤15 若A包含无效点，则删除无效点，直至A不包含无效点，否则，输出A。

步骤16 结束。

下面给出本算法中一些专用名词定义和索力多边形性质。

定义1  设p={p_i|i=1, 2, …, n}为索力边界直线交点集合，其中，p_j、p_k(j≠k；j, k=1, 2, …, n)为p内任意两元素，若∃l∈L，使得线段p_jp_k与l共线，则称p_jp_k为索力多边形的实边；若∀l∈L，都不与线段p_jp_k共线，则称p_jp_k为索力多边形的虚边，如图 5所示。

定义2  设p={p_i|i=1, 2, …, n}为索力边界直线交点集合，p_s={p_sⁱ|i=1, 2, …, m}为索力多边形顶点集合，若∀p_k∈p, p_k∉p_s，则称p_k为无效点；若∀p_k∈p, p_k∈p_s，则称p_k为有效点，如图 5所示。

索力多边形具有以下性质。

设p_c={p_cⁱ|i=1, 2, …, n}为索力多边形顶点集合，且n≥3，其中，p_c^j、p_c^k(j≠k; j, k=1, 2, …, n)为p_c为内任意两个元素，那么，p_c^j与p_c^k不相邻是线段p_c^jp_c^k为索力多边形虚边的充要条件，p_c^j与p_c^k相邻是线段p_c^jp_c^k为索力多边形实边的充要条件，如图 6所示。

定义3  设p_r={p_rⁱ|i=1, 2, …, n}为索力边界直线的交点集合，且n≥3，其中，S_r是p_r的凸包，p_d=p_r-S_r，若点集合p_d可形成封闭区域，则称p_r封闭，反之，称p_r非封闭，如图 7所示。

2.3.3 形心直线

点集合p_c是索力多边形Λ的顶点，且p_cⁱ，i≠1按相对p_c¹的极角大小进行排序，其中p_cⁱ=[λ_i1, λ_i2]^T，且n≥3，这n个顶点可由节2.3.2所述算法计算获得。索力多边形Λ的形心C=[λ_c1, λ_c2]计算如下：

$ \left\{ \begin{array}{l} {\lambda _{c1}} = \frac{1}{{6A}}\sum\limits_{i = 1}^{n - 1} {\left( {{\lambda _{i1}} + {\lambda _{\left( {i + 1} \right)1}}} \right)\left( {{\lambda _{i1}}{\lambda _{\left( {i + 1} \right)2}} - {\lambda _{\left( {i + 1} \right)1}}{\lambda _{i2}}} \right)} ,\\ {\lambda _{c2}} = \frac{1}{{6A}}\sum\limits_{i = 1}^{n - 1} {\left( {{\lambda _{i2}} + {\lambda _{\left( {i + 1} \right)2}}} \right)\left( {{\lambda _{i1}}{\lambda _{\left( {i + 1} \right)2}} - {\lambda _{\left( {i + 1} \right)1}}{\lambda _{i2}}} \right)} . \end{array} \right. $

(32)

其中A是索力多边形的面积。A计算如下：

$ A = \frac{1}{2}\sum\limits_{i = 1}^{n - 1} {\left( {{\lambda _{i1}}{\lambda _{\left( {i + 1} \right)2}} - {\lambda _{\left( {i + 1} \right)1}}{\lambda _{i2}}} \right)} . $

(33)

如图 2所示，索力多边形内的圆点即为Λ的形心C。在形心C处，索力向量T_C=J⁺(-F)+Nλ_C=t_p+Nλ_C。

过形心的直线可通过下式计算获得：

$ \left\{ \begin{array}{l} {\lambda _2} - {\lambda _{c2}} = \left( {\tan \alpha } \right)\left( {{\lambda _1} - {\lambda _{c1}}} \right),\;\;\;\;\alpha \in \left( { - \frac{{\rm{ \mathsf{ π} }}}{2},\frac{{\rm{ \mathsf{ π} }}}{2}} \right);\\ {\lambda _1} = {\lambda _{c1}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\alpha = \frac{{\rm{ \mathsf{ π} }}}{2}. \end{array} \right. $

(34)

其中α为直线与λ₁轴的夹角。

例如，α=α₁+kμ，初值α₁=-5π/12，μ=π/12，当k依次取值0, 1, …, 11时，可分别获得八索并联机构的动平台在位姿x=[0.1, 0.1, 0.8, 5, 5, 5]处，承受外载荷F=[50, 50, 50, 10, 10, 10]时的过形心直线，如图 8所示。

3 仿真验证与分析

为了验证本文所研究内容的正确性与有效性，本节设计了刚度衡量方法和仿真实例对其进行验证与分析。

3.1 刚度衡量方法

为了有效衡量系统刚度变化情况，提出采用“外力与动平台位姿变化量关系”的方法来衡量系统刚度变化。显然，某位姿处，在同等外力作用下，动平台位姿变化越大，机构系统刚度越小，反之亦然。

由式(7)得：

$ \Delta \mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{K}}^ - } \cdot \mathit{\boldsymbol{F}}, $

(35)

其中：Δx=[Δx_o′, Δy_o′, Δz_o′, Δα, Δβ, Δγ]^T表示位姿变化量，K^-=(K^TK)^-1·K^T表示K的广义逆矩阵。

由于外力已知，根据式(35)可以求出动平台位姿变化量，以便于对系统刚度进行分析。

3.2 仿真验证与分析

将上文提出的算法及方法转化为MATLAB程序进行仿真验证，以α=π/4，动平台在位姿x=[0.1, 0.1, 0.8, 5, 5, 5]处，承受外载荷F=[50, 50, 50, 10, 10, 10]为例，对本文方法及理论进行验证。通过式(34)可得到过形心的线段如图 9a所示，通过式(28)可得到沿该线段索力变化曲线如图 9b所示。

由图 9b可以看出，从左至右沿图 9a形心线段，各绳索的索力逐渐增加。通过式(35)可得到沿该形心线段的位姿变化曲线如图 10所示, 端点处的位姿变化值，如表 2所示。

通过测试表明，本文提出的算法可以计算出索并联机构的索力多边形，并能够在索力可行域内选择合适的索力来控制机器人的刚度。当λ₁沿图 9a形心线变化时，如图 9b所示，各绳索索力逐渐增大，此时，除Δy_o′外，位姿各分量变化都逐渐减小，说明随着索力增加，这些方向的刚度逐渐增大。其中，在该位姿处，增加最快的x和γ方向，分别将近增加41%和77%，大幅度的改变了系统的刚度，因此，该方法可以有效地控制系统刚度。

Δy_o′随着索力增加而变大，是由于机构索力增加不对称引起的，如图 9b所示，方向3、4、7、8号绳索索力增加明显大于相反方向4根绳索索力增加速度，这就为提高系统整体刚度的同时降低某方向的刚度提供了可能。

4 结论

本文对八索驱动并联机器人位置与静力学进行分析，通过引入线矢量和微分变换等方式建立其静刚度模型，并提出索力多边形计算算法，在索力可行域内选取合适索力来控制机器人刚度，得出结论: 1)提出的索力多边形计算算法，可以计算动平台在某位姿处的索力可行域；2)在索力可行域内，通过选择合适索力，根据刚度需求，可在一定范围内对机器人刚度进行有效控制。本文提出的索力可行域计算和刚度控制方法，为后续研究索机器人索力优化和刚度自适应控制奠定了基础。