0 引言

发动机、叶片和模具等零部件表面复杂、精度要求较高, 目前主要的加工方式需要大量人力完成主要的或辅助的工作。据调查, 一个熟练工人加工一个汽轮机叶片需要2 h才能完成, 但废品率高达15%, 在我国废品率高达30%。对于提高磨削工作效率和质量, 行业内主要有两种解决方案。一种是采用高精度的数控机床, 一种采是工业机器人。两种方案据调查参数如表 1所示。

由表 1可知, 与工业机器人相比, 数控机床在刚度、定位精度等方面有很大的优势。但是工业机器人在工作空间、工件夹持等方面表现出很大的优势, 且价格只有数控机床的10%。随着机器人在打磨抛光等接触式作业的广泛应用, 对机器人接触式的柔顺力控制要求越来越高。经查阅总结, 工业机器人柔顺力跟踪控制的主要问题有如下两个方面: (1)工业机器人主要是串联机构, 其操作精度和重复定位精度都较差。(2)因为待加工表面的复杂性和精度要求, 很难对机器人和作业对象进行精确建模。

应用工业机器人进行打磨的主要方法是通过离线示教方式对机器人示教, 通过机器人编程记录每个点。为了保证磨削精度, 一条路径需要示教500~600个点才能完成, 对于较复杂的表面, 为达到较高的精度可能需要更多的示教点。要完整完成一个零件的打磨, 通常需要数天的示教调试, 并且容易出错, 对工作人员的经验要求较高^[1-3]。

目前, 机械臂末端力的控制方法主要分为两种。一种是机械结构设计。通过加入一些柔顺弹性部件, 例如弹簧、柔顺关节、机械阻尼调节器等。但这类部件和机械臂结构关联甚密, 虽然在实际操作过程中保证了机械臂具有一定的柔性, 但是对机械臂的控制精度和响应速度造成很大影响。另一种是力矩控制。机械臂通过检测作业过程中各个关节和末端位置的受力情况, 通过改进控制算法来提高机械臂的柔顺性能, 同时也保证机械臂的控制精度^[4-6]。Borenstein等人^[7]提出将检测到的电机输出扭矩与给定的参考力矩相比较的碰撞检测方法, 但是在实际过程中机械臂不能实现与外界接触环境的柔顺碰撞。Yamad等人^[8]提出当末端接触力超出可接受的安全范围时, 通过紧急制动保护作业过程的安全。但是当机械臂在达到安全范围边界时, 如果速度、加速度不为零, 可能会对作业目标或机械臂造成损坏。Albus等人^[9]提出一种碰撞路径规划方案保证末端力可以稳定在期望范围内, 但在实际作业过程中, 目标位置的误差和不确定性, 末端力和关节角度不再具有直接关系。Wilson等人^[10-11]提出了通过卡尔曼滤波器、PID等方法的力位混合控制策略, 在平面运动过程中可以保持±5 N的作用力。但实际作业过程中, 绝大多数的作业面不是平面的, 实用性不强。

上述几种机械臂柔顺力控制方法都是建立在作业对象完全已知的条件下。在实际生产过程中, 要做到作业对象完全已知是非常困难的。如果每次作业前进行精确测量会极大降低作业效率。对于流水线作业对象, 由于其本身的生产误差、装夹误差, 对控制系统都会带来未知信息, 这些未知信息在作业过程中会被放大并产生难以预计的结果。

本研究提出了一种基于模型预测控制和曲面预测相结合的方法对机械臂进行力位混合控制。针对未知的连续作业表面, 根据机械臂当前关节角、末端受力等测得的信息对曲面下一时刻的坐标进行预测, 然后根据逆运动学求得期望的位置信息, 运用动态矩阵控制算法(dynamic matrix control algorithm, DMC)进行滚动优化并对机械臂的输入量进行求解^[12-13]。然后返回施加控制输入后的关节角、末端力等信息进行下一个循环控制, 实现了针对未知的连续作业曲面不用进行点点示教就能完成磨削工作, 并具有较好的磨削质量, 大大节省了磨削工作前期的编程示教时间。

