现代高端装备制造业对工件表面加工质量要求越来越高，基于工业机器人打磨已经成为高端制造业中打磨和抛光这类连续接触式作业的有力抓手，是解决人工作业的工作效率低和工作环境危害大问题以及提升工件加工质量的重要工具之一.

针对无先验模型且型面形状复杂多变的曲面工件机器人打磨，通常可以采用离线和在线打磨轨迹规划2种方式. 离线打磨轨迹规划主要通过传感器测量方式，预先获取工件模型相关先验数据，再根据先验数据进行打磨轨迹规划. Yang等^[1]利用触发式探头获取工件表面点云数据，采用B样条方法重构曲面工件进行机器人打磨轨迹规划. Li等^[2]借助结构光扫描仪，提出一种工件表面形状匹配方法，为离线轨迹规划获取先验数据. Ren等^[3]将工件点云数据拟合为曲面，根据曲面的曲率规划机器人打磨轨迹. Princely等^[4]通过相机获取工件二维形状数据，自动生成机器人去毛刺轨迹. Zhen等^[5] 结合3D点云和机器视觉技术，提出一种机器人自适应打磨轨迹规划方法. 针对多种类、小批量且模型未知的复杂曲面工件，由于标定和安装误差的影响，通过测量获取的先验数据与工件真实模型具有一定误差^[6]，对此类工件的离线打磨轨迹规划较为耗时，降低了打磨效率.

在线打磨轨迹规划主要通过传感器实时获取打磨工具与工件的接触状态信息，进而预测曲面的几何特性，在线规划或修正打磨轨迹. Kiguchi等^[7]提出一种打磨力方向在线调整方法，使接触力方向跟踪工件轮廓曲面法线方向. Jeon等^[8]提出一种视觉和力传感器信息相融合的物体轮廓跟踪方法，在线调整机器人轨迹. Li等^[9]提出一种工具与工件接触状态实时估计方法，用于跟踪工件曲面法线. Oba等^[10]在工件先验数据未知且无力反馈的条件下，提出一种姿态和法向抛光力同步控制方法，使抛光工具的姿态自动适应工件曲面形状. 上述在线轨迹规划方法在针对空间复杂曲面工件打磨时，难以较准确地刻画空间曲面工件的几何特性，且具有较大的力跟踪控制误差.

本研究针对工业机器人柔顺打磨作业需求及法向打磨力需要稳定且准确可控的特点，以无先验模型曲面为打磨对象，开展机器人打磨主动自适应在线轨迹预测研究. 首先，设计一种柔顺浮动打磨力控末端执行器，将导纳控制和比例积分微分（PID）控制应用于末端执行器，以获取稳定准确的法向打磨力. 其次，提出一种无先验模型曲面的机器人打磨主动自适应在线轨迹预测方法. 最后，通过虚拟样机联合仿真分析，以及无先验模型圆台曲面工件自适应在线轨迹预测实验，验证所提出方法的有效性.

机器人打磨系统主要由工业机器人与末端执行器构成，柔顺打磨通过采用主动力控制的末端执行器完成，末端执行器的位置和姿态（简称位姿）由工业机器人控制. 根据 Preston 材料去除理论，材料去除率与工具对工件的法向打磨力正相关^[11]，为了实现材料去除率可控、保证打磨质量，要求法向打磨力稳定、准确和可控；同时，为了实现此类受限运动，须对工业机器人运动轨迹进行预测或修正^[12]. 仅通过机器人与刚性末端执行器的运动控制及被动力控制技术开展打磨作业，难以满足上述需求. 带有主动力控制的柔顺浮动末端执行器，可以将机器人打磨系统对打磨工具的力控制和位姿控制解耦，进而增加机器人打磨系统力控制带宽，提高系统动态响应性能^[13].

设计柔顺浮动打磨力控末端执行器，如图1所示. 该末端执行器以端面打磨方式进行打磨，采用主动力控制进行打磨力控制，并结合机器人进行位姿控制，实现柔顺打磨，具有机器人打磨系统的力控制和位姿控制解耦的特点.

所设计的柔顺浮动打磨力控末端执行器的关键组件主要由以下3部分构成. 1) 法向打磨力控制组件：采用主动力控制（由力传感器2和磁栅编码器4作为传感元件），通过控制音圈电机11的输出力，实现法向打磨力的稳定、准确控制，减小由刚性接触引起的过切和欠切. 2) 浮动运动组件：由直线导轨模组（直线导轨动子5、直线导轨6）实现移动平台7的往复直线运动，保证打磨主轴10的浮动运动精度. 3) 打磨作业组件：打磨主轴10带动打磨工具8高速旋转打磨作业，可以灵活适配不同打磨工具，适应多场景的打磨作业.

