作为日常维护的重要内容，船舶清洗不容忽视，目前主要通过人工或爬壁机器人完成^[1]。船体外壁面积较大，作业人员难以实现无死角作业。虽然爬壁机器人在一定程度上可提高作业效率，但是船体表面形状不规则限制了爬壁机器人的作业范围^[2-3]。因此，开发一种新型高效的自动化清洗装备尤显重要。

绳驱动并联机器人(cable-driven parallel robots，CDPRs)属于典型的并联机构机器人，具有工作空间大、运动惯性小及易重构等优点。针对不同的作业场景，研究人员对绳驱动并联机构开展了广泛研究。RoboCrane是最早的绳驱动并联机构之一，由Bostelman等^[4]开发，现已广泛应用于飞机维护、材料搬运等场景。Nan^[5]采用绳驱动并联机构作为望远镜接收机的支撑系统，并进一步提出一种500 m口径的球面射电望远镜(five-hundred-meter aperture spherical radio telescope，FAST)，为开展天文观测研究提供了更精确的工具。为实现自动化码垛、搬运作业，El-Ghazaly等^[6]设计并开发了CoGiRo机器人，该机器人可以帮助提高生产效率并降低人力成本。李政清等^[7]将绳驱动并联机构应用于物流领域，可在大型仓库开展货物运输作业。王晓光等^[8]将绳驱动机器人应用于飞行器测试，为风洞试验提供了一种新的解决方式。在建筑领域，Wu等^[9]开发了一种绳驱动砖结构建筑机器人(CU-brick cable-driven robot, CU-Brick)，用于实现建筑物料的大范围运输。文[10-11]采用平面绳驱动并联机构作为支撑系统，为建筑物外墙自动化清洗作业提供了一种新的解决方案。受到上述应用启发，本文提出一种新型变结构空间绳驱动清洗机器人，旨在解决人工作业效率低、风险高等问题。与FAST相似，清洗机器人的支撑系统也采用绳驱动并联机构，通过在运动平台上安装清洗工具，实现大范围船壁清洗作业的目标。

绳驱动并联机构的动力学模型一直是研究热点。赵雅聪等^[12]采用虚功原理建立变长度索单元的动力学模型，并将其应用于FAST的馈源支撑系统。Caverly等^[13]采用集总质量法和零空间法推导平面绳驱动并联机构的动力学模型，并考虑了运动平台移动对绳索刚度及绞车惯性的影响。Abdallah等^[14]使用绳驱动技术稳定飞艇吊舱，并以Kirchhoff方程建立了系统动力学模型。Kino等^[15]基于Lagrange法建立了三自由度绳驱动机器人的动力学模型，该动力学模型考虑了负载惯性参数的不确定性对控制精度的影响。通过总结现有的绳驱动并联机构的动力学模型可以发现，大多数研究集中于定基座工况，尚未考虑基座运动对模型精度的影响。此外，大多数研究均简化处理了运动平台所受外部扰动力，无法精确了解扰动力对跟踪精度的影响。清洗机器人在作业时会受到船舶的多自由度运动及外部扰动的作用，现有模型并不适用。

迄今，学者们已经提出了许多绳驱动并联机构的鲁棒控制策略。Wang等^[16]提出一种基于时延估计的自适应超扭转非奇异快速终端滑模控制器，该控制器可在系统未知且不确定的情况下，完成高精度轨迹跟踪。Jia等^[17]采用二阶滑模控制器控制绳索张力，改善了绳索间的同步运动关系，并可以预测模型的不确定性及外部干扰。Caverly等^[18]提出一种正时控制器，该控制器可使绳索张力保持正值。韦慧玲等^[19]提出一种比例-微分修正前馈控制器，在运动平台高速运动状态下，该控制器具备良好的稳态跟踪性能。通过总结以上研究可知，大多数控制策略并未考虑基座运动的影响，一般仅适用于定基座工况。清洗机器人的应用场景为船舶，船舶运动将不可避免地影响轨迹跟踪精度，甚至造成系统闭环不稳定。因此，动基座激励下如何实现高精度轨迹跟踪仍是一个开放性问题。

首先，采用Newton-Euler法建立包含船舶运动和外部扰动的动力学模型。其次，在动力学模型基础上提出了一种模糊自适应比例-积分(proportional-integral, PI)滑模控制器(fuzzy adaptive PI sliding mode controller，FAPI-SMC)，该控制器可在船舶运动下实现快速响应及高精度轨迹跟踪。最后，采用Lyapunov理论证明控制系统的稳定性，并通过仿真结果对比说明FAPI-SMC的有效性。

1 清洗机器人结构及工作原理

绳驱动并联机器人按照结构形式可分为平面、空间机器人。平面机器人仅能实现面内作业，而空间机器人的绳索布置复杂，易与其他结构发生干涉。通过改进平面绳驱动并联机构，设计的空间变结构清洗机器人如图 1所示。空间变结构可概括如下：导绳机构可在滑台驱动下移动，进而改变运动平台的空间位置，使平面绳驱动并联机构能够在三维空间实现清洗作业。

驱动单元用于控制运动平台的空间位姿，包括电机、滑台、编码器和张力传感器，其中，电机用于收放绳索，滑台用于改变导绳机构的空间位置，编码器用于测量电机旋转圈数，张力传感器用于测量绳索张力。运动平台用于完成清洗作业，包括电机1与2、水枪、3D摄像头和激光测距仪；其中，电机1与2分别用于实现水枪水平和垂直方向旋转，3D摄像头用于采集作业面的坐标位置，激光测距仪用于测量运动平台与作业面间的距离。为适应船舱内壁或外壁不规则形状的曲率变化，清洗机器人可以通过调整滑台行程，从而实现多曲率舱壁的清洗作业。此外，所选用的滑台行程越大，对应的运动平台可旋转角度及工作空间越大，清洗机器人的适用曲率范围也相应扩大。

清洗机器人的作业原理如图 2所示。在作业开始前，运动平台在作业区域内快速移动，3D摄像头扫描作业面的位置信息。待作业面的位置坐标扫描完毕后，操作人员向控制系统输入运动平台与作业面间的距离(喷射距离)和运动平台的运动轨迹。在常规作业下，通过控制滑台，使运动平台与作业面间的距离保持为设定值。针对常规作业无法到达的区域(死区)，通过控制电机1和2，使水枪小范围摆动，进一步扩大清洗机器人的作业范围。当运动平台到达预定位置后，水枪喷射清洗剂进行清洗作业。在清洗内凹面时，由于绳索可能会与壁面发生干涉，因此可通过调节水射流压力改变清洗剂的喷射距离，从而完成清洗作业。

2 动力学分析 2.1 滑台运动机理

对文中涉及的坐标系及动力学参数进行解释说明，如图 3所示。其中，f_i为第i根绳索张力，i=1, 2, 3, 4；f_dis和w_dis分别为外部扰动力和外部扰动力矩；G=mg为运动平台所受重力。

O-xyz为全局坐标系，O为坐标原点，坐标轴方向与初始时刻的船舶方向一致，以下在符号表示中以坐标原点代表相应的坐标系。S-xyz为船舶坐标系，坐标原点S位于船舶质心处，x轴方向与船体纵向基线一致，y轴方向与船体横向基线一致，z轴方向由右手定则确定。L₀-xyz为滑台坐标系，坐标原点L₀位于滑台4的初始位置点，坐标轴方向与船舶坐标轴一致。M₀-xyz为运动平台坐标系，坐标原点M₀位于运动平台质心处，坐标轴方向与运动平台位姿有关。C-xyz为3D摄像头坐标系，坐标原点C位于3D摄像头设定的参考坐标点，坐标轴方向与M₀-xyz相同。此外，M_i和L_i分别为运动平台连索点和滑轮连索点，i=1, 2, 3, 4。

本文定义坐标系间的位姿转换关系为$\boldsymbol{q}_{p j}^j=$ $\left[\begin{array}{ll}\boldsymbol{X}_{p j}^j & \boldsymbol{\beta}_{p j}^j\end{array}\right]^{\mathrm{T}}$。其中, $\boldsymbol{X}_{p j}^j$表示位置矢量, $\boldsymbol{\beta}_{p j}^j$表示角度矢量, $p$表示当前坐标系, $j$表示参考坐标系。此外, $\boldsymbol{R}_p^j$表示$p$相对于$j$的旋转矩阵, $\boldsymbol{\omega}_{p j}^j$表示$p$相对于$j$的旋转角速度。

