基于磁场非线性的中低速磁悬浮鲁棒自适应控制

摘要

为了保障磁悬浮列车的悬浮稳定性，研究了中低速磁悬浮列车在磁场非线性和轨道不平顺激扰下的悬浮控制问题。首先，基于有限元方法分析了动态和静态磁场特性，建立了考虑磁饱和及涡流效应的悬浮力模型，并以该悬浮力模型为基础建立了单个悬浮单元的数学模型；然后，提出了一种鲁棒自适应控制方法，在广义（比例‒积分）PI控制框架下，设计了自适应律灵活调节控制参数，并采用李雅普诺夫方法证明了闭环系统内所有信号均是最终一致有界的；最后，在整车动力学模型上进行多种工况的仿真分析，验证了所提出控制方法的有效性。结果表明，在正弦和随机激扰下，鲁棒自适应控制的气隙跟踪误差都降低70%以上；在竖曲线工况下，相较于传统（比例‒积分‒微分）PID控制，同一个悬浮模块的前后悬浮点最大气隙跟踪误差的差值分别从1.571 8 mm和1.227 8 mm下降到0.195 2 mm和0.396 2 mm。

关键词

磁悬浮列车; 磁饱和; 涡流效应; 鲁棒自适应控制

磁悬浮列车因具有无接触、噪声小、易爬坡以及过弯能力强等优点而备受关注。随着磁悬浮列车在复杂环境下的长距离运行，其悬浮控制系统也暴露出许多问题，这主要是由轨道不平顺、风扰动以及参数不确定性等复杂干扰所引起的。因此，在复杂运行条件下保证悬浮系统的鲁棒性显得尤为重要。

目前通常采用磁路计算法与有限元分析法建立悬浮力模型^［1-2］。为了简化计算，磁路计算法常建立于多种假设之上，从而导致计算误差过大；有限元分析法可以有效提高建模精度，但其模型过于复杂且计算效率低，不利于悬浮控制器的设计。为了解决这一问题，Ni等^［3］基于有限元分析法与解析计算法建立了非线性电磁力半解析模型，并设计了反馈线性化控制器，有效抑制了气隙磁场中的谐波波动。Schmid等^［4］采用磁路计算法建立了考虑磁饱和的电磁力模型，并在此基础上提出了一种基于非线性模型预测的控制方案，提高了悬浮系统在不平顺激扰下的鲁棒性。除了磁场的静态特性，动态特性也应被考虑到悬浮力模型中，Yang等^［5］从解析计算与有限元分析两方面解释了涡流效应，并提出了一种轨道结构方案以降低涡流效应对悬浮力的影响。

悬浮控制系统作为磁悬浮列车的关键系统，对于维持车辆的悬浮安全性与乘坐舒适性至关重要，近年来受到众多学者的关注。在悬浮控制器设计早期，大多数研究者将悬浮力在平衡位置处线性展开^［6-7］，这样可以降低控制器设计难度，但同时导致电磁铁在远离平衡位置处控制性能恶化。随着智能算法的发展，众多非线性控制算法在悬浮控制系统中得到应用。陈琛等^［8］提出了一种带加速度补偿的滑模控制算法，有效抑制了电磁铁振动。Sun等^［9-11］分别采用径向基函数神经网络、神经‒模糊切换律和深度信念网络完成滑模控制器的自适应调参，在降低控制量抖振的同时提高了悬浮系统的抗干扰能力。为了在复杂扰动下实现悬浮系统的精确位置控制，Wang等^［12］设计了广义比例‒积分（PI）观测器来补偿时变干扰，并采用自适应固定时间控制器保证系统的抗扰能力与稳态性能。

现有的非线性算法都是基于简化的静态电磁力模型建立的，在车辆复杂运行工况下会出现控制精度下降的问题。针对这一问题，本文提出了一种考虑动态与静态磁场效应的电磁力模型，并在此基础上设计自适应控制方法以提高控制器对于不同运行条件的适应性。首先，建立了一种同时考虑磁饱和与涡流效应的电磁力模型；然后，设计了鲁棒自适应控制器，在不需要精确模型参数的情况下提高悬浮系统对磁场非线性、轨道不平顺以及外部干扰等的鲁棒性；最后，将所提出的电磁力模型和控制器应用于整车动力学模型，从悬浮控制性能和整车动态响应两方面反映所设计控制器的可靠性。

1 非线性电磁力建模

磁悬浮列车主要由车体和多个悬浮架组成，车体通过二系悬挂与悬浮架连接，悬浮架将负载传递到每个悬浮模块上，如图1所示。每个悬浮模块由4个线圈组成，每2个线圈组成一组最小悬浮控制单元。单个悬浮模块上的2个控制单元刚性连接，相互间的耦合力可视为外部扰动，故单点悬浮系统可以作为磁悬浮列车悬浮系统的基本控制单元^［13］。

图1  悬浮模型

Fig.1  Levitation model

如图1所示，考虑到轨道不平顺$z_{\mathrm{r}}(t)$和电磁铁位置$z_{\mathrm{m}}(t)$，悬浮间隙可描述为$\delta (t)=z_{\mathrm{m}}(t)-z_{\mathrm{r}}(t)$。此外，鉴于机械耦合、风扰动以及变化载荷等引起的外部干扰，悬浮单元受到一个总外力$F_{\mathrm{e}}(t)$，再结合牛顿第二定律，单个电磁铁的运动方程为