音圈电机选用MOTICONT公司LVCM-051-064-02及810系列驱动器，最大持续输出力为68.2 N，PWM控制频率为50000 Hz；打磨主轴选用NAKANISHI公司BMS-4020电主轴，最高转速为20000 r/min；力传感器选用ATI公司F/T Delta IP65六维力传感器，力采集精度为1/8 N，采样频率为3000 Hz；磁栅编码器选用RLS公司HILIN系列直线式增量磁栅编码器，最高分辨率为0.1 μm；打磨工具采用球形红刚玉磨头.

机器人系统柔顺打磨的关键是依据打磨工况主动自适应进行打磨工具浮动调节. 本设计通过安装在音圈电机支架12下方的音圈电机11，驱动由安装在移动平台7夹板之间的主轴夹具固连的打磨主轴10，在安装于导轨支撑板3内侧的直线导轨6与安装在移动平台外侧的直线导轨动子组成的导轨模组上往复直线运动，带动被打磨主轴前端夹头夹持打磨工具8上下浮动，进而对法向打磨力稳定准确控制，加之打磨主轴10带动打磨工具8高速旋转，实现打磨工具端面打磨的浮动打磨作业；力传感器2用于测量打磨中工件对工具的反作用力，连同磁栅编码器4为主动自适应在线轨迹与姿态预测提供接触和位置状态反馈信息.

在ADAMS软件中构建机器人打磨系统虚拟样机（由柔顺浮动打磨力控末端执行器、六自由度串联工业机器人和曲面工件组成），如图2所示，将打磨系统虚拟样机导出为Matlab/Simulink仿真控制系统.

Matlab/Simulink仿真控制系统虚拟样机的输入分别为机器人关节变量q、音圈电机输入力F_i；输出分别为所设计末端执行器打磨主轴的轴向位置变化量d_L、工具坐标系{T}的位置P、绕基础坐标系{B}的RPY角（Roll-Pitch-Yaw）以及在机器人工具坐标系{T}表示的工件对工具的反作用力 ${}^T{{\boldsymbol{\tau}}_{\text{c}}} $. 反作用力表达式如下：

(1) $ {}^T{{\boldsymbol{\tau}}_{\text{c}}} = {[{}^TF_{\text{c}}^x,{}^TF_{\text{c}}^y,{}^TF_{\text{c}}^z,{}^T\tau _{\text{c}}^x,{}^T\tau _{\text{c}}^y,{}^T\tau _{\text{c}}^z]^{\text{T}}}. $

在机器人打磨系统作业时，法向打磨力稳定性是影响工件打磨质量的重要因素^[14-15]. 工具与工件之间的接触关系一定程度上可以反映出工件的空间几何特性^[7].

针对机器人打磨力控制问题，PID控制具有良好的力跟踪特性^[16-17]，然而难以同时考虑位置控制精度，且在打磨系统颤振时，PID控制器对控制目标的跟踪不具有柔顺性，因此，PID控制难以满足复杂工况机器人柔顺打磨要求. Hogan^[18]提出阻抗控制方法，将机器人末端执行器和环境接触过程的阻抗特性视为阻抗模型进行控制. 基于位置的阻抗控制也称为导纳控制，其对力的跟踪具有一定柔顺性^[19-20]，且其内环包含稳定的位置控制器，对颤振具有一定抑制能力，本研究将其作为法向打磨力控制方式.

将柔顺浮动打磨力控末端执行器和位置控制器合并为位置控制系统，并简化为如图3(a)所示的理想二阶阻抗模型. 图中，M_e、D_e和K_e分别表示理想二阶阻抗模型的惯性系数、阻尼系数和刚度系数，X_e为期望打磨量，F_n为法向打磨力. 将导纳控制器简化为理想二阶导纳模型，如图3(b)所示.

将打磨中法向打磨力F_n与期望力F_d之间的偏差F_e作为导纳控制器的输入，通过导纳模型的传递，计算出对位置控制系统输入的补偿量X_r，其数学模型如下：

(2) $ {F_{\text{e}}} = {F_{\text{d}}} - {F_{\text{n}}}\text{，} $

(3) $ M{\ddot X_{\text{r}}}+D{\dot X_{\text{r}}}+K{X_{\text{r}}} = {F_{\text{e}}} . $

