机器人 2023, Vol. 45 Issue (3): 333-344  
0
引用本文
杜玉红, 刘栋财, 董广宇. 基于广义坐标形式动力学的6-RUS并联机器人零力控制[J]. 机器人, 2023, 45(3): 333-344.  
DU Yuhong, LIU Dongcai, DONG Guangyu. Free-Force Control of 6-RUS Parallel Robot Based on Dynamics of Generalized Coordinate Form[J]. ROBOT, 2023, 45(3): 333-344.  

基于广义坐标形式动力学的6-RUS并联机器人零力控制
杜玉红1,2 , 刘栋财1,2 , 董广宇1,2     
1. 天津工业大学机械工程学院, 天津 300387;
2. 天津市现代机电装备重点实验室, 天津 300387
摘要:针对6-RUS并联机器人的拖动示教零力控制问题, 提出基于广义坐标形式动力学模型的解决方法。首先, 引入动态摩擦模型, 建立广义坐标形式的动力学模型; 利用全局区域法获取电机在工作空间某位置的实际输出转矩, 并采用Savitzky-Golay(SG)算法平滑噪声, 分析机构自身噪声对转矩数据的影响。然后, 结合6-RUS并联机器人动力学方程构建电机的期望输出转矩, 补偿电机的期望转矩与实际转矩之间的误差, 实时跟踪电机转矩实现零力控制。最后, 在6-RUS并联机器人平台上进行实验。结果表明, 在不同负载情况下, 电机的预测与实际输出转矩之间误差均在8.25% 以下; 在完成相同示教过程中, 本文方法所需的末端力矩仅为2.2 N·m, 而基于重力矩与摩擦力矩补偿的零力控制方法为3.4 N·m, 验证了本文零力控制法的有效性。与传统零力控制方法相比, 该方法简化了动力学建模步骤, 提高了构建期望电机转矩的效率, 解决了机器人静止、启动、运动反向阶段期望电机转矩预测不精确的问题, 同时实时补偿了转矩误差, 使得机器人获得更好的拖动效果。
关键词零力控制    动力学模型    动态摩擦模型    并联机器人    拖动示教    
中图分类号:TH112            文献标志码:A            文章编号:1002-0446(2023)-03-0333-12
Free-Force Control of 6-RUS Parallel Robot Based on Dynamics of Generalized Coordinate Form
DU Yuhong1,2 , LIU Dongcai1,2 , DONG Guangyu1,2     
1. School of Mechanical Engineering, Tiangong University, Tianjin 300387, China;
2. Key Laboratory of Advanced Mechatronics Equipment Technology, Tianjin 300387, China
Abstract: A dynamic model based on generalized coordinates is proposed to solve the free-force control problem in drag teaching of 6-RUS parallel robot. Firstly, the dynamic friction model is introduced, and a dynamic model is established in the form of generalized coordinates; the actual output torque of the motor in a workspace position is obtained by the global area method, and the Savitzky-Golay (SG) algorithm is used to smooth the noise and analyze the influence of the mechanism's own noise on the torque data. Then, the expected output torque of the motor is constructed based on the dynamic equation of 6-RUS parallel robot, the error between the expected motor torque and the actual torque is compensated, and the motor torque is tracked in real time to achieve free-force control. Finally, experiments are carried out on the 6-RUS parallel robot platform. The results show that the errors between the predicted and actual output torques of the motor are less than 8.25% in different load conditions. In the same teaching process, the end torque required by the proposed method is only 2.2 N·m, and 3.4 N·m is required by the free-force control method based on the compensation of gravity and friction torques, which verifies the effectiveness of the proposed free-force control method. Comparing with the traditional free-force control method, the proposed method simplifies the dynamic modeling steps, improves the efficiency of constructing the expected motor torque, and solves the problem of inaccurate prediction of the expected motor torque in the stationary, starting, and reverse motion phases of the robot. The torque error is compensated in real time, so that a better drag effect can be obtained.
Keywords: free-force control    dynamic model    dynamic friction model    parallel robot    drag teaching    

1 引言(Introduction)

零力控制的研究是机器人整体拖动示教的基础,在零力控制下,只需较小拖动力就能使机器人产生运动,即人用极小力拖动机器人使其在工作空间中任意运动。并联机器人相对于串联机器人具有刚度大、结构紧凑、位置精度高、承载能力强等优点[1-2],适于在狭小工作空间内生产高精度、高刚度、大载荷的工件,近些年在机器人领域的占比不断增加。因并联机构自身存在闭环结构和运动约束的特点,且动力学建模复杂,所以基于动力学的并联机器人零力控制研究相对较少。目前可以通过改进串联机器人的零力控制方法,来解决并联机器人在进行复杂轨迹规划时存在的定位难、效率低问题,最终实现直接拖动示教,提高人机交互性能[3-4]

目前,零力控制研究对象大多集中在串联、协作机器人上,游有鹏等[5]提出自测量的重力矩及摩擦力矩计算方案,通过补偿重力矩以及摩擦力矩实现串联机器人无传感器零力控制;侯澈等[6]针对串联机器人更换末端执行器以及负载变化的情况,提出了一种变负载自适应零力控制方法;张铁等[7]采用弹性摩擦模型估计关节摩擦力,在关节起动阶段短暂增加摩擦估计值以增加关节驱动力矩,实现免力矩传感器的串联机器人零力控制。李智清等[8]针对零力控制相关的未知碰撞安全问题,提出了一种新型机器人碰撞检测算法,通过实时观测关节输出力矩与动力学估计力矩的偏差实现机器人零力控制与碰撞检测;宋吉来等[9]针对机器人系统对柔顺性、安全性的要求,提出了一种基于电机电流的机器人主动柔顺控制方法,并采用硅片传输机器人进行拖曳式示教和碰撞保护实验。基于动力学模型和运动状态来估计机器人驱动关节所受的外力是实现零力控制的一种方法,吴汶祥等[10]提出一种基于改进傅里叶级数的辨识方法,推导动力学模型的线性形式并准确地重构机器人的关节力矩值,进而进行零力控制研究;张春涛等[11]利用拉格朗日法整理关节重力矩的三角函数回归矩阵,建立机器人快速动力学模型补偿关节力矩,席万强等[12]提出了一种基于人工蚁群算法的动力学参数辨识方法,推导动力学模型的非线性形式,使关节预测力矩与实际力矩有较高的匹配度,为零力控制研究提供新的思路。