$$mg-F_z(t)+F_{\mathrm{e}}(t)=m\frac{{\mathrm{d}}^{\mathrm{2}}{z}_{\mathrm{m}}(t)}{\mathrm{d}{t}^{\mathrm{2}}}$$

（1）

式中：m为电磁铁质量；$F_z(t)$为悬浮力；$t$为时间；g为重力加速度。$F_z(t)$计算式为

$$F_z(t)=k\frac{I^{\mathrm{2}}(t)}{{\delta }^{\mathrm{2}}(t)}$$

（2）

式中：$I(t)$为线圈电流；$k={\mu }_{\mathrm{0}}{N}^{\mathrm{2}}A/\mathrm{4}$，其中${\mu }_{\mathrm{0}}$为空气磁导率，$N$为线圈匝数，$A$为电磁铁有效磁极面积。

将系统状态设为$x_{\mathrm{1}}=z_{\mathrm{m}}(t)$，$x_{\mathrm{2}}={\dot{z}}_{\mathrm{m}}(t)$，并将$I^{\mathrm{2}}$作为系统的输入$u$，结合式（1）、（2）可得单点悬浮单元的状态方程为：

$$\left\{\begin{array}{l}{\dot{x}}_{\mathrm{1}}=x_{\mathrm{2}}\\ {\dot{x}}_{\mathrm{2}}=g-\frac{k}{m(x_{\mathrm{1}}-z_{\mathrm{r}}{)}^{\mathrm{2}}}u+\frac{\mathrm{1}}{m}{F}_{\mathrm{e}}\\ y=x_{\mathrm{1}}\end{array}\right.$$

（3）

式（2）是在忽略磁场强度与磁通密度间非线性关系的假设下计算得到的。然而，实际系统中控制电流随车辆负载增大而增大，在复杂线路条件下控制电流可能会超过标称值而进入磁饱和区，导致模型失配，造成控制精度下降。为了提高电磁力计算精度，建立了悬浮模块的三维有限元模型，公式计算结果与有限元分析结果对比如图2所示。结果表明，电磁力随电流增大而增大，且当电流较大时计算结果和有限元分析结果误差较大。

图2  电磁力计算结果和有限元分析结果对比

Fig.2  Comparison of electromagnetic force between calculation results and FEA results

为了解决大电流作用下的悬浮力模型失配问题，由非线性磁感应强度‒磁场强度（B‒H）曲线引起的气隙磁场饱和效应以修正系数$\xi (\delta ,I)$的形式加入电磁力的计算式中^［14-15］，可得

$$\xi (\delta ,I)=\frac{p_{\mathrm{1}}+p_{\mathrm{3}}\delta +p_{\mathrm{5}}I+p_{\mathrm{7}}{\delta }^{\mathrm{2}}+p_{\mathrm{9}}{I}^{\mathrm{2}}+p_{\mathrm{11}}\delta I}{\mathrm{1}+p_{\mathrm{2}}\delta +p_{\mathrm{4}}I+p_{\mathrm{6}}{\delta }^{\mathrm{2}}+p_{\mathrm{8}}{I}^{\mathrm{2}}+p_{\mathrm{10}}\delta I}$$

（4）

式中，$p_{\mathrm{1}}$~$p_{\mathrm{11}}$为修正系数。修正后的电磁力为

$$F_{z\mathrm{s}}=\xi (\delta ,I)k\frac{I^{\mathrm{2}}(t)}{{\delta }^{\mathrm{2}}(t)}$$

（5）

采用拟牛顿法对式（4）中未知参数进行估计，结果如表1所示。

表1  修正系数的估计参数

Tab.1  Estimated parameters of correction coefficients

参数	值	参数	值
$p_1$	0.040 3	$p_7$	706.901 0
$p_2$	326.074 3	$p_8$	0.041 3
$p_3$	-11.190 0	$p_9$	-0.011 7
$p_4$	0.329 0	$p_{10}$	116.527 0
$p_5$	0.161 0	$p_{11}$	252.005 0
$p_6$	2 927.887 0

不同气隙下电磁力‒电流关系如图2所示。传统的公式计算结果与有限元分析结果在额定电流（30 A）范围内吻合度较高，但在较高电流下误差较大。相较而言，加入磁饱和修正系数后的公式计算结果和有限元分析结果在大电流范围内都吻合良好，在气隙15 mm、电流40 A工况下最大拟合误差为13.15%，额定悬浮间隙在±4 mm范围内波动，拟合误差均在4%以内。此外，在电流较小时拟合电磁力结果和实测结果^［16］的误差在10%以内，因此修正后的电磁力模型是有效的。

涡流效应因子$\varphi (v)$与悬浮模块的速度和结构参数相关^［17］。为了提高悬浮控制器的鲁棒性，对模型加入涡流效应的影响，改进电磁力计算式为

$$F_{z\mathrm{s}}=\xi (\delta ,I)\varphi (v)k\frac{I^{\mathrm{2}}(t)}{{\delta }^{\mathrm{2}}(t)}=\lambda (\cdot )k\frac{I^{\mathrm{2}}(t)}{{\delta }^{\mathrm{2}}(t)}$$

（6）

式中，$\lambda (\cdot )$为磁饱和及涡流效应综合影响因子。受磁饱和及涡流效应修正系数的影响，实际的控制输入为

$$u_{\mathrm{a}}=\lambda (\cdot )u+u_{\mathrm{h}}(v,m)$$

（7）

式中：$\mathrm{0}\le \lambda (\cdot )\le \mathrm{1}$；$u_{\mathrm{h}}$为未知项。将式（7）代入式（3）可得修正后的系统状态方程为

$$\left\{\begin{array}{l}{\dot{x}}_{\mathrm{1}}=x_{\mathrm{2}}\\ {\dot{x}}_{\mathrm{2}}=g(\cdot )\lambda (\cdot )u+g(\cdot )u_{\mathrm{h}}+f(\cdot )\\ y=x_{\mathrm{1}}\end{array}\right.$$

（8）

式中：$g(\cdot )=\frac{-k}{m(x_{\mathrm{1}}-z_{\mathrm{r}}{)}^{\mathrm{2}}}$；$f(\cdot )=g+\frac{\mathrm{1}}{m}{F}_{\mathrm{e}}$。

本文的目标是为悬浮单元设计一种轨道不平顺、非线性磁场以及外部干扰下的鲁棒自适应控制方案，使系统在实现稳定悬浮的同时保证所有内部信号都是连续有界的。为了设计控制策略，提出以下假设条件：

假设1  　控制增益$g$$(\cdot )$是时变的、未知的且有界的，即存在未知常数$\underline{g}$和$\overline{g}$使得$\mathrm{0}<\underline{g}\le g(\cdot )\le \overline{g}<\mathrm{\infty }$。

假设2  　期望气隙$y_{\mathrm{d}}$及其导数${\dot{y}}_{\mathrm{d}}$、${\ddot{y}}_{\mathrm{d}}$是有界的，即存在未知常数$y_{\mathrm{m}}$使得$\left|y_{\mathrm{d},n}\right|\le y_{\mathrm{m}}<\mathrm{\infty }$，$\forall t\ge t_{\mathrm{0}}$，其中$n=\mathrm{0,1},\mathrm{2}$。

假设3  　对于非线性不确定，存在未知常数$\alpha \ge \mathrm{0}$和已知标量函数$\phi (x,t)\ge \mathrm{0}$，使得$\left|f(\cdot )\right|\le \alpha \phi (x,t)$。若$x$是有界的，则$\phi (x,t)$也是有界的。

假设4  　$\lambda$$(\cdot )$和$u_{\mathrm{h}}$是未知的、时变的，但存在未知常数${\lambda }_{\mathrm{m}}$和$\overline{h}$使得$\mathrm{0}<{\lambda }_{\mathrm{m}}\le \lambda (\cdot )\le \mathrm{1}$及$\left|u_{\mathrm{h}}\right|\le \overline{h}$。

在解决系统（8）的跟踪控制问题时，大多数现有工作^［18-20］常采用假设1；对于悬浮系统，期望气隙值是常数，其导数为零，因此假设2是合理的；假设3为提取非线性模型的核心信息，这对于含有未知项的不确定系统是容易的^［21］；假设4中磁饱和与涡流效应系数的乘积$\lambda$$(\cdot )$显然满足位于0到1的区间内。

2 鲁棒自适应控制器设计

为了便于设计控制器，引入一个滤波变量

$$\omega ={\beta }_{\mathrm{1}}{\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{2}}$$

（9）

式中：${\beta }_{\mathrm{1}}>\mathrm{0}$；${\epsilon }_{\mathrm{1}}=x-y_{\mathrm{d}}$为位置跟踪误差；${\epsilon }_{\mathrm{2}}=\dot{x}-{\dot{y}}_{\mathrm{d}}={\dot{\epsilon }}_{\mathrm{1}}$为速度跟踪误差。求解微分方程（9）可得

$${\epsilon }_{\mathrm{1}}(t)={\mathrm{e}}^{-{\beta }_{\mathrm{1}}t}{\epsilon }_{\mathrm{1}}(\mathrm{0})+{\mathrm{e}}^{-{\beta }_{\mathrm{1}}t}{\int }_{\mathrm{0}}^t{\mathrm{e}}^{{\beta }_{\mathrm{1}}\tau }\omega (\tau )\mathrm{d}\tau$$

（10）

对式（10）采用洛必达法则，可得

$$\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{\epsilon }_{\mathrm{1}}(t)=\mathrm{0}+\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{{\int }_{\mathrm{0}}^t{\mathrm{e}}^{{\beta }_{\mathrm{1}}\tau }\omega (\tau )\mathrm{d}\tau }{{\mathrm{e}}^{{\beta }_{\mathrm{1}}t}}=$$

$$\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{{\mathrm{e}}^{{\beta }_{\mathrm{1}}t}\omega (t)}{{\beta }_{\mathrm{1}}{\mathrm{e}}^{{\beta }_{\mathrm{1}}t}}=\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\omega (t)}{{\beta }_{\mathrm{1}}}$$

（11）

由式（11）可知，若$\omega (t)\to \mathrm{0}$，则${\epsilon }_{\mathrm{1}}(t)\to \mathrm{0}$且衰减速率与$\omega (t)$相同，由式（9）可知${\epsilon }_{\mathrm{2}}(t)$亦是如此。因此，若ω最终一致有界，则${\epsilon }_{\mathrm{1}}$和${\epsilon }_{\mathrm{2}}$也最终一致有界。

构造一个误差变量

$$s=\omega +\vartheta {\int }_{\mathrm{0}}^t\omega \mathrm{d}\tau$$

（12）

式中，$\vartheta >\mathrm{0}$。同理可证，若$s\to \mathrm{0}$，则ω和${\int }_{\mathrm{0}}^t\omega \mathrm{d}\tau$同样趋近于零且衰减速率与s相同。因此，若s最终一致有界，则ω和${\int }_{\mathrm{0}}^t\omega \mathrm{d}\tau$也最终一致有界。

使用广义PI的控制形式^［22］

$$u=-(k_{\mathrm{p}}+\mathrm{\Delta }{k}_{\mathrm{p}}(t))\omega (t)-(k_{\mathrm{i}}+\mathrm{\Delta }{k}_{\mathrm{i}}(t)){\int }_{\mathrm{0}}^t\omega (\tau )\mathrm{d}\tau$$

（13）

式中，$k_{\mathrm{p}}>\mathrm{0}$，$k_{\mathrm{i}}=\vartheta k_{\mathrm{p}}>\mathrm{0}$为恒定增益。与传统PI控制下的恒定控制增益不同，广义PI控制下的控制增益由恒定项与时变项组成，时变增益$\mathrm{\Delta }{k}_{\mathrm{p}}$和$\mathrm{\Delta }{k}_{\mathrm{i}}$由以下算法自适应调整：

$$\mathrm{\Delta }{k}_{\mathrm{p}}=\frac{\widehat{b}{\psi }^{\mathrm{2}}(\cdot )}{\psi (\cdot )\left|s\right|+\zeta },\mathrm{\Delta }{k}_{\mathrm{i}}=\vartheta \mathrm{\Delta }{k}_{\mathrm{p}}$$

（14）

$$\dot{\widehat{b}}=-\sigma \gamma \widehat{b}+\frac{\sigma {\psi }^{\mathrm{2}}(\cdot )s^{\mathrm{2}}}{\psi (\cdot )\left|s\right|+\zeta }$$

（15）

式中：$b$为虚拟控制参数，$\widehat{b}$为$b$的估计值；$\psi (\cdot )=\mathrm{1}+\phi (x,t)+\left|{\epsilon }_{\mathrm{1}}\right|+\left|{\epsilon }_{\mathrm{2}}\right|$为一个易计算的标量函数；$\sigma >\mathrm{0}$，$\gamma >\mathrm{0}$为设计参数；$\zeta >\mathrm{0}$为一个小常数。

控制律结构如图3所示。广义PI控制中的恒定增益可根据需求灵活调整，时变增益则由算法自适应调优，并且为了简化控制器的设计，将2个增益通过参数$\vartheta$联系。一旦确定了广义PI控制的增益，即使存在建模不确定性和磁场非线性，也能确保系统的稳定。

图3  控制律结构

Fig.3  Control law framework

定理  　针对非线性悬浮系统（8），基于假设1~4，在控制律（13）的控制下，控制增益由式（14）、（15）实现更新。闭环系统可自适应磁场非线性与建模不确定性，保证控制误差最终一致有界。此外，亦可保证闭环系统内的所有信号均为连续有界。

证明  　结合式（8）、（9）、（12）可得

$$\dot{s}=g(\cdot )\lambda (\cdot )u+g(\cdot )u_{\mathrm{h}}+f(\cdot )-{\dot{y}}_{\mathrm{d}}+\vartheta {\beta }_{\mathrm{1}}{\epsilon }_{\mathrm{1}}+\vartheta {\epsilon }_{\mathrm{2}}$$

（16）

通过恒定增益与时变增益将（13）转化为

$$u=-(k_{\mathrm{p}}+\mathrm{\Delta }{k}_{\mathrm{p}})s$$

（17）

构造李雅普诺夫函数$V=(\mathrm{1}/\mathrm{2})s^{\mathrm{2}}+$

$(\mathrm{1}/\mathrm{2}\sigma \underline{g}){\tilde{b}}^{\mathrm{2}}$，对其求导可得

$$\dot{V}=-k_{\mathrm{p}}g(\cdot )\lambda (\cdot )s^{\mathrm{2}}-g(\cdot )\lambda (\cdot )\mathrm{\Delta }{k}_{\mathrm{p}}{s}^{\mathrm{2}}+s(g(\cdot )u_{\mathrm{h}}+$$

$$f(\cdot )-{\dot{y}}_{\mathrm{d}}+\vartheta {\beta }_{\mathrm{1}}{\epsilon }_{\mathrm{1}}+\vartheta {\epsilon }_{\mathrm{2}})+(\mathrm{1}/\sigma \underline{g})\dot{\tilde{b}}\tilde{b}$$

（18）

由假设1~4可知，$\left|g(\cdot )u_{\mathrm{h}}-{\dot{y}}_{\mathrm{d}}+f(\cdot )+\vartheta {\beta }_{\mathrm{1}}{\epsilon }_{\mathrm{1}}+\vartheta {\epsilon }_{\mathrm{2}}\right|\le \overline{g}\overline{h}+y_{\mathrm{m}}+\alpha \phi (x,t)+\vartheta {\beta }_{\mathrm{1}}\left|{\epsilon }_{\mathrm{1}}\right|+\vartheta \left|{\epsilon }_{\mathrm{2}}\right|\le b\psi$$(\cdot )$，其中$b=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\overline{g}\overline{h}+y_{\mathrm{m}},\alpha ,\vartheta {\beta }_{\mathrm{1}},\vartheta \right\}<\mathrm{\infty }$且$\psi (\cdot )=\mathrm{1}+\phi (x,t)+\left|{\epsilon }_{\mathrm{1}}\right|+\left|{\epsilon }_{\mathrm{2}}\right|$，并定义一个虚拟参数估计误差$\tilde{b}=b-\underline{g}{\lambda }_{\mathrm{m}}\widehat{b}$，再将$\mathrm{\Delta }{k}_{\mathrm{p}}$代入（18）可得

