武诗睿, 吴丹
清华大学 机械工程系, 北京 100084
收稿日期：2022-03-25
基金项目：国家自然科学基金面上项目(51575306)
作者简介：武诗睿(1998—)，女，博士研究生
通讯作者：吴丹，教授，E-mail: wud@tsinghua.edu.cn

摘要：LuGre模型可以较为全面地描述动态摩擦现象, 但引入了不可观测的鬃毛变形, 给LuGre摩擦参数的高效准确辨识带来技术挑战。传统的LuGre摩擦参数辨识法需要运动系统采用力矩控制, 且辨识工作量较大。传统方法不再适用一些系统, 如多自由度机械臂这种常采用位置控制模式的系统。该文提出了一种基于摩擦力矩—速度曲线特定区域形状分析的LuGre摩擦参数辨识法。定义形状因子来定量表示Stribeck峰和迟滞环线的形状特征, 并用其辅助辨识, 可得到比粒子群算法局部最优解更好的辨识结果。仿真分析和机械臂辨识实验结果表明：与纯PSO法相比, 该方法的辨识实验量更小, LuGre摩擦参数辨识更准确, 关节力矩预测精度更高。
关键词：LuGre摩擦模型参数辨识机械臂
Parameter identification for the LuGre friction model based on an area-specific shape analysis of the friction torque-velocity curve
WU Shirui, WU Dan
Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China

Abstract: The LuGre model is an advanced friction model that describes dynamic friction characteristics. However, the imported unobservable bristle deformation complicates accurate, efficient identification of the model parameters. The traditional LuGre parameter identification method requires the motion system to use torque control and requires a significant computational load. The traditional LuGre parameter identification method is not applicable to some systems, such as multiple degree-of-freedom manipulators, that use the position control mode. Therefore, this paper presents a modified LuGre parameter identification method based on an area-specific analysis of the friction torque-velocity curve. Shape factors are defined to quantify the shape features of the Stribeck peak and the hysteresis loop which are then used for the LuGre parameter identification. The identification result is better than the local optimal solution of the PSO method. Simulations and hardware identification tests on a manipulator verify the effectiveness of this method, which requires fewer experiments, has better parameter identification accuracy and more accurately predicts the manipulator joint torque than the pure PSO method.
Key words: LuGre friction modelparameter identificationmanipulator
摩擦是运动系统固有的复杂特性，了解和定量表达摩擦特性对于提高运动控制系统的性能具有重要意义。传统的Coulomb黏滞摩擦模型是控制系统最常用的模型^[1]，但在系统启停时常引起较大力矩预测误差。目前研究人员对摩擦模型的改进主要有3类方法。第1类是增加更多速度相关项。例如，Bompos等^[2]使用改进的Stribeck模型描述Coulomb摩擦，Wernholt等^[3]和Han等^[4]分别使用三角非线性函数模型和更高阶的速度项模型来描述黏滞摩擦的变化。但这些模型仍认为摩擦只与速度相关。第2类是考虑除速度以外的影响摩擦的更多因素，如力矩、温度与角度位置^[5-7]等因素对摩擦的影响，这使得辨识参数增多，实验量与数据量增大，且要求配置多个传感器。第3类是摩擦动态特性建模，目前最受认可的为LuGre模型^[8-10]。它是一种典型的动态摩擦模型，其中加入了不可观测却代表内部摩擦状态的鬃毛变形量，可以描述较完整的动态摩擦效应，如预滑动、Stribeck效应、迟滞效应等。这些效应主要表现在系统启停与速度极低时。因此，用LuGre模型可以更为准确地描述系统启停或极低速运动时的摩擦特性，这是前两类模型无法做到的。
然而，LuGre摩擦模型是一个包含较多参数的非线性方程组，其参数辨识具有一定难度。传统方法分稳态与动态2步辨识。首先在稳态匀速运动时进行多组不同速度的实验，辨识与Stribeck曲线相关的稳态参数。然后使用远小于起步阻力(静摩擦力)的驱动力矩，使整体处于预滑动状态，完成动态参数(即鬃毛刚度与阻尼系数)的辨识^[11-14]。已有工作主要集中于动态参数辨识方法的研究，早期以传递函数参数拟合方法为代表^[11]，近年来智能搜索算法有更多应用，如粒子群算法(PSO)、遗传算法(GA)、蚁群算法^[12-14]等。这种辨识思路应用广泛，但实验较复杂，且PSO等搜索算法容易陷入局部最优。
尽管上述方法可以完成LuGre摩擦模型参数辨识，但前提条件是运动系统必须采用力矩控制模型，这就限制了方法的应用范围。例如目前广泛使用的多自由度机械臂，各关节的运动控制主要采用位置控制模式，难以构建传统方法中LuGre动态参数辨识所需的条件，即将驱动力矩控制在静摩擦力以下。而且，当需要同时辨识多个关节摩擦参数时，传统方法只能每个关节逐一辨识，带来辨识效率低下的问题。若为了提高效率只使用PSO进行全部参数的搜寻，则会由于LuGre方程强烈非凸，使搜索算法极易落入局部最优。
本文提出一种基于摩擦力矩—速度曲线特定区域形状分析的LuGre摩擦参数辨识方法，以解决以机械臂为典型代表的采用位置控制模式的多自由度运动系统摩擦参数高效高精度辨识问题。深入分析LuGre模型在摩擦力矩—速度曲线上Stribeck峰与迟滞环线的形状随参数变化的规律，从中提取有代表性的形状因子辅助辨识，并开展了仿真分析和实验验证。
1 LuGre摩擦模型本文建立的LuGre摩擦模型为：

