工业机器人在生产生活的各个领域有着越来越广泛的应用，日益提高的精度需求对机器人运动控制器的控制性能提出了更高的要求。采用传统的PID控制器控制机器人运动时，为保证较高的轨迹跟踪精度，通常采用较大的PID参数值，但这也放大了噪声与干扰信号，导致控制器控制性能退化、系统抗干扰能力下降、出现驱动器超载与振动等系统失稳现象。为改善机器人运动控制性能，一般考虑在运动控制中应用基于动力学模型的控制器。而基于模型的控制器需要以机器人动力学参数作为先验值，而获取参数最常见的方法即实验辨识。目前动力学参数辨识的难点主要在于：(1)预测精度不足，由于辨识算法对模型的限制，辨识方法所采用的模型精度不足，基于该模型的预测力矩精度也因之受到限制；(2)预测力矩波动，这是由于动力学模型在某些节点出现预测值的较大范围不连续跳动引起的；(3)误差峰值，由于传统摩擦模型对特殊运动状态下动力学特性的描述不精确，导致在这些运动状态附近与实际所需力矩值相比产生较大的偏差。

国内外学者提出了各种有效的辨识机器人参数的方法。Gautier等提出了串联机器人最小惯性参数集的推导方法^[1]。Swevers等提出了采用基于有限项傅里叶级数的激励轨迹进行辨识实验^[2]，以解析地获得各个关节的关节速度与关节加速度，这样做可以有效提高关节速度、加速度的计算精度。Park采用最小二乘法进行辨识^[3]，但是该方法只能辨识线性模型，所用的库仑-黏性摩擦模型对摩擦特性描述不精确，因而辨识精度不高。Grotjahn等提出了一种顺序辨识的方案^[4]，该方法中待辨识参数需要以已辨识参数作为先验值，因此存在较大累积误差。上述方法均基于线性动力学模型展开研究，对关节间摩擦因素一般近似处理为库仑-黏性摩擦模型。该模型对关节间低速动态摩擦特性描述较差，是辨识精度不足与预测力矩误差峰值产生的主要原因。为改善上述问题，应考虑采用低速动态特性较好的非线性摩擦模型描述关节间摩擦。经典的Stribeck、Dahl、LuGre模型在描述机器人关节间摩擦时效果并不理想^[5-6]。Deami等提出了一种针对关节机器人的基于实验数据拟合的摩擦模型^[7-8]，该摩擦模型对关节型机器人在低速情况下的复杂摩擦特性进行了拟合，因而能够更好地反映机器人实际的摩擦特性。

以力矩的预测精度为目标函数，动力学参数辨识问题就转换为典型的多维非线性函数优化问题。文献[9]中提出了最大似然估计法优化未知参数集，该算法容易使优化参数陷入局部最优。人工蜂群(Artificial bee colony, ABC)算法在处理多维非线性函数优化问题上具有显著优势^[10-11]，收敛速度快，搜索能力强，不容易陷入局部最优。Ding等^[12]利用ABC算法成功辨识出空中机器人的动力学模型，并证明了该算法比遗传算法辨识精度高。针对水下机器人动力学参数辨识，高明松等^[13]通过ABC算法优化出模型中9个未知参数。为此，本文采用ABC算法辨识6-DOF工业机器人的动力学参数。

在机器人控制器设计方面，杨钢等在机器人系统中应用了模糊PID控制器^[14-18]，改善了系统的动态性能。苏玉鑫等提出了非线性PID控制等控制算法^{[19, 21]}，提高了系统反应速度。但上述控制算法得出计算力矩无法反映机器人运动所需要的实际力矩值。对此，文献[22, 23]中提出了基于动力学模型的前馈控制器，在动力学模型精确已知的前提下，前馈控制器能够有效提升机器人轨迹跟踪精度，且实现难度较小。文献[24, 25]中提出了计算力矩控制算法，这一方案需要精确、完整的各项动力学参数作为先验值，而实际应用中这些参数很难获得，因此计算力矩控制方案并不实用。文献[26~28]中提出了基于在线辨识的结果自适应控制算法，但在实时性要求较高的工业机器人领域并不实用。

