磁悬浮列车随机非线性最优控制研究

Research on Stochastic Nonlinear Optimal Control of Maglev Trains

摘要

为提高磁悬浮系统的稳定性，考虑气动升力和悬浮力的非线性特征，建立基于Hamilton理论的受控磁悬浮系统随机微分方程。以可靠度最大、平均首次穿越时间最长和最大Lyapunov指数最小为目标，建立最优控制策略的动态规划方程。研究结果表明：考虑PD控制和最优控制联合作用的磁悬浮系统提高了磁悬浮系统的条件可靠性，延长了平均首次穿越时间；经随机稳定化后最大Lyapunov指数始终为负，满足磁悬浮系统平凡解概率为1渐近稳定的条件；经过最优控制后，系统的联合概率密度发生了性态的变化，提高了系统的稳定性；在高斯白噪声强度较低时，以可靠度最大为目标的最优控制策略具有较好的性能指标；以最大Lyapunov指数最小为目标时，仅在某一范围内具有较好的性能指标。通过研究磁悬浮列车的最优控制问题，为提高列车稳定性、延长发生首次穿越失效的时间提供理论依据。

Abstract

To improve the stability of the maglev systems， the stochastic differential equations for the controlled maglev systems were established based on Hamilton's theory， where the nonlinear characteristics of the aerodynamic lift and levitation forces were taken into account. The dynamic planning equations for an optimal control strategy were developed with the objectives of maximizing the reliability， extending the longest average first-passage time， and minimizing the maximum Lyapunov exponent. The results show that the conditional reliability of the maglev systems may be improved and the average first-passage time prolonged by considering the joint action of PD control and optimal control. Moreover， the maximum Lyapunov exponent is always negative， satisfying the conditions for the trivial solution of the maglev systems to be asymptotically stable with a probability of 1. After optimal control， the joint probability density of the systems undergoes a change in behavior， which improves the system's stability. When the intensity of Gaussian white noise is low， the optimal control strategy for maximum reliability has better performance indicators. However， the strategy for minimum the maximum Lyapunov exponent only exhibits good performance within a certain range. The study of the optimal control problem of the maglev trains provides a theoretical basis for improving the train's stability and prolonging the time until the first-passage failure occurs.

Graphical abstract

关键词

磁悬浮列车 / 控制策略 / 随机激励 / 非线性随机动力学

Key words

maglev train / control strategy / random excitation / nonlinear stochastic dynamics

Author summay

刘伟渭^*（通信作者），男，1984年生，副教授。研究方向为轨道车辆和磁悬浮列车系统动力学。 E-mail：[e-mail]。

0 引言

磁悬浮列车是一种新型的交通运输方式^［1］。与传统的轮轨接触式运输方式相比，磁悬浮列车具有爬坡能力强、转弯半径小、运行平稳等优点。由于磁悬浮列车悬浮于轨道之上，因此不会受到车轮与轨道之间黏着效应的影响^［2］。常导电磁悬浮（electromagnetic suspension，EMS）列车存在自激振动的现象，当振动过大时，系统出现首次穿越失效，对列车运行的安全性和可靠性产生影响，并且由于悬浮力的非线性特征，使得磁悬浮系统是开环不稳定的^［3］。由此可知，研究磁悬浮系统的控制策略对提高列车运行的可靠性具有重要意义。

对于磁悬浮系统稳定性和可靠性的分析，已有较多学者展开相关研究。XU等^［4］为提高磁悬浮列车运行时的可靠性，提出了一种基于数据驱动的稳定性能监测指标，能够对运行中的控制器进行重新分配。CHEN等^［5］通过建立简化的悬浮模块横向运动和偏航运动的微分方程来对中低速磁悬浮的横向稳定性进行研究，采用周期性系数的非自治微分方程来表征横向稳定性，其稳定性等同于Poincare map的稳定性。YU等^［6］、薄凯等^［7］采用正交试验设计方法给出了电磁铁浮重比的最优参数组合，提高了磁悬浮结构的可靠性。CHANG等^［8］利用最大Lyapunov指数的分析方法来确定磁悬浮混沌运动的起点，并发现悬浮力的非线性特性会导致周期加倍分岔，从而引发混沌运动。CHEN等^［9］针对EMS列车建立基于柔性轨道梁的单个电磁铁悬浮模型，计算了轨道梁3个关键动力学参数的稳定域以及退化平衡点处的Lyapunov系数。

针对磁悬浮系统的控制理论已有较多学者展开了研究，但主要是采用PID控制策略以系统响应最小为控制目标。FENG等^［10］考虑时间延迟的单电磁铁悬浮系统主动控制，对不同运行速度和控制参数在时间延迟作用下的动态响应进行了分析。SUN等^［11-12］采用模糊自适应调节PID控制技术对磁悬浮系统发生的车辆-导轨耦合振动问题进行了Hopf分岔研究，通过调节模糊PID控制器的3个增益参数，使磁悬浮系统远离Hopf分岔点，从而抑制车辆与导轨相互作用的振动。HU等^［13］对一种由高温超导电磁铁和普通导电电磁铁组成的新型混合磁悬浮系统进行了PID控制和无约束鲁棒模型预测控制对比分析实验。XU等^［14］针对两级悬挂的磁悬浮车辆模型，将轨道不平顺产生的随机激励作为振源，采用电压反馈的PID控制分析了不同运行速度下悬浮间隙和悬浮电流的变化。刘伟渭等^［15］建立了弹性约束轮对的随机几何模型，并考虑了轨道随机不平顺性以及橡胶弹簧材质不均匀性等多种非线性因素对轮对系统的影响。总体来看，目前对磁悬浮系统的最优控制研究尚不充分，可以通过研究磁悬浮系统的控制策略来提高列车运行的首次穿越失效可靠性和稳定性。

基于此，本文考虑气动升力和悬浮电磁力的非线性特征，建立随机参激和随机外激作用下的随机动力学方程，基于Hamilton理论探究PD控制和最优控制相结合的磁悬浮系统控制策略，为提高磁悬浮系统的稳定性和可靠性提供参考。

1 受控磁悬浮系统建模与求解

1.1 随机动力学模型建立

本文研究的磁悬浮系统基于常导型磁悬浮列车。磁浮力学模型参考一种基于磁悬浮两级悬浮控制系统的耦合动力学模型。

本文研究基于以下假设条件：①忽略电磁铁错位布置、电磁铁模块内部的动力学耦合、空气弹簧对悬浮系统动态特性的影响等非必要因素，使复杂的系统得以简化，以提高求解效率；②将外部激励、磁浮架与磁浮车质量的波动一致近似为高斯白噪声激励。磁悬浮系统的动力学模型^［16-17］ 如图1所示。车体与悬浮架之间通过空气弹簧传递作用力，悬浮电磁铁固定连接在悬浮架上。通过控制悬浮电磁铁上线圈的电流大小，使导轨与悬浮电磁铁间具有一定的悬浮间隙，从而保证磁悬浮列车的正常运行。

模型考虑空气动力载荷对车体的影响和悬浮电磁力的非线性特征。本文主要研究磁悬浮列车的垂向动力学，忽略列车表面压力和气动阻力的影响。EMS型磁悬浮列车的参数如表1所示。由于磁悬浮列车的迎风速度远大于侧风速度，因此忽略侧风，只考虑列车迎风气流的影响。作用在车体上的气动升力F_L可表示为

F_{\mathrm{L}}=F_{\mathrm{L}\mathrm{0}}+F_{\mathrm{L}\mathrm{\Delta }}（1）

\begin{array}{l}{F}_{\mathrm{L}\mathrm{0}}=\frac{\mathrm{1}}{\mathrm{2}}{C}_{\mathrm{L}}\rho A_{\mathrm{v}}{v}^{\mathrm{2}}\\ F_{\mathrm{L}\mathrm{\Delta }}=A_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\mathrm{\pi }ft)\end{array}