$F(v)=\sigma_{0} z+\sigma_{1}(v) \frac{\mathrm{d} z}{\mathrm{~d} t}+f(v),$

(1)

$\frac{\mathrm{d} z}{\mathrm{~d} t}=v-\sigma_{0} \frac{|v|}{g(v)} z, $

(2)

$g(v)=\alpha_{0}+\alpha_{1} \mathrm{e}^{-\left(v / v_{0}\right)^{2}}, $

(3)

$\begin{equation}\sigma_{1}(v)=\sigma_{1} \mathrm{e}^{-\left(v / v_{\mathrm{d}}\right)^{2}}, \end{equation}$

(4)

$f(v)=\alpha_{0} v+\alpha_{3} \operatorname{sgn}(v) v^{2} .$

(5)

其中：自变量v表示运动速度，m·s^-1；z为鬃毛变形量，可以类比Dahl模型中的应力。可变参数σ₀为鬃毛刚度，N·m^-1；σ₁为鬃毛阻尼，N·s·m^-1；α₀为Coulomb摩擦力，N；α₁为静摩擦力与Coulomb摩擦力之差(本文称其为tribeck峰高)，N；v₀为Stribeck速度参数，m·s^-1；v_d为阻尼速度参数，m·s^-1；α₂和α₃为黏滞摩擦速度参数，N·s²·m^-2。
式(1)中，F(v)为模型计算出的摩擦力，σ₀和σ₁决定鬃毛变形带来的摩擦力，f(v)为与速度相关的黏滞摩擦(见式(5))。为了更准确地描述较大速度范围系统的黏滞摩擦，本文改进了f(v)的表达式，在常用的线性项α₂v的基础上增加了速度二次项α₃sgn(v)v²。
式(3)中，g(v)描述Stribeck效应，即在速度极小的一段区间内，摩擦力随速度减小而增大。式中α₀为速度很高时g(v)的稳定值，由于当v＞2v₀时，g(v)=α₀+0.018α₁ $\simeq $α₀，因此在|v| < v₀内Stribeck效应明显。
在许多LuGre模型应用中，认为σ₁为常数^[11-13]，本文采用了文[10]的表述，认为σ₁为变量且随速度增大而衰减。在文[10]中v_d只是为了保证模型的摩擦力计算值在|v|?v₀后的稳定性而存在，可以认为|v_d|?v₀，对摩擦力—速度曲线影响很小。
鬃毛变形量是一个宏观上的小量^[9]，式(1)等号右侧前2项σ₀z+σ₁(v)dz/dt描述了鬃毛在预滑动下类似弹簧的动态效应，且有

$\begin{array}{c}\sigma_{0} z+\sigma_{1}(v) \mathrm{d} z / \mathrm{d} t \simeq \alpha_{0} \operatorname{sgn}(v), \\|v| \gg v_{0} \text { 且 }|v| \gg v_{\mathrm{d}} .\end{array}$

这是因为由式(3)和(4)明显可得：g(v)在|v|?v₀成立时有g(v)$\simeq $α₀，σ₁(v)在|v|?v_d成立时有σ₁(v)$\simeq $0。式(2)中dz/dt在|v|增大过程中由于z的增大而下降，在|v|?v₀后(即g(v)$\simeq $α₀)z趋于稳定，即dz/dt$\simeq $0，从而由式(2)变形后推出|g(v)|$\simeq $σ₀z。综合考虑关系式g(v)$\simeq $α₀和|g(v)|$\simeq $σ₀z，得到σ₀z$\simeq $α₀sgn(v)。此时式(1)退化为

$\begin{equation}F(v)=\alpha_{0} \operatorname{sgn}(v)+\alpha_{2} v+\alpha_{3} \operatorname{sgn}(v) v^{2} .\end{equation}$

(6)