本文根据中科极控(天津)科技有限公司对复杂曲面小型工件喷涂的要求,选用6-RUS并联机器人作为喷涂机器人的主体结构,建立基于广义坐标形式牛顿—欧拉法的动力学模型,引入动态关节摩擦模型,实时构建机器人运动过程中电机的预测转矩,通过导入电机预测转矩实现并联机器人末端零力平衡,并在自研的实验平台上验证方法的有效性。

2 喷涂机器人构形(Configuration of the spraying robot)

中科极控(天津)科技有限公司对不同复杂曲面形状的汽车LED日行灯进行喷涂,日行灯外观如图 1所示,尺寸为400 mm$ \times $300 mm,上表面近似为椭圆形,下表面为凹凸面,内表面有复杂的曲面纹理。为满足喷涂需求,采用6-RUS并联机器人,结构示意图如图 2所示。

图 1 LED日行灯 Fig.1 LED of daytime running light
图 2 并联喷涂机器人结构图 Fig.2 Structure diagram of the parallel spraying robot

相对6-PUS、6-PSS并联机器人,6-RUS并联机器人工作空间较大、结构对称,具有稳定性强、驱动副固连基座能减轻结构载荷等优点。下平台为动平台,动平台实现3个方向的平动和3个方向的转动;上平台为静平台,由6个结构相同的驱动分支RUS连接2个平台,机构简图如图 3所示,具体参数如表 1所示。

图 3 6-RUS并联机构简图 Fig.3 Sketch of the 6-RUS parallel mechanism 注:$A_1$~ $A_6$分别是6个转动副位置,$B_ 1$~ $B_6$为驱动杆与从动杆之间的万向节,$P_1$~ $P_6$为动平台上6个球副。
表 1 机构基本参数表 Tab. 1 Basic parameter table of the mechanism
3 基于广义坐标形式牛顿—欧拉法的6-RUS并联机器人动力学建模(Dynamic modeling of the 6-RUS parallel robot based on Newton-Euler method of generalized coordinates)

目前,构建6-RUS并联机器人驱动铰链约束力和末端约束力之间的关系时,采用牛顿—欧拉公式[13]的动力学模型容易求解运动副约束反力,但需推导186个微分方程、180个约束方程,代数微分方程矩阵维数达到366,计算过于复杂;采用拉格朗日方程[14]的动力学模型可得到较为简洁的闭环动力学方程,但涉及矢量之间的偏微分求导,难于计算。为简化6-RUS并联机器人动力学建模,对牛顿—欧拉方法[15]进行改进,采用广义坐标形式进行分析,简化动力学建模步骤,提高计算效率。

3.1 动态摩擦力模型

摩擦力是并联机器人动态特性关键指标,考虑关节电机、机构截面材质等,取图 3$ A_{1} $点为机器人关节1进行摩擦力分析,建立精确的6-RUS并联机器人动态摩擦力模型。

(1) 预滑动阶段摩擦力模型

6-RUS并联机器人静止启动时力矩的突然变化会产生振荡,为避免振荡带来的影响,在预滑动阶段[16]用S型函数作为摩擦模型,使其与动态连续阶段模型保持连续,如图 4所示。当关节电机与连杆速度为0时,对应的时刻为$ t_{0} $,令$ t_{\rm a} $为宏观上机构进入动态连续阶段的时间间隔,那么$ t_{0} $~$ t_{\rm a} $为静止启动阶段,此阶段宏观上可认为电机速度为0,因此摩擦力矩仅为时间函数,如图 5所示。

图 4 预滑动摩擦模型 Fig.4 Pre-sliding friction model
图 5 正向预滑动摩擦模型 Fig.5 Forward pre-sliding friction model

根据文[16] 提出的构造方程,当$ t_{0} = $ 0时,根据S型函数表达式,正向预滑动摩擦模型表达式为

$ \begin{align} T_{\rm a} =\frac{k_{1} T_{+}} {1+{\rm e}^{-\alpha t}}+k_{2} \end{align} $ (1)

其中:$ T_{\rm a} $为预滑动阶段总摩擦力矩,$ T_{+} $是正向最大静摩擦力矩,$ \alpha $为曲线斜率,$ k_{1} $$ k_{2} $为调整模型的比例因子与偏差。

对于给定一个常数$ f_{\rm s} $$ 0.9<f_{\rm s} <1 $),以及函数$ f_{t} =1/(1+{\rm e}^{-\alpha t}) $,可找到一个$ \alpha $使得$ f(t_{\rm b})=f_{\rm s} $,那么S型函数的曲线斜率为

$ \begin{align} \alpha =\frac{2\ln (f_{\rm s} /(1-f_{\rm s}))}{\hat{t}_{\rm b}} \end{align} $ (2)

其中:$ \hat{t}_{\rm b} $为过渡时间估计值。

为保证预滑动阶段的连续性,选用参数$ k_{1} $$ k_{2} $使得$ T_{\rm a} (0)=0 $$ T_{\rm a} (t_{\rm b})=0 $,则

$ \begin{equation} \begin{split} k_{1} & =\frac{2}{2f_{\rm s} -1} \\ k_{2} & =-T_{+} \frac{1}{2f_{\rm s} -1} \end{split} \end{equation} $ (3)

$ f_{\rm s} $为控制器设置的常量,$ k_{1} $$ k_{2} $取决于$ T_{+} $$ f_{\rm s} $

(2) 动态阶段的摩擦力矩模型

当机器人运动时,由于并联机器人存在单向折返情况,摩擦函数为不连续函数,所以采用连续可微的摩擦力模型[17]建立并联机器人关于速度函数的摩擦力矩模型。摩擦力矩$ T(\dot{q}) $

$ \begin{align} T(\dot{q})=\, & \delta_{1} (\tanh (\alpha_{1} \dot{q})-\tanh (\alpha {}_{2}\dot{q}))+ \\ & \delta_{2} \tanh (\alpha_{3} \dot{q})+\delta_{3} \dot{q} \end{align} $ (4)