式中：F_L0为常定气动升力；F_LΔ为非常定气动升力；C_L为气动升力系数；ρ为空气密度；A_v为车辆参考面积；v为车速；A₀为非常定气动升力的幅值；f为非常定气动升力的周期频率；t为时间。

EMS型磁悬浮列车依靠固定于悬浮架上的电磁铁与导轨之间的相互吸引来实现悬浮和导向。电磁力具有强非线性特征，根据麦克斯韦方程可以得到电磁力F_m的计算表达式为

F_{\mathrm{m}}=\frac{{\mu }_{\mathrm{0}}{N}^{\mathrm{2}}A}{\mathrm{4}}{(\frac{I(t)}{c(t)})}^{\mathrm{2}}=\frac{{\mu }_{\mathrm{0}}{N}^{\mathrm{2}}A}{\mathrm{4}}{(\frac{I_{\mathrm{0}}+\mathrm{\Delta }I}{c_{\mathrm{0}}+\mathrm{\Delta }c})}^{\mathrm{2}}（2）

式中：{\mu }_{\mathrm{0}}为真空磁导率；N为线圈匝数；A为电磁铁有效面积；I(t)为t时刻下的瞬时悬浮电流；I_{\mathrm{0}}为稳定悬浮电流；c(t)表示t时刻下的瞬时悬浮间隙；c_{\mathrm{0}}为标准悬浮间隙；\mathrm{\Delta }I为相对于稳定悬浮的电流扰动；\mathrm{\Delta }c为悬浮间隙扰动。

由于悬浮力与悬浮间隙的平方成反比，故磁悬浮系统是开环不稳定的。本文通过比例微分（PD）控制器对悬浮电流进行反馈控制，控制规则如下：

\mathrm{\Delta }I=k_{\mathrm{p}}\mathrm{\Delta }c+k_{\mathrm{d}}\mathrm{\Delta }\dot{c}（3）

式中：k_{\mathrm{p}}为悬浮间隙增益反馈系数；k_{\mathrm{d}}为速度反馈系数。

为提高电磁力计算的精度，将式（2）进行泰勒级数展开^［17］，并保留二阶精度，则有：

F_{\mathrm{m}}=F_{\mathrm{m}\mathrm{0}}+\mathrm{\Delta }{F}_{\mathrm{m}}（4）

\begin{array}{c}\mathrm{\Delta }{F}_{\mathrm{m}}=k_{\mathrm{f}}(\mathrm{2}{I}_{\mathrm{0}}{k}_{\mathrm{p}}-\frac{\mathrm{2}{I}_{\mathrm{0}}^{\mathrm{2}}}{c_{\mathrm{0}}})\mathrm{\Delta }c+\mathrm{2}{I}_{\mathrm{0}}{k}_{\mathrm{f}}{k}_{\mathrm{d}}\mathrm{\Delta }\dot{c}+\\ k_{\mathrm{f}}(k_{\mathrm{p}}^{\mathrm{2}}-\frac{\mathrm{4}{I}_{\mathrm{0}}{k}_{\mathrm{p}}}{c_{\mathrm{0}}}+\frac{\mathrm{3}{I}_{\mathrm{0}}^{\mathrm{2}}}{c_{\mathrm{0}}^{\mathrm{2}}})\mathrm{\Delta }{c}^{\mathrm{2}}+k_{\mathrm{f}}(\mathrm{2}{k}_{\mathrm{p}}{k}_{\mathrm{d}}-\frac{\mathrm{4}{I}_{\mathrm{0}}{k}_{\mathrm{d}}}{c_{\mathrm{0}}})·\\ \mathrm{\Delta }c\mathrm{\Delta }\dot{c}+k_{\mathrm{f}}{k}_{\mathrm{d}}^{\mathrm{2}}\mathrm{\Delta }{\dot{c}}^{\mathrm{2}}\\ F_{\mathrm{m}\mathrm{0}}=\frac{{\mu }_{\mathrm{0}}A{N}^{\mathrm{2}}{I}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}{c}_{\mathrm{0}}}\end{array}

k_{\mathrm{f}}=\frac{{\mu }_{\mathrm{0}}A{N}^{\mathrm{2}}}{\mathrm{4}{c}_{\mathrm{0}}^{\mathrm{2}}}{I}_{\mathrm{0}}=\frac{\mathrm{2}{c}_{\mathrm{0}}}{N}\sqrt[]{\frac{(m_{\mathrm{1}}+m_{\mathrm{2}})g-F_{\mathrm{L}\mathrm{0}}}{{\mu }_{\mathrm{0}}A}}

式中：m_{\mathrm{1}}为悬浮架的总质量；m_{\mathrm{2}}为车体的质量；g为重力加速度。

对于简化的磁悬浮系统，将空气弹簧简化为弹簧阻尼系统，将导轨考虑为刚性^［18］，则得到如下受控随机激励的磁悬浮系统运动方程^［19］：

\left.\begin{array}{l}{m}_{\mathrm{1}}{\ddot{z}}_{\mathrm{1}}=k_{\mathrm{2}}(z_{\mathrm{2}}-z_{\mathrm{1}})+c_{\mathrm{2}}({\dot{z}}_{\mathrm{2}}-{\dot{z}}_{\mathrm{1}})-\mathrm{\Delta }{F}_{\mathrm{m}}+\\ {\beta }_{\mathrm{11}}{\xi }_{\mathrm{1}}(t)z_{\mathrm{1}}+{\beta }_{\mathrm{12}}{\xi }_{\mathrm{2}}(t)+u\\ m_{\mathrm{2}}{\ddot{z}}_{\mathrm{2}}=-k_{\mathrm{2}}(z_{\mathrm{2}}-z_{\mathrm{1}})-c_{\mathrm{2}}({\dot{z}}_{\mathrm{2}}-{\dot{z}}_{\mathrm{1}})-F_{\mathrm{L}\mathrm{0}}+\\ {\beta }_{\mathrm{21}}{\xi }_{\mathrm{3}}(t)z_{\mathrm{2}}+{\beta }_{\mathrm{22}}{\xi }_{\mathrm{4}}(t)\end{array}\right\}（5）

式中：k_{\mathrm{2}}为空气弹簧垂向刚度；c_{\mathrm{2}}为空气弹簧垂向阻尼；z_{\mathrm{1}}为悬浮架垂向位移；z_{\mathrm{2}}为车体垂向位移；u为控制力；{\beta }_{\mathrm{11}}为悬浮架随机参激控制项；{\beta }_{\mathrm{21}}为车体随机参激控制项；{\beta }_{\mathrm{12}}、{\beta }_{\mathrm{22}}为随机外激，分别表示非常定气动升力和轨道随机不平顺^［20］；{\xi }_k(t)(k=\mathrm{1,2},\mathrm{3,4})表示4个独立的高斯白噪声。

本文将系统结构属性的随机变化等效为刚度的随机变化，{\xi }_k(t)(k=\mathrm{1,2},\mathrm{3,4})分别对应随机刚度对磁轮的作用、随机刚度对车体的作用、空气动力不平顺过程和轨道不平顺随机过程，是4个均值为0、强度为2D的高斯白噪声过程（D为白噪声强度），其自相关函数是\delta函数，可表示为

E({\xi }_k(t){\xi }_l(t+\Delta ))=\mathrm{2}{D}_{kl}\delta (\Delta )（6）

k,l=\mathrm{1,2},\mathrm{3,4}

式中：\xi (t)表示高斯白噪声；\Delta为时间差；D_kl 为高斯白噪声强度矩阵 D 中的元素，它揭示了不同噪声之间的强度关联；\delta (\Delta )为Dirac \delta函数；E（·）表示数学期望。

1.2 随机平均法求解

采用基于Hamilton理论的随机平均法对随机系统进行求解^［21］。式（5）所示为一维非线性随机激励动力学过程。为对该过程进行进一步分析，需将此随机运动方程转化为Hamilton体系。首先取广义位移q_{\mathrm{1}}=z_{\mathrm{1}}、q_{\mathrm{2}}=z_{\mathrm{2}}和广义动量p_{\mathrm{1}}=m_{\mathrm{1}}{\dot{z}}_{\mathrm{1}}、p_{\mathrm{2}}=m_{\mathrm{2}}{\dot{z}}_{\mathrm{2}}，则式（5）可改写为

