CN114800519B - 一种考虑摩擦的六自由度工业机器人动力学参数辨识方法 - Google Patents

Abstract

本发明公开了一种考虑摩擦的六自由度工业机器人动力学参数辨识方法，属于机器人动力学领域。该方法主要包括以下步骤：基于拉格朗日法构建考虑摩擦的六自由度工业机器人的动力学方程；对六自由度工业机器人学方程进行线性化和参数重组；设计并优化基于傅里叶级数展开的Freudenstein 1‑3‑5六自由度工业机器人的激励轨迹；采集六自由度工业机器人的各个关节的位置信息和各个关节的驱动电流并平滑滤波，计算六自由度工业机器人的各个关节的速度、加速度和采样力矩；使用改进的哈里斯鹰优化算法求解动力学参数。本发明提高了机器人动力学参数辨识的精度和效率，为实现基于动力学模型的机器人高精度控制奠定基础。

Description

一种考虑摩擦的六自由度工业机器人动力学参数辨识方法

技术领域

本发明涉及机器人动力学领域，特别涉及一种考虑摩擦的六自由度工业机器人动力学参数辨识方法。

背景技术

随着工业机器人向高速度化、高精度化、智能化、数字化方向发展，提升机器人的运动性能和控制品质显得尤为重要。机器人系统具有多输入多输出、强耦合、非线性以及时变等特点，并且，工业机器人负载惯量大、存在非线性摩擦、不确定性扰动以及机械装配误差等因素影响其运动性能，基于动力学的控制更需要精确的模型，而构建准确的摩擦力模型对提升动力学模型精度具有重要作用。获取机器人动力学模型参数的有效方法主要是进行机器人动力学辨识实验。

动力学参数辨识的主要步骤是建立动力学模型、进行参数线性化与重组、设计激励轨迹及优化参数、采集与处理实验数据、设计参数辨识算法、进行模型验证，具体流程如图1所示。建立动力学模型的主要方法有牛顿-欧拉法、拉格朗日法、凯恩方法和旋量理论等。现有的摩擦模型包含静态摩擦模型和动态摩擦模型。静态摩擦模型参数的辨识需在高精度低速跟踪条件下进行。而动态摩擦模型虽能较好描述摩擦在低速时的特性，但往往是不连续的或者分段连续的。其次，激励轨迹的设计影响辨识参数的可信度。目前工业机器人动力学参数辨识大都采用傅里叶级数型的轨迹，但是，傅里叶级数不满足速度、加速度的边界条件，容易在轨迹起点和终点产生加速度突变，引起机器人振颤，不利于机器人平稳跟踪，影响辨识精度。在辨识算法上，传统的最小二乘法进行参数估计精度不高，且部分参数不符合物理可行性。

发明内容

本发明的目的在于克服上述不足，提供一种可以提高辨识精度、改善辨识效果的考虑摩擦的六自由度工业机器人动力学参数辨识方法。

本发明的一种考虑摩擦的六自由度工业机器人动力学参数辨识方法，包括以下步骤：

步骤一、建立考虑摩擦的六自由度工业机器人动力学模型，动力学方程表示如下：

其中采用Makkar模型,式中，q,/>分别为六自由度工业机器人的各个关节位移向量、速度向量和加速度向量，/>R^b为b维实数集，M(q)∈R^b×b，为六自由度工业机器人正定对称的惯性矩阵，/>为六自由度工业机器人的离心力与哥氏力矩阵，G(q)∈R^b为六自由度工业机器人的重力矩阵，/>为各个关节摩擦力矩向量，τ∈R^b为各个关节驱动力矩向量；

步骤二、对六自由度工业机器人动力学模型进行线性化和参数重组，得到如下公式：

其中，K_min为六自由度工业机器人待辨识的最小动力学参数集，包含六自由度工业机器人各连杆基本惯性参数和各关节摩擦系数，为六自由度工业机器人系统观测矩阵；