式(6)是普通的只考虑速度效应的Coulomb模型。这表明LuGre模型在速度可观(即|v|?v₀)与速度近0时分别由不同的部分主导：速度近0时，鬃毛变形占绝对主导，二次速度项由于速度太小而几乎无贡献；在速度可观时，鬃毛变形稳定，摩擦表达可退化为式(6)。
2 辨识方法本文提出了形状因子的概念以及基于形状因子的辨识方法。形状因子定义为表征摩擦力矩—速度曲线特定区域形状特征的参数。该辨识方法只需让待辨识的运动系统按预设的一条辨识轨迹运动(位置控制)，即可辨识出摩擦参数，且能避免搜索算法落入局部最优。此方法的核心在于分析真实的摩擦力矩—速度曲线和辨识过程中拟合的摩擦力矩—速度曲线的形状特征，将真实曲线与拟合曲线的形状因子之差作为引导与反馈，不断修正初始解，最终使相邻2次迭代中待求参数之差收敛到指定阈值，得到最终的辨识解。由于选择的形状因子与物理模型有直接联系，故可将落入局部最优的初始解引导至全局最优解。
2.1 LuGre模型的摩擦力矩—速度曲线分析式(1)—(5)描述的LuGre模型共有8个参数。其中σ₀、σ₁、α₀、α₁、v₀、v_d与微观鬃毛直接相关，这6个参数决定了极低速(预滑动)区中摩擦力矩—速度曲线的形状。剩下的2个参数α₂、α₃的作用体现在|v|?v₀的速度可观区域，因预滑动区速度极小，可忽略这2个参数在预滑动区对摩擦的影响。
为了辅助分析LuGre低速形状参数与摩擦力矩—速度曲线形状的对应关系及影响，本文以多自由度机械臂为对象，建立了仿真分析模型。对机械臂而言，各关节LuGre摩擦模型中的速度v= $ \dot{q}$，位移x=q，其中q表示关节角位移。
辨识采用工业上典型的S型加减速轨迹^[15]，一条速度曲线示例如图 1所示。选择S型加减速轨迹的原因一是它在工业上较为常用，不用人工专门编程；二是它具有较大的速度区间，且在速度为0时加速度也为0(即在启停时有较长时间的速度近零)，便于LuGre模型中摩擦动态特征的体现与激励。仿真所用参数示例见表 1，其中仿真数据量级的选择与最初提出LuGre模型的文[9]一致。

图 1 S型加减速轨迹速度曲线示例

图选项

表 1 一组LuGre低速形状参数

σ0		σ1		α0		α1		v0		vd
N·m-1		N·s·m-1		N		N		m·s-1		m·s-1
350 000		800		8		3		0.01		0.08

表选项






图 2给出了该参数下图 1轨迹的摩擦力矩—速度曲线(低速段)，可见到速度从零增加时的迟滞环线(hysteresis loop)与速度减至零时的Stribeck峰^[8]。

图 2 摩擦力矩—速度曲线在不同参数下的对比

图选项

图 2中迟滞环线的形成类似应力应变曲线，z越大，变形率dz/dt就越小，最终z达到最大值，这过程中对应的摩擦力矩—速度曲线就是迟滞环线。dF/dv越大代表着迟滞环线斜率和上升速度越大，即在越小的速度区间内达到相同的摩擦，在图中的直观表现为迟滞环线越“窄”。
事实上，在相同加速情况下，迟滞环线宽度与鬃毛硬度和阻尼系数相关，下面通过推导进行解释。由式(2)变形得

$\begin{equation}\begin{array}{c}\frac{\mathrm{d} z}{\mathrm{~d} x} \frac{\mathrm{d} x}{\mathrm{~d} t}=\frac{\mathrm{d} z}{\mathrm{~d} x} v=v-\frac{\sigma_{0}|v|}{g(v)} z \Rightarrow \\\frac{\mathrm{d} z}{\mathrm{~d} x}=1-\frac{\sigma_{0} \operatorname{sgn}(v)}{g(v)} z.\end{array}\end{equation}$

(7)

迟滞环区速度极低，摩擦力主要源于鬃毛变形，且由于|vd|?v₀，在这个小区间内σ₁变化很小，f(v)可忽略，因此有

$\begin{equation}F(v)=\sigma_{0} z+\sigma_{1} \frac{\mathrm{d} z}{\mathrm{~d} t}\end{equation}.$

(8)

将式(8)对x求导，并代入式(7)可得

$\begin{array}{c}\frac{\mathrm{d} F}{\mathrm{~d} x}=\sigma_{0} \frac{\mathrm{d} z}{\mathrm{~d} x}+\sigma_{1} \frac{\mathrm{d}(\mathrm{d} z / \mathrm{d} x)}{\mathrm{d} t}= \\\sigma_{0}\left(1-\frac{F(v) \operatorname{sgn}(v)}{g(v)}\right)+\sigma_{1} \sigma_{0} \frac{\mathrm{d}(1 /|g(v)|)}{\mathrm{d} t} .\end{array}$

(9)