\left.\begin{array}{l}{\dot{q}}_{\mathrm{1}}=\frac{p_{\mathrm{1}}}{m_{\mathrm{1}}}\\ {\dot{p}}_{\mathrm{1}}=-k_{\mathrm{11}}{q}_{\mathrm{1}}+k_{\mathrm{12}}{q}_{\mathrm{2}}-k_{\mathrm{13}}{p}_{\mathrm{1}}+k_{\mathrm{14}}{p}_{\mathrm{2}}-k_{\mathrm{15}}{q}_{\mathrm{1}}^{\mathrm{2}}-\\ k_{\mathrm{16}}{p}_{\mathrm{1}}^{\mathrm{2}}-k_{\mathrm{17}}{q}_{\mathrm{1}}{p}_{\mathrm{1}}+k_{\mathrm{18}}\xi (t)q_{\mathrm{1}}+k_{\mathrm{19}}\xi (t)+u\\ {\dot{q}}_{\mathrm{2}}=\frac{p_{\mathrm{1}}}{m_{\mathrm{2}}}\\ {\dot{p}}_{\mathrm{2}}=k_{\mathrm{21}}{q}_{\mathrm{1}}-k_{\mathrm{21}}{q}_{\mathrm{2}}+k_{\mathrm{22}}{p}_{\mathrm{1}}-k_{\mathrm{23}}{p}_{\mathrm{2}}-\\ k_{\mathrm{24}}+k_{\mathrm{25}}\xi (t)q_{\mathrm{2}}+k_{\mathrm{26}}\xi (t)\end{array}\right\}（7）

\begin{array}{c}{k}_{\mathrm{11}}=k_{\mathrm{2}}+k_{\mathrm{f}}(\mathrm{2}{I}_{\mathrm{0}}{k}_{\mathrm{p}}-\frac{\mathrm{2}{I}_{\mathrm{0}}^{\mathrm{2}}}{c_{\mathrm{0}}})k_{\mathrm{12}}=k_{\mathrm{2}}\\ k_{\mathrm{13}}=\frac{c_{\mathrm{2}}+\mathrm{2}{I}_{\mathrm{0}}{k}_{\mathrm{f}}{k}_{\mathrm{d}}}{m_{\mathrm{1}}}{k}_{\mathrm{14}}=\frac{c_{\mathrm{2}}}{m_{\mathrm{2}}}\\ k_{\mathrm{15}}=k_{\mathrm{f}}(k_{\mathrm{p}}^{\mathrm{2}}-\frac{\mathrm{4}{I}_{\mathrm{0}}{k}_{\mathrm{p}}}{c_{\mathrm{0}}}+\frac{\mathrm{3}{I}_{\mathrm{0}}^{\mathrm{2}}}{c_{\mathrm{0}}^{\mathrm{2}}})k_{\mathrm{16}}=\frac{k_{\mathrm{f}}{k}_{\mathrm{d}}^{\mathrm{2}}}{m_{\mathrm{1}}^{\mathrm{2}}}\\ k_{\mathrm{17}}=\frac{k_{\mathrm{f}}(\mathrm{2}{k}_{\mathrm{p}}{k}_{\mathrm{d}}-\frac{\mathrm{4}{I}_{\mathrm{0}}{k}_{\mathrm{d}}}{c_{\mathrm{0}}})}{m_{\mathrm{1}}}{k}_{\mathrm{18}}={\beta }_{\mathrm{11}}{k}_{\mathrm{19}}={\beta }_{\mathrm{12}}\\ k_{\mathrm{21}}=k_{\mathrm{2}}{k}_{\mathrm{22}}=\frac{c_{\mathrm{2}}}{m_{\mathrm{1}}}{k}_{\mathrm{23}}=\frac{c_{\mathrm{2}}}{m_{\mathrm{2}}}\\ k_{\mathrm{24}}=F_{\mathrm{L}}{k}_{\mathrm{25}}={\beta }_{\mathrm{21}}{k}_{\mathrm{26}}={\beta }_{\mathrm{22}}\end{array}

然后可得该磁悬浮系统的Hamilton函数为

H=\frac{p_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{1}}}+\frac{p_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{2}}}+U(q_{\mathrm{1}},q_{\mathrm{2}})=\frac{p_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{1}}}+\frac{p_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{2}}}+

\frac{k_{\mathrm{11}}}{\mathrm{2}}{q}_{\mathrm{1}}^{\mathrm{2}}+\frac{k_{\mathrm{21}}}{\mathrm{2}}{q}_{\mathrm{2}}^{\mathrm{2}}-k_{\mathrm{12}}{q}_{\mathrm{1}}{q}_{\mathrm{2}}+\frac{k_{\mathrm{15}}}{\mathrm{3}}{q}_{\mathrm{1}}^{\mathrm{3}}+k_{\mathrm{24}}{q}_{\mathrm{2}}（8）

其中，U(q_{\mathrm{1}},q_{\mathrm{2}})表示系统广义势能。U(q_{\mathrm{1}},q_{\mathrm{2}})不能继续分离变量，并且不存在与Hamilton函数H独立、对合的首次积分，则磁悬浮系统（式（5））为拟不可积Hamilton系统。由于随机激励只与位移有关，因此Wong-Zaikai修正项为0。将式（7）转化为如下Itô随机微分方程：

\left.\begin{array}{l}\mathrm{d}{q}_{\mathrm{1}}=\frac{\partial H}{\partial p_{\mathrm{1}}}\mathrm{d}t\\ \mathrm{d}{p}_{\mathrm{1}}=-(\frac{\partial H}{\partial q_{\mathrm{1}}}+m_{\mathrm{11}}\frac{\partial H}{\partial p_{\mathrm{1}}}+m_{\mathrm{12}}\frac{\partial H}{\partial p_{\mathrm{2}}}-u)\mathrm{d}t+\\ {\sigma }_{\mathrm{11}}\mathrm{d}{B}_{\mathrm{1}}(t)+{\sigma }_{\mathrm{12}}\mathrm{d}{B}_{\mathrm{2}}(t)\\ \mathrm{d}{q}_{\mathrm{2}}=\frac{\partial H}{\partial p_{\mathrm{2}}}\mathrm{d}t\\ \mathrm{d}{p}_{\mathrm{2}}=-(\frac{\partial H}{\partial q_{\mathrm{2}}}+m_{\mathrm{21}}\frac{\partial H}{\partial p_{\mathrm{1}}}+m_{\mathrm{22}}\frac{\partial H}{\partial p_{\mathrm{2}}})\mathrm{d}t+\\ {\sigma }_{\mathrm{21}}\mathrm{d}{B}_{\mathrm{3}}(t)+{\sigma }_{\mathrm{22}}\mathrm{d}{B}_{\mathrm{4}}(t)\end{array}\right\}（9）

其中，B_k(t)(k=\mathrm{1,2},\mathrm{3,4})表示相互独立的标准Wiener过程，变量代换如下：

\begin{array}{c}{m}_{\mathrm{11}}=c_{\mathrm{2}}+\mathrm{2}{I}_{\mathrm{0}}{k}_{\mathrm{f}}{k}_{\mathrm{d}}+\frac{k_{\mathrm{f}}{k}_{\mathrm{d}}^{\mathrm{2}}}{m_{\mathrm{1}}}{p}_{\mathrm{1}}+k_{\mathrm{f}}(\mathrm{2}{k}_{\mathrm{p}}{k}_{\mathrm{d}}-\frac{\mathrm{4}{I}_{\mathrm{0}}{k}_{\mathrm{d}}}{c_{\mathrm{0}}})q_{\mathrm{1}}\\ m_{\mathrm{12}}=-c_{\mathrm{2}}{m}_{\mathrm{21}}=-c_{\mathrm{2}}{m}_{\mathrm{22}}=c_{\mathrm{2}}\\ {\sigma }_{\mathrm{11}}=k_{\mathrm{18}}{q}_{\mathrm{1}}\sqrt[]{\mathrm{2}D}{\sigma }_{\mathrm{12}}=k_{\mathrm{19}}\sqrt[]{\mathrm{2}D}\\ {\sigma }_{\mathrm{21}}=k_{\mathrm{25}}{q}_{\mathrm{2}}\sqrt[]{\mathrm{2}D}{\sigma }_{\mathrm{22}}=k_{\mathrm{26}}\sqrt[]{\mathrm{2}D}\end{array}