步骤三、采用基于傅里叶级数展开的Freudenstein 1-3-5模型

设计并优化六自由度工业机器人各个关节的激励轨迹，Freudenstein 1-3-5模型中边界条件设定为：

为观测矩阵最小条件数；q_min、q_max分别为六自由度工业机器人各个关节可转动位置的下限、上限，/>分别为六自由度工业机器人各个关节的速度上限和加速度上限，q(0)、/>分别为0时刻六自由度工业机器人关节位置、速度和加速度，q(t_f)、/>分别为时刻t_f时的六自由度工业机器人各个关节的位置、速度和加速度；

步骤四、采集六自由度工业机器人的各个关节的位置信息q和各个关节的驱动电流i_m，计算六自由度工业机器人的各个关节的速度加速度/>和采样力矩τ_c，以及观测矩阵/>

第一步，在机器人控制器中完成对步骤三获取的六自由度工业机器人各个关节的激励轨迹速度规划与插补，将插补后的六自由度工业机器人各个关节的激励轨迹输入到机器人驱动器中，驱动机器人关节转动；

第二步，在六自由度工业机器人按激励轨迹运行过程中采集六自由度工业机器人的各个关节的位置信息q和各个关节的驱动电流i_m，生成各个关节的位置曲线和驱动电流曲线，然后对位置曲线和驱动电流进行平滑滤波；

第三步，对滤波后的关节位置曲线拟合为改进的傅里叶级数，然后，对改进的傅里叶级数进行微分获得六自由度工业机器人的各个关节的速度和加速度/>

第四步，各个关节的采样力矩τ_c用下式估计得到：

τ_c＝k_ai_m

k_a为驱动电机转矩常数；

第五步，将本步骤第三步获得的六自由度工业机器人的各个关节的位置q、速度和加速度/>代入观测矩阵/>的具体表达式中，求解得到观测矩阵/>的结果；

步骤五、选用基于改进的哈里斯鹰优化算法进行动力学参数辨识，对公式进行求解得到K_min，其中理论力矩τ采用采样力矩τ_c代替。

本发明使用连续可微的摩擦模型、利用基于傅里叶级数展开的弗洛伊德斯坦(Freudenstein)1-3-5轨迹、使用改进的哈里斯鹰(MHHO)智能优化算法确定动力学参数，能一次辨识所有动力学参数，提高了辨识效率和精度。

附图说明

图1六自由度工业机器人动力学模型参数辨识流程图；

图2 Makkar摩擦模型曲线图；

图3哈里斯鹰优化算法流程图。

具体实施方式

下面结合附图和实施例对本发明加以详细说明。

如附图所示的本发明的一种考虑摩擦的六自由度工业机器人动力学参数辨识方法，包括以下步骤：

步骤一、建立考虑摩擦的六自由度工业机器人动力学模型，动力学方程表示如下：

其中，分别为六自由度工业机器人的各个关节位移向量、速度向量和加速度向量，/>R^b为b维实数集，M(q)∈R^b×b，为六自由度工业机器人正定对称的惯性矩阵，/>为六自由度工业机器人的离心力与哥氏力矩阵，G(q)∈R^b为六自由度工业机器人的重力矩阵，/>为各个关节摩擦力矩向量，τ∈R^b为各个关节驱动力矩向量。

具体地，

q＝[q₁,q₂,……,q_n]^T

G(q)＝[g₁,g₂,……,g_n]^T

式中，q_n为第n个关节的位移，m_nn为第n个关节的加速度产生对第n关节力矩的力矩分量，c_nn为第n个关节在速度/>时产生的离心力矩或科氏力矩，g_n为第n个关节受到的重力矩，/>为第n个关节受到的摩擦力矩。