1 机械臂模型以及工作环境模型 1.1 受限机械臂动力学分析

因为工业过程领域较难建立控制对象的模型, 预测控制算法通常通过模型辨识建立预测模型。而机械臂物理参数在生产加工时是确定的, 因此可通过运动学、动力学分析对机械臂建立预测模型。常用的三自由度机械臂模型如图 1所示。

机械臂的控制和其动力学问题是密切相关的, 与外部环境相接触的三自由度机械臂动力学方程可表示为:

$ \mathit{\boldsymbol{D}}\left( \mathit{\boldsymbol{q}} \right)\mathit{\boldsymbol{\ddot q}} + \mathit{\boldsymbol{C}}\left( {\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right)\mathit{\boldsymbol{\dot q}} + \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{q}} \right)\mathit{\boldsymbol{q}} - {\mathit{\boldsymbol{J}}^{\rm{T}}}\mathit{\boldsymbol{F}} = \mathit{\boldsymbol{\tau }}, $

(1)

式中: q, $ \boldsymbol{\dot q}$, $\boldsymbol{ \ddot q}$∈Rⁿ, 分别为机械臂关节的角度、角速度矢量和角加速度矢量; τ∈Rⁿ为各个关节控制力矩的输入矢量; J为机器人学中的雅克比矩阵; F为机械臂末端受力矢量; D(q)∈R^n×n, C(q, $\boldsymbol{\dot q} $)∈R^n×n, G(q)∈R^n×n分别为机器人动力学中的惯量矩阵、哥氏力和离心力矢量、重力矢量。由式(1)可以看出, 机械臂动力学方程是一个强耦合、非线性的时变方程, 不能采用基于线性叠加原理的模型预测控制算法直接进行控制, 因此采取小区间线性化处理再运用模型预测控制方法。为了保证整个控制过程的正确性, 需要对控制模型进行实时优化, 保证算法建立在有效的模型基础上。

根据机械臂动力学方程式(1), 选定x=[x₁, x₂]^T=[q, $\boldsymbol{\dot q} $]^T作为系统的状态矢量, 则系统的状态方程可表示为:

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot x}}}_1} = {\mathit{\boldsymbol{x}}_2}\\ {{\mathit{\boldsymbol{\dot x}}}_2} = - \mathit{\boldsymbol{D}}{\left( \mathit{\boldsymbol{q}} \right)^{ - 1}}\mathit{\boldsymbol{C}}\left( {\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right) + \mathit{\boldsymbol{B\tau }}\left( t \right) - \mathit{\boldsymbol{D}}{\left( \mathit{\boldsymbol{q}} \right)^{ - 1}}\mathit{\boldsymbol{G}}。\\ \mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{x}}_1} \end{array} \right. $

(2)

利用泰勒公式对式(2)进行离散化, 采样周期记为T。离散化后的状态矢量

$ \mathit{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_1}}\\ {{\mathit{\boldsymbol{x}}_2}} \end{array}} \right] = \left[ \begin{array}{l} {\mathit{\boldsymbol{x}}_1} + \mathit{\boldsymbol{T}}{{\mathit{\boldsymbol{\dot x}}}_1} + \frac{{{\mathit{\boldsymbol{T}}^2}}}{2}{{\mathit{\boldsymbol{\ddot x}}}_1}\\ {\mathit{\boldsymbol{x}}_2} + \mathit{\boldsymbol{T}}{{\mathit{\boldsymbol{\dot x}}}_2} \end{array} \right]。$

(3)

