柔性关节机械臂的轨迹跟踪控制问题在过去几年中得到了高度重视^[1-9]。工业机器人的谐波传动系统都不可避免的存在柔性环节，当系统需要快速响应时，就会造成明显的扭转震颤。同时，弹性环节存在于机械臂系统的减速箱、传送带、连接腱等诸多构件当中，适用于刚性机械臂的成熟控制算法应用于含有弹性环节的动态系统中时，其响应速度和动态精度都会受到很大限制。因此，机械臂关节系统中弹性环节的动态特性在机械臂的建模与控制问题中都是不可或缺的^[1-6]。柔性关节的引入增加了机械臂动力学系统的复杂性，系统阶次变成了相同连杆数量机械臂系统的两倍，系统输出自由度的数量超过了控制输入的数量，给控制系统的设计带来了新的挑战^[7-10]。

鉴于此，国内外学者提出了各种控制策略来解决柔性关节为机械臂控制带来的问题。例如反馈线性化^{[3, 11-12]}、基于无源性理论的控制方法^[13-14]、奇异扰动控制策略^[15]、反步法^[15-17]、自适应控制^[18-20]等。以上文献中提出的方法存在如下几个方面的问题。首先，系统对于广义干扰的鲁棒性是以干扰上界已知为前提，控制策略具有较高的保守性。其次，控制策略要求柔性关节机械臂系统的全部状态已知。众所周知，工业机械臂传动系统机构设计非常紧凑，在狭小的空间内为伺服电机传动轴和关节末端同时安装角度传感器和速度传感器，不仅提高了设计成本，而且增加了工程设计难度^{[1-2, 6, 10]}。同时，一些速度传感器由于本身设计机理的原因，其抗噪性不能满足复杂工业现场的要求^{[11, 20-21]}。因此，从经济适用性和物理约束两个方面综合考虑，无速度传感器的柔性关节机械臂控制系统具有更高的研究价值。

近年来，一些学者研究了基于观测器的输出跟踪控制问题^{[3, 12, 22-25]}。文献[17]针对一类非线性动态系统，结合反步法与小增益定理，提出了一种自适应鲁棒模糊输出跟踪控制策略；Lightcap等^[26]设计了扩展卡尔曼滤波器对柔性关节机械臂进行状态估计；文献[27]设计了柔性关节系统的反步控制器，利用自适应神经网络对柔性关节机械臂进行状态估计；文献[28]针对含有参数不确定性和外界干扰的柔性关节机械臂系统设计了滑模鲁棒观测器，并且可以扩展到一般形式的系统状态方程中去。上述文献所提出的方法虽解决了部分状态可测条件下的柔性关节系统的输出跟踪问题，但在控制器设计方面大多仅设计了基于反馈回路的鲁棒控制策略来解决模型不确定性抑或外界干扰对柔性关节机械臂轨迹跟踪问题带来的影响。

此外，对于机械臂系统的跟踪控制而言，伺服电机的减速器以及机械臂的外部传动系统等部件会不可避免地将非线性摩擦环节引入控制系统当中^[29-32]，为柔性关节机械臂的高精度控制带来更大的挑战。在控制器策略的设计中，如果不能精确地、恰当地对系统的非线性摩擦进行补偿，那么在系统控制中，尤其对机械臂的低速控制问题，将不可避免的出现较大的跟踪误差和粘滑运动^[29-30]。因此，为了对潜在的非线性摩擦进行补偿，学者们建立了多种摩擦环节的数学模型^[31-32]。其中，LuGre模型由于其结构简易同时能对多种摩擦特性进行模拟而被广泛使用。但是，LuGre摩擦模型在工程实践当中也存在着诸多不足。首先，在控制器设计当中使用LuGre模型，要单独设计估计器对鬃毛形变状态变量进行估计。当同时考虑模型不确定性时，控制器的一体化设计目的很难达到。其次，在高速情况下，LuGre模型急剧增加的摩擦刚性会为控制器的数字实现产生不利影响，因而不得不对控制器进行进一步改进。再次，基于LuGre模型的摩擦补偿器设计会不可避免地为系统引入符号函数, 导致控制输出存在高频抖振，降低了控制器对高频噪声的鲁棒性。实验表明，优良的静态摩擦模型能够对实际系统存在的非线性摩擦进行高于90%的补偿^[29-30]。同时文献[31]表明，与动态摩擦模型相比，静态摩擦模型可以对包含极限环在内的摩擦现象进行几乎同等程度的估计。综上考虑，为了便于工程实现，基于静态摩擦模型的补偿器在柔性关节机械臂的跟踪控制中具有更强的现实意义^[33]。