其中采用Makkar模型，该模型的数学表达式具体参见Makkar C,Dixon W E,Sawyer W G,et al.A new continuously differentiable friction model for controlsystems design[C].Monterey:Proceedings of the IEEE/ASME InternationalConference on Advanced Intelligent Mechatronics,2005.(马卡,迪克逊,索耶,等.一种新的连续可微摩擦控制系统设计模型[C]，蒙特雷:IEEE/ASME高级智能机电一体化国际会议论文集,2005.)公式如下：

其中，γ_i∈R，i＝1,2…6为未知正常数。描述机器人关节的库伦摩擦特性，描述机器人关节的粘滞摩擦特性，γ₁+γ₄描述机器人关节的静摩擦特性，/>描述机器人关节的斯特里贝克效应，其摩擦曲线如图2所示，该模型为简易的、连续可微的、能集中摩擦主要特性的摩擦模型。而现有的库伦模型、粘滞模型、静摩擦模型、斯特里贝克(Stribeck)模型及其组合模型不连续或者分段连续，在高性能控制器的开发中存在某些问题，不利于高精度控制。

本动力学方程的推导过程如下：

第一步，给定六自由度工业机器人型号及参数，根据改进的D-H参数法建立六自由度工业机器人运动学模型，得到六自由度工业机器人的雅可比矩阵J和旋转矩阵T_R；

第二步，使用拉格朗日法对六自由度工业机器人进行动力学建模，过程为：首先，按照《机器人学导论》(机械工业出版社,2006)中关于绘制机器人的连杆模型原则建立六自由度工业机器人连杆模型，在连杆模型上定义六自由度工业机器人的广义坐标系，在连杆模型上选定独立的广义关节变量q及对应的广义力τ；

第三步，将六自由度工业机器人系统的拉格朗日函数定义为整个六自由度工业机器人系统总动能/>减去系统总势能p(q)，即：

第四步，对拉格朗日函数利用下式进行求导，

整理可得机器人动力学方程：

步骤二、对六自由度工业机器人动力学模型进行线性化和参数重组，得到如下公式：

其中，K_min为六自由度工业机器人待辨识的最小动力学参数集，包含六自由度工业机器人各连杆基本惯性参数和各关节摩擦系数，为六自由度工业机器人系统观测矩阵。

本步骤具体过程如下：

第一步，将六自由度工业机器人动力学模型导入数学计算软件(如可以采用Matlab符号计算工具箱或者软件Symoro中编程获得)中；

第二步，对六自由度工业机器人动力学模型中除之外的各项参数进行QR分解(具体参见《矩阵分析与计算，国防工业出版社，2010年》)得到六自由度工业机器人动力学模型的一次线性形式，即：

其中，为六自由度工业机器人系统回归矩阵，K∈R^10b×1为六自由度工业机器人惯性参数集，且K＝[k₁,k₂,…,k_i,…,k_n]。连杆惯性参数可表述为：

k_i＝[m_i,m_ir_xi,m_ir_yi,m_ir_zi,I_xxi,I_xyi,I_xzi,I_yyi,I_yzi,I_zzi]^T

式中，k_i∈R¹⁰，为六自由度工业机器人惯性参数，m_i为第i根连杆的质量，r_xi、r_yi、r_zi分别为第i根连杆的质心在六自由度工业机器人的广义坐标系下的x、y、z轴分坐标，m_ir_xi、m_ir_yi、m_ir_zi分别为第i根连杆关于x、y、z轴的一阶质量矩。I_xxi、I_xyi、I_xzi为第i根连杆绕x轴对x轴、y轴、z轴产生的惯性矩，I_yyi、I_yzi为第i根连杆绕y轴对y轴、z轴产生的惯性矩，I_zzi为第i根连杆绕z轴对z轴产生的惯性矩。

第三步，对进行线性化，得到六自由度工业机器人动力学模型的二次线性公式为：/>其中，/>为六自由度工业机器人为六自由度工业机器人系统回归矩阵，K′为六自由度工业机器系统动力学参数矩阵。

对进行线性化的过程为：