根据Khasminskii定理，利用一维扩散过程来近似得到H，故支配该一维扩散过程中未受最优控制的平均Itô型随机微分方程为

\mathrm{d}H=\overline{m}(H)\mathrm{d}t+\overline{\sigma }(H)\mathrm{d}B(t)（10）

其中，\overline{m}(H)、\overline{\sigma }(H)分别为一维扩散过程的漂移系数和扩散系数^［22］，其表达式分别如下：

\begin{array}{c}\overline{m}(H)=\frac{\mathrm{1}}{T(H)}{\int }_{\Omega }(-m_{ij}\frac{\partial H}{\partial p_i}\frac{\partial H}{\partial p_j}+\frac{\mathrm{1}}{\mathrm{2}}{\sigma }_{ik}{\sigma }_{jk}\frac{{\partial }^{\mathrm{2}}H}{\partial p_i\partial p_j})/\\ \frac{\partial H}{\partial p_{\mathrm{1}}}\mathrm{d}{p}_{\mathrm{2}}\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{q}_{\mathrm{2}}=-(\frac{c_{\mathrm{2}}}{m_{\mathrm{1}}}+\frac{c_{\mathrm{2}}}{m_{\mathrm{2}}}+\frac{\mathrm{2}{k}_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{m_{\mathrm{1}}})H+\\ \frac{D{k}_{\mathrm{26}}^{\mathrm{2}}}{m_{\mathrm{2}}}+\frac{D{k}_{\mathrm{19}}^{\mathrm{2}}}{m_{\mathrm{1}}}+R^{\mathrm{2}}(\frac{c_{\mathrm{2}}}{\mathrm{4}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{4}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{k_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{\mathrm{2}{m}_{\mathrm{1}}}+\end{array}

\frac{c_{\mathrm{2}}}{\mathrm{4}{m}_{\mathrm{2}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{4}{k}_{\mathrm{21}}{m}_{\mathrm{2}}})+R^{\mathrm{4}}(\frac{k_{\mathrm{15}}{k}_{\mathrm{d}}{k}_{\mathrm{f}}{k}_{\mathrm{p}}}{\mathrm{12}{k}_{\mathrm{11}}^{\mathrm{2}}{m}_{\mathrm{1}}}-\frac{k_{\mathrm{15}}{k}_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{\mathrm{6}{c}_{\mathrm{0}}{k}_{\mathrm{11}}^{\mathrm{2}}{m}_{\mathrm{1}}})（11）

\begin{array}{l}{\overline{\sigma }}^{\mathrm{2}}(H)=\frac{\mathrm{1}}{T(H)}{\int }_{{\Omega }_{}}{\sigma }_{ik}{\sigma }_{jk}\frac{\partial H}{\partial p_i}\frac{\partial H}{\partial p_j}/\frac{\partial H}{\partial p_{\mathrm{1}}}\mathrm{d}{p}_{\mathrm{2}}\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{q}_{\mathrm{2}}=\\ (\frac{\mathrm{2}D{k}_{\mathrm{19}}^{\mathrm{2}}}{m_{\mathrm{1}}}+\frac{\mathrm{2}D{k}_{\mathrm{26}}^{\mathrm{2}}}{m_{\mathrm{2}}})H-R^{\mathrm{4}}(\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{6}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{6}{k}_{\mathrm{21}}{m}_{\mathrm{2}}})+\end{array}

R^{\mathrm{2}}(\frac{DH{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{DH{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{21}}{m}_{\mathrm{2}}}-\frac{D{k}_{\mathrm{19}}^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{1}}}-\frac{D{k}_{\mathrm{26}}^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{2}}})（12）

T(H)={\int }_{{\Omega }_{}}\frac{m_{\mathrm{1}}}{p_{\mathrm{1}}}\mathrm{d}{p}_{\mathrm{2}}\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{q}_{\mathrm{2}}=\frac{\mathrm{2}{\mathrm{\pi }}^{\mathrm{2}}{R}^{\mathrm{2}}\sqrt[]{m_{\mathrm{1}}{m}_{\mathrm{2}}}}{\sqrt[]{k_{\mathrm{11}}{k}_{\mathrm{21}}}}（13）

\frac{R^{\mathrm{2}}}{\mathrm{2}}-\frac{\sqrt[]{k_{\mathrm{12}}}{R}^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\phi \mathrm{s}\mathrm{i}\mathrm{n}\phi }{\sqrt[]{k_{\mathrm{11}}}}+\frac{k_{\mathrm{15}}{R}^{\mathrm{3}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{3}}\phi }{\mathrm{3}{k}_{\mathrm{11}}\sqrt[]{k_{\mathrm{11}}}}+

\frac{R\mathrm{c}\mathrm{o}\mathrm{s}\phi k_{\mathrm{18}}}{\sqrt[]{k_{\mathrm{11}}}}+\frac{R\mathrm{s}\mathrm{i}\mathrm{n}\phi k_{\mathrm{24}}}{\sqrt[]{k_{\mathrm{21}}}}-H=\mathrm{0}（14）

\Omega =\{(q_{\mathrm{1}},q_{\mathrm{2}},p_{\mathrm{2}})|H(q_{\mathrm{1}},q_{\mathrm{2}},\mathrm{0},p_{\mathrm{2}})\le H\}　　 \mathrm{0}\le \phi \le \mathrm{2}\mathrm{\pi }

式中：R为式（14）的最大正根。

2 受控磁悬浮系统随机最优控制

随机最优控制主要研究Markov扩散过程中的反馈控制，其目的是在控制方案中找到一个最优的控制方案，使受控的动力学系统或运动过程由某个初始状态转移到指定目标状态的同时性能指标最优。对于该磁悬浮系统，虽然悬浮控制器已采用PD控制来对电磁力进行控制，但是仍然可能存在系统不稳定而发生首次穿越，因此可考虑随机最优控制与PD控制相结合，将PD控制作为初始条件进行性能指标的最优控制。随机振动系统可通过可靠性函数和首次穿越时间等指标来对系统的可靠性进行描述^［23］。

2.1 可靠度最大化

定义磁悬浮系统的安全域为[\mathrm{0},H_{\mathrm{c}})，初始值为H_{\mathrm{0}}，将条件可靠性作为系统的性能指标，则值函数V(H,t)的表达式为

V(H,t)=\underset{u\in U}{\mathrm{s}\mathrm{u}\mathrm{p}}(P(H(s,u)\in [\mathrm{0},H_{\mathrm{c}}),

（15）

式中：U为所有控制过程的样本集；\tau为首次穿越时间；t_{\mathrm{f}}为终止时间，；sup（·）为上确界，表示从所有可能的控制中找到一个最优的使可靠性最大化的控制；P（·）为条件可靠性函数，表示系统停留在安全域的概率；s为瞬时时间。

将可靠度最大化作为目标的随机最优控制，其动态规划方程为

\frac{\partial V}{\partial t}=-\underset{u\in U}{\mathrm{s}\mathrm{u}\mathrm{p}}(\frac{\mathrm{1}}{\mathrm{2}}{\overline{\sigma }}^{\mathrm{2}}(H)\frac{{\partial }^{\mathrm{2}}V}{\partial H^{\mathrm{2}}}+(\overline{m}(H)+\mathrm{〈}u\frac{\partial H}{\partial p_{\mathrm{1}}}\mathrm{〉})\frac{\partial V}{\partial H})（16）

\mathrm{0}\le t\le t_{\mathrm{f}}H\in [\mathrm{0},H_{\mathrm{c}})

〈·〉=\frac{\mathrm{1}}{T(H)}{\int }_{\Omega }(·/\frac{\partial H}{\partial p_{\mathrm{1}}})\mathrm{d}{q}_{\mathrm{1}}\cdots \mathrm{d}{q}_n\mathrm{d}{p}_{\mathrm{2}}\cdots \mathrm{d}{p}_n

其边界条件为

V(H_{\mathrm{c}},t)=\mathrm{0}V(\mathrm{0},t)\in (\mathrm{0,1}]（17）

其终值条件为

V(H,t_{\mathrm{f}})=\mathrm{1}H\in [\mathrm{0},H_{\mathrm{c}})（18）

对于受控的磁悬浮拟不可积Hamilton系统，其平均Itô随机微分方程满足：

\mathrm{d}H=(\overline{m}(H)+\mathrm{〈}u\frac{\partial H}{\partial p_{\mathrm{1}}}\mathrm{〉})\mathrm{d}t+\overline{\sigma }(H)\mathrm{d}B(t)（19）

对式（16）等号右边取极大值，可得到此时的最优控制律为

u^{\mathrm{*}}=-b\mathrm{s}\mathrm{g}\mathrm{n}({\dot{q}}_{\mathrm{1}})（20）

此时最优控制律为一种开关控制，控制力幅值为b。将式（20）表征的最优控制律u^{\mathrm{*}}代替控制力，并代入到受控系统的平均Itô随机微分方程（式（19））中。为保持漂移系数和扩散系数的连续性，并避免在H\to \mathrm{0}时发生计算精度的溢出，将R^{\mathrm{2}}=\mathrm{2}H代入到式（11）和式（12）中，得到在H\to \mathrm{0}时漂移系数\overline{\overline{m}}(H)和扩散系数{\overline{\sigma }}^{\mathrm{2}}(H)的渐进表达式分别如下：

\overline{\overline{m}}(H)=\overline{m}(H)+\mathrm{〈}u\frac{\partial H}{\partial p_{\mathrm{1}}}\mathrm{〉}=

H(-\frac{c_{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{1}}}-\frac{c_{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{2}}}+\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{21}}{m}_{\mathrm{2}}}-\frac{k_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{m_{\mathrm{1}}})+

H^{\mathrm{2}}(\frac{k_{\mathrm{15}}{k}_{\mathrm{d}}{k}_{\mathrm{f}}{k}_{\mathrm{p}}}{\mathrm{3}{k}_{\mathrm{11}}^{\mathrm{2}}{m}_{\mathrm{1}}}-\frac{\mathrm{2}{k}_{\mathrm{15}}{k}_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{\mathrm{3}{c}_{\mathrm{0}}{k}_{\mathrm{11}}^{\mathrm{2}}{m}_{\mathrm{1}}})+\frac{D{k}_{\mathrm{19}}^{\mathrm{2}}}{m_{\mathrm{1}}}+\frac{D{k}_{\mathrm{26}}^{\mathrm{2}}}{m_{\mathrm{2}}}-\frac{\mathrm{4}\sqrt[]{\mathrm{2}}b}{\mathrm{3}\mathrm{\pi }\sqrt[]{m_{\mathrm{1}}}}{H}^{\frac{\mathrm{1}}{\mathrm{2}}}（21）

{\overline{\sigma }}^{\mathrm{2}}(H)=H^{\mathrm{2}}(\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{3}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{3}{k}_{\mathrm{21}}{m}_{\mathrm{2}}})+H(\frac{D{k}_{\mathrm{19}}^{\mathrm{2}}}{m_{\mathrm{1}}}+\frac{D{k}_{\mathrm{26}}^{\mathrm{2}}}{m_{\mathrm{2}}})（22）

此时受控系统的条件可靠性函数R_{\mathrm{c}}(t_{\mathrm{1}}|H_{\mathrm{0}})可通过差分法进行数值计算。将边界条件（式（17））和终值条件（式（18））设为已知，可求解式（16）所示的动态规划方程得到值函数V(H,t)。由于值函数V(H,t)与条件可靠性函数对两者进行时间上的变换，即t_{\mathrm{1}}=t_{\mathrm{f}}-t，因此两者存在如下关系：

R_{\mathrm{c}}(t_{\mathrm{f}}|H_{\mathrm{0}})=V(H_{\mathrm{0}},\mathrm{0})（23）

设置控制力幅值b={\mathrm{10}}^{\mathrm{4}}，绘制出受控和未受控时的条件可靠性函数如图2所示。对比图2a所示的仅受PD控制作用下的条件可靠性函数和图2b所示的受到PD控制和随机最优控制联合作用的条件可靠性函数发现，与不受控系统相比，添加随机最优控制后系统的条件可靠性下降速度更缓慢，即系统达到失效的时间延长，磁悬浮系统的使用寿命更久，从而提高了系统的可靠性。

2.2 平均首次穿越时间最长

随机振动系统首次越出安全域的时间即可表示为系统的寿命。在实践中，人们对平均首次穿越时间更为感兴趣，因此建立平均首次穿越时间最长的随机最优控制问题^［24］。取性能指标

J(u)=E(\tau (H,u))（24）

定义值函数

V_{\mathrm{1}}(H)=\underset{u\in U}{\mathrm{s}\mathrm{u}\mathrm{p}}(E(\tau (H,u)))（25）

此时的动态规划方程为

\underset{u\in U}{\mathrm{s}\mathrm{u}\mathrm{p}}(\frac{\mathrm{1}}{\mathrm{2}}{\overline{\sigma }}^{\mathrm{2}}(H)\frac{{\partial }^{\mathrm{2}}{V}_{\mathrm{1}}}{\partial H^{\mathrm{2}}}+(\overline{m}(H)+\mathrm{〈}u\frac{\partial H}{\partial p_{\mathrm{1}}}\mathrm{〉})\frac{\partial V_{\mathrm{1}}}{\partial H})=-\mathrm{1}（26）

其边界条件为

\left.\begin{array}{l}{V}_{\mathrm{1}}(H_{\mathrm{c}})=\mathrm{0}\\ V_{\mathrm{1}}(\mathrm{0})\in (\mathrm{0,1}]\end{array}\right\}（27）

类似地，对式（27）等号左侧取极大值，同样取最优控制律u^{\mathrm{*}}=-b\mathrm{s}\mathrm{g}\mathrm{n}({\dot{q}}_{\mathrm{1}})，形式与式（20）相同。通过Grank-Nicolson差分法和四阶Runge-Kutta分别结合边界条件（式（27））对动态规划方程（式（26））进行数值计算。设置控制力幅值b={\mathrm{10}}^{\mathrm{4}}，得到受控和未受控系统的平均首次穿越时间{\mu }_{\mathrm{1}}，如图3所示，其中用线条表示Grank-Nicolson差分法结果，用各图形表示四阶Runge-Kutta计算的数值结果。

由图3可知，对于受到PD控制和随机最优控制联合作用的磁悬浮系统，其平均首次穿越时间相比于仅受PD控制作用得到了延长。并且随着高斯白噪声强度D的增大，受控系统和未受控系统的平均首次穿越时间在逐渐缩短，系统的稳定性也在降低。采用差分法和四阶Runge-Kutta计算得到的结果几乎一致，从而验证了计算结果的准确性。

针对本文受控磁悬浮Hamilton系统的平均Itô随机微分方程（式19），对应推导出的Fokker-Planck-Kolmogorov（FPK）方程如下：

\frac{\partial p}{\partial t}=-\frac{\partial }{\partial H}(\overline{\overline{m}}(H)p(H))+\frac{\mathrm{1}}{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}}{\partial H^{\mathrm{2}}}({\sigma }^{\mathrm{2}}(H)p(H))（28）

式中：p(H)表示平稳概率密度。

求解FPK方程可得

p(H)=C\mathrm{e}\mathrm{x}\mathrm{p}(-\lambda (H))（29）

\begin{array}{c}\lambda (H)={\int }_{\mathrm{0}}^H(\frac{\mathrm{d}{\sigma }^{\mathrm{2}}(u)}{\mathrm{d}u}-\mathrm{2}m(u))/{\sigma }^{\mathrm{2}}(u)\mathrm{d}u\\ C=\frac{\mathrm{1}}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{e}\mathrm{x}\mathrm{p}(-\lambda (H))\mathrm{d}H}\end{array}

其中，C为归一化常数。将H代入式（33）中，得到磁悬浮系统车体的平稳联合概率密度p(p_{\mathrm{2}},q_{\mathrm{2}})。

通过数值模拟计算，得到图4所示的受控和未受控磁悬浮系统联合概率密度图。设置控制力幅值b=\mathrm{0}、高斯白噪声强度D=\mathrm{0.02}时，得到仅受PD控制作用的联合概率密度。此时联合概率密度呈现“火山口”形状，磁悬浮系统以较大概率收敛于极限环位置，磁悬浮系统的平凡解概率不稳定。设置控制力幅值b={\mathrm{10}}^{\mathrm{4}}时，此时磁悬浮系统受到PD控制和随机最优控制联合作用。磁悬浮系统的联合概率密度从“火山口”形状变为具有极大值的“单峰”状，此时磁悬浮系统以较大概率稳定于p_{\mathrm{2}}=\mathrm{0}\mathrm{、}{q}_{\mathrm{2}}=\mathrm{0}，磁悬浮系统的平凡解概率渐近稳定。

3 随机稳定化

在随机参数激励的作用下，随机振动系统往往不稳定，即随机振动系统的解无法收敛于平衡位置或平稳状态。随机稳定化，即通过反馈控制将原本不稳定的系统变成稳定的，或提高原系统的稳定性。本文采用基于Hamilton系统遍历控制及最大Lyapunov指数的随机稳定化控制策略^［25-26］，去掉外激k_{\mathrm{19}}=k_{\mathrm{26}}=\mathrm{0}，考虑式（19）所示的遍历控制，以最大Lyapunov指数最小为准则，取性能指标如下：

J=\underset{T\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{1}}{T}{\int }_{\mathrm{0}}^T(f_{\mathrm{1}}(H(s))+\mathrm{〈}{\mathit{u}}^{\mathrm{T}}(s)\mathit{R}\mathit{u}(s)\mathrm{〉})\mathrm{d}s（30）

式中： u 为控制力矩阵；f_{\mathrm{1}}(H)表示待定的函数；\mathit{R}为权重矩阵；T为时间参数，用于定义计算平均性能的时间窗口长度。

设置最优控制律为

u^{\mathrm{*}}=-\frac{\mathrm{1}}{\mathrm{2}R}\frac{\mathrm{d}V}{\mathrm{d}H}{\dot{q}}_{\mathrm{1}}（31）

式（31）满足条件\mathrm{d}V/\mathrm{d}H>\mathrm{0}，则动态规划方程可表示为

\gamma =\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}(H)\frac{{\mathrm{d}}^{\mathrm{2}}V}{\mathrm{d}{H}^{\mathrm{2}}}+\overline{m}(H)\frac{\mathrm{d}V}{\mathrm{d}H}-

\frac{\mathrm{1}}{\mathrm{4}{R}_{\mathrm{c}}}\mathrm{〈}{(\frac{\partial H}{\partial p_i})}^{\mathrm{2}}\mathrm{〉}{(\frac{\mathrm{d}V}{\mathrm{d}H})}^{\mathrm{2}}+f_{\mathrm{1}}(H)（32）

式中：γ为最优平均成本；p_i（i=1，2）为广义动量；R_c为权衡控制效果与控制能量耗散的权重系数。

根据式（19）支配的磁悬浮系统平均Itô随机微分方程，将R²=2H代入式（31）得到在H→0时的渐进表达式为

\mathrm{〈}{u}^{\mathrm{*}}\frac{\partial H}{\partial P_{\mathrm{1}}}\mathrm{〉}=-\frac{\mathrm{1}}{\mathrm{2}{R}_{\mathrm{c}}{m}_{\mathrm{1}}}\frac{\mathrm{d}V}{\mathrm{d}H}\frac{\mathrm{1}}{T(H)}{\int }_{\Omega }^{}(p_{\mathrm{1}})\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{q}_{\mathrm{2}}\mathrm{d}{p}_{\mathrm{2}}=

-\frac{\mathrm{1}}{\mathrm{2}{R}_{\mathrm{c}}{m}_{\mathrm{1}}}\frac{\mathrm{d}V}{\mathrm{d}H}\frac{H}{\mathrm{2}}（33）

对待定函数f₁（H）作如下规定：

f_{\mathrm{1}}(H)=KH+\gamma \frac{\mathrm{d}V}{\mathrm{d}H}=C_{\mathrm{c}}（34）

式中：K为任意常数；C_{\mathrm{c}}为控制系数。

将式（34）代入到式（32）所示的动态规划方程中，得到受最优控制律u^*作用的最终动态规划方程为

-\frac{\mathrm{1}}{\mathrm{8}{R}_{\mathrm{c}}{m}_{\mathrm{1}}}H{C}_{\mathrm{c}}^{\mathrm{2}}+{\overline{m}}_{\mathrm{0}}H{C}_{\mathrm{c}}+KH=\mathrm{0}（35）

对式（35）进行求解，可得

C_{\mathrm{c}}=\mathrm{4}{R}_{\mathrm{c}}{m}_{\mathrm{1}}({\overline{m}}_{\mathrm{0}}+\sqrt[]{{\overline{m}}_{\mathrm{0}}^{\mathrm{2}}+\frac{K}{\mathrm{2}{R}_{\mathrm{c}}{m}_{\mathrm{1}}}})（36）

忽略漂移系数和扩散系数的高阶项，得到未受控磁悬浮系统在H→0时漂移系数和扩散系数的渐进表达式为

\left.\begin{array}{l}{\overline{m}}_{\mathrm{0}}=-\frac{c_{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{1}}}-\frac{c_{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{2}}}+\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{21}}{m}_{\mathrm{2}}}-\frac{k_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{m_{\mathrm{1}}}\\ {\overline{\sigma }}_{\mathrm{0}}^{\mathrm{2}}=\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{3}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{3}{k}_{\mathrm{21}}{m}_{\mathrm{2}}}\end{array}\right\}（37）