综上所述，本文针对柔性关节机械臂的输出跟踪问题，提出了基于迭代学习观测器(iterative learning observer，ILO)的滤波反步控制策略，解决柔性关节系统状态估计问题的同时兼顾了控制系统的干扰抑制性能。针对柔性关节的动态特性，分析了系统不可测状态和模型不确定性对跟踪控制精度带来的影响。ILO作为系统内环，对系统的内部状态和由模型不确定性引起的复合扰动进行在线估计；利用指令滤波器解决了传统反步法控制器中存在的微分膨胀问题；同时，基于一种新型连续静态摩擦模型，设计了非线性摩擦补偿器，进一步提高了柔性关节机械臂输出跟踪的精度。通过Lyapunov稳定性理论分析了闭环控制系统的全局稳定性。

1 柔性关节机械臂控制系统设计 1.1 柔性关节机械臂模型建立

本文将单臂柔性关节机械臂作为研究对象，转动副由直流电机直接驱动，将关节视为一个具有线性弹性系数的扭转弹簧系统，如图 1所示。

单臂柔性机械臂的动力学方程可表示为^{[9, 18-19]}

$ \left\{ \begin{array}{l} I{{\ddot q}_1} + MgL\sin \left( {{q_1}} \right) + K\left( {{q_1} - {q_2}} \right) = 0\\ J{{\ddot q}_2} - K\left( {{q_1} - {q_2}} \right) = {\tau _u} + {\tau _d} + f\left( {{{\dot q}_2}} \right) \end{array} \right. $

(1)

式中：q₁为机械臂转角，q₂为电机转角，I为转动臂的转动惯量，J为电机的惯性矩，K为关节的弹性刚度系数，τ_u是输入转矩，τ_d是外加干扰，$ f\left( {{{\dot q}_2}} \right)$为系统摩擦，M和L为转动臂的质量和长度，g为重力加速度。

根据文献[30-33]，本文建立如下连续可微的静态非线性摩擦模型，为后续的补偿器设计打下基础：

$ \begin{array}{l} f\left( {{{\dot q}_2}} \right) = - {a_1}\tanh \left( {{c_1}{{\dot q}_2}} \right) - \\ {a_2}\left[ {\tanh \left( {{c_2}{{\dot q}_2}} \right) - \tanh \left( {{c_3}{{\dot q}_2}} \right)} \right] \end{array} $

(2)

式中：a₁与a₂代表不同等级的摩擦系数，c₁、c₂和c₃为拟合不同摩擦效应所需的成型参数。

为了便于后续的控制系统设计，定义如下状态变量：x₁=q₁, $ {x_2} = {{\dot q}_1} $, x₃=q₂, $ {x_4} = {{\dot q}_2} $，将式(1)改写为

$ \left\{ \begin{array}{l} {{\dot x}_1} = {x_2}\\ {{\dot x}_2} = {a_{21}}{x_1} + {a_{22}}\sin {x_1} + {a_{23}}{x_3}\\ {{\dot x}_3} = {x_4}\\ {{\dot x}_4} = {a_{41}}{x_1} + {a_{42}}{x_3} + {b_{41}}\left( {{\tau _u} + {\tau _d} + f\left( {{x_4}} \right)} \right) \end{array} \right. $

(3)

式中：a₂₁=-K/I，a₂₂=-MgL/I，a₂₃=K/I，a₄₁=K/J，a₄₂=-K/J，b₄₁=1/J。

考虑系统模型存在的参数不确定性，则有

$ \left\{ \begin{array}{l} {{\dot x}_1} = {x_2}\\ {{\dot x}_2} = \left( {{a_{21}} + \Delta {a_{21}}} \right){x_1} + \\ \;\;\;\left( {{a_{22}} + \Delta {a_{22}}} \right)\sin {x_1} + \left( {{a_{23}} + \Delta {a_{23}}} \right){x_3}\\ {{\dot x}_3} = {x_4}\\ {{\dot x}_4} = \left( {{a_{41}} + \Delta {a_{41}}} \right){x_1} + \left( {{a_{42}} + \Delta {a_{42}}} \right){x_3} + \\ \left( {{b_{41}} + \Delta {b_{41}}} \right){\tau _u} + {b_{41}}{\tau _d} + {b_{41}}f\left( {{x_4}} \right) \end{array} \right. $

(4)

式中Δ(·)表示参数的摄动值。

将(4)整理成一般形式，则有：

$ \left\{ \begin{array}{l} {{\dot x}_1} = {x_2}\\ {{\dot x}_2} = {a_{21}}{x_1} + {a_{22}}\sin {x_1} + {a_{23}}{x_3} + {d_1}\\ {{\dot x}_3} = {x_4}\\ {{\dot x}_4} = {a_{41}}{x_1} + {a_{42}}{x_3} + {b_{41}}\left( {{\tau _u}} \right) + {F_f}\left( {{x_4}} \right) + {d_2} \end{array} \right. $

(5)

式中：F_f(x₄)=b₄₁f(x₄)；d₁、d₂为由模型不确定性和外加干扰造成的复合干扰，假设复合扰动为有界变量，定义为

$ {d_1} = \Delta {a_{21}}{x_1} + \Delta {a_{22}}\sin {x_1} + \Delta {a_{23}}{x_3} $

$ {d_2} = \Delta {a_{41}}{x_1} + \Delta {a_{42}}{x_3} + \Delta {b_{41}}{\tau _u} + {b_{41}}{\tau _d} $

1.2 系统控制目标

对于跟踪控制任务而言，输出反馈控制结构通常无法获得理想的跟踪性能，且非线性系统存在的模型不确定性及外界扰动会对控制性能产生较大影响。因此，本文将整个控制系统分为内环观测器和外环控制器分别进行设计。对于内环系统，考虑到系统模型不确定性和外界扰动的影响，设计一个ILO对未知状态和系统不确定性进行在线估计，并在控制系统中对复合扰动进行补偿，将系统转化为标称形式。进而利用系统状态估计值将外环控制器设计为状态反馈形式，从而获得理想的跟踪性能。控制系统结构如图 2所示。

1.3 迭代学习观测器设计

本节首先根据非线性系统(5)设计ILO，实现对于系统内部状态和复合扰动的在线估计。

由于系统状态x₁可测，可以将式(5)所示的非线性系统转化为如下标准形式:

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{A}}x + \mathit{\boldsymbol{\xi }}\left( \mathit{\boldsymbol{x}} \right) + \mathit{\boldsymbol{Bu}} + {\mathit{\boldsymbol{d}}_t}\\ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{C}}x \end{array} \right. $

(6)

其中：$ \mathit{\boldsymbol{A = }}\left[{\begin{array}{*{20}{c}} 0&1&0&0\\ {{a_{21}}}&0&{{a_{23}}}&0\\ 0&0&0&1\\ {{a_{41}}}&0&{{a_{42}}}&0 \end{array}} \right] $，$ \mathit{\boldsymbol{\xi }}\left( \mathit{\boldsymbol{x}} \right) = \left[{\begin{array}{*{20}{c}} 0\\ {{a_{22}}\sin {x_1}}\\ 0\\ 0 \end{array}} \right] $，$\mathit{\boldsymbol{B}}\left[{\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ {{b_{41}}} \end{array}} \right], \mathit{\boldsymbol{C = }}{\left[\begin{array}{l} 1\\ 0\\ 0\\ 0 \end{array} \right]^{\rm{T}}}, {\mathit{\boldsymbol{d}}_t}\left[{\begin{array}{*{20}{c}} 0\\ {{d_1}}\\ 0\\ {{F_f}\left( {{x_4}} \right) + {d_2}} \end{array}} \right] $。

此时，根据式(6)所述系统模型，可以设计如下ILO：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\hat {\dot x}}} = \mathit{\boldsymbol{A}}\hat x + \mathit{\boldsymbol{\xi }}\left( {\mathit{\boldsymbol{\hat x}}} \right) + \mathit{\boldsymbol{Bu}} + \mathit{\boldsymbol{L}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\hat y}}} \right) + \mathit{\boldsymbol{\gamma }}\left( t \right)\\ \mathit{\boldsymbol{\gamma }}\left( t \right) = {\mathit{\boldsymbol{\mu }}_1}\mathit{\boldsymbol{\gamma }}\left( {t - \tau } \right) + {\mathit{\boldsymbol{\mu }}_2}\left[ {\mathit{\boldsymbol{y}}\left( t \right) - \mathit{\boldsymbol{\hat y}}\left( t \right)} \right]\\ \mathit{\boldsymbol{\hat y}} = \mathit{\boldsymbol{C\hat x}} \end{array} \right. $