$$\dot{V}\le -k_{\mathrm{p}}\underline{g}{\lambda }_{\mathrm{m}}{s}^{\mathrm{2}}+\left(-\frac{g(\cdot )\lambda (\cdot )\widehat{b}{\psi }^{\mathrm{2}}(\cdot )s^{\mathrm{2}}}{\psi (\cdot )\left|s\right|+\zeta }+\left|s\right|b\psi (\cdot )\right)+(\mathrm{1}/\sigma \underline{g})\dot{\tilde{b}}\tilde{b}\le$$

$$\begin{array}{l}-k_{\mathrm{p}}\underline{g}{\lambda }_{\mathrm{m}}{s}^{\mathrm{2}}+\left(-\underline{g}{\lambda }_{\mathrm{m}}\frac{\widehat{b}{\psi }^{\mathrm{2}}(\cdot )s^{\mathrm{2}}}{\psi (\cdot )\left|s\right|+\zeta }+\frac{b{\psi }^{\mathrm{2}}(\cdot )s^{\mathrm{2}}+b\psi (\cdot )\left|s\right|\zeta }{\psi (\cdot )\left|s\right|+\zeta }\right)+(\mathrm{1}/\sigma \underline{g})\dot{\tilde{b}}\tilde{b}\le \\ -k_{\mathrm{p}}\underline{g}{\lambda }_{\mathrm{m}}{s}^{\mathrm{2}}+\left((b-\underline{g}{\lambda }_{\mathrm{m}}\widehat{b})\frac{{\psi }^{\mathrm{2}}(\cdot )s^{\mathrm{2}}}{\psi (\cdot )\left|s\right|+\zeta }+b\zeta \right)+(\mathrm{1}/\sigma \underline{g})\dot{\tilde{b}}\tilde{b}\le \\ -k_{\mathrm{p}}\underline{g}{\lambda }_{\mathrm{m}}{s}^{\mathrm{2}}+\left((b-\underline{g}{\lambda }_{\mathrm{m}}\widehat{b})\left(\frac{{\psi }^{\mathrm{2}}(\cdot )s^{\mathrm{2}}}{\psi (\cdot )\left|s\right|+\zeta }-\frac{\dot{\widehat{b}}}{\sigma }\right)+b\zeta \right)\le \end{array}$$