本文首先建立机器人的刚体动力学可辨识模型，并用描述关节型机器人关节间摩擦特性效果更好的Deami-Heimann模型描述关节摩擦，替代了传统库仑-黏性摩擦模型，得到完整的机器人动力学模型。设计了辨识实验与激励轨迹，并利用ABC算法、基于实验数据辨识机器人动力学参数。与传统最小二乘法(Least square，LS)得出的辨识结果进行对比，分别计算总残差均方根并比较，以验证本文所述方法的先进性。最后，设计基于模型的前馈控制器，通过比较传统PID控制器与前馈控制器的控制性能，实验验证控制器的优越性。

1 关节型机器人动力学建模

对n自由度关节型机器人，已知关节运动，采用Newton-Euler法计算驱动器力矩，得到机器人刚体动力学方程

$ {\mathit{\boldsymbol{\tau }}_r} = \mathit{\boldsymbol{M}}\left( \mathit{\boldsymbol{q}} \right)\mathit{\boldsymbol{\ddot q}} + \mathit{\boldsymbol{C}}\left( {\mathit{\boldsymbol{q}}, {\rm{ }}\mathit{\boldsymbol{\dot q}}} \right) + \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{q}} \right) $

(1)

式中：${\mathit{\boldsymbol{\tau }}_r} \in {\mathit{\boldsymbol{R}}^n}, \mathit{\boldsymbol{q}} \in {\mathit{\boldsymbol{R}}^n}, \mathit{\boldsymbol{\dot q}} \in {\mathit{\boldsymbol{R}}^n}, \mathit{\boldsymbol{\ddot q}} \in {\mathit{\boldsymbol{R}}^n}$分别为仅考虑刚体动力学时的关节力矩向量、关节转角向量、关节角速度向量及关节角加速度向量；M(q)∈R^n×n为质量矩阵；$\mathit{\boldsymbol{C}}\left( {\mathit{\boldsymbol{q}}, {\rm{ }}\mathit{\boldsymbol{\dot q}}} \right) \in {\mathit{\boldsymbol{R}}^n}$为科氏力与离心力项；G(q)∈Rⁿ为重力项。

将刚体动力学方程线性化并消去系数矩阵中的线性相关项，得出一组最小惯性参数集，进而得到可辨识模型

$ {\mathit{\boldsymbol{\tau }}_{{\rm{rgd}}}} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_r}\left( {\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}, \mathit{\boldsymbol{\ddot q}}} \right){\mathit{\boldsymbol{\theta }}^r} $

(2)

式中：θ^r∈R^m为最小惯性参数集，其中m为最小惯性参数集参数个数；Φ_r∈R^n×m为对应的系数矩阵。

为使动力学模型能够正确反映关节型机器人驱动器力矩与关节运动的关系，需要对关节间摩擦进行补偿。采用文献[7, 8]提出的Deami-Heimann模型，这一模型具有如下形式

$ {\mathit{\boldsymbol{\tau }}_f} = {\mathit{\boldsymbol{f}}_0} + {\mathit{\boldsymbol{f}}_1}{\rm{sign}}(\mathit{\boldsymbol{\dot q}}) + {\mathit{\boldsymbol{f}}_2}\mathit{\boldsymbol{\dot q}} + {\mathit{\boldsymbol{f}}_3}{\rm{arctan}}({\mathit{\boldsymbol{f}}_4}\mathit{\boldsymbol{\dot q}}) $

(3)

式中：τ_f∈Rⁿ为摩擦力矩向量；f₀∈Rⁿ为摩擦力矩零移；f₁∈Rⁿ为库仑摩擦系数；f₂∈Rⁿ为黏性摩擦系数；等式右边第四项是Deami等提出的基于实验数据的补偿项，其中f₃∈Rⁿ和f₄∈Rⁿ为基于实验数据拟合而得出的经验系数，$\mathit{\boldsymbol{\dot q}} \in {\mathit{\boldsymbol{R}}^n}$为角速度向量，该模型是未知参数的非线性函数。本文将采用该模型描述机器人关节间摩擦。由式(2，3)可知，机器人刚体动力学力矩与机器人关节摩擦力矩都是关节位置、关节速度、关节加速度与可辨识参数的函数。因此，在刚体动力学方程基础上添加摩擦力矩向量，得到如动力学方程