(7)

式中：$ \mathit{\boldsymbol{\hat x}}$为待估计的系统内部状态，τ是采样周期，y(t-τ)为上一个采样周期的观测输出，L∈R⁴，μ₁ ∈R^4×4, μ₂∈R⁴为待确定的观测器增益矩阵。

由(6)、(7)可得ILO观测误差动态方程：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\dot {\tilde x} = }}\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{LC}}} \right)\mathit{\boldsymbol{\tilde x}} + \left[ {\mathit{\boldsymbol{\xi }}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{\xi }}\left( {\mathit{\boldsymbol{\hat x}}} \right)} \right] + {\mathit{\boldsymbol{d}}_t} - \mathit{\boldsymbol{\gamma }}\left( t \right)\\ \mathit{\boldsymbol{\tilde y}} = \mathit{\boldsymbol{C\tilde x}} \end{array} \right. $

(8)

式中$\mathit{\boldsymbol{\tilde x}} \buildrel \Delta \over = \mathit{\boldsymbol{x-\hat x}} $为观测器的估计误差。

设$ \mathit{\boldsymbol{\tilde d}}\left( t \right) = \mathit{\boldsymbol{d}}\left( t \right)-\mathit{\boldsymbol{\gamma }}\left( t \right) $, $ {\mathit{\boldsymbol{d}}_\tau } = \mathit{\boldsymbol{d}}\left( t \right)-{\mathit{\boldsymbol{\mu }}_1}\mathit{\boldsymbol{d}}\left( {t-\tau } \right) $。则有