因为机械臂在实际跟踪过程中位移和速度都较小, 因此假定机械臂在短区间内为线性。根据式(2)(3)可得三自由度机械臂离散化后的状态空间模型为:

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{x}}\left( {k + 1} \right) = \boldsymbol{x}\left( {t = \left( {k + 1} \right)T} \right) = \mathit{\boldsymbol{Ax}}\left( k \right) + \mathit{\boldsymbol{B\tau }}\left( k \right) + \mathit{\boldsymbol{G}}\\ \mathit{\boldsymbol{y}}\left( k \right) = \mathit{\boldsymbol{Hx}}\left( k \right),k = 0, \cdots ,N - 1 \end{array} \right., $

(4)

式中:$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_3}}&{\mathit{\boldsymbol{T}}{\mathit{\boldsymbol{I}}_3}}\\ 0_{3\times3}&{{\mathit{\boldsymbol{I}}_3}} \end{array}} \right],\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {\frac{{{\mathit{\boldsymbol{T}}^2}}}{2}\mathit{\boldsymbol{M}}{{\left( \mathit{\boldsymbol{q}} \right)}^{ - 1}}}\\ {\mathit{\boldsymbol{TM}}{{\left( \mathit{\boldsymbol{q}} \right)}^{ - 1}}} \end{array}} \right],$$ \mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{{\mathit{\boldsymbol{T}}^2}}}{2}\mathit{\boldsymbol{D}}{{\left( \mathit{\boldsymbol{q}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{C}}\left( {\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right)\mathit{\boldsymbol{\dot q}} + \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{q}} \right)} \right)}\\ { - \mathit{\boldsymbol{TD}}{{\left( \mathit{\boldsymbol{q}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{C}}\left( {\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right)\mathit{\boldsymbol{\dot q}} + \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{q}} \right)} \right)} \end{array}} \right],$ $\mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_3}}&{{0_{3 \times 3}}} \end{array}} \right]。$

式中:A、B、G在对每个期望点的跟踪过程中为恒值, 即系统为时不变系统。因为跟踪层的预测模型是在模型(4)的基础上进行求解, 为保证预测模型的准确性, 在对每个跟踪点跟踪过程结束之后, 通过机械臂当前时刻的状态对模型进行优化。

1.2 未知作业曲面预测方程

在机械臂作业过程中, 如果作业面不是平面, 传统的控制方法很难保证对机械臂下个时刻的状态进行很好的控制, 尤其是机械臂这种高精度的刚体对另一个高刚度对象作业时, 轻微的误差可能造成很大的应力, 严重时可能导致机械臂本体或零件损坏。传统的控制方式大多根据反馈力进行控制调整, 当曲面斜率较大, 其控制精度较差, 因此, 在作业过程中对作业对象下一时刻的状态进行预测估计十分必要。

机械臂末端和作业对象接触时, 机械臂和力传感器的刚度均远远大于作业对象的刚度, 接触力产生的微形变主要发生在工作曲面上, 因此主要考虑工作曲面产生的微形变。机械臂末端和外部环境的接触简图如图 2所示。机械臂末端受力F_E=0时, 末端和作业表面的接触位置为x_E0。那么作业对象受力后的微小形变Δx_E和末端受力为零时的接触位置x_E0之间的关系可以表示为:

$ \Delta {x_{\rm{E}}} = {x_{\rm{E}}} - {x_{{\rm{E0}}}} = {F_{\rm{E}}}/{K_{\rm{E}}}, $

式中:x_E为末端受力后的位置; F_E为末端所受的外力；K_E为作业对象的刚度。由此可得出机械臂在曲面上的局部放大图如图 3所示。

根据力传感器测得的力和机械臂返回的末端位置计算曲面预测模型, 并引入曲面预测误差e_XE进行反馈调节。曲面模型预测可表示为:

$ \left\{ \begin{array}{l} {X_{{\rm{E}},k + 1}} = {X_{{\rm{E}},k}} + {K_k}\Delta y + {e_{{X_{\rm{E}}}}}\\ {K_k} = {K_{k - 1}} + \Delta {K_{k - 1}}\\ \Delta {K_{k - 1}} = {K_{k - 1}} - {K_{k - 2}}\\ {e_{{X_{\rm{E}}}}} = {X_{\rm{R}}} - {F_{k - 1}}/{K_{\rm{E}}} - {X_{{\rm{E}},k}} \end{array} \right., $

(5)

式中:X_{E, k+1}为工作曲面k+1时刻的预测位置; K_k为当前时刻的预测斜率; Δy表示机械臂末端位置控制空间中每个时刻移动的距离; e_XE为过去k-2时刻预测的曲面信息和过去k-1时刻经过传感器求得的实际曲面信息的差值; ΔK为工作曲面k-1时刻位置处的斜率变化率; X_R为机械臂末端的实际位置。根据当前时刻的位置、斜率和误差计算下个时刻的位置, 然后输入到控制器中进行求解并控制机械臂的移动。当机械臂移动至下个位置时, 根据力传感器、末端位置等状态信息对斜率、斜率的变化率、误差等参数进行更新。

2 控制器设计 2.1 轨迹规划层控制器设计

机械臂是一个高刚度、高精度的集合体, 在实际作业过程中, 在力空间内轻微的位置变化都会产生很大的作用力。过大的作用力极有可能影响生产的质量甚至损坏零件和设备。因此, 要保证控制的快速性和精准性, 需要严苛限制控制过程中的超调问题。本研究基于模型预测思想, 对机械臂力控制空间的控制问题设计如下控制算法^[14]。

根据DMC算法设计预测模型和反馈校正环节, 其稳态模型为:

$ \Delta {\mathit{\boldsymbol{F}}_{{\rm{ss\_TTSS}},k + 1}} = {\mathit{\boldsymbol{K}}_{{\rm{E}},k}}\Delta {\mathit{\boldsymbol{x}}_{{\rm{w}},k + 1}},k = 1. \cdots ,N, $

式中:ΔF_{ss_TTSS, k+1}为机械臂末端预测输出力; Δx_{w, k+1}为机械臂末端k+1时刻在力空间上的微位移变化量; K_{E, k}为1.2节求得的作业曲面的刚度。根据上一时刻的稳态输出力、输入的控制序列以及当前时刻的控制输入检测值求解当前时刻的稳态输出力

$ {\mathit{\boldsymbol{F}}_{{\rm{ss}},k}} = {\mathit{\boldsymbol{F}}_{{\rm{ss}},k - 1}} + {\mathit{\boldsymbol{K}}_{{\rm{E}},k}}\left( {{\mathit{\boldsymbol{x}}_{{\rm{ss\_measure}},k}} - {\mathit{\boldsymbol{x}}_{{\rm{ss}},k - 1}}} \right), $

式中: F_{ss, k-1}为系统k-1时刻的稳态输出力; x_{ss_measure, k}为k时刻机械臂末端在力空间中的测量位置。根据力传感器检测到的当前时刻的输出力F(k), 与k时刻的输出预测值F_{ss_TTSS, k}相减得出误差

$ \mathit{\boldsymbol{e}}\left( k \right) = \mathit{\boldsymbol{F}}\left( k \right) - {\mathit{\boldsymbol{F}}_{{\rm{ss\_TTSS}},k}}, $

使用误差对机械臂下一时刻的预测输出力进行校正:

$ {\mathit{\boldsymbol{F}}_{{\rm{ss\_cor}},k + 1}} = {\mathit{\boldsymbol{F}}_{{\rm{ss\_TTSS}},k + 1}} + \mathit{\boldsymbol{e}}\left( k \right), $

然后根据期望值和校正值之间的差值, 通过稳态模型计算出对应的最优控制输入Δx_{ss_calc}。根据控制响应速度的需求, 可以对控制的节拍参数N进行调整^[15-17]。

机械臂当前时刻的控制增量

$ \Delta {\mathit{\boldsymbol{x}}_{\rm{u}}} = \Delta {\mathit{\boldsymbol{x}}_{{\rm{ss}}\_{\rm{calc}}}}/N, $

力空间内实际期望轨迹为

$ {\mathit{\boldsymbol{x}}_{\rm{z}}}\left( {k + 1} \right) = {\mathit{\boldsymbol{x}}_{\rm{z}}}\left( k \right) + \Delta {\mathit{\boldsymbol{x}}_{\rm{u}}} + \Delta {\mathit{\boldsymbol{x}}_{\rm{E}}}, $

式中: x_z(k+1)为k+1时刻末端的期望位置; x_z(k)为k时刻机械臂末端的实际位置; Δx_u为节拍参数控制后的控制增量; Δx_E为工作曲面在力空间中k时刻至k+1时刻的位置增量。

根据求得的实际期望轨迹, 通过机械臂求逆解得出机械臂3个关节角期望轨迹。求逆解方法如下:

对于图 1所示的机械臂, 其结构简图如图 4所示, OA、OB分别代表其两个连杆。假定已知机械臂末端位置为(x, y, z), 根据机器人学定义θ₁、θ₂、θ₃分别表示机械臂三个关节的关节角。则根据图 4可求得:

$ OC = \sqrt {{x^2} + {y^2}} , $

$ {\theta ''_2} = \arctan \left( {\frac{z}{{OC}}} \right), $

$ OB = \sqrt {O{C^2} + {z^2}} 。$

$ {\theta '_2} = \arccos \left( {\frac{{OB}}{{2OA}}} \right), $

$ {\theta _2} = {\theta '_2} + {\theta ''_2}, $

$ {\theta _3} = - 2{\theta _2}, $

$ {\theta _1} = \arctan \left( {\frac{y}{x}} \right), $

将求得的机械臂关节角作为跟踪层的目标输出设定值进行跟踪算法的设计。

2.2 轨迹跟踪层控制器设计

采用局部线性化简化模型, 并对模型失配的问题进行调整^[18]。因此对控制过程进行了如下设计:

(1) 预测模型

根据1.1节求出的机械臂状态空间模型可得:

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{x}}\left( 0 \right) = {\mathit{\boldsymbol{x}}_{\rm{s}}}\\ \mathit{\boldsymbol{x}}\left( {k + 1} \right) = \mathit{\boldsymbol{Ax}}\left( k \right) + \mathit{\boldsymbol{B\tau }}\left( k \right) + \mathit{\boldsymbol{G}}\\ \mathit{\boldsymbol{y}}\left( {k + 1} \right) = \mathit{\boldsymbol{Hx}}\left( {k + 1} \right)\\ \mathit{\boldsymbol{\tau }}\left( k \right) = \mathit{\boldsymbol{\tau }}\left( {k - 1} \right) + \Delta \mathit{\boldsymbol{\tau }}\left( k \right)\\ \mathit{\boldsymbol{\tau }}\left( 0 \right) = {\mathit{\boldsymbol{\tau }}_{\rm{s}}} \end{array} \right., $

式中: x_s为系统控制过程当前时刻的状态; τ_s为系统上次控制过程的输入值。

记$ {\boldsymbol{y}_{\rm PM}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{x}\left( 1 \right)}\\ \vdots \\ {\boldsymbol{x}\left( N \right)} \end{array}} \right]$, ${\boldsymbol{u}_\tau } = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{\tau} \left( 0 \right)}\\ \vdots \\ {\boldsymbol{\tau} \left( {N - 1} \right)} \end{array}} \right] $, 则系统的预测模型可表示为:

$ {\mathit{\boldsymbol{y}}_{{\rm{PM}}}}\left( k \right) = \mathit{\boldsymbol{L}}{\mathit{\boldsymbol{x}}_{\rm{s}}}\left( k \right) + \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{u}}_\tau }\left( k \right) + \mathit{\boldsymbol{O}}\left( k \right), $

式中: $\mathit{\boldsymbol{L}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{A}}\\ \vdots \\ {{\mathit{\boldsymbol{A}}^N}} \end{array}} \right],\mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{B}}& \cdots &0\\ \vdots &{}& \vdots \\ {{\mathit{\boldsymbol{A}}^{N - 1}}\mathit{\boldsymbol{B}}}& \cdots &\mathit{\boldsymbol{B}} \end{array}} \right], $ $\mathit{\boldsymbol{O}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{G}}\\ {\mathit{\boldsymbol{AG}} + \mathit{\boldsymbol{G}}}\\ {{\mathit{\boldsymbol{A}}^{N - 1}}\mathit{\boldsymbol{G}} + \cdots + \mathit{\boldsymbol{G}}} \end{array}} \right]。$

(2) 滚动优化

因为矩阵M不一定为方阵, 因此M不一定可逆, 那么输入的求解可采用极小化如下性能指标的方式逼近设定值

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{J}}\left( {{\mathit{\boldsymbol{y}}_{{\rm{PM}}}}\left( t \right),{\mathit{\boldsymbol{u}}_\tau }\left( {t - 1} \right)} \right) = \sum\limits_{j = 1}^N {\left\| {{\mathit{\boldsymbol{y}}_{{\rm{PM}}}}\left( {t + i\left| t \right.} \right) - } \right.} }\\ {\left. {{\mathit{\boldsymbol{y}}_{{\rm{ref}}}}\left( {t + i\left| t \right.} \right)} \right\|_Q^2 + \sum\limits_{j = 1}^{{N_{\rm{c}}}} {\left\| {\Delta {\mathit{\boldsymbol{u}}_\tau }\left( {t + i\left| t \right.} \right)} \right\|_\mathit{\boldsymbol{R}}^2} ,} \end{array} $