$ \begin{array}{l} {\mathit{\boldsymbol{\tau }}_{{\rm{dyn}}}} = {\mathit{\boldsymbol{\tau }}_r} + {\mathit{\boldsymbol{\tau }}_f} = {\rm{ }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_r}\left( {\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}, \mathit{\boldsymbol{\ddot q}}} \right){\mathit{\boldsymbol{\theta }}^r} + {\mathit{\boldsymbol{f}}_0} + \\ {\mathit{\boldsymbol{f}}_1}{\rm{sign}}(\mathit{\boldsymbol{\dot q}}) + {\mathit{\boldsymbol{f}}_2}\mathit{\boldsymbol{\dot q}} + {\mathit{\boldsymbol{f}}_3}{\rm{arctan}}({\mathit{\boldsymbol{f}}_4}\mathit{\boldsymbol{\dot q}}) = \\ \mathit{\boldsymbol{F}}\left( {\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}} \right), \mathit{\boldsymbol{\ddot q}}, {\mathit{\boldsymbol{\theta }}^{{\rm{dyn}}}} \end{array} $

(4)

式中：τ_dyn∈Rⁿ为考虑关节间摩擦时的关节力矩向量；θ^dyn∈R^m+5n为待辨识参数集；$\mathit{\boldsymbol{F}}\left( {\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}, \mathit{\boldsymbol{\ddot q, }}{\mathit{\boldsymbol{\theta }}^{{\rm{dyn}}}}} \right) \in {\mathit{\boldsymbol{R}}^n}$，其中每个分量都是关节力矩以动力学参数、关节运动为自变量的非线性函数，由刚体动力学力矩分量与关节摩擦力矩分量通过简单叠加得到。要使方程有唯一最优解，应当令未知量的个数小于方程个数。因此，在N个时间节点采集N组关节转角、关节速度、关节加速度及关节力矩信号(nN>m+5n)，均为(n×1)的列向量。组合这N组信号形成如下非线性可辨识方程

$ \begin{array}{l} {\mathit{\boldsymbol{\tau }}_N} = \left[\begin{array}{l} {\tau _{{t_1}}}\\ {\tau _{{t_2}}}\\ \vdots \\ {\tau _{{t_N}}} \end{array} \right] = \left[\begin{array}{l} \mathit{\boldsymbol{F}}\left( {{\mathit{\boldsymbol{q}}_{{t_1}}}, {{\mathit{\boldsymbol{\dot q}}}_{{t_1}}}, {{\mathit{\boldsymbol{\ddot q}}}_{{t_1}}}\mathit{\boldsymbol{, }}{\mathit{\boldsymbol{\theta }}^{{\rm{dyn}}}}} \right)\\ \mathit{\boldsymbol{F}}{\left( {{\mathit{\boldsymbol{q}}_{{t_2}}}, \mathit{\boldsymbol{\dot q}}} \right)_{{t_2}}}, \left( {{{\mathit{\boldsymbol{\ddot q}}}_{{t_2}}}\mathit{\boldsymbol{, }}{\mathit{\boldsymbol{\theta }}^{{\rm{dyn}}}}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\ \mathit{\boldsymbol{F}}{\left( {{\mathit{\boldsymbol{q}}_{{t_N}}}, \mathit{\boldsymbol{\dot q}}} \right)_{{t_N}}}, \left( {{{\mathit{\boldsymbol{\ddot q}}}_{{t_N}}}\mathit{\boldsymbol{, }}{\mathit{\boldsymbol{\theta }}^{{\rm{dyn}}}}} \right) \end{array} \right] = \\ {\mathit{\boldsymbol{F}}_N}\left( {\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}, \mathit{\boldsymbol{\ddot q, }}{\mathit{\boldsymbol{\theta }}^{{\rm{dyn}}}}} \right) \end{array} $

(5)

式中：τ_t1, τ_t2, …, τ_tN∈Rⁿ为各时间节点由伺服电机电流间接得出的关节力矩；τ_N∈R^Nn×1为N个时间点关节力矩向量的组合；$\mathit{\boldsymbol{F}}\left( {\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}, {{\mathit{\boldsymbol{\ddot q}}}_{{t_i}}}\mathit{\boldsymbol{, }}{\mathit{\boldsymbol{\theta }}^{{\rm{dyn}}}}} \right) \in {\mathit{\boldsymbol{R}}^n}$, i=1, 2, …, N为t_i时刻关节力矩以动力学参数、关节运动为因变量的非线性函数，由式(4)获得；$\mathit{\boldsymbol{q}}, \mathit{\boldsymbol{\dot q}}, \mathit{\boldsymbol{\ddot q}} \in {\mathit{\boldsymbol{R}}^n}$分别为关节转角、速度、加速度向量；θ^dyn∈R^m+5n为动力学参数向量，通过求解式(5)就可以得到θ^dyn的辨识值。