$ \mathit{\boldsymbol{\tilde d}}\left( t \right) = {\mathit{\boldsymbol{\mu }}_1}\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) - {\mathit{\boldsymbol{\mu }}_2}\mathit{\boldsymbol{C\tilde x}}\left( t \right) + {\mathit{\boldsymbol{d}}_\tau } $

(9)

此时，观测误差动态方程可改写为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot {\tilde x} = }}\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{LC}}} \right)\mathit{\boldsymbol{\tilde x + }}{\mathit{\boldsymbol{\mu }}_1}\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) - }\\ {{\mathit{\boldsymbol{\mu }}_2}\mathit{\boldsymbol{C\tilde x}}\left( t \right){\mathit{\boldsymbol{d}}_\tau } + \left[ {\mathit{\boldsymbol{\xi }}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{\xi }}\left( {\mathit{\boldsymbol{\hat x}}} \right)} \right]} \end{array} $

(10)

定理1  考虑式(4)所描述的柔性关节机械臂动力学系统和式(7)所述ILO，若以下条件成立，则最终观测误差是有界收敛的。

$ {\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{LC}}} \right)^{\rm{T}}}\mathit{\boldsymbol{P}} + \mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{LC}}} \right) + \mathit{\boldsymbol{R}} + {\mathit{\boldsymbol{P}}^{\rm{T}}}\mathit{\boldsymbol{P}} = - \mathit{\boldsymbol{Q}} $

(11)

$ \mathit{\boldsymbol{P}} = \mathit{\boldsymbol{\rho }}{\left( {{\mathit{\boldsymbol{\mu }}_2}\mathit{\boldsymbol{C}}} \right)^{\rm{T}}} $

(12)

$ 0 < \left( {\rho + {m^2}} \right)\mathit{\boldsymbol{\mu }}_1^{\rm{T}}{\mathit{\boldsymbol{\mu }}_1} \le \mathit{\boldsymbol{I}} $

(13)

证明  定义如下Lyapunov函数，其中P、R为正定对称阵。

$ {\mathit{\boldsymbol{V}}_1} = {{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\mathit{\boldsymbol{P\tilde x}} + \int_{t - \to \tau }^t {{{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( \alpha \right)\mathit{\boldsymbol{\tilde d}}\left( \alpha \right){\rm{d}}\alpha } $

(14)

将V对时间求导，并结合ILO的观测误差动态方程(10)可得

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\dot V}}}_1} = {{\mathit{\boldsymbol{\dot {\tilde x}}}}^{\rm{T}}}\mathit{\boldsymbol{P}}{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P\dot {\tilde x}}}\left( t \right) + 2{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{\mu }}_1}\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) - }\\ {2\mathit{\boldsymbol{\tilde x}}\left( t \right)\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{\mu }}_2}\mathit{\boldsymbol{C}}{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{d}}_\tau }\left( t \right) - \varepsilon {{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{\tilde d}}\left( t \right) + }\\ {\rho {{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{\tilde d}}\left( t \right) - {{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( {t - \tau } \right)\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) + }\\ {2{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P}}\left[ {\mathit{\boldsymbol{\xi }}\left( \mathit{{x}} \right) - \mathit{\boldsymbol{\xi }}\left( {\mathit{{\hat x}}} \right)} \right]} \end{array} $

(15)

式中ρ-ε=1。

结合式(9)，可进一步得出

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\dot V}}}_1} = {{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left[ {{{\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{LC}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{P}} + \mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{LC}}} \right)} \right]\mathit{\boldsymbol{\tilde x + }}}\\ {2{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{\mu }}_1}\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) + 2{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left( t \right)\left[ {\mathit{\boldsymbol{\xi }}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{\xi }}\left( {\mathit{\boldsymbol{\hat x}}} \right)} \right] - }\\ {2\mathit{\boldsymbol{\tilde x}}\left( t \right)\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{\mu }}_2}\mathit{\boldsymbol{C\tilde x}}\left( t \right) + 2{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{d}}_\tau }\left( t \right) - \varepsilon {{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{\tilde d}}\left( t \right) + }\\ {\rho \mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right)\mu _1^{\rm{T}}{\mu _1}{{\tilde d}_\tau }\left( {t - \tau } \right) + \rho \mathit{\boldsymbol{d}}_\tau ^{\rm{T}}{\mathit{\boldsymbol{d}}_\tau } + }\\ {2\mathit{\boldsymbol{p\tilde d}}\left( {t - \tau } \right)\mathit{\boldsymbol{\mu }}_1^{\rm{T}}{\mathit{\boldsymbol{d}}_\tau }\left( t \right) - 2\rho {{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}}{{\left( {{\mathit{\boldsymbol{\mu }}_2}\mathit{\boldsymbol{C}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{d}}_\tau }\left( t \right) - }\\ {{{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( {t - \tau } \right)\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right)} \end{array} $

(16)

由于d(t)有界，则辅助变量d_τ也为有界变量，即‖d_τ‖≤d_τ。

由式(11)、(12)有

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\dot V}}}_1} \le - {\lambda _{\min }}\left( \mathit{\boldsymbol{Q}} \right){{\left\| {\mathit{\boldsymbol{\tilde x}}} \right\|}^2} - \varepsilon {{\left\| {\mathit{\boldsymbol{\tilde d}}\left( t \right)} \right\|}^2} - }\\ {{{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( {t - \tau } \right)\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) + - \rho \mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right)\mathit{\boldsymbol{\mu }}_1^{\rm{T}}{\mathit{\boldsymbol{\mu }}_1}\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) + }\\ {2\rho \mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right)\mathit{\boldsymbol{\mu }}_1^{\rm{T}}{\mathit{\boldsymbol{d}}_\tau }\left( t \right) + \rho \mathit{\boldsymbol{\tilde d}}_\tau ^2 + 2\eta \left\| \mathit{\boldsymbol{P}} \right\|{{\left\| {\mathit{\boldsymbol{\tilde x}}} \right\|}^2}} \end{array} $