式中：M、D和K分别为导纳模型的惯性系数、阻尼系数和刚度系数.

对式(3)进行拉普拉斯变换，得到位置补偿量X_r和力偏差F_e的关系为

(4) $ {X_{\text{r}}} = \frac{{{F_{\text{e}}}}}{{M{s^2}+Ds+K}}. $

式中：s为拉普拉斯算子.

将位置补偿X_r与期望位置X₀及位置变化量d_L叠加，得到位置控制器的输入量X_i：

(5) $ {X_{\text{i}}} = {X_{\text{r}}}+{X_0} - {d_{\text{L}}}. $

将X_i作用于末端执行器与位置控制器构成的位置控制系统，与导纳控制器一同实现对期望力F_d的跟踪控制，法向打磨力控制框图如图4所示.

定义打磨主轴的轴线方向从主轴后端沿主轴轴线指向打磨工具. 在对曲面工件打磨时，工业机器人实时调整末端执行器的位姿. 当主轴轴线方向与重力加速度方向不同向时，如图5所示，须计算当前状态下对末端执行器音圈电机的输出力进行重力补偿，用于平衡重力作用，从而实现控制器对法向打磨力的准确控制.

在虚拟样机中实时获取工具坐标系{T}在基础坐标系{B}中的RPY角，即横滚角α、偏航角β以及俯仰角γ，根据RPY角与姿态矩阵的转换关系，求得工具坐标系{T}的姿态矩阵为

(6) $ {\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {{\text{c}}\alpha \;{\text{c}}\beta }&\;{{\text{c}}\alpha \;{\text{s}}\beta \;{\text{s}}\gamma - \;\;{\text{s}}\alpha \;{\text{c}}\gamma }&\;{{\text{c}}\alpha \;{\text{s}}\beta \;{\text{c}}\gamma +\;\;{\text{s}}\alpha \;{\text{s}}\gamma } \\ {{\text{s}}\alpha \;{\text{c}}\beta }&{\;{\text{s}}\alpha \;{\text{s}}\beta \;{\text{s}}\gamma +\;\;{\text{c}}\alpha \;{\text{c}}\gamma }&{\;{\text{s}}\alpha \;{\text{s}}\beta \;{\text{c}}\gamma - \;\;{\text{c}}\alpha \;{\text{s}}\gamma } \\ { - {\text{s}}\beta }&\;{{\text{c}}\beta \;{\text{s}}\gamma }&\;{{\text{c}}\beta \;{\text{c}}\gamma } \end{array}} \right]. $

式中： ${\text{c}}\;\omega $表示 $\cos\; \omega $， ${\text{s}}\;\omega $表示 $\sin\; \omega $， $\omega $可取 $\alpha $、 $\; \beta$或 $\gamma $. 工具坐标系{T}的姿态矩阵在基础坐标系{B}中可表示为

(7) $ {\boldsymbol{R}}_T^B = \left[ {\begin{array}{*{20}{c}} {{n_x}}&{{o_x}}&{{a_x}} \\ {{n_y}}&{{o_y}}&{{a_y}} \\ {{n_z}}&{{o_z}}&{{a_z}} \end{array}} \right] . $

式中：n_x、n_y和n_z表示坐标系{T}的x轴单位方向向量在坐标系{B}下的方向余弦，o_x、o_y和o_z表示坐标系{T}的y轴单位方向向量在坐标系{B}下的方向余弦，a_x、a_y和a_z表示坐标系{T}的z轴单位方向向量在坐标系{B}下的方向余弦.

由打磨主轴的轴线方向与工具坐标系{T}的z轴重合，可知主轴轴线方向单位矢量 $ {}^B{\boldsymbol{a}} = {\left[ {{a_x},{a_y},{a_z}} \right]^{\text{T}}} $，再根据基础坐标系{B}中的单位重力矢量 $ {\boldsymbol{g}}={\left[0，0，-1\right]}^{\text{T}} $，可以得出单位重力矢量 $ {\boldsymbol{g}} $与主轴轴线方向单位矢量 $ {}^B{\boldsymbol{a}} $之间的夹角 $\varphi $为

(8) $ \varphi = {\text{arc}}\cos\; \left( {\frac{{{}^B{\boldsymbol{a}} \cdot {\boldsymbol{g}}}}{{\left\| {{}^B{\boldsymbol{a}}} \right\| \cdot \left\| {\boldsymbol{g}} \right\|}}} \right). $