2 辨识实验 2.1 实验设计

如图 1所示是动力学参数辨识实验所用的某型号6自由度关节型机器人。本文仅辨识该关节型机器人前三关节的动力学参数。基于文献[29]的结论，在辨识前三关节的动力学参数时可以忽略后三关节的影响，如此辨识更具针对性。实验使用基频为0.04 Hz的形如式(6)所示的五项傅里叶级数的激励轨迹，如图 2所示。图 1, 2令机器人连续跟踪20个周期的激励轨迹，分别利用伺服电机编码器和电机电流间接获取关节角位移与力矩信息

$ \begin{array}{l} {{\mathit{\boldsymbol{\dot q}}}_i}\left( t \right) = {\mathit{\boldsymbol{q}}_{i, 0}} + \sum\limits_{k = 1}^5 {\left( {{\mathit{\boldsymbol{a}}_{i, k}}{\rm{sin}}\left( {k{\omega _f}t} \right) + {\mathit{\boldsymbol{b}}_{i, k}}{\rm{cos}}\left( {k{\omega _f}t} \right)} \right)} \\ i = 1, 2, 3 \end{array} $

(6)

式中:q_i(t)∈Rⁿ为关节i在关节空间中的轨迹关于时间t的函数；q_{i, 0}∈Rⁿ为各关节在关节空间中t=0时刻的位置；a_{i, k}∈Rⁿ和b_{i, k}∈Rⁿ为关节i的k阶傅里叶级数因子；ω_f=0.04为离散傅里叶级数的基频。

2.2 数据预处理

为降低噪声信号对辨识精度的影响，需要对测量信号作预处理。对采集得到的20个周期的信号进行平均化滤波处理，得到关节角位移与关节力矩在一个周期内的采样平均值^[30]。为进一步降低测量噪声的影响，将关节角位移信号重新拟合为有限项傅里叶级数，并采用式(7，8)解析地计算关节速度与关节加速度

$ \begin{array}{l} {{\mathit{\boldsymbol{\dot q}}}_i}\left( t \right) = \sum\limits_{k = 1}^N {\left( {{\mathit{\boldsymbol{a}}_{i, k}}k{\omega _f}{\rm{cos}}(k{\omega _f}t) - {\mathit{\boldsymbol{b}}_{i, k}}k{\omega _f}{\rm{sin}}(k{\omega _f}t)} \right)} \\ i = 1, 2, \ldots, n \end{array} $

(7)

$ \begin{array}{l} {{\mathit{\boldsymbol{\ddot q}}}_i}\left( t \right) = \\ - {\rm{ }}\sum\limits_{k{\rm{ }} = {\rm{ }}1}^N {\left( {{\mathit{\boldsymbol{a}}_{i, k}}{k^2}\omega _f^2{\rm{sin}}(k{\omega _{f}}t) + {\mathit{\boldsymbol{b}}_{i, k}}{k^2}\omega _f^2{\rm{cos}}(k{\omega _{f}}t)} \right)} \\ i{\rm{ }} = {\rm{ }}1, 2, \ldots, n \end{array} $