(17)

由三角不等式可得

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot V}} \le - \left( {{\lambda _{\min }}\left( \mathit{\boldsymbol{Q}} \right) - 2\eta {\lambda _{\max }}\left( \mathit{\boldsymbol{P}} \right)} \right){{\left\| {\mathit{\boldsymbol{\tilde x}}} \right\|}^2} - }\\ {\varepsilon {{\left\| {\mathit{\boldsymbol{\tilde d}}\left( t \right)} \right\|}^2} + {{\mathit{\boldsymbol{\tilde d}}}^{\rm{T}}}\left( {t - \tau } \right)\left[ {{{\left( {\rho + m} \right)}^2}\mathit{\boldsymbol{\mu }}_1^{\rm{T}}{\mathit{\boldsymbol{\mu }}_1} - } \right.}\\ {\left. \mathit{\boldsymbol{I}} \right]\mathit{\boldsymbol{\tilde d}}\left( {t - \tau } \right) + \left( {\rho + \frac{{{\rho ^2}}}{{{m^2}}}} \right)\mathit{\boldsymbol{\tilde d}}_\tau ^2} \end{array} $

(18)

式中m为正实数。此时，由式(13)可得

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot V}} \le - \left( {{\lambda _{\min }}\left( \mathit{\boldsymbol{Q}} \right) - 2\eta {\lambda _{\max }}\left( \mathit{\boldsymbol{P}} \right)} \right){{\left\| {\mathit{\boldsymbol{\tilde x}}} \right\|}^2} - }\\ {\varepsilon {{\left\| {\mathit{\boldsymbol{\tilde d}}\left( t \right)} \right\|}^2} + \left( {\rho + \frac{{{\rho ^2}}}{{{m^2}}}} \right)\mathit{\boldsymbol{\tilde d}}_\tau ^2} \end{array} $

(19)

由于d_τ是有界的，因此观测器对系统状态和等效扰动的估计误差是一致最终有界的。

1.4 控制器设计

为了获得期望的跟踪性能，利用ILO估计得到的系统内部状态，设计滤波反步控制器，使系统状态实现对于期望轨迹的跟踪。控制器设计包括以下几个步骤：

1) 定义闭环系统的跟踪误差为: e₁=x₁-x_1d。此时，设计虚拟控制量为

$ {\beta _1} = - {k_1}{e_1} + {{\dot x}_{1{\rm{d}}}} $

(20)

式中k₁为正实数。

为了避免下一步中对虚拟控制量进行求导运算，引入如下的滤波器：