式中： $\varphi $为0 °~180 °. 由式(8)可知，夹角 $\varphi $仅与主轴轴向方向矢量和重力单位矢量有关，与主轴以及移动平台的质心位置无关，因此，依据夹角 $\varphi $和主轴以及移动平台质量M_p，可以得到音圈电机输出力的补偿量：

(9) $ {F_{\text{g}}} = {M_{\text{p}}} \cos\; \varphi .$

叠加补偿量F_g与位置控制器的输出力F_co，得到音圈电机的输出力F_i：

(10) $ {F_{\text{i}}} = {F_{{\text{co}}}}+{F_{\text{g}}} .$

进而实现法向打磨力控制的重力补偿，以消除主轴和移动平台自身重力对控制器性能的影响. 叠加重力补偿的法向打磨力控制框图如图6所示.

针对无先验模型的曲面工件，一般通过获取工件曲面的几何信息，求解曲面与基准平面的截交线，将截交线作为机器人打磨的待跟踪曲线，进而进行打磨轨迹规划^[21]. 为了提高对无先验模型曲面工件打磨的主动自适应性，本研究提出一种无先验模型曲面的机器人打磨主动自适应在线轨迹预测方法.

进行如下设定.

1) 通过在机器人工作空间内设置一个起始点A、一个终止点B和一个法向量n_p，构建如图7所示的有向进给平面. 图中，点A和点B分别为预规划打磨路径的起始点和终止点，法向量n_p为过点A和点B平面的法线，一般情况下，法向量n_p平行于坐标系{B}平面xoy. 平面进给方向可以根据工件待加工曲面的加工条件设定.

2) 设定TCP（tool central point）为音圈电机处于零位时，固连在打磨工具的末端点，并随工业机器人与打磨主轴浮动改变自身位姿；规定工具坐标系{T}的原点为TCP，工具坐标系{T}的z轴为打磨主轴的轴线方向，z轴正方向为工具接近工件的方向，x轴始终沿机器人末端法兰轴线向内.

3) 机器人系统打磨时工件对工具的反作用力F_c被分解到工具坐标系{T}的x轴、y轴和z轴上，均由安装在工业机器人法兰上的六维力传感器测得，如图8所示.

所提出在线轨迹预测方法的控制框图如图9所示. 该方法主要包括两部分：1) 基于反作用力分解的在线姿态预测（曲面法线跟踪预测）；2) 基于有向进给平面截交线的在线位置预测. 主要思路如下：通过分解工件对工具的反作用力，对曲面法线进行主动自适应预测，使末端执行器主轴的轴线方向主动跟踪预测的曲面法线，并根据预测的曲面法线求解TCP处曲面切平面与有向进给平面的截交线，进一步预测出下一时刻TCP应达到的目标位姿；根据机器人逆向运动学，解得机器人下一时刻的目标关节变量，将目标关节变量作用于机器人关节驱动，从而使机器人主动自适应调整自身及末端执行器的位姿，实现无先验模型曲面的机器人打磨主动自适应在线轨迹预测. 所提出方法对于无先验模型的复杂多变曲面或平面均具有适用性.

根据Hertz接触理论，打磨工具与曲面工件的接触应力在曲面表面法线方向上最大^[22]，而打磨时工具对工件的切削力和摩擦力相对较小^[23]，可以近似认为工件对工具的反作用力方向主要沿着曲面的法线方向. 为了消除轻微打磨颤振对轨迹预测的影响，通过对力矢量信号做低通滤波，将其分解到工具坐标系{T}的坐标轴上，即

(11) $ {}^T{{\boldsymbol{f}}_{\text{c}}} = {[{F_x},{F_y},{F_z}]^{\text{T}}} . $

再将该力矢量变换到基础坐标系{B}上，坐标变换为

(12) $ {}^B{{\boldsymbol{f}}_{\text{c}}} = {\boldsymbol{R}}_T^B \;{}^T{{\boldsymbol{f}}_{\text{c}}} . $

式中： $ {\boldsymbol{R}}_T^B $为工具坐标系{T}相对于基础坐标系{B}的旋转变换矩阵，可以通过式(6)、(7)求得.

根据式(12)得到工件对工具的反作用力矢量在基础坐标系{B}下的表示 ${}^B{{\boldsymbol{f}}_{\text{c}}}$，其方向由自身单位向量描述，将 ${}^B{{\boldsymbol{f}}_{\text{c}}}$单位化后，得到基础坐标系{B}下该反作用力方向的单位向量 ${}^B{{\boldsymbol{f}}_{\text{e}}}$，下文称其为反作用力方向向量.