其中:$ \delta_{1}, \delta_{2}, \delta_{3}, \alpha_{1}, \alpha_{2}, \alpha_{3} $为未知参数,$ \tanh (\alpha_{1} \dot{q})-\tanh (\alpha _{2}\dot{q}) $是Stribeck现象[18]$ \delta_{2} \tanh (\alpha_{3} \dot{q}) $代表库仑摩擦,$ \delta_{3} \dot{q} $为黏性项。

3.2 动力学建模

选择6-RUS动平台相对静平台的位置与姿态坐标建立广义坐标:

$ \begin{align} \mathit{\boldsymbol{q}} =(x, y, z, \alpha, \beta, \gamma)^{\rm T} \end{align} $ (5)

各构件的平动与转动加速度方程为

$ \begin{align} \begin{pmatrix} \ddot{\mathit{\boldsymbol{u}}}_{\rm p} \\ \dot{\mathit{\boldsymbol{w}}}_{\rm p} \end{pmatrix}= \begin{pmatrix} {\mathit{\boldsymbol{J}}}_{\rm Tp} \\ {\mathit{\boldsymbol{J}}}_{\rm Rp} \end{pmatrix} \ddot{\mathit{\boldsymbol{q}}}+ \begin{pmatrix} {\mathit{\boldsymbol{\delta}}}_{\rm p} \\ {\mathit{\boldsymbol{\gamma}}}_{\rm p} \end{pmatrix} \end{align} $ (6)

其中:$ {\mathit{\boldsymbol{\delta}}}_{\rm p} =(\dot{\mathit{\boldsymbol{J}}}_{\rm Tp1}, \cdots, \dot{\mathit{\boldsymbol{J}}}_{\rm Tp6})^{\rm T} $, $ {\mathit{\boldsymbol{\gamma}}}_{\rm p} =(\dot{\mathit{\boldsymbol{J}}}_{\rm Rp1}, \cdots, \dot{\mathit{\boldsymbol{J}}}_{\rm Rp6})^{\rm T} $, $ {\mathit{\boldsymbol{J}}}_{\rm Tp} $$ {\mathit{\boldsymbol{J}}}_{\rm Rp} $分别表示构件相对于广义坐标的平动和转动雅可比矩阵。

建立6-RUS并联机器人动平台与驱动杆的牛顿—欧拉方程:

$ \begin{align} & \begin{pmatrix} {m_{\rm p}} & 0 & 0 & 0 \\ 0 & {{\mathit{\boldsymbol{I}}}_{\rm p}} & 0 & 0 \\ 0 & 0 & {m_{{\rm c}i}} & 0 \\ 0 & 0 & 0 & {{\mathit{\boldsymbol{I}}}_{{\rm c}i}} \end{pmatrix} \begin{pmatrix} {\ddot{\mathit{\boldsymbol{u}}}_{\rm p}} \\ {\dot{w}_{i}} \\ {\ddot{\mathit{\boldsymbol{u}}}_{{\rm c}i}} \\ {\dot{w}_{i}} \end{pmatrix}+\begin{pmatrix} 0 \\ {{\mathit{\boldsymbol{I}}}_{\rm p}} \\ 0 \\ {{\mathit{\boldsymbol{I}}}_{{\rm c}i}} \end{pmatrix}w_{i}^{2} \\ = & \begin{pmatrix} {{\mathit{\boldsymbol{F}}}_{{\rm p}i}^{\rm n}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm p}i}^{\rm n}} \\ {{\mathit{\boldsymbol{F}}}_{{\rm c}i}^{\rm n}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm c}i}^{\rm n}} \end{pmatrix}+\begin{pmatrix} {{\mathit{\boldsymbol{F}}}_{{\rm p}i}^{\rm a}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm p}i}^{\rm a}} \\ {{\mathit{\boldsymbol{F}}}_{{\rm c}i}^{\rm a}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm c}i}^{\rm a}} \end{pmatrix} \end{align} $ (7)

其中:$ {\mathit{\boldsymbol{F}}}_{{\rm c}i}^{\rm n} $$ {\mathit{\boldsymbol{M}}}_{{\rm c}i}^{\rm n} $为驱动铰链处约束力以及电机转矩,$ {\mathit{\boldsymbol{F}}}_{{\rm c}i}^{\rm a} $$ {\mathit{\boldsymbol{M}}}_{{\rm c}i}^{\rm a} $为中间铰链约束力以及力矩,$ {\mathit{\boldsymbol{M}}}_{{\rm c}i}^{\rm n} = $ $ {\mathit{\boldsymbol{M}}}_{i}^{\rm n} -{\mathit{\boldsymbol{T}}}_{i}^{\rm n} $$ {\mathit{\boldsymbol{T}}}_{i}^{\rm n} $为电机摩擦力矩,$ {\mathit{\boldsymbol{F}}}_{i}^{\rm n} $$ {\mathit{\boldsymbol{M}}}_{{\rm p}i}^{\rm n} $为动平台受到的铰链约束力及质心的转矩,$ {\mathit{\boldsymbol{F}}}_{i}^{\rm a} $$ {\mathit{\boldsymbol{M}}}_{{\rm p}i}^{\rm a} $为外部负载与力矩,$ m $$ \mathit{\boldsymbol{I}} $分别为构件的质量及转动惯量。

合并6-RUS并联机器人各构件的牛顿—欧拉方程:

$ \begin{align} \mathit{\boldsymbol{m\ddot{u}}}& ={\mathit{\boldsymbol{F}}}^{\rm n}+{\mathit{\boldsymbol{F}}}^{\rm a} \end{align} $ (8)
$ \begin{align} \mathit{\boldsymbol{I\dot{w}}}+{\mathit{\boldsymbol{w}}}\times \mathit{\boldsymbol{Iw}} & = \mathit{\boldsymbol{M}}^{{{\rm n}}}+{\mathit{\boldsymbol{M}}}^{{{\rm a}}} \end{align} $ (9)

其中:$ i =1, 2, 3, 4, 5, 6 $