$ \left\{ \begin{array}{l} \dot \beta _{1,1}^c = {\omega _1}\beta _{1,2}^c\\ \dot \beta _{1,2}^c = - 2{\zeta _1}{\omega _1}\beta _{1,2}^c - \omega _1^2\left( {\beta _{1,1}^c - {\beta _1}} \right) \end{array} \right. $

(21)

式中:滤波器的初值设置为β_{1, 1}^c(t₀)=β₁(t₀)，β_{1, 2}^c(t₀)=0, t₀为系统初始时刻。此时，β_{1, 1}^c和β_{1, 2}^c分别为对于β₁及其一阶导数的估计值。

定义一个补偿跟踪误差系统：

$ \left\{ \begin{array}{l} {v_1} = {e_1} - {\xi _1}\\ {{\dot \xi }_1} = - {k_1}{\xi _1} + \left( {\beta _{1,1}^c - {\beta _1}} \right) + {\xi _2} \end{array} \right. $

(22)

初值选取为ξ(t₀)=0，ξ₂的定义将在步骤2给出。

2) 定义系统的跟踪误差为: e₂=x₂-β_{1, 1}^c。此时，设计虚拟控制量为

$ {\beta _2} = \frac{1}{{{a_{23}}}}\left[ { - {k_2}{{\hat e}_2} + \beta _{1,2}^c - {a_{21}}{{\hat x}_1} - {a_{22}}\sin \left( {{{\hat x}_1}} \right) - {v_1} - {{\hat d}_1}} \right] $

(23)

式中: k₂为正实数，$ {\hat e_2} \buildrel \Delta \over = {\hat x_2}-\beta _{1, 1}^c $为使用ILO状态估计值得到的跟踪误差, $ {\hat d_1}$为系统复合扰动d₁的估计值。

为了避免下一步中对虚拟控制量进行求导运算，同样引入如下滤波器：

$ \left\{ \begin{array}{l} \dot \beta _{2,1}^c = {\omega _2}\beta _{2,2}^c\\ \dot \beta _{2,2}^c = - 2{\zeta _2}{\omega _2}\beta _{2,2}^c - \omega _2^2\left( {\beta _{2,1}^c - {\beta _2}} \right) \end{array} \right. $

(24)

并引入一个跟踪误差补偿系统：

$ \left\{ \begin{array}{l} {v_2} = {e_2} - {\xi _2}\\ {{\dot \xi }_2} = - {k_2}{\xi _2} + {a_{23}}\left( {\beta _{2,1}^c - {\beta _2}} \right) + {a_{23}}{\xi _3} \end{array} \right. $

(25)

其中，滤波器和跟踪误差补偿系统的初值设置与1)相同，ξ₃的定义将在步骤3给出。

3) 定义系统的跟踪误差为: e₃=x₃-β_{2, 1}^c。此时，设计虚拟控制量为

$ {\beta _3} = - {k_3}{{\hat e}_3} + \beta _{2,2}^c - {a_{23}}{v_2} $

(26)

式中：k₃为正实数，$ {\hat e_3} \buildrel \Delta \over = {\hat x_3}-\beta _{2, 1}^c $为使用ILO状态估计值得到的跟踪误差。

为了避免步骤3中对虚拟控制量进行求导运算，同样引入如下的滤波器：

$ \left\{ \begin{array}{l} \dot \beta _{3,1}^c = {\omega _3}\beta _{3,2}^c\\ \dot \beta _{3,2}^c = - 2{\zeta _3}{\omega _3}\beta _{3,2}^c - \omega _3^2\left( {\beta _{3,1}^c - {\beta _3}} \right) \end{array} \right. $

(27)

并引入一个跟踪误差补偿系统：

$ \left\{ \begin{array}{l} {v_3} = {e_3} - {\xi _3}\\ {{\dot \xi }_3} = - {k_3}{\xi _3} + \left( {\beta _{3,1}^c - {\beta _3}} \right) \end{array} \right. $

(28)

滤波器和跟踪误差补偿系统的初值设置与步骤1)相同。

4) 定义系统的跟踪误差为：e₄=x₄-β_{3, 1}^c。此时，设计得到最终的控制量为

$ \begin{array}{*{20}{c}} {{\beta _4} = \frac{1}{{{b_{41}}}}\left[ { - {k_4}{{\hat e}_4} + \beta _{3,2}^c - {a_{41}}{{\mathit{\boldsymbol{\hat x}}}_1} - } \right.}\\ {\left. {{a_{42}}{{\mathit{\boldsymbol{\hat x}}}_3} - {\mathit{\boldsymbol{v}}_3} - {F_f}\left( {{{\mathit{\boldsymbol{\hat x}}}_4}} \right) - {{\mathit{\boldsymbol{\hat d}}}_2}} \right]} \end{array} $

(29)

式中：k₄为正实数，$ {\hat e_4} \buildrel \Delta \over = {\hat x_4}-\beta _{3, 1}^c$为使用ILO状态估计值得到的跟踪误差，${{\mathit{\boldsymbol{\hat d}}}_2}$为系统复合扰动d₂的估计值。

此时，取v₄=e₄，得到补偿跟踪误差系统的动态方程为

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot v}}}_1} = - {k_1}{\mathit{\boldsymbol{v}}_1} + {\mathit{\boldsymbol{v}}_2}\\ {{\mathit{\boldsymbol{\dot v}}}_2} = - {k_2}{\mathit{\boldsymbol{v}}_2} - {\mathit{\boldsymbol{v}}_1} + {a_{23}}{\mathit{\boldsymbol{v}}_3} + {k_2}{{\mathit{\boldsymbol{\tilde x}}}_2} + {a_{21}}{{\mathit{\boldsymbol{\tilde x}}}_1} + \\ \;\;\;\;\left. {{a_{22}}\left( {\sin {\mathit{\boldsymbol{x}}_1} - \sin {{\mathit{\boldsymbol{\hat x}}}_1}} \right)} \right) + \left( {{\mathit{\boldsymbol{d}}_1} - {{\mathit{\boldsymbol{\hat d}}}_1}} \right)\\ {{\mathit{\boldsymbol{\dot v}}}_3} = - {k_3}{v_3} + {v_3} - {a_{23}}{v_2} + {k_3}{{\tilde x}_3}\\ {{\mathit{\boldsymbol{\dot v}}}_4} = - {k_4}{\mathit{\boldsymbol{v}}_4} - {\mathit{\boldsymbol{v}}_3} + {k_4}{{\mathit{\boldsymbol{\tilde x}}}_4} + {a_{41}}{{\mathit{\boldsymbol{\tilde x}}}_1} + \\ \;\;\;\;{a_{42}}{{\mathit{\boldsymbol{\tilde x}}}_3} + \left( {{F_f}\left( {{\mathit{\boldsymbol{x}}_4}} \right) - {f_f}\left( {{{\mathit{\boldsymbol{\hat x}}}_4}} \right)} \right) + \left( {{\mathit{\boldsymbol{d}}_2} - {{\mathit{\boldsymbol{\hat d}}}_2}} \right) \end{array} \right. $