根据打磨工具的进给方向和工具与工件的接触状态，构建进给方向与曲面法线关系，如图10所示. 图中， $ {S_{\text{e}}} $为工件曲面； $ {{\boldsymbol{V}}_{\text{n}}} $为曲面真实法线单位向量，由曲面自身特性决定；Σ_⊥为过点TCP且与V_t−1垂直的平面； ${{\boldsymbol{V}}_{t - 1}}$为上一时刻打磨进给速度单位向量，可以通过计算当前时刻TCP位置 $ {}^B{{\boldsymbol{P}}_t} $（可以通过虚拟样机输出的位置矢量P获得）与上一时刻TCP位置 $ {}^B{{\boldsymbol{P}}_{t - 1}} $之差，并将其单位化获得，即

(13) $ {{\boldsymbol{V}}_{t - 1}} = \frac{{{}^B{{\boldsymbol{P}}_t} - {}^B{{\boldsymbol{P}}_{t - 1}}}}{{\left\| {{}^B{{\boldsymbol{P}}_t} - {}^B{{\boldsymbol{P}}_{t - 1}}} \right\|}}. $

在工具旋转作业或进给运动时，若直接利用工件对工具的反作用力方向向量 ${}^B{{\boldsymbol{f}}_{\text{e}}}$近似表示曲面法线方向向量 ${{\boldsymbol{V}}_{\text{n}}}$，则切削力和摩擦力会导致轨迹预测的曲面法线方向与真曲面法线方向产生偏差，如图10所示. 将反作用力方向向量 ${}^B{{\boldsymbol{f}}_{\text{e}}}$分解为沿工具进给反方向的进给阻力向量 ${}^B{{\boldsymbol{f}}_{\text{v}}}$和与进给方向垂直的预测曲面法线向量 ${}^B{{\boldsymbol{u}}^{\text{n}}}$，并使用 ${}^B{{\boldsymbol{u}}^{\text{n}}}$近似表示真实曲面法线方向 ${{\boldsymbol{V}}_{\text{n}}}$，从而减小摩擦和切削作用对曲面法线预测的影响. 由图10可以得到

(14) $ {}^B{{\boldsymbol{f}}_{\text{v}}} = \left( {{}^B{{\boldsymbol{f}}_{\text{e}}} \cdot {{\boldsymbol{V}}_{t - 1}}} \right){{\boldsymbol{V}}_{t - 1}} , $

(15) $ {}^B{{\boldsymbol{u}}^{\text{n}}} = {}^B{{\boldsymbol{f}}_{\text{e}}} - {}^B{{\boldsymbol{f}}_{\text{v}}} . $

将下一时刻打磨系统工具坐标系{T}的z轴取为近似曲面法线的负方向 $- {}^B{{\boldsymbol{u}}^{\text{n}}}$，即

(16) $ {}^B{\boldsymbol{a}} = - {}^B{{\boldsymbol{u}}^{\text{n}}}. $

为了实现姿态的稳定连续预测，使基础坐标系{B}下工具坐标系{T}的y轴正方向 ${}^B{\boldsymbol{o}}$始终与基础坐标系{B}的yoz面平行，即 $ {}^B{o_x} \equiv 0 $，通过工具坐标系{T}的z轴与y轴的垂直关系 ${}^B{\boldsymbol{o}} \cdot {}^B{\boldsymbol{a}} = 0$，可以得到

(17) $ \left. {\begin{array}{*{20}{c}} {{}^B{a_y}{}^B{o_y}+{}^B{a_z}{}^B{o_z} = 0}, \\ {{}^Bo_y^2+{}^Bo_z^2 = 1}. \end{array}} \right\} $

进一步可以得到

(18) $ \left. {\begin{array}{*{20}{c}} {{}^B{o_y} = \pm \sqrt {\dfrac{{{}^Ba_y^2}}{{{}^Ba_z^2+{}^Ba_y^2}}} } \\ {{}^B{o_z} = \mp \sqrt {\dfrac{{{}^Ba_z^2}}{{{}^Ba_z^2+{}^Ba_y^2}}} } \end{array}} \right\} $

即有2个解

(19) $ {{}^B}{\boldsymbol{o}} = {\left[ {0,\sqrt {\frac{{{{}^B}a_y^2}}{{{{}^B}a_z^2+{{}^B}a_y^2}}} , - \sqrt {\frac{{{{}^B}a_z^2}}{{{{}^B}a_z^2+{{}^B}a_y^2}}} } \right]^{\text{T}}} , \quad$