设定$\boldsymbol{q}_{M_0 L_0}^{L_0}=\left[x_{M_0 L_0}^{L_0} \quad y_{M_0 L_0}^{L_0} \quad z_{M_0 L_0}^{L_0} \quad \theta_{M_0 L_0}^{L_0}\right.$  $\left.\phi_{M_0 L_0}^{L_0} \quad \varphi_{M_0 L_0}^{L_0}\right]^{\mathrm{T}}$为$M_0-x y z$相对于$L_0-x y z$的位姿。其中: $x_{M_0 L_0}^{L_0} 、y_{M_0 L_0}^{L_0}$和$z_{M_0 L_0}^{L_0}$表示$M_0-x y z$相对于$L_0-x y z$的位置, $\theta_{M_0 L_0}^{L_0} 、\phi_{M_0 L_0}^{L_0}$和$\varphi_{M_0 L_0}^{L_0}$表示$M_0-x y z$相对于$L_0-x y z$的转动角度。如前所述, 可以通过$3 \mathrm{D}$摄像头实时采集作业面的坐标位置, 设采集到的作业面的坐标信息在$C-x y z$下表示为$\boldsymbol{X}_{\mathrm{ope}}^C=$ $\left[\begin{array}{lll}x_{\mathrm{ope}}^C & y_{\mathrm{ope}}^C & z_{\mathrm{ope}}^C\end{array}\right]^{\mathrm{T}}$, 其在$L_0-x y z$下表示如下:

$ \left\{\begin{array}{l} \boldsymbol{X}_{\mathrm{ope}}^{M_0}=\boldsymbol{X}_{\mathrm{ope}}^C+\boldsymbol{X}_{C M_0}^{M_0}, \\ \boldsymbol{X}_{\mathrm{ope}}^{L_0}=\boldsymbol{X}_{M_0 L_0}^{L_0}+\boldsymbol{R}_{M_0}^{L_0} \boldsymbol{X}_{\mathrm{ope}}^{M_0} . \end{array}\right. $

(1)

其中$\boldsymbol{X}_{\mathrm{ope}}^{L_0}=\left[\begin{array}{lll}x_{\text {ope }}^{L_0} & y_{\text {ope }}^{L_0} & z_{\mathrm{ope}}^{L_0}\end{array}\right]^{\mathrm{T}}$。

在实际作业中, 将滑台按照上下或左右分为两组, 组内同步运动。如图 4所示, 以滑台分为上下两组为例, 作业面在$L_0-x y z$下的坐标$\boldsymbol{X}_{\mathrm{ope}}^{L_0}$与运动平台位姿$y_{M_0 L_0}^{L_0} 、\theta_{M_0 L_0}^{L_0}$的关系表示如下：

$ \left\{\begin{array}{l} y_{M_0 L_0}^{L_0}=y_{\mathrm{ope}}^{L_0}+d_{\mathrm{spr}}, \\ \theta_{M_0 L_0}^{L_0}=\arctan \left(\frac{y_{\mathrm{ope}}^{L_0}(t)-y_{\mathrm{ope}}^{L_0}(t-1)}{z_{\mathrm{ope}}^{L_0}(t)-z_{\mathrm{ope}}^{L_0}(t-1)}\right) . \end{array}\right. $

(2)

其中: $d_{\mathrm{spr}}$为设定的喷射距离, $z_{\mathrm{ope}}^{L_0}(t)$及$y_{\mathrm{ope}}^{L_0}(t)$为$t$时刻的坐标信息。在绳索张紧情况下, 滑台行程$\left(h_1\right.$和$\left.h_2\right)$与$y_{M_0 L_0}^{L_0} 、z_{M_0 L_0}^{L_0}$和$\theta_{M_0 L_0}^{L_0}$有关, 由三角形相似定理得

$ \left\{\begin{array}{l} h_1=y_{M_0 L_0}^{L_0}+z_{M_0 L_0}^{L_0} \tan \left(\theta_{M_0 L_0}^{L_0}\right)-\tan \left(\theta_{M_0 L_0}^{L_0}\right) H, \\ h_2=y_{M_0 L_0}^{L_0}+z_{M_0 L_0}^{L_0} \tan \left(\theta_{M_0 L_0}^{L_0}\right) . \end{array}\right. $

(3)

其中：H为作业面高度，为已知量。

根据式(1)—(3)，在X_ope^C已知的情况下，可求得滑台组的运动行程。

2.2 运动平台动力学模型

在进行作业时，船舶会受到海浪激励的作用而发生摇晃，进而对机器人的作业精度产生影响。因此，建模过程需考虑海浪激励的作用。

$M_0-x y z$相对于$O-x y z$的位置$\boldsymbol{X}_{M_0 O}^O$表示如下：

$ \begin{gathered} \boldsymbol{X}_{M_0 O}^O=\boldsymbol{X}_{S O}^O+\boldsymbol{R}_S^O \boldsymbol{X}_{M_0 S}^S, \\ \boldsymbol{X}_{M_0 S}^S=\boldsymbol{X}_{M_0 L_0}^{L_0}+\boldsymbol{X}_{L_0 S}^S . \end{gathered} $

(4)

对式(4)求导, 可得运动平台的速度$\dot{\boldsymbol{X}}_{M_0 O}^O$及加速度$\ddot{\boldsymbol{X}}_{M_0 O}^O$, 表示如下:

$ \left\{\begin{array}{l} \dot{\boldsymbol{X}}_{M_0 O}^O=\dot{\boldsymbol{X}}_{S O}^O+\left[\boldsymbol{\omega}_{S O}^O \times\right] \boldsymbol{R}_S^O \boldsymbol{X}_{M_0 S}^S+\boldsymbol{R}_S^O \dot{\boldsymbol{X}}_{M_0 S}^S, \\ \ddot{\boldsymbol{X}}_{M_0 O}^O=\ddot{\boldsymbol{X}}_{S O}^O+\left[\dot{\boldsymbol{\omega}}_{S O}^O \times\right] \boldsymbol{R}_S^O \boldsymbol{X}_{M_0 S}^S+ \\ {\left[\boldsymbol{\omega}_{S O}^O \times\right]\left[\boldsymbol{\omega}_{S O}^O \times\right] \boldsymbol{R}_S^O \boldsymbol{X}_{M_0 S}^S+} \\ 2\left[\boldsymbol{\omega}_{S O}^O \times\right] \boldsymbol{R}_S^O \dot{\boldsymbol{X}}_{M_0 S}^S+\boldsymbol{R}_S^O \ddot{\boldsymbol{X}}_{M_0 S}^S . \end{array}\right. $

(5)

其中$\left[\boldsymbol{\omega}_{S O}^O \times\right]$为$\boldsymbol{\omega}_{S O}^O$的斜对称矩阵。与速度及加速度相似, 进一步求得角速度$\boldsymbol{\omega}_{M_0 O}^O$及角加速度$\dot{\boldsymbol{\omega}}_{M_0 O}^O$, 表示如下:

$ \left\{\begin{array}{l} \boldsymbol{\omega}_{M_0 O}^O=\boldsymbol{R}_S^O \boldsymbol{\omega}_{S O}^S+\boldsymbol{R}_S^O \boldsymbol{\omega}_{M_0 S}^S, \\ \dot{\boldsymbol{\omega}}_{M_0 O}^O=\boldsymbol{R}_S^O \dot{\boldsymbol{\omega}}_{S O}^S+\left[\boldsymbol{\omega}_{S O}^O \times\right] \boldsymbol{R}_S^O \boldsymbol{\omega}_{M_0 S}^S+\boldsymbol{R}_S^O \dot{\boldsymbol{\omega}}_{M_0 S}^S . \end{array}\right. $

(6)

根据图 3所示的运动平台受力情况，采用Newton-Euler法建立运动平台的动力学模型。运动平台所受的合力及合力矩分别表示如下：

$ \left\{\begin{array}{l} \sum\limits_{i=1}^4 \boldsymbol{u}_i^O f_i+m \boldsymbol{G}+\boldsymbol{f}_{\text {dis }}=m \ddot{\boldsymbol{X}}_{M_0 O}^O, \\ \sum\limits_{i=1}^4 f_i\left(\boldsymbol{a}_{M_i}^O \times \boldsymbol{u}_i^O\right)+\boldsymbol{w}_{\mathrm{dis}}=\boldsymbol{I} \dot{\boldsymbol{\omega}}_{M_0 O}^O+\boldsymbol{\omega}_{M_0 O}^O \times\left(\boldsymbol{I} \boldsymbol{\omega}_{M_0 O}^O\right), \\ \quad \boldsymbol{I}=\boldsymbol{R}_S^O \boldsymbol{R}_{M_0}^{L_0} \boldsymbol{I}_{M_0 M_0}\left(\boldsymbol{R}_{M_0}^{L_0}\right)^{\mathrm{T}}\left(\boldsymbol{R}_S^O\right)^{\mathrm{T}}, \end{array}\right. $

$ \boldsymbol{I}_{M_0 M_0}=\left[\begin{array}{ccc} \ell_{x x} & -\ell_{x y} & -\ell_{x z} \\ -\ell_{y x} & \ell_{y y} & -\ell_{y z} \\ -\ell_{z x} & -\ell_{z y} & \ell_{z z} \end{array}\right] \text {. } $

(7)

其中: $\boldsymbol{u}_i^O$为绳长单位向量, $i=1, 2, 3, 4 ; \boldsymbol{a}_{M_i}^O$为运动平台质心与其连索点间的向量, $i=1, 2, 3$, $4 ; \boldsymbol{I}$为全局坐标系下的惯性矩阵。

进一步总结式(7)，可得：

$ \begin{array}{c} & \boldsymbol{A} \ddot{\boldsymbol{q}}_{M_0 O}^O+\boldsymbol{B} \boldsymbol{q}_{M_0 O}^O-\boldsymbol{G}-\boldsymbol{F}_{\mathrm{dis}}=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{F}, \\ & \boldsymbol{A}=\left[\begin{array}{cc} m \boldsymbol{\mu}_{3 \times 3} & \mathbf{0}_{3 \times 3} \\ \mathbf{0}_{3 \times 3} & \boldsymbol{I} \end{array}\right], \\ & \boldsymbol{B}=\left[\begin{array}{cc} \boldsymbol{0}_{3 \times 3} & \boldsymbol{0}_{3 \times 3} \\ \boldsymbol{0}_{3 \times 3} & {\left[\boldsymbol{\omega}_{M_0 O}^O \times\right] \boldsymbol{I}} \end{array}\right], \\ & \boldsymbol{F}_{\text {dis }}=\left[\begin{array}{ll} \boldsymbol{f}_{\text {dis }} & \boldsymbol{w}_{\text {dis }} \end{array}\right]^{\mathrm{T}}, \\ & \boldsymbol{G}=\left[\begin{array}{llllll} 0 & 0 & m \boldsymbol{g} & 0 & 0 & 0 \end{array}\right]^{\mathrm{T}}, \\ & \boldsymbol{J}=\left[\begin{array}{ll} \left(\boldsymbol{u}_1^O\right)^{\mathrm{T}} & \left(\boldsymbol{a}_{M_1}^O \times \boldsymbol{u}_1^O\right)^{\mathrm{T}} \\ \left(\boldsymbol{u}_2^O\right)^{\mathrm{T}} & \left(\boldsymbol{a}_{M_2}^O \times \boldsymbol{u}_2^O\right)^{\mathrm{T}} \\ \left(\boldsymbol{u}_3^O\right)^{\mathrm{T}} & \left(\boldsymbol{a}_{M_3}^O \times \boldsymbol{u}_3^O\right)^{\mathrm{T}} \\ \left(\boldsymbol{u}_4^O\right)^{\mathrm{T}} & \left(\boldsymbol{a}_{M_4}^O \times \boldsymbol{u}_4^O\right)^{\mathrm{T}} \end{array}\right] \\ & \end{array} $

(8)

其中: $\boldsymbol{q}_{M_0 O}^O$和$\ddot{\boldsymbol{q}}_{M_0 O}^O$分别为运动平台的速度和加速度矩阵; $\boldsymbol{A}$为质量项, $\boldsymbol{\mu}_{3 \times 3}$为单位矩阵; $\boldsymbol{B}$为动力学模型的Coriolis项; $\boldsymbol{G}$为重力项; $\boldsymbol{F}_{\text {dis }}$为扰动项; $\boldsymbol{J}$为动力学Jacobi矩阵; $\boldsymbol{F}$为绳索张力项。

2.3 电机动力学模型

考虑到电机的迟滞效应，建立电机动力学模型，表示如下^[20]：

$ \begin{array}{c} & \boldsymbol{A}_{\text {inc }} \ddot{\boldsymbol{\alpha}}+\boldsymbol{B}_{\text {vis }} \dot{\boldsymbol{\alpha}}+\boldsymbol{r}_{\text {win }} \boldsymbol{F}=\tau, \\ & \boldsymbol{A}_{\text {ine }}=\operatorname{diag}\left[A_{\text {ine, } 1}, A_{\text {ine, } 2}, \cdots, A_{\text {ine, } 4}\right] \text {, } \\ & \boldsymbol{B}_{\text {vis }}=\operatorname{diag}\left[B_{\text {vis }, 1}, B_{\text {vis }, 2}, \cdots, B_{\text {vis }, 4}\right] \text {, } \\ & \boldsymbol{r}_{\text {win }}=\operatorname{diag}\left[r_{\text {win }, 1}, r_{\text {win }, 2}, \cdots, r_{\text {win }, 4}\right], \\ & \boldsymbol{\alpha}=\left[\alpha_1, \alpha_2, \cdots, \alpha_4\right]^{\mathrm{T}}, \\ & \boldsymbol{\tau}=\left[\tau_1, \tau_2, \cdots, \tau_4\right]^{\mathrm{T}} . \\ & \end{array} $

(9)

其中: $A_{\text {ine, } i}$为电机转动惯量系数; $B_{\text {vis, } i}$为电机黏性摩擦系数; $r_{\text {win }, i}$为绞盘半径; $\alpha_i$为电机转角, $\dot{\boldsymbol{\alpha}}$和$\ddot{\boldsymbol{\alpha}}$分别为电机转角的速度矩阵和加速度矩阵; $\tau_i$为电机转矩; $i=1, 2, 3, 4$。

2.4 系统动力学模型

在2.2和2.3节的基础上，进一步建立系统动力学模型。绳长与电机转角矩阵的关系式表示如下：

$ \boldsymbol{r}_{\text {win }} \boldsymbol{\alpha}=\Delta l=(l(t)-l(t-1)), $

$ \boldsymbol{\alpha}=\boldsymbol{r}_{\text {win }}^{-1}(l(t)-l(t-1)) . $

(10)

其中：Δl为绳长变化量，l(t)为时刻t的绳索长度。联立式(8)—(10)，获得的系统动力学模型表示如下：

$ \begin{gathered} \boldsymbol{A}_{\mathrm{e}} \ddot{\boldsymbol{q}}_{M_0 O}^O+\boldsymbol{B}_{\mathrm{e}} \dot{\boldsymbol{q}}_{M_0 O}^O-\boldsymbol{G}_{\mathrm{e}}-\boldsymbol{F}_{\mathrm{de}}=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{\tau}, \\ \boldsymbol{A}_{\mathrm{e}}=\boldsymbol{r}_{\text {win }} \boldsymbol{A}+\boldsymbol{r}_{\text {win }}^{-1} \boldsymbol{J}^{\mathrm{T}} \boldsymbol{A}_{\text {ine }} \boldsymbol{J}, \end{gathered} $

(11)

$ \begin{gathered} \boldsymbol{B}_{\mathrm{e}}=\boldsymbol{r}_{\text {win }} \boldsymbol{B}+\boldsymbol{r}_{\text {win }}^{-1} \boldsymbol{J}^{\mathrm{T}} \boldsymbol{A}_{\text {inc }} \dot{\boldsymbol{J}}+\boldsymbol{r}_{\text {win }}^{-1} \boldsymbol{J}^{\mathrm{T}} \boldsymbol{B}_{\text {vis }} \boldsymbol{J}, \\ \boldsymbol{G}_{\mathrm{e}}=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{r}_{\text {win }}\left(\boldsymbol{J}^{\mathrm{T}}\right)^{-1} \boldsymbol{G}, \\ \boldsymbol{F}_{\mathrm{de}}=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{r}_{\text {win }}\left(\boldsymbol{J}^{\mathrm{T}}\right)^{-1} \boldsymbol{F}_{\mathrm{dis}} . \end{gathered} $

(12)

其中: $\boldsymbol{A}_{\mathrm{e}} 、\boldsymbol{B}_{\mathrm{e}} 、\boldsymbol{G}_{\mathrm{e}}$和$\boldsymbol{F}_{\mathrm{de}}$分别为引人电机动力学模型后的质量项、Coriolis项、重力项和扰动项。

2.5 运动平台受力分析

平面绳驱动并联机构无法承受法向力。运动平台受力如图 5所示, $\boldsymbol{G}_{\mathrm{ms}}$和$\boldsymbol{G}_{\mathrm{mc}}$为运动平台重力$\boldsymbol{G}$的分力; $f_{\text {up }}$和$f_{\text {down }}$分别为上绳索、下绳索的面内合力; $f_{\text {uc }}$和$f_{\text {us }}$为$f_{\text {up }}$的分力; $f_{\text {dc }}$和$f_{\text {ds }}$为$f_{\text {down }}$的分力。在斜向作业状态下, 重力及外部扰动力将产生垂直于运动平台表面的分力, 导致平台沿着分力方向移动。

为抵消重力分力及外部扰动力对运动平台的影响，只需在运动平台发生移动的初始时刻增大绳索张力，使其满足如下关系：

$ \left\{\begin{array}{l} \boldsymbol{f}_{\mathrm{us}}+\boldsymbol{f}_{\mathrm{ds}} \geqslant \boldsymbol{G}_{\mathrm{ms}}+\boldsymbol{f}_{\mathrm{dis}}, \\ \boldsymbol{f}_{\text {uc }}-\boldsymbol{f}_{\mathrm{de}}=\boldsymbol{G}_{\mathrm{mc}} . \end{array}\right. $

(13)

2.6 绳索张力优化

清洗机器人属于欠约束类型，无法采用常规的方法求解绳索张力。为了获得绳索张力的最优解，必须进行优化处理。

在作业时，绳索张力F应满足如下约束条件：

$ \begin{gathered} \boldsymbol{F}_{\text {min }} \leqslant \boldsymbol{F} \leqslant \boldsymbol{F}_{\text {max }}, \\ \boldsymbol{F}_{\text {min }}=\left[f_{\text {min }, 1}, f_{\text {min }, 2}, \cdots, f_{\text {min }, 4}\right]^{\mathrm{T}}, \end{gathered} $

(14)

其中：f_{min, i}为第i根绳索的最小允许张力；f_{max, i}为第i根绳索的最大允许张力，由电机的最大转矩和绳索的破断张力决定；i=1, 2, 3, 4。

进一步，将绳索张力分解，

$ \begin{gathered} \boldsymbol{F}=\boldsymbol{F}_{\text {par }}+\boldsymbol{F}_{\mathrm{gcn}}= \\ \boldsymbol{J}^{\mathrm{T}+}\left(\boldsymbol{A} \ddot{\boldsymbol{q}}_{M_0 O}^O+\boldsymbol{B} \boldsymbol{q}_{M_0 O}^O-\boldsymbol{G}-\boldsymbol{F}_{\mathrm{dis}}\right)+\boldsymbol{N}(\boldsymbol{J}) \lambda, \\ \boldsymbol{J}^{\mathrm{T}+}=\left(\boldsymbol{J} \boldsymbol{J}^{\mathrm{T}}\right)^{-1} \boldsymbol{J}, \\ \boldsymbol{F}_{\mathrm{par}}=\left[F_{\mathrm{par}, 1}, F_{\mathrm{par}, 2}, \cdots, F_{\mathrm{par}, 4}\right]^{\mathrm{T}} . \end{gathered} $

(15)

其中：F_par为张力特解项，F_gen为张力通解项；J^T+为伪逆矩阵；N (J)为J的一维零空间基底；λ为张力系数。

本文将相关力的最小二范数作为目标函数，优化模型表示如下^[21]：

$ \left\{\begin{array}{l} \min \left[\sum\limits_{i=1}^4\left(f_i-f_{\mathrm{ave}, i}\right)^2\right]^{1 /2}, \\ f_{\mathrm{ave}, i}=\left(f_{\min , i}+f_{\max , i}\right) /2, i=1, 2, 3, 4 ; \end{array}\right. $

(16)

$ \left\{\begin{array}{l} \text { s. t. } \boldsymbol{A} \ddot{\boldsymbol{q}}_{M_0 O}^O+\boldsymbol{B} \dot{\boldsymbol{q}}_{M_0 O}^O-\boldsymbol{G}-\boldsymbol{F}_{\mathrm{dis}}=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{F}, \\ \underline{\lambda} \leqslant \lambda \leqslant \bar{\lambda}, \\ \underline{\lambda}=\min\limits_{1 \leqslant i \leqslant 4}\left(\max\limits_{1 \leqslant i \leqslant 4}\left(\frac{f_{\min , i}-F_{\mathrm{par}, i}}{\boldsymbol{N}(\boldsymbol{J})_i}, \frac{f_{\max , i}-F_{\mathrm{par}, i}}{\boldsymbol{N}(\boldsymbol{J})_i}\right)\right), \\ \bar{\lambda}=\max\limits_{1 \leqslant i \leqslant 4}\left(\min\limits_{1 \leqslant i \leqslant 4}\left(\frac{f_{\min , i}-F_{\mathrm{par}, i}}{\mathbf{N}(\boldsymbol{J})_i}, \frac{f_{\max , i}-F_{\mathrm{par}, i}}{\mathbf{N}(\boldsymbol{J})_i}\right)\right) . \end{array}\right. $

(17)

其中：f_{ave, i}为张力均值，λ为λ的下限，λ为λ的上限。

3 机器人控制策略

通过设计合适的控制策略，控制电机收放绳索和滑台线性移动，实现运动平台的位姿控制。为便于分析，将机器人控制系统分为电机控制系统和滑台控制系统2部分。

3.1 电机控制策略

1) PI滑模控制策略。

清洗机器人的动力学模型存在不确定性，这将导致模型误差增大和控制系统闭环不稳定。由于滑模控制具有响应快、鲁棒性强等特点，因此被广泛用于不确定系统的控制^[16]。

将清洗机器人的动力学模型分解为已建模部分和建模误差，并分别设计控制律。在此基础上，引入比例项和积分项对滑模面进行修正，进而得到PI滑模控制器(PI sliding mode controller，PI-SMC)，有效提高系统性能和鲁棒性。以下是PI-SMC的设计过程。

根据式(11)可得系统动力学模型，表示如下：

$ \boldsymbol{A}_{\mathrm{s} y} \ddot{\boldsymbol{q}}_{\mathrm{M}_0 O}^O+\boldsymbol{B}_{\mathrm{sys}} \boldsymbol{\boldsymbol { q }}_{M_0 O}^O-\boldsymbol{G}_{\mathrm{sys}}-\boldsymbol{F}_{\mathrm{de}}=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{\tau} . $

(18)

其中A_sys、B_sys和G_sys为考虑建模误差后的动力学模型参数。然而，实际工程往往无法精准建立机器人真实的动力学模型，一般只能建立名义模型。A_sys、B_sys和G_sys表示如下：

$ \left\{\begin{array}{l} \boldsymbol{A}_{\mathrm{sys}}=\boldsymbol{A}_{\mathrm{e}}+\boldsymbol{E}_1, \\ \boldsymbol{B}_{\mathrm{sys}}=\boldsymbol{B}_{\mathrm{e}}+\boldsymbol{E}_2, \\ \boldsymbol{G}_{\mathrm{sys}}=\boldsymbol{G}_{\mathrm{e}}+\boldsymbol{E}_3 . \end{array}\right. $

(19)

其中：E₁、E₂和E₃分别为A_sys、B_sys和G_sys的建模误差，与由绳索的单向受力和悬链线效应等因素导致的模型不确定性相关。

运动平台的轨迹误差矢量E_tra表示如下：

$ \boldsymbol{E}_{\mathrm{tra}}=\boldsymbol{q}_{\mathrm{des}}-\boldsymbol{q} . $

(20)

其中：q_des为期望轨迹，q为实际轨迹。

进一步定义滑模函数，表示如下：

$ \begin{gathered} \boldsymbol{s}=\dot{\boldsymbol{E}}_{\text {tra }}+\boldsymbol{\varLambda} \boldsymbol{E}_{\mathrm{tra}}, \\ \boldsymbol{\varLambda}=\operatorname{diag}\left[\varLambda_1, \varLambda_2, \cdots, \varLambda_6\right] . \end{gathered} $

(21)

其中: $\dot{\boldsymbol{E}}_{\mathrm{tra}}$为速度误差矢量, $\boldsymbol{\varLambda}$为与滑模面相关的正定对角矩阵。

定义参考矢量为q_ref，表示如下：

$ \left\{\begin{array}{l} \boldsymbol{q}_{\text {ref }}=\int_0^t \boldsymbol{s} \mathrm{d} t+\boldsymbol{q}, \\ \dot{\boldsymbol{q}}_{\text {ref }}=\boldsymbol{s}+\dot{\boldsymbol{q}}, \\ \ddot{\boldsymbol{q}}_{\text {ref }}=\dot{\boldsymbol{s}}+\ddot{\boldsymbol{q}} . \end{array}\right. $

(22)

其中: $\dot{\boldsymbol{q}}_{\text {ref }}$和$\ddot{\boldsymbol{q}}_{\text {ref }}$分别为$\boldsymbol{q}_{\text {ref }}$的一阶和二阶导数。联立式(20)-(22) 得

$ \left\{\begin{array}{l} \boldsymbol{q}_{\text {ref }}=\boldsymbol{q}_{\text {des }}+\boldsymbol{\varLambda} \boldsymbol{E}_{\text {tra }}, \\ \ddot{\boldsymbol{q}}_{\text {ref }}=\ddot{\boldsymbol{q}}_{\mathrm{des}}+\boldsymbol{\varLambda} \dot{\boldsymbol{E}}_{\mathrm{tra}} . \end{array}\right. $

(23)

PI滑模控制矢量Q_con表示如下：

$ \begin{array}{c} \boldsymbol{Q}_{\mathrm{con}}=\boldsymbol{J}^{\mathrm{T}} \boldsymbol{\tau}=\boldsymbol{A}_{\mathrm{sys}} \ddot{\boldsymbol{q}}+\boldsymbol{B}_{\mathrm{sys}} \boldsymbol{q}-\boldsymbol{G}_{\mathrm{sys}}-\boldsymbol{F}_{\mathrm{de}}= \\ \boldsymbol{A}_{\mathrm{sys}}\left(\ddot{\boldsymbol{q}}_{\mathrm{des}}-\ddot{\boldsymbol{E}}_{\mathrm{tra}}\right)+\boldsymbol{B}_{\mathrm{sys}}\left(\dot{\boldsymbol{q}}_{\mathrm{des}}-\dot{\boldsymbol{E}}_{\mathrm{tra}}\right)-\boldsymbol{G}_{\mathrm{sys}}-\boldsymbol{F}_{\mathrm{de}}= \\ \boldsymbol{A}_{\mathrm{e}} \ddot{\boldsymbol{q}}_{\mathrm{ref}}+\boldsymbol{B}_{\mathrm{e}} \boldsymbol{q}_{\mathrm{ref}}-\boldsymbol{G}_{\mathrm{e}}+\boldsymbol{E}_{\mathrm{tot}}-\boldsymbol{A}_{\mathrm{sys}} \dot{\boldsymbol{s}}-\boldsymbol{B}_{\mathrm{sys}} \boldsymbol{s}-\boldsymbol{F}_{\mathrm{de}} . \end{array} $

(24)

其中: $\boldsymbol{E}_{\text {tot }}=\boldsymbol{E}_1 \ddot{\boldsymbol{q}}_{\text {ref }}+\boldsymbol{E}_2 \boldsymbol{q}_{\text {ref }}-\boldsymbol{E}_3$为动力学建模误差总和。$Q_{\text {con }}$的组成表示如下：

$ \begin{array}{c} & \boldsymbol{Q}_{\text {con }}=\boldsymbol{Q}_{\text {e }}+\boldsymbol{Q}_{\text {PI }}+Q_{\text {rob }}, \\ & \boldsymbol{Q}_{\mathrm{e}}=\boldsymbol{A}_{\mathrm{e}} \ddot{\boldsymbol{q}}_{\text {ref }}+\boldsymbol{B}_{\mathrm{e}} \dot{\boldsymbol{q}}_{\text {ref }}-\boldsymbol{G}_{\mathrm{e}}, \\ & \boldsymbol{Q}_{\mathrm{Pl}}=\boldsymbol{K}_{\mathrm{pro}} \boldsymbol{s}+\boldsymbol{K}_{\text {int }} \int_0^t \boldsymbol{s} \mathrm{d} t, \\ & \boldsymbol{Q}_{\text {rob }}=\boldsymbol{K}_{\text {rob }} \operatorname{sgn}(s), \\ & \boldsymbol{K}_{\text {pro }}=\operatorname{diag}\left[K_{\text {pro. } 1}, K_{\text {pro, } 2}, \cdots, K_{\text {pro, } 6}\right] \text {, } \\ & \boldsymbol{K}_{\text {int }}=\operatorname{diag}\left[K_{\text {int }, 1}, K_{\text {int }, 2}, \cdots, K_{\text {int, } 6}\right] \text {, } \\ & \boldsymbol{K}_{\text {rob }}=\operatorname{diag}\left[K_{\text {rob, } 1}, K_{\text {rob, } 2}, \cdots, K_{\text {rob, } 6}\right] \text {, } \\ & \end{array} $

$ \boldsymbol{D}=\boldsymbol{E}_{\mathrm{tot}}-\boldsymbol{F}_{\mathrm{de}} . $

(25)

其中: $\boldsymbol{Q}_{\mathrm{e}}$为名义模型控制率; $\boldsymbol{Q}_{\mathrm{PI}}$为比例-积分项; $\boldsymbol{Q}_{\text {rob }}$为鲁棒项, 用于补偿外界扰动对系统的影响; $\operatorname{sgn}(\boldsymbol{s})$为$\boldsymbol{s}$的符号函数; $\boldsymbol{K}_{\mathrm{pro}}$为比例增益矩阵, $\boldsymbol{K}_{\mathrm{int}}$为积分增益矩阵; $\boldsymbol{K}_{\mathrm{rob}}$为鲁棒增益矩阵, 系数选取应满足$K_{\mathrm{rob}, i} \geqslant\left|\boldsymbol{D}_i\right|, i=1, 2, \cdots, 6, \boldsymbol{D}$为建模误差与外界扰动项总和。$\boldsymbol{K}_{\mathrm{pro}} 、\boldsymbol{K}_{\mathrm{int}}$及$\boldsymbol{K}_{\mathrm{rob}}$均为正定对角矩阵。

以下是证明系统稳定性的过程。Lyapunov函数为

$ V_1=\frac{1}{2} \boldsymbol{s}^{\mathrm{T}} \boldsymbol{A}_{\mathrm{sys}} \boldsymbol{s}+\frac{1}{2}\left(\int_0^t \boldsymbol{s} \mathrm{d} t\right)^{\mathrm{T}} \boldsymbol{K}_{\text {int }} \int_0^t \boldsymbol{s} \mathrm{d} t . $

(26)

对式(26)求导，得到

$ \begin{gathered} \dot{V}_1=\boldsymbol{s}^{\mathrm{T}} \boldsymbol{A}_{\text {sys }} \dot{\boldsymbol{s}}+\frac{1}{2} \boldsymbol{s}^{\mathrm{T}} \dot{\boldsymbol{A}}_{\text {sys }} \boldsymbol{s}+\boldsymbol{s}^{\mathrm{T}} \boldsymbol{K}_{\text {int }} \int_0^t \boldsymbol{s} \mathrm{d} t= \\ \boldsymbol{s}^{\mathrm{T}}\left(\boldsymbol{A}_{\text {sys }} \dot{\boldsymbol{s}}+\boldsymbol{B}_{\text {sys }} \boldsymbol{s}+\boldsymbol{K}_{\text {int }} \int_0^t \boldsymbol{s} t\right)+\frac{1}{2} \boldsymbol{s}^{\mathrm{T}}\left(\dot{\boldsymbol{A}}_{\text {sys }}-2 \boldsymbol{B}_{\text {sys }}\right) \boldsymbol{s} . \end{gathered} $

(27)

根据动力学方程的斜对称性质，式(27)可简化为

$ \dot{V}_1=\boldsymbol{s}^{\mathrm{T}}\left(\boldsymbol{A}_{\text {sys }} \dot{\boldsymbol{s}}+\boldsymbol{B}_{\text {sys }} \boldsymbol{s}+\boldsymbol{K}_{\text {int }} \int_0^t \boldsymbol{s} \mathrm{d} t\right) . $

(28)

结合式(24)和(25)得

$ \begin{aligned} & \boldsymbol{A}_{\mathrm{sys}} \dot{\boldsymbol{s}}+\boldsymbol{B}_{\mathrm{sys}} \boldsymbol{s}+\boldsymbol{K}_{\mathrm{int}} \int_0^t \boldsymbol{s} \mathrm{d} t= \\ & -\boldsymbol{K}_{\mathrm{pro}} \boldsymbol{s}-\boldsymbol{K}_{\mathrm{rob}} \operatorname{sgn}(\boldsymbol{s})+\boldsymbol{D} . \end{aligned} $

(29)

将式(29)代入式(28)可得

$ \dot{V}_1=-\boldsymbol{s}^{\mathrm{T}} \boldsymbol{K}_{\text {pro }} \boldsymbol{s}-\boldsymbol{s}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{rob}} \operatorname{sgn}(\boldsymbol{s})+\boldsymbol{s}^{\mathrm{T}} \boldsymbol{D} . $

(30)

因为$K_{\mathrm{rob}, i} \geqslant\left|\boldsymbol{D}_i\right|$, 且$\boldsymbol{K}_{\mathrm{pro}}$为正定对角矩阵, 所以可得

$ \dot{V}_1 \leqslant-s^{\mathrm{T}} \boldsymbol{K}_{\text {pro }} \boldsymbol{s} \leqslant 0 . $

(31)

由LaSalle不变集定理^[22]可知, 当$t$趋于无穷时, $s$趋于0, 即$\boldsymbol{E}_{\mathrm{tra}}$和$\dot{\boldsymbol{E}}_{\mathrm{tra}}$也趋于0。因此, 控制系统是稳定的。

2) 模糊自适应PI滑模控制策略。

为提高系统性能，进一步引入自适应律和模糊控制规则进行优化。自适应律可自动调整PI-SMC的参数，从而适应外部环境和内部变化。同时，模糊控制规则可对PI-SMC的增益系数进行实时调整，进而提高系统的响应速度和控制性能。在此基础上，FAPI-SMC的控制策略如下。

假定D是有界矩阵，且满足如下关系^[23]：

$ \|\boldsymbol{D}\|<b_0+b_1\|\boldsymbol{q}\|+b_2\|\dot{\boldsymbol{q}}\|^2 . $

(32)

其中: $b_0 、b_1$和$b_2$为正常数, $\|\cdot\|$为向量的Euclid范数。

为提高FAPI-SMC的抗干扰特性，采用自适应项Q_{rob, a}替换Q_rob项，表示如下：

$ \begin{gathered} \boldsymbol{Q}_{\text {rob, a }}=\boldsymbol{\varsigma} \boldsymbol{s}+\left(\hat{b}_0+\hat{b}_1\|\boldsymbol{q}\|+\hat{b}_2\|\dot{\boldsymbol{q}}\|^2\right) \operatorname{sgn}(\boldsymbol{s}), \\ \boldsymbol{\varsigma}=\operatorname{diag}\left[\varsigma_1, \varsigma_2, \cdots, \varsigma_6\right] . \end{gathered} $

(33)

其中: $\boldsymbol{\varsigma}$为与自适应项相关的正定对角矩阵; $\hat{b}_i$是$b_i$的估计值, $\dot{\hat{b}}_i$为$\hat{b}_i$的一阶导数, 满足$\dot{\hat{b}}_0=$ $k_0\|\boldsymbol{s}\|, \dot{\hat{b}}_1=k_1\|\boldsymbol{s}\|\|\boldsymbol{q}\|, \quad \dot{\hat{b}}_2=k_2\|\boldsymbol{s}\|$ $\|\dot{\boldsymbol{q}}\|^2, k_i$为正常数, $i=0, 1, 2$。引人自适应项后, 控制矢量表示如下：

$ \boldsymbol{Q}_{\text {con, a }}=\boldsymbol{Q}_{\mathrm{e}}+\boldsymbol{Q}_{\text {PI }}+\boldsymbol{Q}_{\text {rob, a }} . $

(34)

进一步，将K_pro及K_int采用模糊控制规则确定：

$ \boldsymbol{Q}_{\text {P1, fuz }}=\boldsymbol{K}_{\text {pro, fuz }} \boldsymbol{s}+\boldsymbol{K}_{\text {int, fuz }} \int_0^t \boldsymbol{s} \mathrm{d} t . $

(35)

模糊控制器的相关参数如下：

(1) 为使运动平台快速平稳地完成轨迹跟踪，将控制器设计为双输入-双输出的结构。输入设定为E_tra和$\dot{\boldsymbol{E}}_{\mathrm{tra}}$，输出设定为$\boldsymbol{K}_{\mathrm{pro}, \mathrm{fuz}}$和$\boldsymbol{K}_{\mathrm{int} \text {, fuz }}$。

(2) 模糊控制器输入和输出语言选定为正大(positive big，PB)、正中(positive medium，PM)、正小(positive small，PS)、零(zero，ZO)、负小(negative small，NS)、负中(negative medium，NM)、负大(negative big，NB)。设定输入变量$\boldsymbol{E}_{\mathrm{tra}} \in(-0.01, 0.01), \dot{\boldsymbol{E}}_{\mathrm{tra}} \in(-1, 1)$; 输出变量$\boldsymbol{K}_{\text {pro, fuz }} \in(0, 300), \boldsymbol{K}_{\text {int }, \mathrm{fuz}} \in(0, 300)$。此外，本文选用灵敏度较强的三角形隶属度函数。

(3) 参考专家经验及整定原则^[24]，建立模糊控制规则如表 1和表 2所示。在模糊控制规则的基础上，采用重心法进行解模糊化。

针对该控制系统的稳定性进行证明，选取的Lyapunov函数表示如下：

$ \begin{gathered} V_2=\frac{1}{2} \boldsymbol{s}^{\mathrm{T}} \boldsymbol{A}_{\text {sys }} \boldsymbol{s}+ \\ \frac{1}{2}\left(\int_0^t \boldsymbol{s} \mathrm{d} t\right)^{\mathrm{T}} \boldsymbol{K}_{\text {int, fuz }} \int_0^t \boldsymbol{s} \mathrm{d} t+\sum\limits_{i=0}^2 \frac{1}{2 k_i}\left(b_i-\hat{b}_i\right)^2 . \end{gathered} $

(36)

结合式(27)，对式(36)求导可得：

$ \begin{gathered} \dot{V}_2=\boldsymbol{s}^{\mathrm{T}}\left(-\boldsymbol{K}_{\text {pro, fuz }} \boldsymbol{s}-\boldsymbol{Q}_{\text {rob, a }}+\boldsymbol{D}\right)+ \\ \sum\limits_{i=0}^2 \frac{1}{k_i}\left(\hat{b}_i-b_i\right) \dot{\hat{b}}_i= \\ \boldsymbol{s}^{\mathrm{T}}\left(-\boldsymbol{K}_{\text {pro, fuz }} \boldsymbol{s}-\left(\boldsymbol{\varsigma} \boldsymbol{s}+\left(\hat{b}_0+\hat{b}_1\|\boldsymbol{q}\|+\right.\right.\right. \\ \left.\left.\left.\hat{b}_2\|\dot{\boldsymbol{q}}\|^2\right) \operatorname{sgn}(\boldsymbol{s})\right)+\boldsymbol{D}\right)+ \\ \sum\limits_{i=0}^2 \frac{1}{k_i}\left(\hat{b}_i-b_i\right) \dot{\hat{b}}_i \leqslant \boldsymbol{s}^{\mathrm{T}}\left(-\left(\hat{b}_0+\hat{b}_1\|\boldsymbol{q}\|+\right.\right. \\ \left.\left.\hat{b}_2\|\dot{\boldsymbol{q}}\|^2\right) \operatorname{sgn}(\boldsymbol{s})+\boldsymbol{D}\right)+ \\ \sum\limits_{i=0}^2 \frac{1}{k_i}\left(\hat{b}_i-b_i\right) \dot{\hat{b}}_i \leqslant\|\boldsymbol{D}\|\|\boldsymbol{s}\|- \\ \left.\frac{1}{k_0}+\hat{b}_1\|\boldsymbol{q}\|+\hat{b}_2\|\boldsymbol{q}\|^2\right)\|\boldsymbol{s}\|+ \\ \left(\hat{b}_0-b_0\right) k_0\|\boldsymbol{s}\|+\frac{1}{k_1}\left(\hat{b}_1-b_1\right) k_1\|\boldsymbol{s}\|\|\boldsymbol{q}\|+ \\ \frac{1}{k_2}\left(\hat{b}_2-b_2\right) k_2\|\boldsymbol{s}\|\|\dot{\boldsymbol{q}}\|^2 \leqslant-\delta\|\boldsymbol{s}\| \leqslant 0 . \end{gathered} $

(37)

其中$\delta=b_0+b_1\|\boldsymbol{q}\|+b_2\|\dot{\boldsymbol{q}}\|^2-\|\boldsymbol{D}\|$。由式(32) 可得, $\delta>0$。根据有限时间稳定性理论, 轨迹误差可以收敛至 ^[25]。

(1) 扰动的不确定性会导致轨迹误差无法精准收敛至0，自适应参数也会随之变化，使系统变得更复杂。因此，本文采用死区技术改善自适应律^[26]：

$ \dot{\hat{b}}_0= \begin{cases}k_0\|{\boldsymbol{s}}\|, & \|{\boldsymbol{s}}\| \geqslant \varepsilon, \\ 0, & \|s\|<\varepsilon ;\end{cases} $

(38)

$ \dot{\hat{b}}_1= \begin{cases}k_1\|\boldsymbol{s}\|\|\boldsymbol{q}\|, & \|\boldsymbol{s}\| \geqslant \varepsilon, \\ 0, & \|\boldsymbol{s}\|<\varepsilon ;\end{cases} $

(39)

$ \dot{\hat{b}}_2= \begin{cases}k_2\|\boldsymbol{s}\|\|\dot{\boldsymbol{q}}\|^2, & \|\boldsymbol{s}\| \geqslant \varepsilon, \\ 0, & \|\boldsymbol{s}\|<\varepsilon .\end{cases} $

(40)

其中ε为死区系数。

(2) 由于符号函数$\operatorname{sgn}(\boldsymbol{s})$具有不连续性, 因此系统在运行过程中容易产生抖振。本文采用连续函数$\operatorname{sigmoid}(\boldsymbol{s})$代替$\operatorname{sgn}(\boldsymbol{s})$, 表示如下

$ \begin{gathered} \operatorname{sigmoid}\left(s_i\right)=\frac{2}{1+E_{\text {tra }, 1}^{-\varGamma s_i}}-1, \\ i=1, 2, \cdots, 6 . \end{gathered} $

(41)

其中：E_{tra, i}为E_tra的元素，s_i为s的元素，Γ为抖振改善系数。

在此基础上，电机控制逻辑如图 6所示。

3.2 滑台的控制策略

根据式(2)和(3)，滑台的期望行程表示如下：

$ \left\{\begin{array}{l} h_{\text {des, } 1}=y_{\mathrm{ope}}^{L_0}+d_{\mathrm{spr}}+z_{M_0 L_0}^{L_0} \tan \left(\theta_{M_0 L_0}^{L_0}\right)-\tan \left(\theta_{M_0 L_0}^{L_0}\right) H, \\ h_{\text {des, } 2}=y_{\mathrm{ope}}^{L_0}+d_{\mathrm{spr}}+z_{M_0 L_0}^{L_0} \tan \left(\theta_{M_0 L_0}^{L_0}\right) . \end{array}\right. $

(42)

通过滑台端部的激光测距仪实时测量实际行程h₁、h₂，进一步可得滑台的行程误差表示如下：

$ E_{h, i}=h_{\text {des }, i}-h_i, i=1, 2 . $

(43)

基于式(43)，采用传统的比例-积分-微分(proportional-integral-derivative，PID)控制器控制策略，相应的控制器输出表示如下：

$ \eta_i=\boldsymbol{K}_{\text {pro }, h} \boldsymbol{E}_{h, i}+\boldsymbol{K}_{\text {int }, h} \int_0^t \boldsymbol{E}_{h, i} \mathrm{~d} t+\boldsymbol{K}_{\text {der }, h} \dot{\boldsymbol{E}}_{h, i} . $

(44)

其中：η_i为第i个滑台的电机转速，i=1, 2。

4 仿真环境设定 4.1 船舶六自由度运动

本文采用软件AQWA模拟船舶在设定海况下的运动情况，所选船舶参数如表 3所示^[27]。

为验证清洗机器人适应不同海况的能力，分别采用规则波和不规则波作为船舶运动激励，相关参数分别如表 4和5所示。不规则波选用JONSWAP谱，其被广泛用于南海海域的海浪模拟，参数选取参照5级海况^[27]。

将船体模型导入软件AQWA并划分网格，如图 7所示。

水动力分析结果如图 8所示。将船舶运动数据作为运动激励输入至系统动力学模型，验证海浪激励下机器人的作业稳定性。

4.2 水枪反作用力

当水枪向作业面喷射清洗剂时，运动平台会受到M₀-xyz下的y方向水枪反作用力F_gun。本文采用Fluent进行流体仿真，模拟水射流冲击作业面的过程。在给定水射流速度和喷射距离的条件下，探究水射流是否可达作业面及F_gun大小，仿真参数如表 6所示。

将三维模型导入Fluent，进一步建立仿真环境并定义相关参数，流场模型如图 9a所示。图 9b为水枪流速迹线仿真，由图可知，在设定的条件下水射流可以喷射至作业面。水枪反作用力变化趋势如图 9c所示，仿真开始后，F_gun迅速增大，说明清洗剂可喷洒至作业面。将F_gun的数据输入清洗机器人动力学模型。

4.3 风载荷

挪威船级社(Det Norske Veritas，DNV)提出采用风压投影法计算风载荷^[28]，表示如下：

$ \left\{\begin{array}{l} F_{\text {wind }}=\varPi_{\text {sha }} P \vartheta_{\text {arc }} \sin \kappa_{\text {wind }}, \\ P=\frac{1}{2} \rho \varUpsilon_{\text {wind }}^2, \\ \varUpsilon_{\text {wind }}\left(z_{M O}^O\right)=\varUpsilon_0\left(1+5.73 \times 10^{-2} \times\right. \\ \left.\sqrt{1+0.148 \varUpsilon_0} \times \ln \frac{z_{M O}^O}{10}\right) . \end{array}\right. $

(45)

其中: $\vartheta_{\mathrm{arc}}$为结构物的迎风面面积; $\kappa_{\mathrm{wind}}$为风向与结构物迎风面的夹角; $\varPi_{\mathrm{sha}}$为形状系数, 与结构物形状有关; $P$为风压; $\rho$为空气密度, 一般采用$15{ }^{\circ} \mathrm{C}$时的空气密度, $\rho=1.226 \mathrm{~kg} /\mathrm{m}^3 ; \varUpsilon_{\text {wind }}$为设计风速; $r_0$为基本风速, 与海况及风力等级相关。

5 数值仿真

先采用MATLAB/Simulink建立机器人的动力学模型，再将运动平台的位姿参数输入控制器。滑台的控制方法较成熟，在仿真时不需要建立滑台的控制策略，可以直接将滑台的运动参数输入控制器。通过将不同控制器应用至相同作业场景，比较控制效果的优劣。

5.1 仿真参数设定

清洗机器人仿真所需参数如下：$H=15.0 \mathrm{~m}$, 作业面宽$W=20.0 \mathrm{~m}$; 运动平台高为$0.7 \mathrm{~m}$, 运动平台宽为$0.5 \mathrm{~m}$, 运动平台质量$m=50 \mathrm{~kg}$; 转动惯量$\ell_{x x}=1.052 \mathrm{~kg} \cdot \mathrm{m}^2, \ell_{y y}=2.052 \mathrm{~kg} \cdot \mathrm{m}^2, \ell_{z z}=$ $3.083 \mathrm{~kg} \cdot \mathrm{m}^2$; 重力加速度$g=-9.8 \mathrm{~N} /\mathrm{kg}$。电机参数如下: $r_{\text {win, } i}=0.1 \mathrm{~m}, A_{\text {ine }, i}=0.01 \mathrm{~kg} \cdot \mathrm{m}^2$, $B_{\text {vis }, i}=0.01 \mathrm{Nms} f_{\min , i}=20 \mathrm{~N}, f_{\max , i}=2000 \mathrm{~N}$。风载荷参数如下: $\kappa_{\text {wind }}=90^{\circ}, \varPi_{\mathrm{sha}}=1$; 基本风速的选取参照Beaufort风级表中的$9-10$级风速, $\varUpsilon_0=$ $25.72 \mathrm{~m} /\mathrm{s}$。$L_0-x y z$相对于$S-x y z$的位置矢量设定为$\boldsymbol{X}_{L_0 S}^S=\left[\begin{array}{lll}0 & 0 & 5\end{array}\right]^{\mathrm{T}} \; \mathrm{m}$。控制器参数的选取需要综合考虑电机转矩及轨迹跟踪精度，通过反复整定控制器参数，选定的参数如表 7所示。PI-SMC及FAPI-SMC的参数的选定原则如下：在满足张力限值的情况下，使轨迹跟踪精度达到最优。为与FAPI-SMC的响应速度进行对比，PID控制器参数的选定原则如下：参照专家经验及整定原则，在提高系统响应速度及满足张力限值的情况下，使轨迹跟踪精度达到最优^[24]。

船舶相对于$O-x y z$的初始位姿设定为$\boldsymbol{X}_{S O}^O(0)=\left[\begin{array}{lll}0 & 0 & 0\end{array}\right]^{\mathrm{T}}, \boldsymbol{\beta}_{S O}^O(0)=\left[\begin{array}{lll}0 & 0 & 0\end{array}\right]^{\mathrm{T}}$; 运动平台相对于$O-x y z$的初始位姿设定为$\boldsymbol{X}_{M_0 O}^O(0)=$ $\left[\begin{array}{lll}7.2 & 2.7 & 12.2\end{array}\right]^{\mathrm{T}}, \boldsymbol{\beta}_{M_0 O}^O(0)=\left[\begin{array}{lll}0 & 0 & 0\end{array}\right]^{\mathrm{T}}$; 运动平台在$O-x y z$下的初始速度及角速度为0。

运动平台在$L_0-x y z$下的运动轨迹设定为花瓣形状, 表示如下：

$ \left\{\begin{array}{l} x_{M_0 L_0}^{L_0}=\left(n_{\text {cir }}-\gamma_{\text {amp }} \sin \left(2 {\rm{\mathsf{π}}} v n_{\text {cir }} \tilde{\omega} t\right)\right) \cos (\tilde{\omega} t), \\ z_{M_0 L_0}^{L_0}=\left(n_{\text {cir }}-\gamma_{\text {amp }} \sin \left(2 {\rm{\mathsf{π}}} v n_{\text {cir }} \tilde{\omega} t\right)\right) \sin (\tilde{\omega} t), \\ \phi_{M_0 L_0}^{L_0}=0 . \end{array}\right. $

(46)

其中: 基圆半径$n_{\text {cir }}=2.0 \mathrm{~m}$, 正弦幅度$\gamma_{\text {amp }}=0.5 \mathrm{~m}$, 运动频率$v=0.478 \mathrm{~Hz}$, 运动角速度$\tilde{\omega}=2 {\rm{\mathsf{π}}} /T$, 运动周期$T=30 \mathrm{~s}$。进一步, 滑台的行程被设定如下：

$ \left\{\begin{array}{l} h_1=2(\sin t+1), \\ h_2=2(\cos t+1) . \end{array}\right. $

(47)

在实际作业中, 扰动不确定且不可忽略。为进一步验证控制器的鲁棒性, 对运动平台施加扰动力。由图 1可知, 运动平台的侧面及顶面为空心结构, 迎风面面积较小, $M_0-x y z$下的$x$和$z$方向海风扰动可以忽略。运动平台主要受到$y$方向扰动力, 由$F_{\mathrm{gun}}$和海风作用力$F_{\mathrm{wind}}$组成。同时, 为了验证控制器的抗扰性能, 在第$4 \mathrm{~s}$对其他方向施加扰动力, 表示如下：

$ \begin{aligned} & \left\{\begin{array}{l} f_{\text {dis, } x}=60 \cos (0.6(t-4)), \\ f_{\text {dis, } y}=F_{\text {gun }}+F_{\text {wind }}, \\ f_{\text {dis, } z}=60 \cos (0.6(t-4)) ; \end{array}\right. \\ & \left\{\begin{array}{l} w_{\text {dis, } x}=0.4 \cos (0.6(t-4)), \\ w_{\text {dis, } y}=0.6 \cos (0.6(t-4)), \\ w_{\text {dis, } z}=0.4 \cos (0.6(t-4)) . \end{array}\right. \end{aligned} $

(48)

其中: $f_{\text {dis, } x} 、f_{\text {dis, } y}$和$f_{\text {dis, } z}$分别为运动平台在$M_0-$ $x y z$下的$x 、y 、z$方向扰动力大小, $w_{\text {dis, } x} 、w_{\text {dis, } y}$和$w_{\mathrm{dis}, z}$分别为运动平台在$M_0-x y z$下的$x 、y 、z$方向扰动力矩大小。

外部扰动如图 10所示，y方向扰动力在120~160 N波动。

5.2 规则波激励下的曲面工况仿真分析

为进一步凸显FAPI-SMC的性能，将FAPI-SMC与PI-SMC和PID控制器进行对比，3种控制器的相关参数见表 7，规则波激励下的轨迹误差对比如图 11所示。

由图 11可知，当绳索张力在预设范围时，FAPI-SMC均可完成轨迹跟踪，验证清洗机器人动力学模型有效。与PID相比，FAPI-SMC的稳态误差明显减小；PID控制器虽然在一定程度上提高了系统响应速度，但是超调量较大，导致系统振荡较严重。与PI-SMC相比，FAPI-SMC在引入自适应律后，其轨迹跟踪的超调量及稳态误差均得到大幅改善。进一步可发现，所提出的FAPI-SMC具有更好的控制性能，该控制器具有更快的收敛速度和更高的轨迹跟踪精度。其中，位置稳态误差约±0.02 m，角度稳态误差约±0.02°。

针对y方向的轨迹跟踪情况进行对比分析，分别比较不同控制器的最大误差、稳态误差和上升时间，控制器性能比较结果如表 8所示。最大误差指的是误差首次减小至0后出现的最大值。FAPI-SMC的最大误差较PI-SMC和PID分别降低了0.171、0.154 m。稳态误差指系统趋于稳定时输出值与期望值之差。FAPI-SMC的稳态误差较PI-SMC和PID分别降低了0.073、0.101 m。上升时间指误差首次减小至0所需的时间，是系统响应速度的一种度量指标。FAPI-SMC的上升时间较PI-SMC及PID分别减少了1.44、5.34 s。通过比较控制器的综合性能可以发现，FAPI-SMC的最大误差减少量、响应速度和稳态性能均优于其他控制器。FAPI-SMC的最大误差较PID减少6%，响应速度较PID提升57%，稳态性能较PID提升3%。

图 12为FAPI-SMC控制下的速度及角速度误差。由图 12可知，在施加扰动后，速度及角速度误差瞬间增大，在FAPI-SMC作用下又趋于平稳，这表明FAPI-SMC具备良好的抗扰性能。其中，速度稳态误差约±0.1 m/s，角速度稳态误差约±8.7×10^-3 rad/s。

图 13为FAPI-SMC控制下的绳索张力及电机转矩。由图 13a可知，优化前绳索张力出现负值，与绳索单向受力特性不符，故需要进一步优化。图 13b—13c分别为优化后绳索张力及电机力矩，优化后绳索张力在设定范围内连续均匀变化。由于在第4 s施加外部扰动，因此FAPI-SMC需要调节输出张力以克服扰动对系统的影响，导致张力出现轻微波动。

FAPI-SMC的自适应参数变化趋势如图 14所示。由图 14可知，自适应参数在起始时刻迅速增大，表明FAPI-SMC具备快速响应能力。此外，当在第4 s施加了外部扰动时，自适应参数突然增大。在FAPI-SMC的作用下，自适应参数得以保持稳定，表明系统趋于稳定。

5.3 不规则波激励下的曲面工况仿真分析

针对不规则波激励下的清洗作业进行仿真分析，运动平台的初始位置矢量设定为X_M0O^O(0)=[7.30 3.05 11.85]^T，采用不规则波激励，其余参数保持不变，仿真结果如图 15所示。

仿真结果表明：FAPI-SMC的响应速度及控制精度均优于PI-SMC及PID，PI-SMC的控制精度略高于PID。与规则波激励下的仿真结果相似，FAPI-SMC的位置稳态误差约±0.02 m，角度稳态误差约±0.01°，这表明所提出的控制策略在不同形式的波浪激励下均具备良好的轨迹跟踪性能。

在曲面工况下，FAPI-SMC控制下的电机转矩如图 16所示。因为运动轨迹没有发生变化，所以电机转矩变化趋势与图 13c相似。此外，由图 16可知，在运动平台的运动轨迹保持不变的情况下，不同形式的外部激励对电机转矩影响较小，这表明FAPI-SMC具有良好的输出稳定性。

5.4 不规则波激励下的平面工况仿真分析

平面清洗是机器人日常作业的主要工况，为验证清洗机器人的变工况适应能力与平面作业精度，开展了平面作业仿真分析。仿真开始前，将滑台行程设置为0并保持不变，仿真采用不规则波激励，运动平台的初始位置矢量设定为X_M0O^O(0)=[7.30 0 11.85]^T，其他参数保持不变，仿真结果如图 17所示。

由图 17可知，在平面工况下，FAPI-SMC的轨迹跟踪精度和响应速度优于其他控制器。其中，FAPI-SMC控制下的位置稳态误差约±0.01 m，角度稳态误差约±0.01°。由图 17b可知，与曲面工况相比，PID的稳态误差降低约0.16 m，PI-SMC的稳态误差降低约0.19 m，FAPI-SMC的稳态误差变化较小，约0.01 m。结果表明工况变化将影响轨迹跟踪精度，同时也说明FAPI-SMC相较PI-SMC和PID可有效降低工况变化对轨迹跟踪精度造成的影响，具备良好的作业稳定性。

在平面工况下，FAPI-SMC控制下的电机转矩如图 18所示。由图 18可知，平面工况下的电机转矩大小和变化趋势与图 16相似，表明滑台运动对电机转矩的影响较小。

6 结论

针对船舶清洗作业的自动化、高效率需求，提出一种变结构空间绳驱动清洗机器人。在考虑船舶运动及外部扰动的作用下，基于Newton-Euler法推导了清洗机器人的系统动力学模型。在此基础上，分别提出控制器PI-SMC和FAPI-SMC。在不同的作业场景下，FAPI-SMC控制下的运动平台位置稳态误差为±0.02 m，角度稳态误差为±0.02°。与PI-SMC和PID相比，FAPI-SMC的最大误差分别减少8%、6%，响应速度提升18%、57%，稳态性能提升2%、3%。上述结果表明，本文所提出的控制策略在船舶运动下具备高精度、快响应性能，可为绳驱动机器人的实船应用提供理论支撑。在未来的研究工作中，将通过搭建实验样机进一步验证理论研究的正确性。