根据对机器人摩擦特性的研究结果可获知γ₂,γ₃,γ₅取值，为已知值，γ₂,γ₃,γ₅分别为105.8、10.0、93.6。因此，摩擦力矩向量可线性表示为：

第四步，动力学方程全部线性化为:

式中，a_i、b_i、c_i为系统回归矩阵/>中的元素，/>

γ_i1、γ_i2、γ_i3、γ_i4、γ_i5、γ_i6为第i关节受到的摩擦力矩中的常系数，i＝1,2,……,n。

第五步，由于回归矩阵中含有线性相关列，通常不满秩。因此，对回归矩阵/>中各列线性重组相消和去除全为0列，得到满秩的观测矩阵/>对六自由度工业机器系统动力学参数矩阵K′进行参数重组和线性化得到最小动力学参数集K_min。此步中涉及到对观测矩阵/>和最小动力学参数集K_min的求解可在软件symbolicmath toolbox中编程实现，得到观测矩阵/>和最小动力学参数集K_min的具体表达式，具体推导参见《矩阵分析与计算，国防工业出版社，2010年》。从而得到如下表达式：

步骤三、采用基于傅里叶级数展开的弗洛伊德斯坦Freudenstein 1-3-5模型(具体参见Biagiotti,L.,&Melchiorri,C.(2008).Trajectory planning for automaticmachines and robots.比亚乔蒂，L.，&梅尔奇奥里，C.(2008)。《自动机械和机器人的轨迹规划》)设计并优化六自由度工业机器人各个关节的激励轨迹，Freudenstein 1-3-5模型中边界条件设定为：

为观测矩阵最小条件数，可在Matlab中编程求值。q_min、q_max分别为六自由度工业机器人各个关节可转动位置的下限、上限，/>分别为六自由度工业机器人各个关节的速度上限和加速度上限，q(0)、/>分别为0时刻六自由度工业机器人关节位置、速度和加速度，q(t_f)、/>分别为时刻t_f时的六自由度工业机器人各个关节的位置、速度和加速度。

步骤四、采集六自由度工业机器人的各个关节的位置信息q和各个关节的驱动电流i_m，计算六自由度工业机器人的各个关节的速度加速度/>和采样力矩τ_c，以及观测矩阵

第一步，在机器人控制器中完成对步骤三获取的六自由度工业机器人各个关节的激励轨迹速度规划与插补，将插补后的六自由度工业机器人各个关节的激励轨迹输入到机器人驱动器中，驱动机器人关节转动。

第二步，在六自由度工业机器人按激励轨迹运行过程中采集六自由度工业机器人的各个关节的位置信息q和各个关节的驱动电流i_m，生成各个关节的位置曲线和驱动电流曲线，然后对位置曲线和驱动电流进行平滑滤波；

第三步，对滤波后的关节位置曲线拟合为改进的傅里叶级数，然后，对改进的傅里叶级数进行微分获得六自由度工业机器人的各个关节的速度和加速度/>

本步骤中对收集到的位置曲线和电流曲线可以采用Savitzky-Golay滤波器进行平滑滤波。

第四步，各个关节的采样力矩τ_c可用下式估计得到：

τ_c＝k_ai_m

k_a为驱动电机转矩常数，通过电机手册查找可知。

第五步，将本步骤第三步获得的六自由度工业机器人的各个关节的位置q、速度和加速度/>代入观测矩阵/>的具体表达式中，求解得到观测矩阵/>的结果。

步骤五、选用基于改进的哈里斯鹰优化算法(MHHO)(具体参见刘骏鹏.哈里斯鹰算法的改进及应用研究[D].浙江大学,2021.)进行动力学参数辨识，对公式进行求解得到K_min，由于无法直接获得理论力矩τ，采样力矩τ_c代替理论力矩τ来进行计算，采样力矩τ_c以及观测矩阵/>采用步骤四中获得的数值；

该算法可以在Matlab中编程实现。如图3所示，作为本发明的一种实施方式，K_min的求解过程如下：