式中:Q、R为设定的权重矩阵; y_ref为2.1节求得的2个关节的期望轨迹; 等式右边第一项是针对系统的实际输出和1.2节中位置空间中求解出的参考输出进行约束; 等式右边第二项是针对系统保证完成任务目标的情况下优化控制增量, 保证其可以平稳运行, 不会出现突变。

通过极值必要条件$ \frac{{{\rm d}\boldsymbol{J}\left( k \right)}}{{{\rm d}\Delta {\boldsymbol{u}_\tau }\left( k \right)}} = 0$求得:

$ \Delta {\mathit{\boldsymbol{u}}_\tau }\left( k \right) = {\left( {{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{QM}} + \mathit{\boldsymbol{R}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{Q}}\left( {{\mathit{\boldsymbol{y}}_{{\rm{ref}}}} - {\mathit{\boldsymbol{y}}_{{\rm{p0}}}}} \right), $

式中: y_ref为期望轨迹; y_p0=Lx_s+N。在跟踪过程中, 把控制量Δu(k)=S^TΔu_τ(k)作为每次实际的输入量输入到控制系统中, 其中S^T=[1 0 … 0]。

(3) 反馈校正

把求得的控制增量作为实际的输入量输入到系统中, 得到未来时刻的输出位置预测值

$ {\mathit{\boldsymbol{y}}_{{\rm{PM}}}}\left( k \right) = \mathit{\boldsymbol{L}}{\mathit{\boldsymbol{x}}_{\rm{s}}}\left( k \right) + \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{u}}_\tau }\left( k \right) + \mathit{\boldsymbol{N}}\left( k \right)。$

因为实际作业过程中可能存在模型适配问题、环境干扰等未知因素, 模型预测位置输出值会偏离实际值, 如果不对模型进行纠正, 后面的优化将会偏离真实值。因此在下一个时刻把检测到的实际位置输出值y(k+1)和预测位置值y_PM(k+1|k)相减作为预测误差

$ \mathit{\boldsymbol{e}}\left( {k + 1} \right) = \mathit{\boldsymbol{y}}\left( {k + 1} \right) - {\mathit{\boldsymbol{y}}_{{\rm{PM}}}}\left( {k + 1\left| k \right.} \right)。$

因为未来的预测输出误差和现在的预测误差缺乏因果性描述, 采用时间序列法对未来的预测误差进行校正。

校正后的输出位置预测值

$ {{\mathit{\boldsymbol{\tilde y}}}_{{\rm{cor}}}}\left( {k + 1} \right) = {\mathit{\boldsymbol{y}}_{{\rm{PM}}}}\left( k \right) + \mathit{\boldsymbol{he}}\left( {k + 1} \right), $

式中: $ {\boldsymbol{\tilde y}_{\rm cor}}\left( {k + 1} \right) = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\tilde y}_{\rm cor}}\left( {k + 1\left| {k + 1} \right.} \right)}\\ \vdots \\ {{\boldsymbol{\tilde y}_{\rm cor}}\left( {k + N\left| {k + 1} \right.} \right)} \end{array}} \right]$; h是N维校正向量, h=[h₁, …, h_N]^T, 取h全为1。然后取校正值作为下个时刻的初始值进行计算^[19-20]。