$ \begin{align*} {\mathit{\boldsymbol{m}}} & = \begin{pmatrix} {{{m}}_{{\rm c}i} {\mathit{\boldsymbol{E}}}_{3}} & \cdots & 0 \\ \vdots & {{{m}}_{{\rm d}i} {\mathit{\boldsymbol{E}}}_{3}} & \vdots \\ 0 & \cdots & {{{m}}_{\rm p} {\mathit{\boldsymbol{E}}}_{3}} \end{pmatrix} \\ {\mathit{\boldsymbol{I}}} & = \begin{pmatrix} {{\mathit{\boldsymbol{I}}}_{{\rm c}i}} & \cdots & 0 \\ \vdots & {{\mathit{\boldsymbol{I}}}_{{\rm d}i}} & \vdots \\ 0 & \cdots & {{\mathit{\boldsymbol{I}}}_{\rm p}} \end{pmatrix}, \; {\mathit{\boldsymbol{F}}}^{\rm n}= \begin{pmatrix} {{\mathit{\boldsymbol{F}}}_{{\rm c}i}^{\rm n}} \\ {{\mathit{\boldsymbol{F}}}_{{\rm d}i}^{\rm n}} \\ {{\mathit{\boldsymbol{F}}}_{{\rm p}i}^{\rm n}} \end{pmatrix}, \; {\mathit{\boldsymbol{F}}}^{\rm a}= \begin{pmatrix} {{\mathit{\boldsymbol{F}}}_{{\rm c}i}^{\rm a}} \\ {{\mathit{\boldsymbol{F}}}_{{\rm d}i}^{\rm a}} \\ {{\mathit{\boldsymbol{F}}}_{{\rm p}i}^{\rm a}} \end{pmatrix}\\ {\mathit{\boldsymbol{M}}}^{\rm a}& = \begin{pmatrix} {{\mathit{\boldsymbol{M}}}_{{\rm c}i}^{\rm a}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm d}i}^{\rm a}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm p}i}^{\rm a}} \end{pmatrix}, {\mathit{\boldsymbol{M}}}^{\rm n}= \begin{pmatrix} {{\mathit{\boldsymbol{M}}}_{{\rm c}i}^{\rm n}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm d}i}^{\rm n}} \\ {{\mathit{\boldsymbol{M}}}_{{\rm p}i}^{\rm n}} \end{pmatrix} \end{align*} $

$ {\mathit{\boldsymbol{\delta}}} ={\mathit{\boldsymbol{w}}}+{\mathit{\boldsymbol{Iw}}} $$ \dot{\mathit{\boldsymbol{w}}}=( {\mathit{\boldsymbol{J}}}_{\rm T}, {{\mathit{\boldsymbol{J}}}_{\rm R}})^{\rm T} $,代入式(8)(9) 得出:

$ \begin{align} \mathit{\boldsymbol{\varPhi H\ddot{q}}}+{\mathit{\boldsymbol{K}}}({\mathit{\boldsymbol{q}}}, \dot{\mathit{\boldsymbol{q}}})=\overline{\mathit{\boldsymbol{F}}}^{\rm a}+\mathit{\boldsymbol{W\eta}} \end{align} $ (10)

其中:$ \bm\eta $为活动铰链约束力$ {\mathit{\boldsymbol{\eta}}} =(\eta_{\rm c}, \eta_{\rm d}, \eta_{\rm p})^{\rm T} $$ \overline{\mathit{\boldsymbol{F}}}^{\rm a} $为构件受到的力与力矩组成的矩阵,$ \mathit{\boldsymbol{W}} $为铰链处的传递矩阵,$ {\mathit{\boldsymbol{\varPhi}}} = \begin{pmatrix} {{\mathit{\boldsymbol{m}}}} & \\ & {{\mathit{\boldsymbol{I}}}} \end{pmatrix} $$ {\mathit{\boldsymbol{H}}}= \begin{pmatrix} {\mathit{\boldsymbol{J}}}_{\rm T} \\ {\mathit{\boldsymbol{J}}}_{\rm R} \end{pmatrix} $$ {\mathit{\boldsymbol{K}}}(\mathit{\boldsymbol{q, \dot{q}}})= \begin{pmatrix} {\mathit{\boldsymbol{mJ}}}_{\rm T} \dot{\mathit{\boldsymbol{q}}} \\ {\mathit{\boldsymbol{IJ}}}_{\rm R} \dot{\mathit{\boldsymbol{q}}}+{\mathit{\boldsymbol{w}}}({\mathit{\boldsymbol{Iw}}}) \end{pmatrix} $

消去铰链约束反力项,基于广义坐标的6-RUS并联机器人运动微分方程为

$ \begin{align} {\mathit{\boldsymbol{H}}}^{\rm T}\mathit{\boldsymbol{\varPhi H\ddot{q}}}+{\mathit{\boldsymbol{H}}}^{\rm T}{\mathit{\boldsymbol{K}}}({\mathit{\boldsymbol{q}}}, \dot{\mathit{\boldsymbol{q}}})={\mathit{\boldsymbol{H}}}^{\rm T}\overline{\mathit{\boldsymbol{F}}}^{\rm a} \end{align} $ (11)

求解6个运动微分方程,再求解180个约束代数方程即可得到动力学响应以及运动副约束力,通过运动学反解,计算获得各构件的笛卡儿坐标。

3.3 零力控制方案

目前,基于力矩模式的零力控制策略大多采用电机提供驱动力补偿连杆的重力矩与摩擦力矩,但未考虑机器人启动阶段以及运动方向改变时的关节摩擦变化和基于动力学模型的期望力矩与实际转矩之间误差的影响。

采用基于广义坐标形式的6-RUS并联机器人动力学模型计算预测转矩,估计连杆重力矩以及摩擦力矩,进一步估计预测转矩与实际转矩的误差,从而构建输出转矩$ T(\dot{q}) $,零力控制策略如图 6所示。

$ \begin{align} T(\dot{q})=T_{\rm f} (\dot{q})+G(\dot{q})+\hat{M} \end{align} $ (12)
图 6 零力控制策略框图 Fig.6 Diagram of the free-force control strategy

其中:$ T_{\rm f} (\dot{q}) $为动态摩擦模型对关节摩擦力矩的估计值,$ G(\dot{q}) $为重力矩估计值,$ \hat{M} $为实际转矩$ M_{\rm d} $与预测转矩$ M_{{\rm c}i}^{\rm n} $之间的误差值。