经过随机稳定化后，扩散系数的表达式与式（37）一致，漂移系数表达式如下：

{\overline{\overline{m}}}_{\mathrm{0}}=-\frac{c_{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{1}}}-\frac{c_{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{2}}}+\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{2}{k}_{\mathrm{21}}{m}_{\mathrm{2}}}-

\frac{k_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{m_{\mathrm{1}}}-\frac{C_{\mathrm{c}}}{\mathrm{4}{R}_{\mathrm{c}}{m}_{\mathrm{1}}}（38）

根据式（37）和式（38），得到未受最优控制和受最优控制系统的最大Lyapunov指数分别为{\lambda }_{\mathrm{u}}和{\lambda }_{\mathrm{c}}，其表达式如下：

\left.\begin{array}{l}{\lambda }_{\mathrm{c}}=-\frac{\sqrt[]{{\overline{\overline{m}}}_{\mathrm{0}}^{\mathrm{2}}+\frac{K}{\mathrm{2}{R}_{\mathrm{c}}{m}_{\mathrm{1}}}}}{\mathrm{2}}-\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{12}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}-\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{12}{k}_{\mathrm{21}}{m}_{\mathrm{2}}}\\ {\lambda }_{\mathrm{u}}=-\frac{c_{\mathrm{2}}}{\mathrm{4}{m}_{\mathrm{1}}}-\frac{c_{\mathrm{2}}}{\mathrm{4}{m}_{\mathrm{2}}}+\frac{D{k}_{\mathrm{18}}^{\mathrm{2}}}{\mathrm{6}{k}_{\mathrm{11}}{m}_{\mathrm{1}}}+\frac{D{k}_{\mathrm{25}}^{\mathrm{2}}}{\mathrm{6}{k}_{\mathrm{21}}{m}_{\mathrm{2}}}-\frac{k_{\mathrm{d}}{k}_{\mathrm{f}}{I}_{\mathrm{0}}}{\mathrm{2}{m}_{\mathrm{1}}}\end{array}\right\}（39）