(30)

可以看出，v_i是跟踪误差e_i的微分同胚，镇定e_i等价于镇定虚拟跟踪误差v_i，进而，由于在定义v_i的过程中引入了指令滤波跟踪误差ξ_i，因此，针对v_i进行稳定性分析能够在镇定跟踪误差e_i的同时，镇定指令滤波器跟踪误差ξ_i。

1.5 稳定性分析

定理2  考虑式(4)所述的柔性关节机械臂模型，考虑系统的外部扰动和内部不确定性，设计式(6)所述的ILO和式(28)所述的滤波反步控制器，则系统的跟踪误差、状态估计误差和扰动估计误差是一致最终有界的。

证明  根据ILO和滤波反步法设计步骤，定义如下的Lyapunov函数：

$ \mathit{\boldsymbol{V}} = \sum\limits_{i = 1}^4 {v_i^2 + {\mathit{\boldsymbol{V}}_1}} $

(31)

考虑式(29)所述的跟踪误差补偿的动态方程，及式(18)所述的观测器稳定性结论，得到上述Lyapunov函数对时间的导数为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot V}} \le - \sum\limits_{i = 1}^4 {{k_1}\mathit{\boldsymbol{v}}_i^2} + {{\mathit{\boldsymbol{\dot V}}}_1} + {k_2}{\mathit{\boldsymbol{v}}_2}{{\mathit{\boldsymbol{\tilde x}}}_2} + {k_3}{\mathit{\boldsymbol{v}}_3}{{\mathit{\boldsymbol{\tilde x}}}_3} + {k_4}{\mathit{\boldsymbol{v}}_4}{{\mathit{\boldsymbol{\tilde x}}}_4} + }\\ {\left( {{a_{21}} + {a_{22}}} \right){\mathit{\boldsymbol{v}}_2}{{\mathit{\boldsymbol{\tilde x}}}_1} + {a_{41}}{\mathit{\boldsymbol{v}}_4}{{\mathit{\boldsymbol{\tilde x}}}_1} + {a_{42}}{\mathit{\boldsymbol{v}}_4}{{\mathit{\boldsymbol{\tilde x}}}_3} + }\\ {{\mathit{\boldsymbol{v}}_2}\left( {{\mathit{\boldsymbol{d}}_1} - {{\mathit{\boldsymbol{\hat d}}}_1}} \right) + {v_4}\left( {{\mathit{\boldsymbol{F}}_f}\left( {{\mathit{\boldsymbol{x}}_4}} \right) - {\mathit{\boldsymbol{F}}_f}\left( {{{\mathit{\boldsymbol{\hat x}}}_4}} \right)} \right) + {\mathit{\boldsymbol{v}}_4}\left( {{\mathit{\boldsymbol{d}}_2} - {{\mathit{\boldsymbol{\hat d}}}_2}} \right)} \end{array} $

(32)

由于F_f(x₄)满足Lipschitz条件，即

$ \left\| {{\mathit{\boldsymbol{F}}_f}\left( {{\mathit{\boldsymbol{x}}_4}} \right) - {\mathit{\boldsymbol{F}}_f}\left( {{{\mathit{\boldsymbol{\hat x}}}_4}} \right)} \right\| \le \xi \left\| {{\mathit{\boldsymbol{x}}_4} - {{\mathit{\boldsymbol{\hat x}}}_4}} \right\|,\xi > 0 $

(33)

