货运列车的车体运动方程和车轮动态模型描述为^[

式中：M为整车质量；$v_{\mathrm{t}}$为车体速度；$v_{\mathrm{d}}$为车轮线速度；$J$为电机转动惯量；$\omega$为车轮角速度；$r$为车轮半径；$T_{\mathrm{m}}$为牵引电机的输出转矩；$T_{\mathrm{z}}$为期望的牵引电机转矩；$T_{\mathrm{L}}$为货运列车的负载转矩；$F_{\mathrm{a}}$为车辆的黏着力；$F_{\mathrm{z}}$为车辆的行驶阻力；$B$为车轮黏滞摩擦因数；$\alpha$, $\epsilon$为大于0的牵引参数。

行驶阻力的表达式为^[

式中：$a_{\mathrm{1}},a_{\mathrm{2}},a_{\mathrm{3}}$为行驶阻力参数，在不同行车环境下行驶阻力参数取不同值，并不能对其直接给定^[

车辆的黏着力$F_{\mathrm{a}}$与负载转矩$T_{\mathrm{L}}$表达式为^[

式中：$\mu$为黏着因数；$W$为车辆轴重；$g$为重力加速度；$R_{\mathrm{g}}$为齿轮箱传动比。

黏着因数$\mu$与蠕滑速度$v_{\mathrm{s}}$之间存在高度的非线性关系^[

式中：$a,b,c$为轨面参数。

对于不同的轨面条件，其参数如表1所示^[

表1  不同轨面黏着参数

Table 1  Adhesion parameters on different rail surface

轨面条件	轨面参数			最优黏着因数
	a	b	c
干燥	0.54	1.00	1.200	0.286
潮湿	0.54	0.72	0.206	0.206
雨雪	0.54	0.42	1.200	0.100

Download:  CSV Download:  Table Images

联立式(1)、式(2)、式(4)和式(5)，可得

货运列车的黏着过程有着明显的非线性特征，以干燥轨面为例（如图1所示），在黏着区内，随着蠕滑速度的增加，黏着因数逐步增大；在空转区内，随着蠕滑速度的增加，黏着因数则快速减小。其余轨面的黏着特性曲线变化规律相似。黏着控制的约束问题是将列车的黏着点尽可能控制在黏着特性曲线峰值的小范围邻域内，即货运列车在行驶过程中，若能将蠕滑速度控制在黏着特性曲线峰值对应点的附近，可以避免车轮出现打滑或空转。

图1  干燥轨面黏着特性曲线

Fig. 1  Characteristics curve of adhesion on dry rail surface

Download:  Full-size image | High-res image | Low-res image

本文所提的黏着控制约束的数学表达为

式中：$v_{\mathrm{s}\mathrm{o}\mathrm{p}\mathrm{t}}$为列车行驶轨面的最优蠕滑速度。

针对式(2)设计黏着因数全维状态观测器^[

式中：${\mathit{A}}_{\mathrm{1}}=\left[\begin{array}{cc}-\frac{B}{J}& -\frac{\mathrm{1}}{J}\\ \mathrm{0}& \mathrm{0}\end{array}\right]$；${\mathit{B}}_{\mathrm{1}}=\left[\begin{array}{c}\frac{\mathrm{1}}{J}\\ \mathrm{0}\end{array}\right]$；$\mathit{C}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right]$；$u=T_{\mathrm{m}}$；$\mathit{x}={\left[\begin{array}{cc}\omega & T_{\mathrm{L}}\end{array}\right]}^{\mathrm{T}}$；$y=\omega$；L为观测器增益矩阵，$\mathit{L}={\left[\begin{array}{cc}{L}_{\mathrm{1}}& L_{\mathrm{2}}\end{array}\right]}^{\mathrm{T}}$。

利用极点配置，可得到的观测结果为：

货运列车的车体运动方程中有2个力是无法直接测量的，一个是轮轨之间的黏着力$F_{\mathrm{a}}$^[

选择车体速度$v_{\mathrm{t}}$和车轮角速度$\omega$作为状态变量，将式(1)~式(4)改写如下：

令$x_{\mathrm{1}}=v_{\mathrm{t}},x_{\mathrm{2}}=\omega ,u=T_{\mathrm{m}}$，则有

式(15)可改写为：

在车辆的行驶阻力中，需要辨识的参数为$a_{\mathrm{1}}\mathrm{、}{a}_{\mathrm{2}}\mathrm{、}$

$a_{\mathrm{3}}$。根据最小二乘法原理，可得

设计辨识误差指标为

式中：$N$为输入信号的维数；${\tau }_i$为模型第$i$个样本的输入。

为减小PSO陷入局部最优的可能性，设计动态权重为$w(i)=\frac{w_{\mathrm{m}\mathrm{a}\mathrm{x}}-w_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}}$${\mathrm{e}}^{-\frac{i}{G}{\left(w_{\mathrm{m}\mathrm{a}\mathrm{x}}-w_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)}^{\mathrm{2}}}$，其中$w(i)$为第$i$次迭代的权重，G为PSO的最大迭代次数，$w_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{0.1},w_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{0.8}$。设计非平衡学习因子$c_{\mathrm{1}}=\mathrm{0.4},c_{\mathrm{2}}=\mathrm{0.1}$。

货运列车黏着控制的前提是针对蠕滑速度在线寻优^[

图3  蠕滑速度滑模极值搜索结构

Fig. 3  Sliding mode extremum seeking structure at slip velocity

Download:  Full-size image | High-res image | Low-res image

在图3中，$g(t)$为极值搜索的辅助函数，其一阶导数$\rho$决定了极值搜索的收敛速度；滑模面为$\delta (t)=\mu -g(t)$，极值搜索增益$k$决定了稳态振荡的幅值。若要收敛速度快，则$\rho$值必须较大，但又增加了系统脱离滑模面的可能，这反而导致了SMESC无法持续锁定极值；若要削弱稳态振荡，则需要$k$值足够小。但SMESC的滑模可达条件为$\frac{\rho }{k}$足够小，这又与快速收敛和削弱稳态振荡的参数要求相悖。

针对上述问题，本文设计的改进SMESC思路为：闭环系统在极值搜索的稳态阶段，黏着观测器的观测误差是逐步减小的，以$e_{\mathrm{o}}$表示观测误差，考虑将蠕滑速度导数估计值改为$k{e}_{\mathrm{o}}\mathrm{s}\mathrm{g}\mathrm{n}\left(\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{\pi }\delta }{\beta }\right)$，以实现搜索增益在到达稳态的阶段逐步减小；而$\left|g{(t)}^{\text{'}}\right|$值需要逐步增大，因此考虑将辅助函数的斜率改为时变项$\rho +\xi t$，其中$\xi$是待设计且大于0的数，改进SMESC结构示意图如图4所示。

图4  改进滑模极值搜索结构

Fig. 4  Improved sliding mode extremum seeking structure

Download:  Full-size image | High-res image | Low-res image

改进SMESC有3个阶段：到达段、滑动段和稳态段，本小节将从上述3个阶段进行收敛性分析。为便于分析，将相关变量关系定义如下：

定义滑模函数为：

切换函数为

那么，子切换函数对应的滑模面为

以$\sigma \to \mathrm{2}n\beta$为例展开收敛性分析：

①到达段收敛性分析。

取Lyapunov函数为$V=\frac{\mathrm{1}}{\mathrm{2}}{\eta }^{\mathrm{2}}$，结合式(23)~式(25)可推出

其中，$\mathrm{\Theta }=\frac{\partial \mu (v_{\mathrm{s}})}{\partial v_{\mathrm{s}}}k\cdot \mathrm{s}\mathrm{g}\mathrm{n}\left(\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{\pi }\sigma }{\beta }\right)+\rho +\xi t$，根据Young不等式，式(26)可进一步写为

设初始状态有$V=V_{\mathrm{0}}$，滑模面到达状态有$V=V_{\mathrm{e}\mathrm{n}\mathrm{d}}=\mathrm{0}$，对式(27)进行处理可得

求解式(28)，可得

因此，系统变量可在有限时间收敛至滑模面。值得注意的是，本文提出的改进策略中，$\rho \mathrm{和}\xi$都是待设计的正数，可设计使得$\mathrm{\Theta }>\mathrm{0}$且随时间递增，则可在一定程度上缩短系统状态收敛到滑模面的时间。

②滑动段收敛性分析。

系统状态到达滑模面后，极值搜索进入滑动段，该阶段仍保持由式(26)所示的稳定性，系统状态将在滑模面$\sigma =\mu (v_{\mathrm{s}})-g(t)$上滑动，并趋近于极值。滑模面的动态方程为

当系统状态达到滑模面后有$\sigma ={\sigma }^{\text{'}}=\mathrm{0}$，整理式(30)移项可得

式中：$\varphi (v_{\mathrm{s}})=\frac{\mathrm{d}\mu (v_{\mathrm{s}})}{\mathrm{d}{v}_{\mathrm{s}}}$。

设蠕滑速度极值点为$v_{\mathrm{s}\mathrm{o}\mathrm{p}\mathrm{t}}$，则有以下方程成立：

式(32)意味着系统变量$v_{\mathrm{s}}$被迫以速度$\rho +\xi t$向极大值$v_{\mathrm{s}\mathrm{o}\mathrm{p}\mathrm{t}}$收敛。

③稳态段收敛性分析。

在传统的SMESC中（如图3所示），当系统状态收敛到极值以后，滑模面的动态方程为${\sigma }^{\text{'}}=\varphi (v_{\mathrm{s}})\cdot k\cdot \mathrm{s}\mathrm{g}\mathrm{n}(\eta )-\rho$，系统不再满足滑模条件。根据黏着特性可知，一旦速度越过黏着峰值点，则黏着曲线呈现下降趋势，此时$\varphi (v_{\mathrm{s}})<\mathrm{0}$，结合式(31)可知

此阶段$\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{\pi }\sigma }{\beta }<\mathrm{0}$，那么$\mathrm{s}\mathrm{g}\mathrm{n}\left(\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{\pi }\sigma }{\beta }\right)=-\mathrm{1}$，故式(33)本质上是$v_{\mathrm{s}}^{\text{'}}=-k$；而极值搜索在系统状态脱离滑模面后会迫使其再度进入到达阶段，此过程有$v_{\mathrm{s}}^{\text{'}}=k$。因此，在稳态段滑模极值搜索的输出会在$(-k,k)$区间振荡。在上文设计的观测器，定义其观测误差为$e_{\mathrm{o}}=\widehat{\mu }(v_{\mathrm{s}})-\mu (v_{\mathrm{s}})$，经过极点配置后，观测器的误差会随时间的增大趋近于0。因此，将蠕滑速度导数设计为$v_{\mathrm{s}}^{\text{'}}=k\cdot e_{\mathrm{o}}\mathrm{s}\mathrm{g}\mathrm{n}\left(\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{\pi }\sigma }{\beta }\right)$，可使得在稳态段有$v_{\mathrm{s}}^{\text{'}}=\pm k\cdot e_{\mathrm{o}}\to \mathrm{0}$，有望消除或大幅减弱SMESC的稳态振荡。

黏着控制是一种非典型的输出受限问题，具体表现为控制过程以限制轮轨黏着状态在最佳区域邻域内为目标。本文在实现无稳态振荡最优蠕滑速度搜索后，以锁定的最优蠕滑速度为跟踪目标，基于障碍Lyapunov函数的黏着控制器过程如下：

步骤1：定义蠕滑速度的跟踪误差为$e_{\mathrm{s}}$，其动态方程如下：

由于黏着是动态过程，构造时变障碍参数$k_{\mathrm{a}}=v_{\mathrm{s}\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathrm{*}}+k_{\mathrm{s}}$，其中$v_{\mathrm{s}\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathrm{*}}$为极值搜索的实时结果，$k_{\mathrm{s}}$为$V_{\mathrm{2}}$中的参数，取Lyapunov函数$V_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{4}}\left[\mathrm{1}-\mathrm{s}\mathrm{g}\mathrm{n}(e)\right]\mathrm{l}\mathrm{g}\frac{k_{\mathrm{a}}^{\mathrm{2}}}{k_{\mathrm{a}}^{\mathrm{2}}-e^{\mathrm{T}}e}$，以限定系统的整体障碍界限，对其求导可得

为保障误差$e_{\mathrm{s}}$稳定收敛，可设计虚拟控制律为

步骤2：定义转矩跟踪误差为$e_{\mathrm{T}}=T_{\mathrm{m}}-T_{\mathrm{m}}^d$，其动态方程可表达为

构造Lyapunov函数$V_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{4}}\left[\mathrm{1}+\mathrm{s}\mathrm{g}\mathrm{n}\left(e_T\right)\right]\mathrm{l}\mathrm{g}\frac{k_{\mathrm{s}}^{\mathrm{2}}}{k_{\mathrm{s}}^{\mathrm{2}}-e_{\mathrm{T}}}$，对其求导可得

为保障$V_{\mathrm{2}}^{\text{'}}\le \mathrm{0}$，控制律设计如下：

取整体系统Lyapunov函数为$V=\frac{\mathrm{1}}{\mathrm{4}}\left[\mathrm{1}-\mathrm{s}\mathrm{g}\mathrm{n}(e_{\mathrm{s}})\right]V_{\mathrm{1}}+\frac{\mathrm{1}}{\mathrm{4}}\left[\mathrm{1}+\mathrm{s}\mathrm{g}\mathrm{n}(e_{\mathrm{s}})\right]V_{\mathrm{2}}$，对其求导可得

情况1：当$e_{\mathrm{s}}<\mathrm{0}$时，$\mathrm{s}\mathrm{g}\mathrm{n}(e_{\mathrm{s}})=-\mathrm{1}$，障碍Lyapunov函数表达式为

情况2：当$e_{\mathrm{s}}>\mathrm{0}$时，$\mathrm{s}\mathrm{g}\mathrm{n}(e_{\mathrm{s}})=\mathrm{1}$，障碍Lyapunov函数表达式为

综上所述，控制律在闭环系统下能够使系统稳定。