绳牵引并联机器人是一类将驱动器的运动和力通过绳索传递给末端执行器的并联机器人^[1]。由于其具有质量轻、可拆装重组、结构简单、运动速度快、承载能力大、运动链惯性小、干涉发生几率小、累积误差小以及工作空间大等优点^[2]，一直都是研究的热点。目前，绳牵引并联机器人在物料搬运^[3]、天文观测^[4]、患者康复^[5]、运动模拟^[6]、风洞实验^[7]、三维打印^[8-10]等领域得到了广泛的应用，成功的工程应用范例包括运用于体育赛事转播的绳牵引并联摄像机器人Skycam^[11]、应用于大型构件吊装、起重工作的RoboCrane^[12-13]，应用在2015年米兰世博会德国馆的绳索牵引移动显示屏^[14]，以及应用于500 m口径球面射电望远镜(five-hundred meter aperture spherical telescope，FAST) 的绳牵引馈源支撑系统^[15-22]。

目前，大多数绳牵引并联机器人的研究把绳索理想化为无质量直线模型进行，该方法容易忽略绳索质量对运动性能的影响，造成较大的模型误差。当绳牵引并联机构的跨度比较大时，绳索质量和下垂对机器人的影响是不可忽略的^[23]。因此，采用可真实反映绳索特性的悬链线模型能够推导出更精确的动力学模型，进而设计出精度更高的控制策略。目前，研究人员通过绳索悬链线模型对绳牵引并联机器人进行了研究^{[17, 23-25]}，但尚未对悬链线模型超越方程解析解的高精度求解和演化规律、精确数值解及其分布规律进行研究。Irvine^[25]用悬链线模型描述了绳索在自重作用下的形状。国内研究学者在不考虑绳索惯性力前提下根据非线性悬链线模型来描述绳索形状^[26-29]，对500 m口径球面射电望远镜(FAST)的绳牵引馈源支撑系统的动力学进行研究。Du等^[30]根据非线性悬链线模型通过有限元离散化的方法研究了大跨度绳索收放运动时机构的动力学行为。该模型尚未考虑绳索惯性力对动力学的影响。Su等^[31]利用非线性悬链线模型构建了快速时变绳索动力学模型，建立了包含绳索惯性力的索力优化求解模型和优化迭代算法，但尚未对悬链线模型的解析解和数值解进行研究。

综上可知，目前对绳牵引并联机器人绳索非线性悬链线模型超越方程的降维解算研究较少。因此，本文针对绳索非线性悬链线模型高效快速降维解算问题进行研究，为实现动力学精确建模和稳定运动实时控制奠定基础。首先，建立绳索非线性悬链线模型超越方程并确定其边界条件；其次，建立非线性悬链线模型超越方程的系数矩阵；进而，根据系数矩阵行列式为零的方法对超越方程进行降维，从而根据Taylor展开方法求出悬链线模型的解析解；最后，利用Newton迭代法对悬链线模型精确数值解的关系特征和变化规律进行研究，确定数值解的正确符号。

1 悬链线模型的差分方程和边界条件

考虑一根两端悬挂的长度为L_i的柔索，如图 1所示。其中，OXY是全局固定参考坐标系，O是原点。柔索i=(1，2，…，m)与滑轮相切的点状铰接点A_i假设固定在全局固定框架上。此外，第i根柔索连接在末端执行器上的点状铰接点B_i处，并且假设该连接点相对于移动平台固定。g是重力加速度。

根据上图可得柔索的悬链线方程如下^[32]:

$ y(x) = {\kappa _0} + {\kappa _1}\cosh \left( {\frac{x}{{{\kappa _1}}} + {\kappa _2}} \right). $

(1)

其中，x和y表示绳索上任意一点的横坐标和纵坐标；κ₀、κ₁和κ₂表示系数，即需要求解的未知量。式(1)代表的数学悬链线方程描述在重力作用下柔索的形状。给定曲线上任意点的坐标，可以完全确定整个柔索曲线。

本文做了3点假设对理论模型进行简化：1) 绳索与末端执行器、机架的连接均为理想铰接；2) 绳索既不能受压同时也不能受弯^[24]，满足理想柔性的条件；3) 绳索的变形符合Hooke定律^[24]。

根据绳牵引并联机器人的结构特点可知柔索须满足以下3个约束条件：

1) 柔索长度约束条件