将机构同一位置与姿态情况下的电机电流信号转化为实际电机转矩,通过对比$ M_{{\rm c}i}^{\rm n} $$ M_{\rm d} $得到误差值$ \hat{M} $。在示教模式下将$ \hat{M} $实时补偿给模型预测转矩,使机构自身达到平衡,即自然人用极小外力拖动机器人到达任意位置,达到拖动示教的目的。

4 基于动力学的并联机器人零力控制实验(Dynamics-based free-force control experiment of parallel robots) 4.1 实验平台

搭建的实验平台如图 7所示,在工作空间内设定6-RUS并联机器人的运行轨迹,对比分析机构末端负载、电机转矩模型预测值与实际值的关系以及误差,根据实际喷涂工况的需求,实验选用分别悬挂30 N、45 N、60 N三种负载的喷枪进行实验,在同一负载条件下,机器人末端沿上述6种空间轨迹曲线运动,伺服驱动器内置传感器每隔1 s采集机器人电机的转矩,经正向动力学映射为末端的力与力矩进行实验验证。

图 7 6-RUS并联机器人实验平台 Fig.7 The experimental platform of 6-RUS parallel robot
4.2 摩擦模型参数辨识

(1) 预滑动阶段

机器人预滑动阶段摩擦模型中参数$ T_{+} $$ \hat{t}_{\rm b} $$ f_{\rm s} $未知,其中$ f_{\rm s} $为控制器设置的常量,根据文[19] 指令加速度与过渡时间的试验数据关系,推导出静、动摩擦力矩过渡时间的估计值$ \hat{t}_{\rm b} $$ T_{+} $的辨识方法参照文[20],各参数的取值如表 2所示。

表 2 预滑动阶段摩擦参数辨识值 Tab. 2 Identification values of friction parameters in pre-sliding stage

(2) 动态阶段

实验中不安装驱动臂,可近似认为机器人关节系统处于无负载状态,且6-RUS并联机器人各关节的结构与安装相同,选用关节1作为对象进行实验,关节以匀速进行往复运动,匀速运动不受科氏力与离心力的影响,可认为编码器测量的关节转矩$ T $为摩擦力矩$ T(\dot{q}) $。因此,采用图 8的速度曲线,匀速运动时间为2 s、速度范围为1~100 rad/s,采样时间间隔为2 ms,让关节进行正转与反转运动,并提取匀速段力矩数据,可得到不同速度下的摩擦力矩。

图 8 关节电机输入速度 Fig.8 Input speed of the joint motor

对各匀速段采集的1000个关节转矩值求取均方根(RMS)可以得到对应速度下的摩擦力矩值:

$ \begin{align} \gamma_{\:\rm RMS}=\sqrt{\frac{T_{1}^{2} +T_{2}^{2} +\cdots +T_{n}^{2}} {n}} \end{align} $ (13)

机器停止运动后,每个关节均得到50对正向与反向转矩,进一步对动态阶段摩擦模型式(4) 进行辨识,摩擦力矩与速度曲线如图 9所示,参数辨识结果如表 3所示。

图 9 摩擦力矩与速度曲线 Fig.9 Friction torque and speed
表 3 动态阶段摩擦模型参数辨识值 Tab. 3 Parameter identification values of the friction model in dynamic stage

引入拟合优度$ R^{2} $与均方根误差(RMSE)评价模型的拟合程度,计算公式为

$ \begin{align} R^{2} & =1-\frac{ \sum\limits_{i=1}^{n} (y-\hat y) }{ \sum\limits_{i=1}^{n}(y-\bar y)^2 } \end{align} $ (14)
$ \begin{align} \gamma_{\:\rm RMSE} & =\sqrt{\frac{ \sum\limits_{i=1}^{n} (y-\hat y) }{n}} \end{align} $ (15)

其中:$ y $为实际采样数据,$ \bar{y} $为实际采样数据的均值,$ \hat{y} $为模型预测值。

根据计算得到$ R^{2} $为0.992、RMSE为0.452,由图 9可知动态阶段摩擦模型曲线对采样数据有很好的拟合效果,并且拟合优度接近1,由此可知摩擦模型精度高。

4.3 零力控制实验

针对汽车LED日行灯的外观尺寸及其曲面复杂特性,采用全局区域法模拟喷涂机器人工作空间进行零力控制实验。在机器人工作空间中定义一个原点以及6条空间射线$ L_{i}: A_{i}X+ B_{i}Y + C_{i}Z + D= 0 $,并在射线上均匀取480个数据采集点,其中6条空间射线$ L_{i} $O-XY坐标系进行投影获得的6条线段之间的夹角均为60$ ^{\circ} $$ L_{i} $由原点向工作空间边缘延展,根据卡瓦列利原理[21],6条空间射线所覆盖的范围为机器人末端工作空间,如图 10所示。

图 10 机器人工作空间 Fig.10 Robot workspace
4.4 电机转矩数据预处理

在机器人每条运行曲线上采集480组电机电流信号与关节位置信号,由于信号周期性以及机构在运动过程受到自身抖动、刚度变化等噪声污染,需对信号进行预处理。采用多项式进行拟合:

$ \begin{align} T_{(2n+1)\times 1} =t_{(2n+1)\times k} +a_{k\times 1} +E_{(2n+1)\times 1} \end{align} $ (16)

进而采用最小二乘法得到SG模型的滤波值:

$ \begin{align} {\mathit{\boldsymbol{P}}}={\mathit{\boldsymbol{t}}}\times({\mathit{\boldsymbol{t}}}^{\text{trans}}\times {\mathit{\boldsymbol{t}}})^{-1}\times {\mathit{\boldsymbol{t}}}^{\text{trans}}\times {\mathit{\boldsymbol{X}}} \end{align} $ (17)

其中:$ \mathit{\boldsymbol{X}} $$ 2 n +1 $$ t $时刻的矩阵,$ \mathit{\boldsymbol{t}}^\text{trans} $$ t $时刻对应矩阵转置。

将采集的信号数据分别以3、5、8点为一个窗口的形式分解为$ n $个窗口进行多项式拟合,对采集的电机输出转矩进行平滑处理,效果图如图 11所示。

图 11 不同窗口SG滤波对轴1电机测量值的平滑效果 Fig.11 The smoothing effect of SG filtering in different windows on the measured motor value of axis 1