第一步，将改进的哈里斯鹰优化算法(MHHO)编程写入Matlab程序，步骤如下：

设定哈里斯鹰种群的参数，然后对哈里斯鹰种群初始化得到每一只哈里斯鹰的初始位置，参数包括种群数量η、优化问题维度d以及最大迭代次数T；

哈里斯鹰种群E_l＝(X₁,X₂,...,X_η)^T，哈里斯鹰初始位置矩阵为：

其中，d为优化问题维度，X_η,d表示第η只哈里斯鹰第d个维度，最小动力学参数集K_min为d×1维矩阵，将最小动力学参数集K_min中的d各元素依次代入上式得到待优化的初始位置矩阵

第二步，将哈里斯鹰优化算法中适应度函数编程写入Matlab程序；

适用度函数F(X)表述如下：

适用度函数F(X)中的元素f(x)计算式如下：

其中，m为采样次数，τ_c为机器人各关节采样力矩，τ为机器人各关节理论力矩，x为哈里斯鹰初始位置矩阵中的各元素。

第三步，在程序中调整哈里斯鹰的种群位置个体位置并更新兔子位置与逃逸能量。根据调整后种群的新位置，依次遍历每个兔子个体的适应度值(在本例中，适应度值越低意味着适应度值越优)，找出哈里斯鹰最佳适应度值，并且每次迭代都需要更新兔子的逃逸能量，能量是一个在-2到2之间变化的值，随着迭代次数增加而自适应衰减，哈里斯鹰能够根据兔子猎物的能量进行策略选择。

步骤501，计算兔子当前的逃逸能量E。

其中，t为当前迭代次数，E₀为(-1,1)之间的随机数；

步骤502，判断E的大小，若|E|≥1，则执行步骤503，若|E|<1则执行步骤504。

步骤503，哈里斯鹰进行全局搜索，寻找兔子。根据发现或未发现兔子，随机产生数σ，并得到下述搜索公式：

其中，X(t+1)表示下次迭代时哈里斯鹰的位置向量，X_r(t)是兔子的位置向量，X_rand(t)表示鹰群的随机个体的位置向量，UB和LB为优化问题维度的上限、下限，X(t)为该鹰的当前位置向量，r₁,r₂,r₃,r₄和σ为(0,1)之间的随机数，X_m(t)是哈里斯鹰种群的当前平均位置，计算公式如下：

X_ρ(t)表示第t次迭代时每只哈里斯鹰的位置。

步骤504，当0<|E|<1时，哈里斯鹰进行局部搜索，用μ<0.5表示兔子逃脱成功，μ≥0.5表示兔子逃脱失败；然后判断E的绝对值是否小于0.5,μ是否小于0.5。若0.5≤|E|<1、μ≥0.5，则执行步骤505，否则执行步骤506。若0<|E|<0.5、μ≥0.5，则执行步骤507，否则执行步骤508。

步骤505，当0.5≤|E|<1、μ≥0.5时，哈里斯鹰执行盘旋围捕策略，哈里斯鹰位置向量更新如下：

X(t+1)＝△X(t)-E|JX_r(t)-X(t)|

△X(t)＝X_r(t)-X(T)

其中，△X(t)为当前迭代次数为t时哈里斯鹰与兔子的位置之差，随机数J∈(0,2)。

步骤506，当0.5≤|E|<1、μ<0.5时，哈里斯鹰执行盘旋围捕与渐进俯冲攻击策略，哈里斯鹰位置向量更新如下：

其中，S为d维优化问题随机行向量，矩阵Y＝X_r(t)-E|JX_r(t)-X(t)|，矩阵Z＝Y+S×LF(d)，将矩阵Y、Z、X(t)代入步骤五第二步中的适用度函数F(X)中，从而求得F(Y)、F(X(t))、F(Z)，LF(d)为调用的Matlab中飞行函数。