$ \int_{{x_A}}^{{x_B}} {\sqrt {1 + {{\left[ {\frac{{{\rm{d}}y(x)}}{{{\rm{d}}x}}} \right]}^2}} } {\rm{d}}x - L = 0. $

(2)

2) 端点A_i的几何约束条件

$ {A_i}{y_A} = {\kappa _0} + {\kappa _1}\cosh \left( {\frac{{{x_A}}}{{{\kappa _1}}} + {\kappa _2}} \right). $

(3)

3) 端点B_i的几何约束条件

$ {y_B} = {\kappa _0} + {\kappa _1}\cosh \left( {\frac{{{x_B}}}{{{\kappa _1}}} + {\kappa _2}} \right). $

(4)

2 超越方程的降维

把式(1)代入式(2)中并积分得到

$ \begin{array}{*{20}{c}} { - {\kappa _1}\sinh {\kappa _2}\cosh \frac{{{x_A}}}{{{\kappa _1}}} - {\kappa _1}\cosh {\kappa _2}\sinh \frac{{{x_A}}}{{{\kappa _1}}} + }\\ {{\kappa _1}\sinh {\kappa _2}\cosh \frac{{{x_B}}}{{{\kappa _1}}} + {\kappa _1}\cosh {\kappa _2}\sinh \frac{{{x_B}}}{{{\kappa _1}}} - L = 0.} \end{array} $

(5)