图 11所示,截取1轴数据进行分析,可以得到5点窗口相比于3点窗口对测量值的平滑效果相对更好,同时8点窗口相对于5点窗口产生了过渡的平滑效果,导致数据失真,因此采用5点窗口对电机输出转矩值进行平滑处理。在对6个轴的转矩数据进行预处理的过程中发现,5、6轴的电机输出转矩相对预测值抖动较大,电机测量值比对结果如图 12所示。由图 12可知,机构的电机在0 ~8 s时实际转矩值与预测转矩值误差相对较小,8~15 s时实际转矩值相对预测转矩值小,15 s后实际转矩值相对预测转矩值大,同时转矩值波动较大,所以5、6轴的从动杆将产生弯曲,同时铰链处将产生一定程度的松动。机构在200 mm$ \times $100 mm工作范围内转矩误差无影响,当超出范围时,铰链松动和杆弯曲对机构的相互应力影响较大,使得电机预测转矩值与实际转矩值不符,实际转矩值随应力变化且数据波动大。

图 12 轴5、6电机测量值比对结果 Fig.12 Comparison results of measured motor values of axes 5 and 6
5 实验结果分析(Analysis on experimental results) 5.1 6-RUS并联机器人动力学模型验证

在6-RUS并联机器人设备上,机构末端预测负载$ F_{\text{w}x}=F_{\text{w}y}= $ 0,$ F_{\text{w}z}= $ 30 N$ , $ 45 N$ , $ 60 N,实际转矩$ M_{\text{w}x}= $ $ M_{\text{w}y}=M_{\text{w}z}= $ 0,机构姿态角$ \alpha = $ 10$ ^{\circ} $$ \beta = $ 15$ ^{\circ} $$ \gamma = $ 0$ ^{\circ} $,对比模型预测负载值与实际值之间的误差,机构末端实际负载如图 13所示。

图 13 机构末端实际负载 Fig.13 The actual load at the end of the mechanism

根据图 13可知,机构末端$ z $方向数值计算的预测负载$ F_{ z} $与实际负载$ F_{\text{w}z} $之间最大误差为8.12%,$ x $$ y $方向数值计算的预测负载$ F_{ x} $$ F_{ y} $与实际负载$ F_{\text{w}x} $$ F_{\text{w}y} $之间误差都在8.5% 以内,预测转矩$ M_{ x} $$ M_{ y} $$ M_{ z} $与实际转矩$ M_{\text{w}x} $$ M_{\text{w}y} $$ M_{\text{w}z} $之间误差都为5%。不考虑机构自身加工精度误差时,结果可验证基于广义坐标形式下6-RUS并联机器人的零力控制动力学模型的有效性。

5.2 不同负载情况下电机转矩分析

根据实际喷涂工况的需求,机器人末端分别悬挂30 N、45 N、60 N三种预测负载的喷枪,本文选用$ F_{\text{w}x}=F_{\text{w}y}= $ 0,$ F_{\text{w}z}= $ 30 N$ , $ 45 N$ , $ 60 N,$ M_{\text{w}x}=M_{\text{w}y}=M_{\text{w}z}= $ 0,姿态角$ \alpha = $ 10$ ^{\circ} $$ \beta = $ 15$ ^{\circ} $$ \gamma = $ 0$ ^{\circ} $进行分析,机构末端运行轨迹曲线为$ L_{1} $~$ L_{6} $,目前串联机器人常用的零力控制算法有基于重力矩与摩擦力矩补偿的零力控制方法[5],该方法采用基于自测量的重力矩与摩擦力矩计算方案进行力矩补偿且效果相对较好,将本文模型的零力控制算法与文[5] 算法进行对比,分析机器人6个电机实际转矩与本文模型预测转矩、文[5] 预测转矩之间的误差,将6个电机转矩数据分为3组,实验结果如图 14~图 16所示。

图 14 第1组(轴1、6)电机转矩 Fig.14 The first group motor torque (axes 1 and 6)
图 15 第2组(轴2、3)电机转矩 Fig.15 The second group motor torque (axes 2 and 3)
图 16 第3组(轴4、5)电机转矩 Fig.16 The third group motor torque (axes 4 and 5)

根据图 14~图 16可知,末端预测负载为30 N时电机转矩的范围为3.5~15.8 N$ \cdot $m,机构末端预测负载为45 N时电机转矩范围为4~18 N$ \cdot $m,机构末端预测负载为60 N时电机转矩范围为4.5~21 N$ \cdot $m,同时本文基于广义坐标形式的动力学模型的预测转矩与基于自测量重力矩与摩擦力矩补偿的零力控制算法的预测转矩均对实际力矩有较好的跟随性;当喷涂路线从工件侧曲面转向内曲面(图 14~图 16中曲线$ L_{1} $~$ L_{6} $之间过渡区域)时,基于广义坐标形式的动力学模型得到的预测转矩相对平滑且与实际转矩更接近,因此采用精确的动态摩擦模型可有效

提高动力学模型的准确性。末端负载增加时电机转矩随之增加,同时本文模型的零力控制算法相比文[5] 算法,预测转矩更接近实际转矩值。在曲线$ L_{5} $$ L_{6} $的过渡阶段,实际转矩值与模型预测转矩值相差比较大,分析可知该位置是运行曲线极小夹角往返点,杆件弯曲导致转矩误差大。

图 14~16中,曲线$ L_{1} $$ L_{2} $处于$ (0, 186, 180) $$ (0, 0, 0) $$ (0, -186, 180) $三点围成的区域内,曲线$ L_{3} $$ L_{4} $处于$ (-130, -80, 70) $$ (0, 0, 0) $$ (130, -80, 70) $三点围成的区域内,曲线$ L_{5} $$ L_{6} $处于$ (150, 80, $ $ 150) $$ (0, 0, 0) $$ (-150, 80, 150) $三点围成的区域内。

为了进一步阐述基于广义坐标形式的动力学模型预测转矩的精确性,图 17给出末端不同负载下,本文模型、文[5] 方法的预测转矩与实际电机转矩误差曲线。

图 17 实际与预测电机转矩误差 Fig.17 Error between actual and theoretical motor torques