由式(9)可知在相同加速情况下，σ₀和σ₁与dF/dx成正比。在摩擦力矩—速度曲线中，dF/dx与dF/dv类似，都代表了摩擦随速度变化的速率，斜率越大，迟滞环线越窄。更大的σ₀与σ₁意味着更大的斜率(dF/dx值更大)与更窄的迟滞环线。从图 2的对比图也可得到验证。
速度方向改变前出现的Stribeck峰的形状主要由g(v)中的参数描绘。α₁决定“峰高”，v₀决定峰的“陡峭度”。由图 3中摩擦力矩—速度曲线局部放大图可见，不同参数组合下(只改变α₁和v₀)Stribeck峰的对比，可直观看到α₁越大则峰越高，v₀越小则峰越陡。

图 3 不同参数组合下Stribeck峰对比

图选项

以上分析说明了参数对曲线特征的主要影响区，作为提出形状因子概念与辨识法的基础。
2.2 形状因子提取与定义2.1节分析表明，摩擦力矩—速度曲线描述的静摩擦存在2个区域：一个是Stribeck峰区，其形状与α₁、v₀有关；另一个是迟滞环线区，其形状与σ₀、σ₁有关。因此，可以用参数定量描述2个区域摩擦特征。
2.2.1 Stribeck峰区形状因子Stribeck峰出现在速度趋近于零的轨迹段，其数值f_Stribeck等于总摩擦力减去只考虑速度效应的Coulomb摩擦力。Stribeck峰摩擦力可表示为

$f_{\text {Stribeck }}=F(v)-\left(\alpha_{0} \operatorname{sgn}(v)+\alpha_{2} v+\alpha_{3} \operatorname{sgn}(v) v^{2}\right) .$

(10)

在Stribeck峰区关注的是峰的高度和陡峭度，分别对应α₁和v₀。为了定量描述Stribeck峰，可对它与速度的关系曲线进行函数拟合。由于式(3)描述的Stribeck效应与速度为指数关系，故定义拟合函数为

$y=a \mathrm{e}^{-\left(\frac{v}{b}\right)^{2}} .$

(11)

其中y为f_Stribeck的拟合值。Stribeck峰在|v|?v₀后几乎消失，所以选择3(°)/s以下的区段进行拟合(该分界值是根据文[8-10]的系统参数而选定的)。
将拟合结果中a与b作为Stribeck区的形状因子，量纲与α₁、v₀相同。从图 3的实验结果中可直观看到，a和b分别用于摩擦力矩—速度曲线中“峰高”和“陡峭度”的简化描述，可以作为α₁和v₀的粗略估计。
2.2.2 迟滞环线区形状因子迟滞环线出现在速度从零逐渐增大的预滑动区，2.1节表明迟滞环线上力矩变化速率与σ₀和σ₁直接相关，这启发可以将F设法写成含有未知量σ₀、σ₁与已知量v、x的表达式，并尽量利用极低速的条件将其线性化。由式(3)可知，在v=0.3v₀时g(v)=α₀+0.91α₁，只比v=0时的g(v)减小了0.09α₁，所以为方便分析中的线性化，将迟滞环线区形状因子的求取区间设置在|v|?v₀内。本文通过初步实验和参考文[13]，将此区间设定在10^-3 rad/s内。
将式(2)变形得

$\begin{equation}\frac{1}{1-\frac{\sigma_{0}}{g(v)} z} \frac{\mathrm{d} z}{\mathrm{~d} t}=\frac{\mathrm{d} x}{\mathrm{~d} t} .\end{equation}$

(12)

可知在|v|?v₀小区间内可认为g(v)几乎不变，所以式(12)两边对时间t积分后得

$\frac{1}{1-\frac{\sigma_{0}}{g(v)} z}=\mathrm{e}^{\frac{\sigma_{0}}{g(v)} x+c} .$

(13)

由于存在未知的初始鬃毛小变形z₀，为简化表达，令c=e^z₀，式(13)可重写为

$\sigma_{0} z=g(v)\left(1-\mathrm{e}^{-\left(\frac{\sigma_{0}}{g(v)} x+c\right)}\right) .$

(14)

将式(14)代入式(8)，可得

$F(v)=C_{0}-g(v) \mathrm{e}^{-\sigma_{0} x / g(v)}+v \sigma_{1} \mathrm{e}^{-\sigma_{0} x / g(v)+c} .$

(15)

其中C₀为与变量无关的常数项。为便于计算，将式(15)的指数项做一阶Taylor展开，经过整理可得

$\Delta F \simeq\left[\begin{array}{llll}1 & x & v & v x\end{array}\right]\left[\begin{array}{c}C_{01} \\\sigma_{0} \mathrm{e}^{c} \\\sigma_{1} \mathrm{e}^{c} \\-\sigma_{1} \sigma_{0} \mathrm{e}^{c} / g(v)\end{array}\right] .$

(16)