$$-k_{\mathrm{p}}\underline{g}{\lambda }_{\mathrm{m}}{s}^{\mathrm{2}}+b\zeta +\gamma \tilde{b}\widehat{b}$$

（19）

由于

$$\tilde{b}\widehat{b}=\tilde{b}(\mathrm{1}/\underline{g}{\lambda }_{\mathrm{m}})(b-\tilde{b})\le (\mathrm{1}/\mathrm{2}\underline{g}{\lambda }_{\mathrm{m}})(b^{\mathrm{2}}-{\tilde{b}}^{\mathrm{2}})$$

则可得

$$\dot{V}\le -k_{\mathrm{p}}\underline{g}{\lambda }_{\mathrm{m}}{s}^{\mathrm{2}}+b\zeta -\frac{\gamma }{\mathrm{2}\underline{g}{\lambda }_{\mathrm{m}}}{\tilde{b}}^{\mathrm{2}}+\frac{\gamma }{\mathrm{2}\underline{g}{\lambda }_{\mathrm{m}}}{b}^{\mathrm{2}}\le -\eta V+\theta$$

（20）

式中：$\eta =\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{2}{k}_{\mathrm{p}}\underline{g}{\lambda }_{\mathrm{m}},\gamma \sigma \right\}$；$\theta =\gamma b^{\mathrm{2}}/\mathrm{2}\underline{g}{\lambda }_{\mathrm{m}}+b\zeta$。由式（20）可确定一个集合${\Omega }_{\mathrm{1}}=\left\{(s,\tilde{b})|V\le \theta /\eta \right\}$，当$(s,\tilde{b})\notin {\Omega }_{\mathrm{1}}$时$\dot{V}<\mathrm{0}$，又$(\mathrm{1}/\mathrm{2})s^{\mathrm{2}}\le V(t)$，因此存在一个确定的时间使得有$\left|s\right|\le \sqrt[]{\mathrm{2}V}\le \sqrt[]{\mathrm{2}\theta /\eta }$，即s是最终一致有界的，误差${\epsilon }_{\mathrm{1}}$和${\epsilon }_{\mathrm{2}}$亦是最终一致有界的。