3 仿真试验

以空间三维关节型机械臂作为控制对象进行仿真试验, 控制框图如图 5所示。

机械臂单次运动步长设定为1 mm。作业曲面为体现不同斜率变化率对试验结果的影响, 选取曲面模型为y=2平面上的正弦曲线, z=0.1·sin(0.1·2π/2 000)。作业曲面对系统是完全未知的。机械臂从坐标点(0, 0.2, 0)开始对曲面预测跟踪, 仿真点设置为2 000个。最小二乘法对刚度计算参考文献[15], 其中的参数a₁=0.05, a₂=0.1, a₃=2;模型预测控制时域和预测时域选择50;每个预测时域内的跟踪轨迹设置为根据上个时刻的斜率变化率预测到的后面50个点位的坐标。预测曲面和实际曲面曲线图如图 6a)所示, 预测误差如图 6b)所示。

由图 6可见, 预测轨迹基本完全跟踪实际轨迹。随着实际轨迹斜率的变化, 预测轨迹出现相应的误差, 但误差较小, 对末端轨迹控制结果影响可以忽略不计。机械臂末端轨迹和实际曲面的曲线图如图 7(a)所示, 末端轨迹和实际曲面之间的距离如图 7(b)所示。

如图 7所示, 机械臂末端轨迹很好的跟踪了实际轨迹, 图 7中所示的末端轨迹保持在实际轨迹下方, 代表作业过程中作业曲面所产生的微形变。虽然受控制算法、扰动的影响, 形变出现了一定的变化, 但仍在可接受范围之内, 对末端力的影响较小。机械臂关节各个关节的跟踪曲线和误差如图 8所示。

由图 8可知, 3个关节的关节角实现了很好的跟踪特性。关节2和关节3主要承受了末端力的影响, 因此末端力的振幅和关节1相比较大。关节1不受末端力的影响, 跟踪效果最好。机械臂末端力传感器返回值如图 9所示。

由图 9可见, 控制系统对机械臂的控制效果达到了期望目标, 精度保持在±4 N之内。且受曲面斜率的变化影响较小, 实现了对未知连续曲面的实际跟踪效果。

4 结论

本研究针对磨削机器人对曲面作业时前期示教编程耗时问题和末端力控制问题, 基于曲面预测和模型预测控制方法, 提出了机械臂针对连续曲面作业时的末端力控制算法, 并经过机械臂运动学、动力学分析, 求解出满足期望性能指标的输入序列。仿真结果表明, 基于曲面预测和模型预测控制方法的末端力控制算法, 使得机械臂对未知连续曲面作业时可以一个较高的实时性能和较高的精度控制末端力, 并在保证磨削质量的前提下节省大量的前期示教编程时间。因此, 基于曲面预测和模型预测控制的磨削机械臂对连续曲面作业时的控制算法能够有效完成曲面的预测和末端力的控制, 并具有很好的精度和拓展性。