在受控和未受控系统中通过各自最大Lyapunov指数的正负即可以判断系统的局部稳定性。当Lyapunov指数大于0时，该系统概率为1渐进不稳定；当Lyapunov指数小于0时，该系统概率为1渐进稳定。为研究受控和未受控磁悬浮系统的稳定性，设置控制参数K=1000、R_{\mathrm{c}}=\mathrm{0.01}。根据式（39）进行MATLAB数值模拟计算，得到的结果如图5所示。

由图5可知，仅受PD控制作用时，随着高斯白噪声强度D的增大，最大Lyapunov指数从逐渐变为{\lambda }_{\mathrm{u}}>\mathrm{0}，即此时磁悬浮系统的平凡解从概率为1渐进稳定变为概率为1渐进不稳定。当{\lambda }_{\mathrm{u}}>\mathrm{0}时，磁悬浮系统的总能量不会耗散，会随着时间的增加逐渐增加。经过随机稳定化后的随机最优控制磁悬浮系统的最大Lyapunov指数始终满足，即磁悬浮系统的平凡解变为概率为1渐进稳定，系统的总能量会逐渐耗散，最终收敛于平衡位置。

为研究经过随机稳定化后磁悬浮系统联合概率密度的变化，根据式（29）进行数值模拟，受控和未受控磁悬浮系统联合概率密度如图6所示。相比于图6a所示的仅受PD控制的磁悬浮系统，经过随机稳定化后联合概率密度从“火山口”形状变为“单峰”状，磁悬浮系统的定性性态发生了变化。由图6b可知，此时磁悬浮系统的平凡解以较大概率稳定于p_{\mathrm{2}}=\mathrm{0}\mathrm{、}{q}_{\mathrm{2}}=\mathrm{0}，但存在较小的概率收敛于其他位置。