其中C₀₁为常数项。从式(16)可见x与v的回归系数σ₀e^c和σ₁e^c之间无耦合，且它们分别与σ₀、σ₁的真实值呈一个模糊的比例关系，可以作为σ₀、σ₁的简化估计。所以定义σ₀e^c和σ₁e^c为迟滞环区的形状因子，量纲分别与σ₀和σ₁相同。
此节定义的形状因子虽不是真实的参数，但与真实值有直接关系且易于求取，所以可以作为反馈迭代的基础。
2.3 基于形状因子法的LuGre参数辨识2.3.1 速度可观区摩擦参数辨识根据1节中对式(6)的分析，将|v|?v₀的区域称为速度可观区。由于在速度可观后摩擦特性可以用式(6)描述，因此，可先在速度可观区将式(6)中的3个参数α₀、α₂、α₃辨识出来，再进行其他LuGre参数的辨识。
辨识α₀、α₂、α₃的方法如下：首先让系统按S型加减速轨迹运动，以获得摩擦力矩—速度曲线，然后在速度大于3(°)/s的摩擦力矩—速度曲线上均匀选取采样点，将采样点的摩擦力矩F(v)做[sgn(v) ??v ??sgn(v)v²]上的最小二乘回归，得到回归参数：

$\left[\begin{array}{lll}\alpha_{0} & \alpha_{2} & \alpha_{3}\end{array}\right]^{\mathrm{T}}=\operatorname{inv}\left(\boldsymbol{W}_{v}^{\mathrm{T}} \boldsymbol{W}_{v}\right) \boldsymbol{W}_{v}^{\mathrm{T}} \boldsymbol{F} .$

(17)

其中W_v为回归矩阵，其行数等于采样点个数，第i行的行向量为[sgn(v_i)??v_i??sgn(v_i)v_i²]；F(v)为各采样点的摩擦力矩组成的列向量。
2.3.2 基于形状因子的LuGre参数辨识对于其他的LuGre参数，可采用2.2节定义的形状因子对不精确的初值进行引导迭代。该方法中参数的初值通过简单的PSO法获得^[12]。下文中的真实曲线指实验采集的摩擦力矩—速度曲线，拟合曲线指利用当前辨识参数计算出的摩擦力矩—速度曲线。整个辨识思路如图 4所示。

图 4 形状因子法

图选项

借用自动控制理论的思想，真实曲线的形状因子作为“参考输入”，拟合曲线的形状因子作为“反馈量”，它们的差作为“控制量”，调整当前辨识参数，直至拟合曲线与真实曲线的形状因子趋近相同，则认为当前辨识参数也趋近了参数的真实值。事实上，真实与拟合下的摩擦力矩—速度曲线的差值是最直接的误差反馈，这也是纯PSO法最小化的目标函数，但整条曲线上误差反馈的物理意义比较模糊，没有提取出误差中最具价值的部分。本文的形状因子法通过分析曲线特定区域的形状特征，提取了整体误差中与待辨识参数误差最直接相关的部分，用于修正待辨识参数，所以可取得更好的效果。
在Stribeck峰区，设真实与拟合的摩擦力矩—速度曲线的形状因子向量分别为[a_r, b_r]^T和[a_s, b_s]^T，[a_r, b_r]^T在迭代中不变，[a_s, b_s]^T随着拟合曲线的变化而变化。将这2种形状因子之差乘以一个利于收敛的小于1的收缩因子w后作为a₁与v₀的增量进行迭代。记PSO的初始解为[α₁⁰, v₀⁰]^T，第k步形状因子为[a_s^k, b_s^k]^T，则从第k到k+1步迭代有

$\left[\begin{array}{c}\alpha_{1}^{k+1} \\v_{0}^{k+1}\end{array}\right]=\left[\begin{array}{c}\alpha_{1}^{k} \\v_{0}^{k}\end{array}\right]+w\left[\begin{array}{l}a_{\mathrm{r}}-a_{\mathrm{s}}^{k} \\b_{\mathrm{r}}-b_{\mathrm{s}}^{k}\end{array}\right] .$

当相邻2步迭代的辨识参数结果之差小于特定的ε₁，则认为迭代完成。由此可实现α₁、v₀的辨识。
在迟滞环线区，根据2.2.2节的分析，取速度在0~10^-3 rad/s内的区段进行分析。设真实轨迹和拟合轨迹中的迟滞环初始段摩擦力矩分别为F_r和F_s(ΔF=F_r-F_s)、鬃毛刚度分别为σ_0r和σ_0s、鬃毛阻尼分别为σ_1r和σ_1s，鬃毛变形分别为z_r和z_s，则根据式(8)得