对于$t\ge \mathrm{0}$，求解式（20）可得$V(t)\le {\mathrm{e}}^{-\eta t}V(\mathrm{0})+\theta /\eta \in {\rho }_{\mathrm{\infty }}$，$s\in {\rho }_{\mathrm{\infty }}$且$\widehat{b}\in {\rho }_{\mathrm{\infty }}$，则存在$\omega \in {\rho }_{\mathrm{\infty }}$和${\int }_{\mathrm{0}}^t\omega \mathrm{d}\tau \in {\rho }_{\mathrm{\infty }}$，进而推得${\epsilon }_{\mathrm{1}}\in {\rho }_{\mathrm{\infty }}$及${\epsilon }_{\mathrm{2}}\in {\rho }_{\mathrm{\infty }}$，因此$x\in {\rho }_{\mathrm{\infty }}$也成立（假设2中已假设$y_{\mathrm{d}}$有界）。根据假设3，$\phi (x,t)\in {\rho }_{\mathrm{\infty }}$亦满足，再由$\psi (\cdot )=\mathrm{1}+\phi (x,t)+\left|{\epsilon }_{\mathrm{1}}\right|+\left|{\epsilon }_{\mathrm{2}}\right|$可知$\psi (\cdot )\in {\rho }_{\mathrm{\infty }}$，进而由式（14）得到$\mathrm{\Delta }{k}_{\mathrm{p}}\in {\rho }_{\mathrm{\infty }}$，$\mathrm{\Delta }{k}_{\mathrm{i}}\in {\rho }_{\mathrm{\infty }}$。根据式（13）、（15）可知$u\in {\rho }_{\mathrm{\infty }}$且$\dot{\widehat{b}}\in {\rho }_{\mathrm{\infty }}$，最后由式（16）可知$\dot{s}\in {\rho }_{\mathrm{\infty }}$。至此，系统内所有信号皆连续有界。