(20) $ {}^B{\boldsymbol{o}}' = = {\left[ {0, - \sqrt {\frac{{{}^Ba_y^2}}{{{}^Ba_z^2+{}^Ba_y^2}}} ,\sqrt {\frac{{{}^Ba_z^2}}{{{}^Ba_z^2+{}^Ba_y^2}}} } \right]^{\text{T}}} . $

为了确保姿态稳定连续变化不发生突变，根据式(6)、(7)获取当前时刻工具坐标系{T}的姿态矩阵 ${}^B{{\boldsymbol{o}}_t}$，并分别求解其与式(18)中2个解的欧氏距离，取距离最近的解作为下一时刻的 ${}^B{\boldsymbol{o}}$. 根据式(17)~(20)求得的 ${}^B{\boldsymbol{a}}$、 $ {}^B{\boldsymbol{o}} $，可以进一步求得工具坐标系{T}的x轴正方向为

(21) $ {}^B{\boldsymbol{n}} = {}^B{\boldsymbol{o}} \times {}^B{\boldsymbol{a}} . $

通过上述计算，得到下一时刻工具坐标系{T}在基础坐标系{B}下的姿态矩阵为

(22) $ {}^B{{\boldsymbol{R}}_{t+1}} = \left[ {{}^B{\boldsymbol{n}}},\;{{}^B{\boldsymbol{o}}},\;{{}^B{\boldsymbol{a}}} \right] .$

将式(22)所示的姿态矩阵作为工具坐标系{T}的目标姿态，通过预测的曲面法线向量 ${}^B{{\boldsymbol{u}}^{\text{n}}}$解得目标姿态矩阵 $ {}^B{{\boldsymbol{R}}_{t+1}} $即可完成曲面法线跟踪预测.

设定有向进给平面及方向如图11所示，该平面的一般方程为

(23) $ {\varSigma}_{0}：U\left(x-{x}_{0}\right)+V\left(y-{y}_{0}\right) +W\left(z-{z}_{0}\right)=0 . $

式中： ${\left[ {U,\;V,\;W} \right]^{\text{T}}}$为平面的法向量n_p， $\left( {{x_0},{y_0},{z_0}} \right)$表示平面上一点的坐标. 由式(15)计算出当前时刻基础坐标系{B}下曲面的法线方向 ${}^B{{\boldsymbol{u}}^{\text{n}}}$，根据当前工具的位置 ${}^B{{\boldsymbol{P}}_t}$，构建出当前TCP处曲面的切平面为

(24) $ {\varSigma }_{{t}}：{}^{B}u{}_{x}^{\text{n}}\left(x-{}^{0}P{}_{x}\right)+{}^{B}u{}_{y}^{\text{n}}\left(y-{}^{0}P{}_{y}\right) +{}^{B}u{}_{z}^{\text{n}}\left(z-{}^{0}P{}_{z}\right)=0 . $

进一步可以求得有向进给平面与曲面切平面截交线的方向向量为

(25) $ {}^B{\boldsymbol{l}} = {\left[ {U,V,W} \right]^{\text{T}}} \times {}^B{{\boldsymbol{u}}^{\text{n}}}. $

由于工件轮廓曲线可以通过多个微元直线近似^[24]，在工具当前位置 $ {}^B{{\boldsymbol{P}}_t} $处沿 ${}^B{\boldsymbol{l}}$按步长因子 $\eta $平移，预测得到下一时刻工具的位置向量 $ {}^B{{\boldsymbol{P}}_{t+1}} $为

(26) $ {}^B{{\boldsymbol{P}}_{t+1}} = {}^B{{\boldsymbol{P}}_t}+{}^B{\boldsymbol{l}} \cdot \eta . $

将 $ {}^B{{\boldsymbol{P}}_{t+1}} $作为下一时刻工具的位置向量.

至此，基于反作用力分解预测的姿态矩阵 ${}^B{{\boldsymbol{R}}_{t+1}}$和位置向量 $ {}^B{{\boldsymbol{P}}_{t+1}} $均已获得，将预测的下一时刻目标齐次矩阵构建为

(27) $ {}^B{{\boldsymbol{T}}_{t+1}} = \left[ {\begin{array}{*{20}{c}} {{}^B{{\boldsymbol{R}}_{t+1}}}&{{}^B{{\boldsymbol{P}}_{t+1}}} \\ {\boldsymbol{0}}&{\boldsymbol{1}} \end{array}} \right] . $