$\left\{\begin{array}{l}F_{\mathrm{r}}=\sigma_{0 \mathrm{r}} z_{\mathrm{r}}+\sigma_{1 \mathrm{r}} \dot{z}_{\mathrm{r}}, \\F_{\mathrm{s}}=\sigma_{0 \mathrm{~s}} z_{\mathrm{s}}+\sigma_{1 \mathrm{~s}} \dot{z}_{\mathrm{s}}\end{array}\right.$

(18)

将式(18)按照式(16)展开并做差，经整理可得

$\Delta F \simeq\left[\begin{array}{llll}1 & x & v & v x\end{array}\right]\left[\begin{array}{c}C_{1} \\\sigma_{0 \mathrm{r}} \mathrm{e}^{c_{\mathrm{r}}}-\sigma_{0 \mathrm{~s}} \mathrm{e}^{c_{\mathrm{s}}} \\\sigma_{1 \mathrm{r}} \mathrm{e}^{c_{\mathrm{r}}}-\sigma_{1 \mathrm{~s}} \mathrm{e}^{c_{\mathrm{s}}} \\\frac{\sigma_{1 \mathrm{~s}} \sigma_{0 \mathrm{~s}} \mathrm{e}^{c_{\mathrm{s}}}-\sigma_{1 \mathrm{r}} \sigma_{0 \mathrm{r}} \mathrm{e}^{c_{\mathrm{r}}}}{g}\end{array}\right].$

(19)

其中C₁为变量无关的常数项，z_r0和z_s0为初始鬃毛小变形，c_r=e^z_r0，c_s=e^z_s0。由此可见式(19)中真实与拟合曲线的形状因子之差(σ_0re^c_r-σ_0se^c_s)和(σ_1re^c_r-σ_1se^c_s)分别代表了真实与拟合的鬃毛刚度和阻尼的差值，可以指导修正。
将采样点的ΔF进行[1 x v vx]上的最小二乘回归，即可得到回归参数。记PSO的初始解为[σ₀⁰ σ₁⁰]^T，第k步迭代中形状因子之差(σ₀₁e^c1－σ₀₂e^c2)和(σ₁₁e^c1－σ₁₂e^c2)分别记为$\begin{equation}\Delta \hat{\sigma}_{0}^{k}\end{equation} $和$\Delta \hat{\sigma}_{1}^{k} $，w为小于1的收缩因子，则第k到k+1步有

$\left[\begin{array}{c}\sigma_{0}^{k+1} \\\sigma_{1}^{k+1}\end{array}\right]=\left[\begin{array}{c}\sigma_{0}^{k} \\\sigma_{1}^{k}\end{array}\right]+\omega\left[\begin{array}{c}\Delta \hat{\sigma}_{0}^{k} \\\Delta \hat{\sigma}_{1}^{k}\end{array}\right] .$

同样地，当相邻2步迭代的辨识参数结果之差小于特定的ε₂，则认为迭代完成，由此实现σ₀、σ₁的辨识。
α₁、v₀、σ₀、σ₁辨识完成后，v_d作为最后一个未确定的参数，整个问题退化为单参数的凸优化问题，使用黄金分割法即可求出v_d的最优值。
综上，使用S型轨迹进行系统运动激励以辨识LuGre摩擦参数时，首先采用式(17)初步辨识α₀、α₂、α₃，然后进行其他摩擦参数辨识。此时，给定α₀、α₂、α₃，以1%的微小变化裕度进行小规模PSO初始搜索，获得α₁、v₀、σ₀、σ₁的初值，再采用基于形状因子的辨识法，经迭代收敛获得终值。最后利用黄金分割法优化获得v_d。
3 仿真分析与辨识实验3.1 仿真分析使用表 1的参数、图 1的S型轨迹及LuGre仿真模型生成的摩擦数据，各方法分别进行100次辨识实验。在采用纯PSO法时，使用100个个体经200代搜寻，未落入局部最优(参数辨识精度在1%以内)的实验占比为63%；采用形状因子法时，前期PSO使用50个体经50代搜索然后进入迭代，最终参数辨识精度在1%内的实验占比94%。表 2给出了形状因子法和纯PSO法的辨识结果，可以看到，形状因子法可以很大程度上避免结果落入局部最优，并且比纯PSO法具有更高的辨识精度。
表 2 真实值及形状因子法与纯PSO法的辨识结果对比

数据类型	σ0		σ1		α0		α1		v0		vd
	N·m-1		N·s·m-1		N		N		m·s-1		m·s-1
真实值	350 000		800.00		8.00		3.00		0.01		0.08
形状因子法辨识值	350 010		800.06		8.00		3.00		0.01		0.13
纯PSO辨识值	259 320		741.00		6.74		1.19		0.99		0.47

表选项






图 5a和5b分别给出了Stribeck和迟滞环区域在形状因子法前期PSO的结果，图 5c、5d为该条件下迭代5次后的结果，图 5中原曲线为LuGre仿真模型在真实参数下生成的摩擦数据，拟合曲线则是辨识结果下的数据。可以看出，非最优的前期PSO结果经过形状因子法的迭代引导后能很好地与真实曲线重合。

图 5 2个形状区域迭代前后对比

图选项

3.2 辨识实验3.2.1 辨识实验设计机械臂是一种典型的多轴运动系统，在工业生产中已得到广泛应用。精确辨识关节摩擦不仅可以准确预测关节力矩，而且有助于提高关节位置控制精度和机械臂启停或低速运动时的平稳性。
辨识实验采用实验室自主研制的七自由度机械臂(见图 6)。该机械臂采用RTX和Ethercat的控制架构，可实时采集位置控制下运动轨迹对应的各轴驱动电流，采样频率为1 kHz。根据驱动电流和驱动力矩的关系，可计算出各轴驱动力矩τ。同时在各轴驱动电机输出端安装有关节力矩传感器，用于测量除摩擦力矩之外的机械臂惯性力矩τ_dyn，用于辨识机械臂的惯性最小参数集^[16-17]。摩擦力矩可由F=τ-τ_dyn直接得出，用于各轴LuGre参数的辨识。

图 6 实验室自研机械臂

图选项

使用图 1的S型加减速轨迹作为训练轨迹进行辨识，并使用另一条随机的S型加减速轨迹作为验证轨迹，检验机械臂关节力矩预测的精度。
3.2.2 实验结果与分析作为示例，图 7和8给出了一些轴的摩擦力矩辨识结果。图 7表明，在Stribeck峰区(见图 7a、7c)和迟滞环线(见图 7b、7d)，迭代后结果都可更好地拟合实际曲线。图 8表明，形状因子辨识结果对应的摩擦力矩—速度曲线比与纯PSO法下的结果更接近于真实摩擦曲线。

图 7 形状因子法迭代前后对比

图选项

图 8 真实摩擦、形状因子辨识、PSO辨识对应的摩擦力矩—速度曲线对比

图选项

利用辨识出的各轴LuGre摩擦参数可预测出验证轨迹的摩擦力矩F_predict，再加上对应的τ_dyn，即得到预测的总关节力矩τ_predict。为了评价关节力矩预测精度，本文采用以下2个评价指标：
1) 预测与真实关节力矩的加权平均绝对误差百分比WAPE，计算方法为