在喷涂机器人启动、终止阶段,由图 14(c)(轴6)、图 15(b)(轴2、3)、图 16(c)(轴5)可知,基于自测量重力矩与摩擦力矩补偿算法的转矩出现突变且与实际转矩误差大,而本文引入预滑动摩擦模型的预测转矩均更接近实际转矩;在喷涂机器人变换运动方向时(即图 14~16$ L_{1} $~$ L_{6} $衔接处),由图 15(b) 轴2的$ L_{4} $$ L_{5} $$ L_{6} $处、图 15(c) 轴2的$ L_{4} $$ L_{5} $$ L_{2} $$ L_{3} $处可知,基于自测量重力矩与摩擦力矩补偿算法的转矩与实际转矩偏离大,而本文引入动态连续摩擦模型的预测转矩对实际转矩有很好的跟随性且运动曲线衔接处也更接近实际转矩。同时由图 17可知,在不同负载的喷涂机器人启动、终止阶段文[5] 方法的误差大于本文算法,同时在第20 s、60 s运动方向改变处误差也相对大,由此说明本文的动态摩擦模型相对于文[5] 的摩擦模型在启动、终止、变向阶段更精确、补偿更好,有利于提高算法的精确性。

图 17可以看出,在机器人喷涂工件的工作空间内,基于广义坐标形式动力学模型的预测转矩误差明显小于文[5] 零力控制方法,本文模型转矩误差在5% 左右,文[5] 转矩为7.5%~10%。在运动曲线过渡区域(误差峰值处),本文模型误差相对较小,最大误差为8.25%。以上分析进一步验证了基于广义坐标形式动力学模型的零力控制方法的精确性。

5.3 末端无负载的拖动示教实验

根据喷涂工件实际工况的需求,6-RUS并联机器人的末端需在工作空间$ X $$ Y $$ Z $方向上运动,同时末端姿态可旋转0~35$ ^{\circ} $,以适应复杂曲面的喷涂工况。因此在零力控制模式下,令自然人拖动并联机器人进行示教,模拟喷涂特殊工况。实验过程中的特殊位置如图 18所示,采用本文零力控制算法与文[5] 方法进行示教,实时跟踪拖动示教中末端受到的手部力,如图 1920所示,实验过程中拖动示教与再现轨迹如图 21所示。

图 18 机器人拖动示教实验截图 Fig.18 The screenshots of robot drag teaching experiment
图 19 示教的末端受力(模型) Fig.19 Taught end force (model)
图 20 示教的末端受力(文[5]) Fig.20 Taught end force (literature [5])
图 21 拖动示教与再现轨迹 Fig.21 Drag teaching and teaching-playback trajectories

图 1920可知,操作者拖动机器人进行示教,对比在本文算法与基于自测量重力矩与摩擦力矩的补偿算法下末端的受力情况,在机器人启动阶段,预滑动摩擦模型补偿下的末端力与力矩相对较小,从而能轻松启动机器人;示教前50 s机器人频繁改变运动方向,在动态摩擦模型补偿下末端力矩数值变化均呈现稳定,末端力相对文[5] 方法均抖动小,由此说明动态摩擦模型提高了算法的精确性。同时本文模型下机器人拖动示教过程最大末端力为2 N、最大末端力矩为2.2 N$ \cdot $m,相比文[5] 方法更小,基本符合人手部的感知力,同时拖动示教的轨迹与再现轨迹基本相近,从而进一步验证了基于广义坐标形式动力学模型的零力控制方法对6-RUS并联机器人拖动示教的有效性。

综上,针对工件复杂曲面喷涂工况,在机器人工作空间内、不同负载情况下,基于广义坐标形式动力学模型的6个电机预测转矩值与实际转矩值始终相近且误差约5%,验证了基于广义坐标形式动力学模型的6-RUS并联机器人零力控制方法的正确性,同时机器人启动与运动转向阶段误差小,表明动态摩擦模型有效提高了动力学模型的精确性。在零力控制模式下,自然人拖动6-RUS并联机器人进行示教,本文模型的零力控制算法末端力最大值为2 N、末端力矩最大值为2.2 N$ \cdot $m,基于自测量重力矩与摩擦力矩补偿的算法需末端力矩为3.4 N$ \cdot $m,可以认为本文模型的零力控制方法相比基于自测量重力矩与摩擦力矩的补偿算法具有更好的拖动喷涂效果,即自然人只需很小的力就能拖动机器人到达工作空间内任意位置,解决了并联机器人在位置模式下拖动困难、无法获取复杂轨迹信息的问题,为多自由度并联机器人拖动示教提供了基础。

6 结论(Conclusion)

针对并联机器人拖动示教应用场景,提出了一种基于广义坐标形式动力学模型的6-RUS并联机器人零力控制方法。

(1) 利用广义坐标形式建立6-RUS并联机器人动力学模型,引入动态摩擦模型,通过在动态摩擦模型中加入预滑动摩擦项,提高关节摩擦精度,解决了机器人静止启动、运动反向阶段构建期望电机转矩误差偏大的问题,从而建立精确的6-RUS并联机器人动力学模型。

(2) 采用全局区域法获取各位置电机实际输出转矩,进一步用SG算法平滑输出转矩数据,分析机构自身噪声对输出转矩数据产生剧烈抖动的影响。

(3) 运用6-RUS逆动力学模型获取电机预测输出转矩,在喷涂工件工作空间内,分析电机实际转矩与预测转矩关系;结果表明电机实际转矩与预测转矩的误差为5% 左右。在拖动示教过程中,实时补偿转矩误差能提高示教模式下电机转矩补偿的精确度,使得机器人末端所需接触力矩仅为2.2 N$ \cdot $m。

实验表明,本文提出的基于广义形式动力学模型的6-RUS并联机器人零力控制方法是可行的,与基于自测量重力矩与摩擦力矩的补偿算法相比,采用提出的零力控制方法能实现机器人任意阶段轻松拖动,减小末端接触力矩,达到更好的拖动喷涂效果。该方法为多自由度并联机器人的拖动示教研究提供了新思路,同时为多自由度并联机器人应用于喷涂行业打下基础。