$ {\rm{令}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sinh {\kappa _2} = \frac{{2\tanh \frac{{{\kappa _2}}}{2}}}{{1 - {{\tanh }^2}\frac{{{\kappa _2}}}{2}}}, $

$ \cosh {\kappa _2} = \frac{{1 + {{\tanh }^2}\frac{{{\kappa _2}}}{2}}}{{1 - {{\tanh }^2}\frac{{{\kappa _2}}}{2}}},\tanh \frac{{{\kappa _2}}}{2} = t, $

$\frac{1}{\kappa_{1}}=\lambda>0 $ ，式(5)可写成

$ \begin{array}{*{20}{c}} {{t^2}\left[ {\sinh \left( {\lambda {x_A}} \right) - \sinh \left( {\lambda {x_B}} \right) - \lambda L} \right] + }\\ t\left[ {2\cosh \left( {\lambda {x_A}} \right) - 2\cosh \left( {\lambda {x_B}} \right)} \right] + \\ \sinh (\lambda {x_A}) - \sinh (\lambda {x_B}) + L\lambda = 0. \end{array} $

(6)

式(6)是关于t的二次多项式方程。

因为柔索同时通过端点 $A_{i} \text { 和 } B_{i}, \frac{1}{\kappa_{1}}=\lambda>0 $ ，式(4)减去式(3)整理得

$ \begin{array}{*{20}{c}} {\lambda {y_B} - \lambda {y_A} + \cosh \left( {\lambda {x_A}} \right)\cosh {\kappa _2} + }\\ {\sinh \left( {\lambda {x_A}} \right)\sinh {\kappa _2} - \cosh \left( {\lambda {x_B}} \right)\cosh {\kappa _2} - }\\ {\sinh \left( {\lambda {x_B}} \right)\sinh \left( {{\kappa _2}} \right) = 0.} \end{array} $

(7)

根据 $\tanh \frac{\kappa_{2}}{2}=t $ 整理式(7)得

$ \begin{array}{*{20}{c}} {{t^2}\left[ {\lambda {y_B} - \lambda {y_A} - \cosh \left( {\lambda {x_A}} \right) + \cosh \left( {\lambda {x_B}} \right)} \right] + }\\ {t\left[ {2\sinh \left( {\lambda {x_B}} \right) - 2\sinh \left( {\lambda {x_A}} \right)} \right] + }\\ {\lambda {y_A} - \lambda {y_B} - \cosh \left( {\lambda {x_A}} \right) + \cosh \left( {\lambda {x_B}} \right) = 0.} \end{array} $

(8)

式(8)是关于t的二次多项式方程。

至此，把式(5)和式(7)的超越方程通过换元法降维变为关于t的二次多项式方程。

令式(6)的系数方程如下：

$ \left\{ {\begin{array}{*{20}{l}} {{a_1} = \sinh \left( {\lambda {x_A}} \right) - \sinh \left( {\lambda {x_B}} \right) - \lambda L,}\\ {{b_1} = 2\cosh \left( {\lambda {x_A}} \right) - 2\cosh \left( {\lambda {x_B}} \right),}\\ {{c_1} = \sinh \left( {\lambda {x_A}} \right) - \sinh \left( {\lambda {x_B}} \right) + \lambda L.} \end{array}} \right. $

(9)

令式(8)的系数方程如下：

$ \left\{ {\begin{array}{*{20}{l}} {{a_2} = \lambda {y_B} - \lambda {y_A} - \cosh \left( {\lambda {x_A}} \right) + \cosh \left( {\lambda {x_B}} \right),}\\ {{b_2} = 2\sinh \left( {\lambda {x_B}} \right) - 2\sinh \left( {\lambda {x_A}} \right),}\\ {{c_2} = \lambda {y_A} - \lambda {y_B} - \cosh \left( {\lambda {x_A}} \right) + \cosh \left( {\lambda {x_B}} \right).} \end{array}} \right. $

(10)

则式(6)与(8)可写成：

$ \left\{ {\begin{array}{*{20}{l}} {{a_1}{t^2} + {b_1}t + {c_1} = 0,}\\ {{a_2}{t^2} + {b_2}t + {c_2} = 0.} \end{array}} \right. $

(11)

根据式(11)无法构建一个n阶矩阵进行计算，式(11)的左右两边同时乘以t，得到可以构成如下4×4的系数矩阵的方程组：

$ \left\{ {\begin{array}{*{20}{l}} {{a_1}{t^3} + {b_1}{t^2} + {c_1}t = 0,}\\ {{a_2}{t^3} + {b_2}{t^2} + {c_2}t = 0,}\\ {{a_1}{t^2} + {b_1}t + {c_1} = 0,}\\ {{a_2}{t^2} + {b_2}t + {c_2} = 0.} \end{array}} \right. $

(12)

系数矩阵R表示如下：

$ \mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {{a_1}}&{{b_1}}&{{c_1}}&0\\ {{a_2}}&{{b_2}}&{{b_2}}&0\\ 0&{{a_1}}&{{b_1}}&{{c_1}}\\ 0&{{a_2}}&{{b_2}}&{{b_2}} \end{array}} \right]. $

(13)

当线性齐次方程组有非零解时，其系数矩阵的行列式应取0值。即当式(6)和(8)同时满足等式有非零解时，系数矩阵的行列式取值为0，即有

$ \begin{array}{*{20}{c}} - 2{\cosh ^2}(\lambda {x_A} - \lambda {x_B}) + \\ {\left[ {{L^2}{\lambda ^2} - y_A^2{\lambda ^2} + 2{y_B}{y_A}{\lambda ^2} - y_B^2{\lambda ^2} + 4} \right] \cdot }\\ {\cosh \left( {\lambda {x_A} - \lambda {x_B}} \right) - {L^2}{\lambda ^2} + y_A^2{\lambda ^2} - }\\ {2{y_B}{y_A}{\lambda ^2} + y_B^2{\lambda ^2} - 2 = 0.} \end{array} $

(14)

3 悬链线模型的解析解

式(14)是关于 $ \cosh \left(\lambda x_{A}-\lambda x_{B}\right)$ 的一元二次方程，求解得出 $ \cosh \left(\lambda x_{A}-\lambda x_{B}\right)$ 的值如下：

$ \begin{array}{*{20}{c}} {\cosh \left( {\lambda {x_A} - \lambda {x_B}} \right) = }\\ {\frac{1}{2}{L^2}{\lambda ^2} - \frac{1}{2}y_A^2{\lambda ^2} + {y_B}{y_A}{\lambda ^2} - \frac{1}{2}y_B^2{\lambda ^2} + 1.} \end{array} $

(15)

因为 $ \lambda>0, x_{B}>x_{A}, \text { 令 } \lambda=\frac{\nu}{x_{B}-x_{A}}>0$ ，整理式(15)得

$ \cosh \nu = \frac{{{\nu ^2}}}{{2{{\left( {{x_B} - {x_A}} \right)}^2}}}\left[ {{L^2} - {{\left( {{y_B} - {y_A}} \right)}^2}} \right] + 1. $

(16)

令 $ C=\frac{\left[L^{2}-\left(y_{B}-y_{A}\right)^{2}\right]}{\left(x_{B}-x_{A}\right)^{2}}$ ，根据三角形勾股定理可知C≥1，式(16)整理为

$ \cosh \nu = 1 + \frac{1}{2}C{\nu ^2}. $

(17)

对cosh ν进行Taylor展开，式(17)可写为

$ 1 + \frac{1}{2}{\nu ^2} + \frac{1}{{24}}{\nu ^4} = 1 + \frac{1}{2}C{\nu ^2}. $

(18)

对式(18)中的ν进行求解得出如下3个值：

$ 0,2{\kern 1pt} {\kern 1pt} \sqrt {3C - 3} , - 2{\kern 1pt} {\kern 1pt} \sqrt {3C - 3} . $

(19)

由于 $\kappa_{1}=\frac{\mu}{T_{0}}=\frac{1}{\lambda}> $0 ^[32]，μ是柔索线密度，T₀是索力，因此，$\lambda=\frac{\nu}{x_{B}-x_{A}}>0 $ ，故ν>0。得到Taylor展开后的解析解为

$ \nu = 2{\kern 1pt} {\kern 1pt} \sqrt {3C - 3} . $

(20)

4 悬链线模型的数值解特性

令 $f\left(\nu_{n}\right)=\cosh (\nu)-1-\frac{1}{2} C \nu^{2} $ ，采用Newton迭代法 $ \nu_{n+1}=\nu_{n}-\frac{f\left(\nu_{n}\right)}{f^{\prime}\left(\nu_{n}\right)}$ ，可以求解C的精确数学值。ν和C的关系如图 2所示，可以看出，C的取值随着ν的增加而增加。

假设μ=1 kg/m，L= 6 m，x_A= 1 m，y_A= 2 m，x_B=4 m，y_B=-1 m，从图上可知，当C=3时，ν= 4.898 979 485 6 (m²·N)/kg，T₀ = 0.783 980 175 4 N，κ₀= [3.633 608 090-5.879 067 841×10^-10i, -2.633 607 826]m，κ₁=1.275 542 458 1 kg/(m·N)，κ₂=[-2.639 549 976+3.141 592 654i, -3.738 162 224] (m²·N)/kg。本研究看到κ₀和κ₂分别有一个虚数和一个实数2个值，因此必须对这2个值进行验证。下面通过例子证明如下：

实例1

假设(x_A，y_A)=(2，3) m和(x_R，y_R)=(5，0) m，L=8 m，μ=1 kg /m。悬链线形状如图 3所示。

显然，正确的解是κ₀和κ₂分别取实数的曲线解2，因为κ₀和κ₂分别取虚数时的曲线解1与重力方向相反。也就是说，曲线解2对应于势能中的最小值，而曲线解1对应于最大值。

函数y(x)的二阶导数的符号能够帮助区分势能中的最小值和最大值：

$ {y^{\prime \prime }}(x) = \frac{{\cosh \left( {\frac{x}{{{\kappa _1}}} + {\kappa _2}} \right)}}{{{\kappa _1}}}. $

(21)

根据反三角函数性质，κ₁应该总是一个实常数，但是κ₂可能实数也可能是虚数。本研究假设κ₂的实部是κ_2r，而其虚部是κ_2i，有

$ \begin{array}{*{20}{c}} {{y^{\prime \prime }}(x) = \frac{{\cosh \frac{x}{{{\kappa _1}}}\cosh {\kappa _{2{\rm{r}}}}\cos {\kappa _{2{\rm{i}}}}}}{{{\kappa _1}}} + }\\ {I\frac{{\cosh \frac{x}{{{\kappa _1}}}\sinh {\kappa _{2{\rm{r}}}}\sin {\kappa _{2{\rm{i}}}}}}{{{\kappa _1}}} + }\\ {\frac{{\sinh \frac{x}{{{\kappa _1}}}\sinh {\kappa _{2{\rm{r}}}}\cos {\kappa _{2{\rm{i}}}}}}{{{\kappa _1}}} + }\\ {\frac{{\sinh \frac{x}{{{\kappa _1}}}\cosh {\kappa _{2{\rm{r}}}}\sin {\kappa _{2{\rm{i}}}}}}{{{\kappa _1}}}.} \end{array} $

(22)

y(x)的二阶导数的实部的符号：

$ \frac{{\cos {\kappa _{2{\rm{i}}}}\cosh \left( {\frac{{{\kappa _{2{\rm{r}}}}{\kappa _1} + x}}{{{\kappa _1}}}} \right)}}{{{\kappa _1}}}. $

(23)

显然，κ₁总是正的实数，否则就没有可行的解，当－ $ \frac{\pi}{2} \leqslant \kappa_{2 \mathrm{i}} \leqslant \frac{\pi}{2} $ 时，后者的表达是非负的。因此，应该保留假想κ_2i对应于|κ_2i|的解 $\kappa_{2} \leqslant \frac{\pi}{2} $ 。因此，例1的等式是

$ \begin{array}{*{20}{c}} {y(x) = - 2.550{\kern 1pt} {\kern 1pt} 561{\kern 1pt} {\kern 1pt} 598 + 0.591{\kern 1pt} {\kern 1pt} 830{\kern 1pt} {\kern 1pt} 280{\kern 1pt} {\kern 1pt} 8 \times }\\ {\cosh \left( {\frac{x}{{0.591{\kern 1pt} {\kern 1pt} 830{\kern 1pt} {\kern 1pt} 280{\kern 1pt} {\kern 1pt} 8}} - 6.308{\kern 1pt} {\kern 1pt} 074{\kern 1pt} {\kern 1pt} 428} \right)}。\end{array} $

总之，κ₁和κ₂的正确解是曲线解2。下面通过求解柔索长度的例子证明κ₀、κ₁和κ₂的值是合理的。

实例2

$ {L^\prime } = \int_{{x_A}}^{{x_B}} {\sqrt {1 + {{\left[ {\frac{{{\rm{d}}y(x)}}{{{\rm{d}}x}}} \right]}^2}} } {\rm{d}}x. $

(24)

求解式(24)得到

$ \begin{array}{*{20}{c}} {{L^\prime } = \frac{1}{2}\left\{ {{\kappa _1}\left[ {\cosh \left( {\frac{{2{\kappa _1}{\kappa _2} + {x_A} + 2{x_B}}}{{{\kappa _1}}}} \right) + } \right.} \right.}\\ {\sinh \left( {\frac{{2{\kappa _1}{\kappa _2} + {x_A} + 2{x_B}}}{{{\kappa _1}}}} \right) - \cosh \left( {\frac{{2{\kappa _1}{\kappa _2} + 2{x_A} + {x_B}}}{{{\kappa _1}}}} \right) - }\\ {\sinh \left( {\frac{{2{\kappa _1}{\kappa _2} + 2{x_A} + {x_B}}}{{{\kappa _1}}}} \right) + \cosh \frac{{{x_B}}}{{{\kappa _1}}} + }\\ {\left. {\sinh \frac{{{x_B}}}{{{\kappa _1}}} - \cosh \frac{{{x_A}}}{{{\kappa _1}}} - \sinh \frac{{{x_A}}}{{{\kappa _1}}}} \right]}\\ {\left. {\left[ {\cosh \left( {\frac{{ - {\kappa _1}{\kappa _2} - {x_A} - {x_B}}}{{{\kappa _1}}}} \right) + \sinh \left( {\frac{{ - {\kappa _1}{\kappa _2} - {x_A} - {x_B}}}{{{\kappa _1}}}} \right)} \right]} \right\}.} \end{array} $

(25)

将x_A= 2 m，x_B= 3 m，κ₀= -2.550 561 598 m，κ₁= 0.591 830 280 8 kg/(m·N)，κ₂= -6.308 074 428 (m²·N)/kg代入方程(23)中，得到L^′=7.999 989 875 m，跟实例1中的L=8 m比较误差仅为0.000 010 125 m，可认为二者是相等的。也就是说κ₀、κ₁和κ₂的取值是合理的，因为求解出的柔索长度和实际柔索长度相等。

5 结论

针对考虑绳索质量的绳牵引并联机器人结构特点，通过柔索微分单元和积分法推导了柔索悬链线模型差分方程。同时，基于柔索长度积分关系和柔索端点坐标确定了悬链线模型的边界条件。提出了一种利用系数矩阵行列式为零的方法对悬链线模型的超越方程进行降维，实现了悬链线模型由超越方程降为一元二次方程。为了研究绳牵引并联机器人悬链线特性的解析解和数值解机理，分别提出了一种基于Taylor展开法的悬链线模型解析解的求解方法和一种基于Newton迭代法的悬链线模型数值解的求解方法。通过Taylor展开方法对悬链线模型解析解进行了求解。利用Newton迭代法对悬链线模型中参数精确解中间变量的关系特征和变化规律进行了研究。最后，结合实例研究悬链线模型数值解的取值范围和特性，确定数值解的正确符号。