$\mathrm{WAPE}=\frac{\sum\limits_{i=1}^{N}\left|T_{\text {predict }}(i)-T_{\text {real }}(i)\right|}{\sum\limits_{i=1}^{N}\left|T_{\text {real }}(i)\right|} \times 100 \% .$

其中：T_predict(i)为轨迹点i的预测力矩，T_real(i)为实验记录的轨迹点i的真实总力矩，N为轨迹总点数。
2) 最大误差百分比maxPE，计算方法为

$\begin{equation}\operatorname{maxPE}=\frac{\max\limits _{i=1: N}\left(\left|T_{\text {predict }}(i)-T_{\text {real }}(i)\right|\right)}{\sum\limits_{i=1}^{N}\left|T_{\text {real }}(i)\right| / N} \times 100 \%\end{equation}.$

WAPE和maxPE分别体现了摩擦影响下的平均辨识精度和最大可能误差。
表 3给出了验证轨迹在不同情况下的各轴关节力矩预测误差，将使用LuGre模型和形状因子法、使用式(10)的Coulomb摩擦模型、使用LuGre模型和纯PSO法辨识分别记为A、B、C。作为示例，图 9给出了验证轨迹在这3种情况下机械臂第1轴的总力矩预测残差。表 3和图 9所示结果表明，采用LuGre摩擦模型大大减小了速度变向处的尖峰误差，且形状因子法可以更好地辨识LuGre参数，有助于减小最大误差。
表 3 验证轨迹的力矩预测误差对比?

%
轴	WAPE				maxPE
	A	B	C		A	B	C
1	8.0	5.7	5.3		115.1	21.4	14.1
2	7.2	6.8	5.5		92.7	46.9	30.0
3	6.6	5.1	4.4		101.4	40.5	21.0
4	7.5	4.8	4.3		116.7	43.1	36.2
5	12.0	6.8	5.2		124.2	31.7	21.5
6	7.2	6.1	5.3		122.8	30.1	21.2

表选项

图 9 各种方法下的关节力矩预测残差

图选项

4 结论本文提出了一种基于形状因子的LuGre摩擦参数辨识方法。利用摩擦力矩—速度曲线上Stribeck峰和迟滞环线的特点，定义形状因子来定量表示Stribeck峰和迟滞环线的形状特征，通过对比当前辨识的摩擦力矩—速度曲线与真实曲线的形状因子，不断迭代优化PSO法的初步辨识结果，使得摩擦的预测结果更精确。得益于使用较大速度范围的辨识轨迹，辨识工作量较小，且不需要在辨识LuGre动态参数时严格将驱动力矩控制在静摩擦以下，该方法可以很好地应用于位置控制模式的运动系统，能同时辨识多个运动副的摩擦参数。在一台机械臂上进行了应用验证，将该方法应用于机械臂关节摩擦辨识，可显著提高关节速度在零附近的关节力矩预测精度。