3 数值仿真分析

3.1　磁悬浮列车磁力耦合模型

本文建立了有23个自由度的中低速磁悬浮列车动力学模型，如图4所示。车体具有浮沉、点头和侧滚3个自由度，悬浮模块具有浮沉和点头2个自由度，相关几何参数与材料参数分别如表2、3所示，其中下标j表示第j个悬浮模块，下标i、k表示与悬浮模块连接的第i个空气弹簧和第k个电磁铁。

图4  中低速磁悬浮列车多体系统

Fig.4  Multi-body system of medium-low speed maglev vehicle

表2  动力学模型几何参数

Tab.2  Geometric parameters of dynamic model

参数	元素数目	数值/m
相邻悬浮架中心纵向距离，$c_{jx}$	j=1，2，3，4，5	5.8，2.9，0，-2.9，-5.8
空气弹簧与悬浮架中心纵向距离$,b_{ix}$	i=1，2	1.3，-1.3
空气弹簧与悬浮架中心横向距离，$b_{(\mathrm{L},\mathrm{R})y}$		0.93
头部电磁铁与悬浮架中心纵向距离，$e_{kx}$	k=1，2，3，4	1.05，0.35，-0.35，-1.05

表3  动力学模型材料参数

Tab.3  Material parameters of dynamic model

参数	值
车体质量，mc/kg	16 770
车体质量惯性矩，Icx，Icy/（kg·m2）	19 700，350 000
悬浮模块质量，mm/kg	1 000
悬浮模块质量惯性矩，Imy/（kg·m2）	400
二系悬挂刚度，kz/（N·m-1）	6$\times {10}^5$
二系悬挂阻尼系数，cz/（N·s·m-1）	3$\times {10}^3$
电磁铁质量，m/kg	$750$
电磁铁线圈匝数，N	270
电磁铁有效面积，A/m2	0.024
空气磁导率，${\mu }_0$	$4\mathrm{\pi }\times {10}^{-7}$

根据达朗贝尔原理，可以建立车体与悬浮模块的运动方程。以平衡位置为参考点，车体的浮沉运动方程为

$$m_{\mathrm{c}}{\ddot{z}}_{\mathrm{c}}=-{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{5}}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}(F_{\mathrm{L}ji}+F_{\mathrm{R}ji})$$

（21）

式中：$F_{(\mathrm{L},\mathrm{R})ji}=k_z{z}_{\mathrm{c}\mathrm{m}(\mathrm{L},\mathrm{R})ji}+c_z{\dot{z}}_{\mathrm{c}\mathrm{m}(\mathrm{L},\mathrm{R})ji}$为二系悬挂力，下标L，R表示左右侧，其中，z_{cm（L，R）ji}和$\dot{z}$_{cm（L，R）ji}表示与第j个悬浮模块相连的第i个空气弹簧的垂向变形量与变形速度，可表示为：