(8)

式中：${{\mathit{\boldsymbol{\dot q}}}_i}\left( t \right) \in {\mathit{\boldsymbol{R}}^n}$和${{\mathit{\boldsymbol{\ddot q}}}_i}\left( t \right) \in {\mathit{\boldsymbol{R}}^n}$为关节i在关节空间中的速度、加速度关于时间t的函数。关节力矩来源于电机电流，其测量噪声较大，简单的平均化滤波过滤高频分量的效果较差，而低通滤波又将引入相位失真。为此，本文采用五点三次平滑法处理关节力矩信号。取平滑次数为100，对关节力矩采样信号进行平滑，三个关节的力矩数据处理前后对比结果分别如图 3~5所示。

由图 3~5可以看出，处理后信号相当完整地保留了处理前信号的特征信息，信号平滑，且未出现相位失真现象。预处理对信号中高频分量的抑制效果显著。

2.3 ABC算法

ABC算法是对蜂群觅食行为的模拟，将机器人待辨识的动力学参数作为食物源，根据优胜劣汰原则，通过采集输入的激励轨迹与输出的关节力矩、角位移等数据作为实验样本，进行优化计算得到最优参数。算法的具体流程如图 6所示^[11]。

图 6中，NP表示蜂群中引领蜂与跟随蜂的总数，limit为侦察蜂搜索的阈值，G为算法执行的最大循环次数。这些参数的取值具有一定的经验性，需要基于特定对象进行调整。具体算法执行流程参考文献[9, 10]。简而言之，辨识以预测力矩与测量力矩的总残差均方根为目标函数，经种群初始化获得初始种群并进行迭代，每次迭代都进行引领蜂搜索、跟随蜂搜索，如有必要，还要进行侦查蜂搜索，每次迭代都记录本次迭代中适应度值。迭代结束时, 当前迭代循环中的最优个体就是算法的全局最优解。

2.4 实验结果与分析

在Matlab环境下，采用ABC算法辨识该工业机器人的动力学参数。令NP=20, limit=15, G=200, 以预测力矩与测量力矩的总残差均方根为目标函数，经过约60次迭代，算法逐步向全局最优收敛，总残差均方根最终收敛于全局最优值28.112 9 N·m。算法主要指标如表 1所示，总残差均方根迭代曲线如图 7所示。

由图 7可以看出，采用人工蜂群算法辨识基于Deami-Heimann摩擦模型的关节型机器人动力学参数时，算法能够以较快的速度定位全局最优所在的区间，且不易陷入局部最优。

为了验证基于本文提出的算法计算得出的动力学参数能够更加有效地描述机器人动力学特性，采用传统最小二乘法结合库仑-黏性摩擦模型，利用相同的测量数据辨识动力学参数，其力矩预测总残差均方根为31.248 0 N·m。而基于本文提出的辨识算法，力矩预测总参差均方根降低至28.112 9 N·m。可见，基于本文所述方法，力矩预测精度提高了约10.03%。因此，本文所述方法相比传统最小二乘法，其辨识精度得到了提升。

为了更形象地展示该辨识摩擦模型的优越性，将两种方法得到的力矩预测结果绘制在同一张图中。以关节2为例，两种摩擦模型的力矩预测对比如图 8所示。观察局部放大图，可以发现基于本文所述方法得出的预测力矩相比传统最小二乘法，其拟合程度较高，且关节预测力矩不存在跳变等现象。这同时也是基于库仑-黏性摩擦模型进行动力学参数辨识的一大缺陷。当关节速度的符号发生变化时，基于简单的库仑-黏性摩擦模型和最小二乘法得出的预测力矩将发生跳变，这对控制器的控制性能有较大损害。由图 8可知，使用基于Deami-Heimann摩擦模型和人工蜂群算法辨识得出的关节型机器人动力学模型进行力矩预测时，力矩预测误差较小，能够更好地拟合测量数据。并且，相比库仑-黏性摩擦模型，基于Deami-Heimann模型的参数辨识还有效地抑制了误差峰值。