参考文献

[1]	DING L, WU H T, YAO Y, et al. Dynamic model identification for 6-DOF industrial robots[J]. Journal of Robotics, 2015, 2015: 471478.
[2]	BOMPOS N A, ARTEMIADIS P K, OIKONOMOPOULOS A S, et al. Modeling, full identification and control of the Mitsubishi PA-10 robot arm[C]//Proceedings of the 2007 IEEE/ASME International Conference on Advanced Intelligent Mechatronics. Zurich, Switzerland: IEEE, 2007: 1-6.
[3]	WERNHOLT E, GUNNARSSON S. Nonlinear identification of a physically parameterized robot model[J]. IFAC Proceedings Volumes, 2006, 39(1): 143-148. DOI:10.3182/20060329-3-AU-2901.00016
[4]	HAN Y, WU J H, LIU C, et al. An iterative approach for accurate dynamic model identification of industrial robots[J]. IEEE Transactions on Robotics, 2020, 36(5): 1577-1594. DOI:10.1109/TRO.2020.2990368
[5]	WAHRBURG A, KLOSE S, CLEVER D, et al. Modeling speed-, load-, and position-dependent friction effects in strain wave gears[C]//2018 IEEE International Conference on Robotics and Automation (ICRA). Brisbane, QLD, Australia: IEEE, 2018: 2095-2102.
[6]	WOLF S, ISKANDAR M. Extending a dynamic friction model with nonlinear viscous and thermal dependency for a motor and harmonic drive gear[C]//2018 IEEE International Conference on Robotics and Automation (ICRA). Brisbane, QLD, Australia: IEEE, 2018: 783-790.
[7]	SIMONI L, BESCHI M, LEGNANI G, et al. On the inclusion of temperature in the friction model of industrial robots[J]. IFAC-PapersOnLine, 2017, 50(1): 3482-3487. DOI:10.1016/j.ifacol.2017.08.933
[8]	JOHANASTROM K, CANUDAS-DE-WIT C. Revisiting the LuGre friction model[J]. IEEE Control Systems, 2008, 28(6): 101-114. DOI:10.1109/MCS.2008.929425
[9]	DE WIT C C, OLSSON H, ASTROM K J, et al. A new model for control of systems with friction[J]. IEEE Transactions on Automatic Control, 1995, 40(3): 419-425. DOI:10.1109/9.376053
[10]	OLSSON H, ?STR?M K J, DE WIT C C, et al. Friction models and friction compensation[J]. European Journal of Control, 1998, 4(3): 176-195. DOI:10.1016/S0947-3580(98)70113-X
[11]	GANDHI P S, GHORBEL F H, DABNEY J. Modeling, identification, and compensation of friction in harmonic drives[C]//Proceedings of the 41st IEEE Conference on Decision and Control. Las Vegas, NV, USA: IEEE, 2002: 160-166.
[12]	IRAKOZE R, YAKOUB K, KADDOURI A. Identification of piezoelectric LuGre model based on particle swarm optimization and real-coded genetic algorithm[C]//2015 IEEE 28th Canadian Conference on Electrical and Computer Engineering (CCECE). Halifax, NS, Canada: IEEE, 2015: 1451-1457.
[13]	陈良波. 2自由度机器人转动关节LuGre摩擦模型参数辨识[D]. 泉州: 华侨大学, 2013. CHEN L B. Parameters identification of LuGre friction model of revolute joints in a 2-degree-of-freedom robot[D]. Quanzhou: Huaqiao University, 2013. (in Chinese)
[14]	DUAN H B, WANG D B, ZHU J Q, et al. Parameter identification of LuGre friction model for flight simulation servosystem based on ant colony algorithm[J]. Transactions of Nanjing University of Aeronautics & Astronau, 2004(3): 179-183.
[15]	李志杰, 蔡力钢, 刘志峰. 加加速度连续的S型加减速规划算法[J]. 计算机集成制造系统, 2019, 25(5): 1192-1201. LI Z J, CAI L G, LIU Z F. S type acceleration & deceleration fast planning algorithm with continuous jerk[J]. Computer Integrated Manufacturing Systems, 2019, 25(5): 1192-1201. (in Chinese)
[16]	SWEVERS J, VERDONCK W, DE SCHUTTER J. Dynamic model identification for industrial robots[J]. IEEE Control Systems, 2007, 27(5): 58-71. DOI:10.1109/MCS.2007.904659
[17]	MAGRINI E, FLACCO F, DE LUCA A. Estimation of contact forces using a virtual force sensor[C]//2014 IEEE/RSJ International Conference on Intelligent Robots and Systems. Chicago, IL, USA: IEEE, 2014: 2126-2133.