$$z_{\mathrm{c}\mathrm{m}(\mathrm{L},\mathrm{R})ji}=(z_{\mathrm{c}}-(c_{jx}+b_{ix}){\theta }_{\mathrm{c}y}+b_{(\mathrm{L},\mathrm{R})y}{\theta }_{\mathrm{c}x})-$$

$$(z_{\mathrm{m}j}-b_{ix}{\theta }_{\mathrm{m}jy})$$

（22）

$${\dot{z}}_{\mathrm{c}\mathrm{m}(\mathrm{L},\mathrm{R})ji}=({\dot{z}}_{\mathrm{c}}-(c_{jx}+b_{ix}){\dot{\theta }}_{\mathrm{c}y}+b_{(\mathrm{L},\mathrm{R})y}{\dot{\theta }}_{\mathrm{c}x})-$$

$$({\dot{z}}_{\mathrm{m}j}-b_{ix}{\dot{\theta }}_{\mathrm{m}jy})$$

（23）

同理，车体的侧滚和点头运动方程可以表示为：

$$I_{\mathrm{c}x}{\ddot{\theta }}_{\mathrm{c}x}=-{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{5}}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}(b_{\mathrm{L}y}(F_{\mathrm{L}ji}-F_{\mathrm{R}ji}))$$

（24）

$$I_{\mathrm{c}y}{\ddot{\theta }}_{\mathrm{c}y}={\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{5}}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}((c_{jx}+b_{ix})(F_{\mathrm{L}ji}+F_{\mathrm{R}ji}))$$

（25）

单个悬浮模块受空气弹簧二系悬挂力与电磁铁悬浮力的共同作用。悬浮模块的浮沉和点头运动方程为：

$$m_{\mathrm{m}j}{\ddot{z}}_{\mathrm{m}(\mathrm{L},\mathrm{R})j}={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}{F}_{(\mathrm{L},\mathrm{R})ji}-{\displaystyle \sum _{k=\mathrm{1}}^{\mathrm{4}}}\Delta F_{\mathrm{e}jkz}$$

（26）

$$I_{\mathrm{m}jy}{\ddot{\theta }}_{\mathrm{m}(\mathrm{L},\mathrm{R})jy}=-{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}{b}_{ix}{F}_{(\mathrm{L},\mathrm{R})ji}+{\displaystyle \sum _{k=\mathrm{1}}^{\mathrm{4}}}{e}_{kx}\mathrm{\Delta }{F}_{\mathrm{e}jkz}$$

（27）

式中，$F_{\mathrm{e}jkz}$为式（6）提出的修正电磁力。

3.2　仿真结果分析

车辆的初始悬浮气隙均设置为20 mm，目标气隙为10 mm。同时，为了验证所提出算法在维持悬浮稳定性方面的优势，基于改进的悬浮力模型（6），将传统PID与鲁棒自适应控制算法（下文简称PIDSAT和RobustPI）应用于所建立的整车动力学模型，并分别在波长23 mm、幅值3 mm的正弦谱与幅值3 mm的随机谱（见图5）2种轨道激扰下进行水平直线与竖曲线2种工况的仿真分析。传统PID控制参数为：$k_{\mathrm{P}}$=6 000，$k_{\mathrm{I}}=\mathrm{0.1}$，$k_{\mathrm{D}}=\mathrm{60}$，鲁棒自适应控制参数为：$\sigma =\mathrm{0.1}$，$\gamma =\mathrm{1}$，$\vartheta =\mathrm{3}$，$k_{\mathrm{p}}$=5 000。由于磁悬浮列车的第1位悬浮模块的第1个悬浮点受涡流效应的影响更强烈^［23］，因此选取第1位悬浮模块的前后悬浮点进行分析。仿真结果如图6~11所示，量化指标如表4所示，包括电磁铁加速度的均方根值$a_{\mathrm{m},\mathrm{R}\mathrm{M}\mathrm{S}}$、车体Sperling平稳性指标W以及气隙跟踪误差的最大绝对值$\mathrm{m}\mathrm{a}\mathrm{x}(\left|e\right|)$和均方根值$e_{\mathrm{R}\mathrm{M}\mathrm{S}}$。

图5  随机轨道不平顺

Fig.5  Random track irregularity

图6  工况1的气隙响应结果

Fig.6  Air gap response results under condition 1

图7  工况1的电流响应结果

Fig.7  Current response results under condition 1

图8  工况1的加速度响应结果

Fig.8  Acceleration response results under condition 1

图9  工况2的气隙响应结果

Fig.9  Air gap response results under condition 2

图10  工况2的电流响应结果

Fig.10  Current response results under condition 2

图11  工况2的加速度响应结果

Fig.11  Acceleration response results under condition 2

表4  2种工况下指标结果

Tab.4  Index results under two working conditions