在机器人打磨系统虚拟样机中，分别在导纳控制和PID控制2种力控制器作用下，将正弦曲面作为模型曲面，开展该无先验曲面机器人打磨主动自适应在线轨迹预测仿真分析，验证所提出方法的有效性及法向打磨力控制的稳定性. 同时开展无先验模型圆台曲面工件自适应在线轨迹预测实验，进一步验证所提出方法的有效性.

在Matlab/Simulink与ADAMS的联合仿真分析中，将法向打磨力的期望值设置为30 N，导纳控制器内部导纳模型的惯性系数、阻尼系数和刚度系数分别设置为23、1300和8.25，在ADAMS中设置工具与工件的接触刚度为1.0×10⁴ N/m，穿透深度为1×10⁻⁵ m，最大阻尼系数3×10⁴ ${\text{N}} \cdot {\text{s/m}}$；法向打磨力控制精度设置为最大力跟踪误差不超过1 N，平移步长因子 $\eta $设置为0.04. 根据打磨工具是否对工件施加稳定的法向打磨力作为各阶段规划的时间划分依据. 将仿真过程的第0~10 s规划为工具接触工件并逐步建立稳定法向打磨力的过程，称之为接触阶段（其中，第0~5 s为工具接触工件过程，第6 ~10 s为建立稳定法向打磨力过程）；将第10~50 s规划为主动自适应在线轨迹与姿态预测过程，称之为在线预测阶段. 在力控制器分别为导纳控制和PID控制作用下，所提出方法对期望位置在x、y和z方向的最大预测误差 $\Delta x$、 $ \Delta y $和 $\Delta z$，以及最大综合位置误差 $ \Delta d $如表1中位置预测所示；对期望RPY角的预测误差 $\Delta \alpha $、 $ \Delta \beta $和 $ \Delta \gamma $，以及最大综合姿态误差（法线跟踪误差） $ \Delta \delta $如表1中姿态预测（法线跟踪预测）所示；法向打磨力的最大跟踪误差 $ \Delta F_{\text{n}}^{\max } $、误差平均值 $ {\Delta {\bar F_{\text{n}}}} $和误差方差 $ \Delta F_{\text{n}}^{\rm{s}} $如表1中法向力跟踪所示.

仿真结果表明，所提出方法在稳定法向打磨力控制器下的位置跟踪较准确，利用位置综合误差^[25]量化评估位置预测的准确性，表达式为

(28) $ \Delta d = \left[{{{\left( {\Delta {x}} \right)}^2}+{{\left( {\Delta {y}} \right)}^2}+{{\left( {\Delta {z}} \right)}^2}}\right]^{1/2} . $

经计算，在采用导纳控制和PID控制时，最大综合位置预测误差 $\Delta d$如表1与图12所示，分别为0.506 、0.677 mm. 因此，所提出方法在稳定力控制器下的位置预测准确性较好.

因曲面工件模型未知，仅通过预测，无法精确跟踪无先验模型曲面，但小于5°的预测角度误差，对打磨质量的影响较小，亦可满足实际应用的需要^[26-27]. 采用机器人打磨仿真过程中虚拟样机输出的RPY角与沿曲面法线规划的期望轨迹RPY角的综合姿态误差，评估所提出方法对曲面法线的跟踪性能. 在力控制器为导纳控制和PID控制时，在线预测阶段的仿真结果如图13、14所示. 图中， $\phi $和 $\Delta \phi $分别表示角度度量和角度误差度量. 仿真结果表明，所提出方法能够较好地跟踪期望轨迹的RPY角，即能够较好地实现姿态预测.

利用姿态综合误差^[25]量化评估曲面法线跟踪性能，定义曲面法线跟踪误差为

(29) $ \Delta \delta = \left[ {{{\left( {\Delta {\alpha }} \right)}^2}+{{\left( {\Delta {\beta }} \right)}^2}+{{\left( {\Delta {\gamma }} \right)}^2}}\right]^{1/2} .$

经计算，所提出方法的曲面法线跟踪误差如表1与图15所示. 在采用导纳控制和PID控制时，最大综合姿态预测误差 $\Delta \delta $分别为4.912°和4.946°. 因此，曲面法线跟踪误差均小于5°，可以满足无先验模型曲面跟踪机器人打磨应用需要^[26-27]，表明所提出方法对无先验模型曲面法线的主动自适应跟踪性能较好，验证了所提出方法对无先验模型曲面法线跟踪的有效性.