因此，将式(32)展开，利用式(33)与柯西不等式可得

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot V}} \le - {k_1}{{\left\| {{\mathit{\boldsymbol{v}}_1}} \right\|}^2} - \frac{{\left( {{k_2} - {a_{21}} - {a_{22}} - 1} \right)}}{2}{{\left\| {{\mathit{\boldsymbol{v}}_2}} \right\|}^2} - }\\ {{k_3}{{\left\| {{\mathit{\boldsymbol{v}}_3}} \right\|}^2} - \left( {{k_4} - {a_{41}} - {a_{42}} - 1} \right){{\left\| {{\mathit{\boldsymbol{v}}_4}} \right\|}^2} - }\\ {\left( {{\lambda _m} - \frac{{{a_{21}} + {a_{22}} + {a_{41}}}}{2}} \right){{\left\| {{{\mathit{\boldsymbol{\tilde x}}}_1}} \right\|}^2} - }\\ {\left( {{\lambda _m} - \frac{{{k_2}}}{2}} \right){{\left\| {{{\mathit{\boldsymbol{\tilde x}}}_2}} \right\|}^2} - \left( {{\lambda _m} - \frac{{{k_3} + {a_{23}}}}{2}} \right){{\left\| {{{\mathit{\boldsymbol{\tilde x}}}_3}} \right\|}^3} - }\\ {\left( {{\lambda _m} - \frac{{{k_4} + {\xi ^2}}}{2}} \right){{\left\| {{{\mathit{\boldsymbol{\tilde x}}}_4}} \right\|}^2} - \left( {\varepsilon - 1} \right){{\left\| {\mathit{\boldsymbol{\tilde d}}\left( t \right)} \right\|}^2} + }\\ {\left( {\rho + \frac{{{\rho ^2}}}{{{m^2}}}} \right)\mathit{\boldsymbol{\tilde d}}_\tau ^2} \end{array} $

(34)

式中λ_m=-(λ_min(Q)-2ηλ_max(P))。

此时，通过选取合适的控制参数，闭环系统的补偿跟踪误差和观测器误差是一致最终有界的。进而，根据滤波反步法中补偿跟踪误差的定义可知，系统的跟踪误差和指令滤波器跟踪误差是一致最终有界的。

2 仿真分析

为了评价本文提出的柔性关节机械臂输出跟踪控制方法的有效性，本节针对柔性关节机械臂的输出跟踪问题进行了数值仿真。

根据文献[18-19]所描述的单臂柔性关节机械臂的数学模型，选取系统参数如下：I=1.35×10^-4 kg·m²，J=2.16×10^-3 kg·m²，K=7.47 N·m/rad，M=0.2 kg，L=0.02 m，g=9.8 m/s²。系统摩擦模型参数为a₁=0.016 3，a₂=0.008 11，c₁=15，c₃=1.5^[30]。在仿真中，选取系统的参考输入信号为:x_d=2πsin t。系统状态初值选取为：x₁=0，x₂=0，x₃=0，x₄=0。系统的参数摄动量为15%，取外界干扰τ_d=10sin t。设置滤波器参数为: ω_n=600 rad/s，ζ=0.9。选取的系统控制器参数为

$ {\mathit{\boldsymbol{\mu }}_1} = {\rm{diag}}\left( {0.45,0.42,0.38,0.39} \right), $

$ {\mathit{\boldsymbol{\mu }}_2} = {\left[ {0.12\;41\;0.39\;293} \right]^{\rm{T}}}, $

$ \mathit{\boldsymbol{L}} = {\left[ {0.81\;983\;12.36\;1079} \right]^{\rm{T}}}, $

$ {k_1} = 53.1,{k_2} = 42.4,{k_3} = 51.7,{k_4} = 67.7。$

分别使用本文所述控制策略与传统滑模反步控制器对柔性关节机械臂的跟踪控制效果进行对比仿真。图 3为使用两种控制策略得到的输出跟踪响应曲线，图 4为上述两种控制方法的控制输入曲线对比。由图 3可以看出，基于ILO的指令滤波反步控制器对柔性关节机械臂动力学中的匹配与非匹配干扰进行了良好的估计，并且对系统存在的摩擦进行了精确补偿，因此能够获得比常规滑模反步法更高的跟踪精度。就控制器输出特性而言，从图 4可以看出，使用滑模反步控制器得到的系统控制输入具有明显的抖振现象，而本文提出的控制策略能够得到平滑的响应曲线，且具有相对较低的功耗。

图 5列出了迭代学习观测器的状态估计曲线。由状态估计曲线可以看出，本文所述的迭代学习观测器能够对柔性关节机械臂两个惯性装置的转角与角速度进行较好的在线估计，且观测误差能够在极短的时间内收敛。柔性关节机械臂是一个典型的非线性四阶系统，ILO对其状态的估计效能表明其对一些具有较高相对阶的非线性系统具有更良好的观测精度。ILO中包含的观测更新律能够同时对系统的复合干扰进行估计，图 6(a)为ILO对非匹配干扰项d₁的估计，图 6(b)为对匹配干扰项d₂的估计。由此可以看出，无论是匹配干扰还是非匹配干扰，ILO都能够对其进行快速跟踪。且对于一个四阶系统，观测器对于扰动和状态估计的超调量较小，使得该方法具有较高的实用价值。

3 结论

1) ILO能够同时对柔性关节机械臂的关节角速度和集总干扰进行高精度估计。

2) 基于ILO的干扰估计功能设计的带有前馈补偿的指令滤波反步控制器解决了反步迭代过程中的项数膨胀问题，实现了对指令角度的跟踪。通过Lyapunov稳定性分析可以得出，闭环系统的状态是一致最终有界的。

本文所述控制策略能够有效实现柔性关节机械臂的鲁棒输出反馈跟踪控制。与传统控制方法相比，本文提出的控制策略具有更高的控制精度和较低的能量损耗。