参考文献(References)
[1]
Mikula A M. Tilting pad thrust bearing. New design cuts power loss[J]. Machine Design, 1987, 59(26): 117-120.
[2]
苑飞虎, 赵铁石, 赵延治, 等. 并联机构承载能力分析[J]. 中国机械工程, 2015, 26(7): 871-877.
Yuan F H, Zhao T S, Zhao Y Z, et al. Analysis of load carrying capacity of parallel mechanism[J]. China Mechanical Engineering, 2015, 26(7): 871-877. DOI:10.3969/j.issn.1004-132X.2015.07.004
[3]
何玉庆, 赵忆文, 韩建达, 等. 与人共融——机器人技术发展的新趋势[J]. 机器人产业, 2015(5): 74-80.
He Y Q, Zhao Y W, Han J D, et al. Co-existence with humans-The new trend of robot technology development[J]. Robot Industry, 2015(5): 74-80.
[4]
Kim Y L, Ahn K H, Song J B. Direct teaching algorithm based on task assistance for machine tending[C]//International Conference on Ubiquitous Robots and Ambient Intelligence. Piscataway, USA: IEEE, 2016: 861-866.
[5]
游有鹏, 张宇, 李成刚. 面向直接示教的机器人零力控制[J]. 机械工程学报, 2014, 50(3): 10-17.
You Y P, Zhang Y, Li C G. Force-free control for the direct teaching of robots[J]. Journal of Mechanical Engineering, 2014, 50(3): 10-17.
[6]
侯澈, 王争, 赵忆文, 等. 面向直接示教的机器人负载自适应零力控制[J]. 机器人, 2017, 39(4): 439-448.
Hou C, Wang Z, Zhao Y W, et al. Load adaptive force-free control for the direct teaching of robots[J]. Robot, 2017, 39(4): 439-448.
[7]
张铁, 洪景东, 刘晓刚, 等. 基于弹性摩擦模型的机器人免力矩传感器拖动示教方法[J]. 农业机械学报, 2019, 50(1): 419-427.
Zhang T, Hong J D, Liu X G, et al. Dragging teaching method without torque sensor for robot based on elastic friction model[J]. Transactions of the Chinese Society of Agricultural Machinery, 2019, 50(1): 419-427.
[8]
李智清, 叶锦华, 吴海斌, 等. 基于卷积力矩观测器与摩擦补偿的机器人碰撞检测[J]. 浙江大学学报(工学版), 2019, 53(3): 427-434.
Li Z Q, Ye J H, Wu H B, et al. Robot collision detection with convolution torque observer and friction compensation[J]. Journal of Zhejiang University: Engineering Science, 2019, 53(3): 427-434.
[9]
宋吉来, 曲道奎, 徐方, 等. 机器人无力传感器主动柔顺控制研究[J]. 电机与控制学报, 2020, 24(8): 159-166.
Song J L, Qu D K, Xu F, et al. Active compliance control research for robots without force sensor[J]. Electric Machines and Control, 2020, 24(8): 159-166.
[10]
吴文祥, 朱世强, 靳兴来. 基于改进傅里叶级数的机器人动力学参数辨识[J]. 浙江大学学报(工学版), 2013, 47(2): 231-237.
Wu W X, Zhu S Q, Jin X L. Dynamic identification for robot manipulator based on modified Fourier series[J]. Journal of Zhejiang University: Engineering Science, 2013, 47(2): 231-237.
[11]
张春涛, 王勇, 穆春阳, 等. 基于快速动力学辨识的机器人力/位混合控制碰撞检测研究[J]. 仪器仪表学报, 2021, 42(6): 161-170.
Zhang C T, Wang Y, Mu C Y, et al. Collision detection for robot force position hybrid control based on fast dynamics identification[J]. Chinese Journal of Scientific Instrument, 2021, 42(6): 161-170.
[12]
席万强, 陈柏, 丁力, 等. 考虑非线性摩擦模型的机器人动力学参数辨识[J]. 农业机械学报, 2017, 48(2): 393-399.
Xi W Q, Chen B, Ding L, et al. Dynamic parameter identification for robot manipulators with nonlinear friction model[J]. Transactions of the Chinese Society of Agricultural Machinery, 2017, 48(2): 393-399.
[13]
Dasgupta B, Choudhury P. A general strategy based on the Newton-Euler approach for the dynamic formulation of parallel manipulators[J]. Mechanism and Machine Theory, 1999, 34(6): 801-824.
[14]
王洪波, 黄真. 六自由度并联式机器人拉格朗日动力方程[J]. 机器人, 1990, 12(1): 23-26.
Wang H B, Huang Z. Lagrange's dynamic equation of six-DOF parallel multi-loop robot manipulator[J]. Robot, 1990, 12(1): 23-26.
[15]
陈根良, 王皓, 来新民, 等. 基于广义坐标形式牛顿-欧拉方法的空间并联机构动力学正问题分析[J]. 机械工程学报, 2009, 45(7): 41-48.
Chen G L, Wang H, Lai X M, et al. Forward dynamics analysis of spatial parallel mechanisms based on the Newton-Euler method with generalized coordinates[J]. Chinese Journal of Mechanical Engineering, 2009, 45(7): 41-48.
[16]
Park E C, Lim H. Position control of X-Y table at velocity reversal using B presliding friction characteristics[J]. IEEE Transactions on Control Systems Technology, 2003, 11(1): 24-31.
[17]
Makkar C, Dixon W E, Sawyer W G, et al. A new continuously differentiable friction model for control system design[C]//IEEE/ASME International Conference on Advanced Intelligent Mechatronics. Piscataway, USA: IEEE, 2005: 600-605.
[18]
Liu L L, Wu Z Y. A new identification method of the Stribeck friction model based on limit cycles[J]. Journal of Mechanical Engineering Science, 2014, 228(15): 2678-2683.
[19]
Xi Z Y, Lu C H, Zhang T. Two methods of friction torque measurement in an alternating current servo worktable[C]//Proceedings of New Advances in Manufacturing Engineering. Xi'an: Northwestern Polytechnical University Press, 2006: 19-25.
[20]
Xi Z Y, Lu C H, Zhang J C, et al. Study on relationships between characteristics of multi-signals related to friction in AC servo feed systems[C]//International Conference on Mechanical Engineering and Mechanics. New York, USA: Science Press USA Inc., 2005: 563-567.
[21]
Makarov B M. Calculation of the volume of multidimensional cones and balls using Cavalieri's principle[J]. American Mathematical Monthly, 2012, 119(8): 685-689.