经过以最大Lyapunov指数最小为目标的随机最优控制后，磁悬浮系统的稳定性得到了提高。通过选取合理的参数K、R_c，可使未受控系统的最大Lyapunov指数减小，从而保证磁悬浮系统满足局部概率为1渐进稳定。

4 最优控制策略性能评价指标

为评价最优控制策略，判断最优控制策略的性能指标，引入如下两个准则^［27］：

K_{\mathrm{h}}=\frac{\sqrt[]{E_{\mathrm{u}}(q_i^{\mathrm{2}})}-\sqrt[]{E_{\mathrm{c}}(q_i^{\mathrm{2}})}}{\sqrt[]{E_{\mathrm{u}}(q_i^{\mathrm{2}})}}（40）

{\mu }_{\mathrm{h}}=\frac{K_{\mathrm{h}}}{\sqrt[]{E(u_i^{\mathrm{2}})/(\mathrm{2}D)}}（41）

\begin{array}{l}{E}_{\mathrm{u}}(q_i^{\mathrm{2}})={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{q}_i^{\mathrm{2}}{p}_{\mathrm{u}}(p_i,q_i)\mathrm{d}{p}_{\mathrm{1}}\mathrm{d}{p}_{\mathrm{2}}\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{q}_{\mathrm{2}}\\ E_{\mathrm{c}}(q_i^{\mathrm{2}})={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{q}_i^{\mathrm{2}}{p}_{\mathrm{c}}(p_i,q_i)\mathrm{d}{p}_{\mathrm{1}}\mathrm{d}{p}_{\mathrm{2}}\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{q}_{\mathrm{2}}\\ E(u_i^{\mathrm{2}})={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{u}_i^{\mathrm{2}}{p}_{\mathrm{c}}(p_i,q_i)\mathrm{d}{p}_{\mathrm{1}}\mathrm{d}{p}_{\mathrm{2}}\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{q}_{\mathrm{2}}\end{array}

式中：K_h为最优控制策略的控制效果值，表示受控系统的均方位移减少量百分比；{\mu }_{\mathrm{h}}为控制效率，表示单位均方控制力下能提供的位移减少量百分比；E_{\mathrm{u}}(q_i^{\mathrm{2}})(i=\mathrm{1,2})为未受控系统的位移均方值；E_{\mathrm{c}}(q_i^{\mathrm{2}})为受控系统的位移均方值；E(u_i^{\mathrm{2}})为控制力的均方值；p_{\mathrm{u}}为不受控系统的广义动量；p_{\mathrm{c}}为受控系统的广义动量。

由于采用可靠度最大化和平均首次穿越时间最长的最优控制策略时，其最优控制律在形式上是等效的，因此针对这两种最优控制策略的性能指标进行同时评价，经过MATLAB的数值模拟结果如图7所示。利用随机稳定化控制策略，以最大Lyapunov指数最小化为优化目标，可得到控制效果和控制效率在不同高斯白噪声强度下的数值结果如图8所示。

对于控制效果K_h和控制效率{\mu }_{\mathrm{h}}，其数值越大，说明控制策略越优。由图7可知，随着高斯白噪声强度D的增大，控制效果K_h减弱，控制效率{\mu }_{\mathrm{h}}逐渐降低，但高斯白噪声强度较低时有较好的控制效果K_h和较高的控制效率{\mu }_{\mathrm{h}}，此时说明控制策略较优。由图8可知，当高斯白噪声强度较低时，控制效果K_h和控制效率{\mu }_{\mathrm{h}}的数值都较小。随着高斯白噪声强度的提高，控制效果K_h在逐渐增强，控制效率{\mu }_{\mathrm{h}}呈现先提高后降低的趋势。由此可知，采用最大Lyapunov指数最小的控制策略在较大的高斯白噪声强度下才具有较好的控制效果，并且控制效率在一定范围内较高。

5 结论

本文针对磁悬浮系统的最优控制问题，基于Hamilton理论的随机平均法建立磁悬浮系统的非线性随机微分方程。建立以可靠度最大化、平均首次穿越时间最长和最大Lyapunov指数最小为控制目标的最优控制策略。主要结论如下：

1）相比于仅受PD控制作用的磁悬浮系统，在PD控制基础上添加最优控制能在一定程度上针对最优控制策略进行优化。选择合理的控制参数能够改善系统的稳定性，提高控制策略的性能。

2）以可靠度最大化为目标和以平均首次穿越时间最长为目标的最优控制策略，两者最优控制律的形式是相同的。添加最优控制后，提高了磁悬浮的可靠性，增加了平均首次穿越时间，即延长了磁悬浮系统的使用寿命。且这两种控制策略在高斯白噪声较低时控制效果均较好、控制效率均较高。

3）以最大Lyapunov指数最小化为目标的最优控制策略使磁悬浮系统的平凡解始终能够满足概率为1渐进稳定的条件。并且磁悬浮系统的联合概率密度从“火山口”形状变为“单峰”状，使磁悬浮系统的平凡解从概率不稳定变为概率稳定，提高了磁悬浮系统的稳定性。但该控制策略对高斯白噪声强度敏感性过大，仅能保证在一定范围内的高斯白噪声强度下系统的控制效果和控制效率较好。

通过对磁悬浮系统进行最优控制，避免磁悬浮列车在运行过程中发生首次穿越失效，否则列车有可能发生跳轨，严重时甚至可能发生导轨下表面与悬浮电磁铁上表面碰撞的现象，从而影响乘坐的舒适性和安全性。本文研究结果对磁悬浮列车的控制策略优化具有理论指导意义。

参考文献

[1] 熊嘉阳， 邓自刚. 高速磁悬浮轨道交通研究进展［J］. 交通运输工程学报， 2021， 21（1）： 177-198.

[2] XIONG Jiayang， DENG Zigang. Research Progress of High-speed Maglev Rail Transit［J］. Journal of Traffic and Transportation Engineering， 2021， 21（01）： 177-198.

[3] 唐文冰， 肖立业， 王粟， 等. 磁悬浮轨道交通中的磁悬浮导向方式研究综述［J］. 电工电能新技术， 2022， 41（5）： 45-60.

[4] TANG Wenbing， XIAO Liye， WANG Su， et al. Summary of Research on Levitation-guidance Modes in maglev Rail Transportation Technology［J］. Advanced Technology of Electrical Engineering and Energy， 2022， 41（5）： 45-60.

[5] LI J， LI J. A Practical Nonlinear Controller for Levitation System with Magnetic Flux Feedback［J］. Journal of Central South University， 2016， 23（7）： 1729-1739.

[6] XU Y， LONG Z， ZHAO Z， et al. Real-time Stabi-lity Performance Monitoring and Evaluation of Maglev Trains' Levitation System： a Data-driven Appro-ach［J］. IEEE Transactions on Intelligent Transportation Systems， 2022， 23（3）： 1912-1923.

[7] CHEN X， MA W， LUO S， et al. Study on Lateral Stability of Levitation Modules for Low and Medium-speed Maglev Trains［J］. Archive of Applied Mechanics， 2020， 90（2）： 437-447.

[8] YU X， WANG Y， BO K， et al. Research on Parameters Optimum Design of Levitation Electromagnet Based on Orthogonal Design for EMS Maglev Train［J］. Lecture Notes in Electrical Engineering， 2022， 890： 50-58.

[9] 薄凯， 陈俊全， 王东， 等. 常导型高速磁浮列车电磁系统参数优化设计［J］. 电机与控制学报， 2022， 26（11）： 49-57.

[10] BO Kai， CHEN Junquan， WANG Dong， et al. Research on Parameters Optimum Design Electrimagnet System for EMS Maglev Train［J］. Electric Machines and Control， 2022， 26（11）： 49-57.

[11] CHANG S， LUE Y. A Study of the Nonlinear Response and Chaos Suppression in a Magnetically Levitated System［J］. Australian Journal of Mechanical Engineering， 2020， 18（1）： 94-105.

[12] CHEN X， MA W， LUO S. Study on Stability and Bifurcation of Electromagnet-track Beam Coupling System for EMS Maglev Vehicle［J］. Nonlinear Dynamics， 2020， 101（4）： 2181-2193.

[13] FENG Y， ZHAO C， WU D， et al. Effect of Levitation Gap Feedback Time Delay on the EMS Maglev Vehicle System Dynamic Response［J］. Nonlinear Dynamics， 2023， 111（8）： 7137-7156.

[14] SUN Y， XU J， QIANG H， et al. Hopf Bifurcation Analysis of Maglev Vehicle-guideway Interaction Vibration System and Stability Control Based on Fuzzy Adaptive Theory［J］. Computers in Industry， 2019， 108： 197-209.

[15] SUN Y， HE Z， XU J， et al. Dynamic Analysis and Vibration Control for a Maglev Vehicle-guideway Coupling System with Experimental Verification［J］. Mechanical Systems and Signal Processing， 2023， 188： 109954.

[16] HU W， ZHOU Y， ZHANG Z， et al. Model Predictive Control for Hybrid Levitation Systems of Maglev Trains with State Constraints［J］. IEEE Transactions on Vehicular Technology， 2021， 70（10）： 9972-9985.

[17] XU J， CHEN C， SUN Y， et al. Nonlinear Dynamic Analysis of Maglev Vehicle Based on Flexible Guideway and Random Irregularity［J］. International Journal of Applied Electromagnetics and Mechanics， 2019， 60（2）： 263-280.

[18] 刘伟渭，戴焕云，曾京.弹性约束轮对系统的随机稳定性研究［J］.中国机械工程，2013，24（6）：799-804.

[19] LIU Weiwei， DAI Huanyun， ZENG Jing. Research on Stochastic Stability of Elastic Constraint Wheelset System［J］. China Mechanical Engineering， 2013，24（6）：799-804.

[20] 吴晗， 曾晓辉. 气动升力下磁浮车辆非线性响应研究［J］. 机械工程学报， 2021， 57（14）： 223-231.

[21] WU Han， ZENG Xiaohui. Nonlinear Dynamics of Malev Vehicle under Aerodynamic Lift［J］. Journal of Mechanical Engineering， 2021， 57（14）： 223-231.

[22] WU H， ZENG X H， GAO D G， et al. Dynamic Stability of an Electromagnetic Suspension Maglev Vehicle under Steady Aerodynamic Load［J］. Applied Mathematical Modelling， 2021， 97： 483-500.

[23] 陈志贤， 李忠继， 杨吉忠， 等. 常导高速电磁悬浮车辆二系悬挂结构对比优化［J］. 机械工程学报， 2022， 58（10）： 160-168.

[24] CHEN Zhixian， LI Zhongji， YANG Jizhong， et al. Comparison and Optimization of Secondary Suspension Structure of High Speed EMS Vehicle［J］. Journal of Mechanical Engineering，2022， 58（10）： 160- 168.

[25] 胡帛茹， 赵春发， 蔡文锋， 等. 两种中低速磁浮车辆动力学性能仿真对比分析［J］. 机械工程学报， 2023， 59（10）： 311-322.

[26] HU Boru， ZHAO Chunfa， CAI Wenfeng， et al. Numerical Simulation and Comparative Analysis on Dynamic Performances of Two Types of Medium-low Speed Maglev Vehicles［J］. Journal of Mechanical Engineering， 2023， 59（10）： 311-322.

[27] 陈琛， 徐俊起， 荣立军， 等. 轨道随机不平顺下磁浮车辆非线性动力学特性［J］. 交通运输工程学报， 2019， 19（4）： 115-124.

[28] CHEN Chen， XU Junqi， RONG Lijun， et al. Nonlinear Dynamics Characteristics of Maglev Vehicle under Track Random Irregularities［J］.Journal of Traffic and Transportation Engineering， 2019， 19（4）： 115-124.

[29] SUN J， XIA L， YING Z， et al. Reliability of Nonlinear Stochastic Controlled Systems Considering the Dynamics of Sensors and Actuators［J］. Journal of Vibration and Control， 2022， 28（15/16）：2052-2060.

[30] SUN J， ZHU W， JIANG W， et al. Reliability of a Class of Nonlinear Systems under Switching Random Excitations［J］. Nonlinear Dynamics， 2020， 99（3）： 2083-2094.

[31] 刘伟渭， 姜瑞金， 刘凤伟， 等. 约束轮对首次穿越失效的随机非线性最优控制［J］. 铁道学报， 2018， 40（6）： 30-35.

[32] LIU Weiwei， JIANG Ruijin， LIU Fengwei， et al. Stochastic Nonlinear Optimal Control for First-passage Failure of Constrained Wheelset［J］. Journal of the China Railway Society， 2018， 40（6）： 30-35.

[33] 朱位秋， 应祖光. 拟哈密顿系统非线性随机最优控制［J］. 力学进展， 2013， 43（1）： 39-55.

[34] ZHU Weiqiu， YING Zuguang. Advances in Research on Nonlinear Stochastic Optimal Control of Quasi-Hamiltonian Systems［J］. Advances in Mechanics， 2013， 43（1）： 39-55.

[35] LIU W， YIN X， GUO Z， et al. Feedback Stabilization of quasi Nonintegrable Hamiltonian Systems under Combined Gaussian and Poisson White Noise Excitations［J］. Probabilistic Engineering Mechanics， 2023， 71： 103-107.

[36] WANG Y， YING Z， ZHU W. Stochastic Minimax Control for Stabilizing Uncertain Quasi-integrable Hamiltonian Systems［J］. Automatica， 2009， 45（8）： 1847-1853.

[37] HUAN R， WU Y， ZHU W. Stochastic Optimal Bounded Control of MDOF Quasi Nonintegrable-Hamiltonian Systems with Actuator Saturation［J］. Archive of Applied Mechanics， 2009， 79（2）： 157-168.

基金资助

国家自然科学基金(51875481)

四川省科技计划重点研发项目(2023YFN0079)