工况	运动方向	激励	控制算法	W	$a_{\mathrm{m},\mathrm{R}\mathrm{M}\mathrm{S}}$/ （$\mathrm{m}\cdot {\mathrm{s}}^{-2}$）	$\mathrm{m}\mathrm{a}\mathrm{x}(\left|e\right|)$/mm		$e_{\mathrm{R}\mathrm{M}\mathrm{S}}$/mm
						前悬浮点	后悬浮点	前悬浮点	后悬浮点
1	水平直线	正弦	PIDSAT	2.01	0.582 4	2.533 7	1.170 4	1.121 5	0.481 1
			RobustPI	1.80	0.428 8	0.527 2	0.295 6	0.190 6	0.107 1
		随机	PIDSAT	2.23	1.092 5	2.425 5	1.235 5	0.899 3	0.365 7
			RobustPI	2.38	2.206 5	0.697 1	0.331 2	0.208 3	0.103 8
2	竖曲线	正弦	PIDSAT	1.78	0.249 1	2.565 3	0.993 5	0.681 5	0.226 2
			RobustPI	1.52	0.233 0	0.462 7	0.267 5	0.137 3	0.072 6
		随机	PIDSAT	2.11	0.806 5	2.202 4	0.974 6	0.811 4	0.369 1
			RobustPI	2.19	0.948 2	0.721 5	0.325 3	0.217 1	0.097 2

3.2.1　水平直线工况

磁悬浮列车的运行速度设置为144 km·h^-1，仿真结果如图6~8所示。结果表明，电磁铁在RobustPI下具有跟随短波激扰、维持悬浮气隙稳定的良好性能。此外，还将前后悬浮点的仿真结果进行对比，以突出涡流效应对悬浮性能的影响。在磁饱和及涡流效应的共同作用下，前悬浮点的气隙跟踪误差指标$\mathrm{m}\mathrm{a}\mathrm{x}(\left|e\right|)$和$a_{\mathrm{m},\mathrm{R}\mathrm{M}\mathrm{S}}$均明显大于后悬浮点，PIDSAT下这一差异更为明显。这是因为在高速运行条件下，涡流效应对悬浮力的衰减作用加剧，导致气隙跟踪误差增大。为了抑制跟踪误差的进一步增大，控制电流会超过标称值（见图7），使得电磁铁进入饱和作用区，从而导致PIDSAT的控制精度明显下降。由图6可知，在RobustPI下，前后悬浮点的气隙偏差仍能控制在很小范围内，但为了精确补偿气隙变化，RobustPI下的电磁铁振动也更为显著，振动响应指标略大。经过二系悬挂的衰减作用，2种控制方法下的车体振动响应指标相差不大，均能满足乘坐舒适性要求。

3.2.2　竖曲线工况

该工况下列车以90 km·h^-1的速度依次进入过渡曲线、斜坡及过渡曲线，过渡曲线选择圆曲线，半径为1 000 m，长为70 m，坡度为7%，仿真结果如图9~11所示。

列车在进入竖曲线过程中悬浮模块会产生俯仰运动，从图9中正弦谱激扰下PIDSAT气隙响应结果可以明显看出，电磁铁因俯仰运动而产生气隙变化不均匀的现象；在RobustPI下整车悬浮性能几乎不受线路条件变化的影响，与PIDSAT相比，在正弦谱与随机谱激扰下同一个悬浮模块前后悬浮点最大气隙跟踪误差的差值分别从1.571 8 mm和1.227 8 mm下降到0.195 2 mm和0.396 2 mm，有效抑制了电磁铁的俯仰运动，保证了悬浮稳定性。同样，为了及时补偿气隙变化，RobustPI下的电流波动幅值也相应增大（见图10），加速了电磁铁的振动，但经过二系悬挂的衰减作用，随机不平顺激扰下车体的舒适性指标为2.19，说明所提出的控制算法在改善气隙响应的同时能够保证良好的乘坐舒适性。

4 结语

为了提高中低速磁悬浮列车在复杂激扰下的悬浮稳定性，首先建立了考虑磁饱和及涡流效应的悬浮力模型，并在悬浮力模型的基础上提出了一种鲁棒自适应控制方法，用以解决在磁场非线性因素及轨道不平顺激扰作用下的悬浮稳定性下降问题。最后，将改进的悬浮力模型与控制方法应用在整车动力学模型上，验证了所提出的控制算法在不同行驶工况及复杂线路条件下的适应性。在周期性激扰下鲁棒自适应控制的气隙跟踪误差降低70%以上，同时电磁铁振动响应及车体平稳性指标均能提高10%左右；在随机激扰下鲁棒自适应控制气隙跟踪误差亦降低70%以上，同时也能保证车体具有较好的平稳性。此外，在竖曲线工况下，相较于传统PID，2种激扰下同一个悬浮模块的前后悬浮点最大气隙跟踪误差的差值分别从1.571 8 mm和1.227 8 mm下降到0.195 2 mm和0.396 2 mm，有效抑制了电磁铁俯仰运动，从悬浮安全性与乘坐舒适性两方面验证了鲁棒自适应控制算法的有效性。

作者贡献声明

张济民：研究构思和设计。

汪杭生：控制器设计与分析，仿真分析及论文撰写。

任 乔：模型构建，算法设计。