3 基于模型的前馈控制

参数辨识的目的是为基于模型的控制器提供参数。动力学前馈控制器在关节机器人控制中已有广泛应用，图 9所示是前馈控制器的控制策略。

该方案就是在原有的PID环节之前增加了一个力矩前馈环节，根据辨识得出的动力学模型与期望轨迹，可以计算得出各机器周期驱动器的输出力矩预测值，并存储于控制器ROM中。在机器人开始运动后，每个机器周期驱动器输出当前力矩预测值，与误差信号一起构成控制信号。如果期望轨迹是先验已知的(这在实际生产活动中是很常见的)，那么力矩预测值将可以离线计算并存储在控制器的存储器中，而不需要实时计算。本文利用辨识得到的动力学参数设计了基于模型的控制器，并应用于机器人的位置控制，通过实验验证基于辨识结果设计的前馈控制器能够有效改善机器人的动力学特性与轨迹跟踪精度。

首先，为3个关节各指定一个周期的激励轨迹，这一个周期的激励轨迹应当在各个关节的关节空间范围内，且关节速度、加速度存在跳变，以验证控制器的性能；使用传统的PID控制器(PID参数在使得驱动器不出现超载现象的前提下经过了优化)，令机器人跟踪如图 10所示的激励轨迹，同时从伺服电机编码器获得各关节的实际位置。持续跟踪50个周期，对得到的实际位置取平均值，结合图 10所示的期望位置，得到跟踪误差如图 11所示；同样地，基于辨识实验得出的机器人动力学参数与期望轨迹计算得出机器人跟踪既定期望轨迹的预测力矩，使用图 9所示的控制器，令研究对象跟踪图 10所示的期望轨迹，同时从伺服电机编码器获得各关节的实际位置，结合图 10所示的期望位置，得到跟踪误差如图 12所示；最后将采用传统PID控制器得出的跟踪误差曲线与使用了前馈控制器时的跟踪误差曲线进行了对比。

对比图 11及图 12可以看出，相比未添加力矩前馈时的情形，添加了力矩前馈环节的运动控制其跟踪精度得到了显著提升。以关节1为例，未添加力矩前馈时其最大跟踪误差绝对值约为0.05 rad，而添加了力矩前馈环节后的最大跟踪误差绝对值减小到约0.01 rad。这证明了：(1)基于机器人动力学模型的前馈控制器能够有效改善机器人轨迹跟踪精度；(2)本文提出的动力学参数辨识算法是有效的，基于辨识结果的模型正确地反映了各关节执行期望轨迹时所需的力矩，基于辨识结果的前馈控制器改善了机器人动力学特性与轨迹跟踪精度。

4 结束语

传统关节机器人辨识通常基于线性摩擦模型，本文采用精度更高、低速动态性能更好的摩擦模型描述关节摩擦，建立动力学模型，并辨识相关参数。辨识过程中，采用有限项傅里叶级数轨迹作为激励轨迹，通过将关节位置输出信号重新拟合为傅里叶级数，有效降低了测量噪声对关节位置数据的影响。提出采用五点三次滤波处理关节力矩信号，噪声抑制效果好。由于最小二乘法只能辨识线性模型，因此为辨识非线性模型中的动力学参数，本文提出采用人工蜂群算法辨识关节型机器人的动力学参数集，并为之设计了实验与激励轨迹。基于实验结果计算得出的预测力矩较之以往的经典模型，其精度得到了提高，误差峰值得以有效抑制，这将显著改善机器人在关节换向附近的动力学特性，特别是抑制大振幅振动现象。基于辨识结果设计基于模型的前馈控制器，并进行运动控制实验，实验结果表明基于模型的控制器相比传统的PID控制显著提高了轨迹跟踪精度，避免了驱动器过载等现象，同时关节换向附近的大振幅振动得到了有效抑制。