为了提高工件打磨质量，在打磨时不仅要使法向打磨力的方向主动自适应跟踪曲面的法向，同时要求其稳定可控. 导纳控制和PID控制作用下对应的反作用力曲线如图16所示. 图中，F_x、F_y和F_z分别表示反作用力在工具坐标系{T}下的x、y和z轴的分量，F_x与法向打磨力F_n相等，因此可以用F_x量化表征法向打磨力.

在导纳控制下，TCP在到达正弦曲面顶端时，会产生一定位置跳动，如图16(a)所示；法向打磨力响应曲线的上升时间为4.24 s；在打磨过程中，系统具有一定打磨力跟踪误差，如表1所示.

在PID控制下，TCP在正弦曲面的任意位置均能实现平稳过渡；但在工具与工件接触初期，因碰撞产生较大的冲击和颤振，法向打磨力瞬间达到58.054 N，如图16(b)所示，而颤振幅值会随PID控制调整过程逐渐减小，最终稳定收敛；法向打磨力响应曲线的上升时间为7.31 s；在法向打磨力控制过程稳定后，法向打磨力可以较准确地跟踪期望力，力跟踪误差如表1所示.

由上述结果可知，所提出方法在导纳控制下，法向打磨力的跟踪误差方差较大，力跟踪的准确性较差，但颤振抑制能力较好，且最大力跟踪误差较小，法向打磨力响应曲线的上升时间更短，快速性更好；而在PID控制下，力跟踪的准确性较好，但颤振抑制能力较差. 因此，所提出方法在法向打磨力控制器为导纳控制和PID控制下均适用，具有较好的适应性.

为了进一步验证所提出在线轨迹预测方法的有效性，搭建如图17所示的机器人打磨系统实验平台，开展实物平台实验研究.

ABB IRB-4600/2.05机器人及IRC控制器和示教器作为机器人打磨系统的位姿控制部件；BECKHOFF公司CX-260354嵌入式控制器和力控末端执行器构成打磨力控制部件；工件为无先验几何模型的圆台曲面铸铁工件，如图18所示. CX-260354控制器通过EtherCAT网络协议获取磁栅编码器和力传感器信号，经过计算后，将预测的TCP位姿及力控制信号发送到工业机器人和音圈电机驱动器，开展主动自适应在线轨迹预测.

有向进给平面如图18所示，法向打磨力采用PID控制，为了防止实验过程中法向打磨力过大引起打磨主轴卡滞，设定期望法向打磨力为15 N. 平台实验如图19所示，工件用夹具固定，在工件表面设定起始点A(1000, 20, 700)和终止点B(1105, −85, 680)，单位为mm；法向量设定为 ${n_{\rm{p}}} = \left( { - {{\sqrt 2 }}/{2}, - {{\sqrt 2 }}/{2},0} \right)$，根据A、B和n_p构建有向进给平面，并设实验时间为30 s，从点A开始，在CX-260354嵌入式控制器运行所提出方法对应程序，使机器人打磨系统完成主动自适应姿态跟踪以及位置预测，记录末端执行器在预测时的RPY预测角轨迹、位置预测轨迹以及法向打磨力数据，并绘制曲线图，分别如图20~22所示.

由实验结果可知，法向打磨力跟踪误差最大值为4.037 N，平均值为−0.086 N. 实验表明，所提出的在线轨迹预测方法能够沿设定的有向进给平面实现该无先验模型工件的RPY角预测和位置预测以及曲面法线跟踪，验证了所提出方法的有效性.

为了提高工业机器人打磨的作业柔顺性及无先验模型曲面工件的机器人打磨自适应性，设计了柔顺浮动打磨力控末端执行器，并提出了无先验模型曲面的机器人打磨主动自适应在线轨迹预测方法. 仿真和实验验证了所提出方法对无先验曲面的在线轨迹预测与法线跟踪的有效性，同时其可以实现对法向打磨力的准确跟踪控制，且方法易实现，实时性高. 所提出方法可以为无先验模型曲面工件的机器人打磨提供一种在线轨迹预测方法，具有一定的参考价值. 后续在所提出方法基础上，可以针对机器人打磨过程中复杂时变非线性动力学耦合与不确定性扰动以及机器人打磨作业颤振对在线轨迹预测的影响开展进一步深入研究，提高在强扰动情况下机器人打磨主动自适应在线轨迹预测的鲁棒性.