步骤507，当0<|E|<0.5、μ≥0.5时，哈里斯鹰执行强势突袭策略，哈里斯鹰位置向量更新如下：

X(t+1)＝X_r(t)-E|△X(t)|

步骤508，当0<|E|<0.5、μ<0.5时，哈里斯鹰执行强势突袭与渐进俯冲攻击策略，哈里斯鹰位置向量更新如下：

第四步，判断当前迭代次数t是否大于或者等于最大迭代次数T或者达到迭代精度，如果是，则输出当前哈里斯鹰位置矩阵，即为所需辨识动力学参数K_min，否则，搜索继续。

以上对本发明的描述仅仅是示意性的，而不是限制性的，所以，本发明的实施方式并不局限于上述的具体实施方式。如果本领域的普通技术人员受其启示，在不脱离本发明宗旨和权利要求所保护范围的情况下，做出其他变化或变型，均属于本发明的保护范围。

Claims (1)

1.一种考虑摩擦的六自由度工业机器人动力学参数辨识方法，其特征在于包括以下步骤：

步骤一、建立考虑摩擦的六自由度工业机器人动力学模型，动力学方程表示如下：

其中采用Makkar模型,式中，/>分别为六自由度工业机器人的各个关节位移向量、速度向量和加速度向量，/>R^b为b维实数集，M(q)∈R^b×b，为六自由度工业机器人正定对称的惯性矩阵，/>为六自由度工业机器人的离心力与哥氏力矩阵，G(q)∈R^b为六自由度工业机器人的重力矩阵，/>为各个关节摩擦力矩向量，τ∈R^b为各个关节驱动力矩向量；

步骤二、对六自由度工业机器人动力学模型进行线性化和参数重组，得到如下公式：

其中，K_min为六自由度工业机器人待辨识的最小动力学参数集，包含六自由度工业机器人各连杆基本惯性参数和各关节摩擦系数，为六自由度工业机器人系统观测矩阵；

步骤三、采用基于傅里叶级数展开的Freudenstein 1-3-5模型

设计并优化六自由度工业机器人各个关节的激励轨迹，Freudenstein 1-3-5模型中边界条件设定为：

为观测矩阵最小条件数；q_min、q_max分别为六自由度工业机器人各个关节可转动位置的下限、上限，/>分别为六自由度工业机器人各个关节的速度上限和加速度上限，/>分别为0时刻六自由度工业机器人关节位置、速度和加速度，/>分别为时刻t_f时的六自由度工业机器人各个关节的位置、速度和加速度；

步骤四、采集六自由度工业机器人的各个关节的位置信息q和各个关节的驱动电流i_m，计算六自由度工业机器人的各个关节的速度加速度/>和采样力矩τ_c，以及观测矩阵

第一步，在机器人控制器中完成对步骤三获取的六自由度工业机器人各个关节的激励轨迹速度规划与插补，将插补后的六自由度工业机器人各个关节的激励轨迹输入到机器人驱动器中，驱动机器人关节转动；

第二步，在六自由度工业机器人按激励轨迹运行过程中采集六自由度工业机器人的各个关节的位置信息q和各个关节的驱动电流i_m，生成各个关节的位置曲线和驱动电流曲线，然后对位置曲线和驱动电流进行平滑滤波；

第三步，对滤波后的关节位置曲线拟合为改进的傅里叶级数，然后，对改进的傅里叶级数进行微分获得六自由度工业机器人的各个关节的速度和加速度/>

第四步，各个关节的采样力矩τ_c用下式估计得到：

τ_c＝k_ai_m

k_a为驱动电机转矩常数；

第五步，将本步骤第三步获得的六自由度工业机器人的各个关节的位置q、速度和加速度/>代入观测矩阵/>的具体表达式中，求解得到观测矩阵/>的结果；

步骤五、选用基于改进的哈里斯鹰优化算法进行动力学参数辨识，对公式进行求解得到K_min，其中理论力矩τ采用采样力矩τ_c代替。