Method of Determining Design Value of Rockfall Block Size Based on the Return Period
-
摘要:
山区交通落石事件频发且具有突发性和随机性，基于概率的防护结构可靠性设计是降低灾害风险的必要和有效手段，然而，作为设计基本参数的落石尺寸设计值的确定方法还不明确，使得可靠性设计不能有效进行. 介绍基于泊松-广义帕累托分布（GPD）复合模型的落石尺寸-重现期预测方法，结合宝成铁路白—上段落石记录数据，进行GPD阈值选取，并分析记录中缺失的最大、最小落石尺寸对阈值选取和对预测结果的影响；与地震作用标准值超越概率对应，提出以“小落石不坏、中落石可修、大落石不塌”三水准设防为标准的落石尺寸设计值确定方法. 研究表明：记录中缺失的最大落石尺寸对阈值选取影响小，但对落石尺寸预测结果影响大，与重现期对应的落石尺寸预测值随设定的最大落石尺寸呈指数型增大趋势，而记录的最小落石尺寸对阈值选取和预测结果的影响都较小；以宝成铁路白—上段落石记录数据为例，根据选取出的合理阈值，设定不同最大落石尺寸值，利用本文方法，预测得到宝成铁路白—上段50、100、
1000 、10000 a重现期落石尺寸区间范围，计算得到50、100 a设计基准期三水准设防标准对应的落石尺寸标准值即设计值，可为类似工程设计提供参考.Abstract:The rockfall events on transport lines in mountainous areas occur frequently and are sudden and random. Probability-based reliability design of protective structures is a necessary and effective means to reduce the hazard risk. However, the method of determining the design value of rockfall size, as a basic parameter of the design, is still unclear, making the reliability design fail to be carried out effectively. The prediction method of rockfall size and the return period was introduced based on the Poisson-generalized Pareto distribution (GPD) composite model. According to the rockfall record data of the Bai−Upper section of Baoji−Chengdu Railway, the GPD threshold was selected, and the impact of the missing maximum and minimum rockfall sizes (vmax and vmin) in the record on the threshold selection and prediction results was analyzed. Corresponding to the exceedance probability of the standard value of earthquake action, a method of determining the design value of rockfall size was proposed, which was based on the three-level protection standard featuring “the structure doesn’t break under small rockfalls, can be repaired under medium rockfalls, and doesn’t collapse under large rockfalls”. The results show that the missing maximum rockfall size in the record has a small impact on the threshold selection but a large impact on the prediction results of the rockfall size. The prediction value of rockfall size corresponding to the return period increases exponentially with the set maximum rockfall size, while the recorded minimum rockfall size has less impact on both the threshold selection and the prediction results. By taking the record data of the Bai−Upper section of Baoji−Chengdu Railway as an example, according to the selected reasonable threshold, different maximum rockfall size values are set. By using the method presented in this paper, the interval range of rockfall size in return periods of 50, 100, 1 000, and 10 000 years on the Bai−Upper section of Baoji−Chengdu Railway is predicted, and the standard value (design value) of rockfall size corresponding to the three-level protection standard in design reference periods of 50 and 100 years is calculated, which can provide a reference for similar engineering designs.
-
Key words:
- protective structure /
- rockfall size /
- return period /
- generalized Pareto distribution /
- design value
-
随着公路交通运输业的发展，传统汽车行业面临新的挑战[1]. 磁悬浮汽车将磁悬浮技术与传统汽车相结合，有望成为未来公路交通运输的发展趋势. 在众多磁浮制式[2-4]中，永磁电动悬浮系统凭借结构简单、系统可靠、悬浮间隙大等优势受到众多学者青睐[5]. 其中，环形Halbach永磁轮能够将阻碍系统运行的磁阻力转化为前进的驱动力，具有悬浮驱动一体化的特点[6]. 自20世纪末，许多学者已将其作为磁悬浮汽车的关键技术，在电磁力解析模型和结构参数优化[7]等方面展开研究.
永磁电动悬浮汽车具有横向欠阻尼的特点[8]，对其进行横向力补偿的横向运行模式研究至关重要. 环形Halbach永磁轮分为轴向轮与径向轮，2000年，Fujii等[9-10]用4个轴向轮构成悬浮系统，实现稳定悬浮，并通过倾斜和部分重叠布置产生横向力以实现横向调节功能. 2016年，Jung[11]也采用轴向轮部分重叠布置方式，设计了无接触运输系统和四轮悬浮驱动样车. 2017年，Jung[12]利用径向轮设计了四轮磁悬浮样机，通过与导体板轨缝作用产生横向恢复力以抵抗横向扰动，实现导向功能. 2023年，刘新等[13]通过将径向轮与导体板倾斜布置，将法向力分解为悬浮力与横向力来补偿系统的横向力. 以上研究均是通过磁轮与导体板间特殊放置被动地产生横向力，只能沿固定轨道行驶，不符合自由灵活的公路运输特点. 2021年，张泽等[14]设计了可在无轨道导体板上主动转向的四轮原理样机模型，并针对垂向和纵向运行模式展开研究，但未考虑横向方向，对横向扰动较为敏感，且在过弯时无法主动调节横向力以克服离心力，只能沿直线行驶.
永磁电动悬浮汽车横向开环系统极不稳定，需要设计合适的控制系统，但相较于较为成熟的电磁悬浮（EMS）制式控制研究[15]，永磁电动悬浮系统（PMEDS）的控制研究仍处于初步探索阶段. 2006年，de Boeij等[16]以电动磁悬浮列车为研究对象，设计一种基于五自由度状态空间模型的具有积分作用的一阶滑模控制器，并通过实验证明该控制器是全局稳定的. 2019年，Wright[8]建立永磁电动悬浮汽车四自由度动力学模型，采用线性二次型调节器、线性二次积分控制等方法进行对比分析，验证了控制器的有效性，但四自由度中没有横向自由度. 2023年，Zhang等[17] 针对磁悬浮汽车偏航自由度控制易超调与抗扰性差问题，设计了双环自抗扰控制器，并通过仿真和实验验证了该控制器的轨迹跟踪性能与抗扰能力. 2024年，Li等[18]分析了永磁电动悬浮汽车浮沉自由度动力学模型并引入温度干扰器，设计了基于径向基函数神经网络滑模控制（RBF-SMC）策略有效抑制导体板温升效应引起的垂向扰动，但并未考虑温升对系统磁阻力的影响，横、纵与偏航方向的温升扰动无法得到有效抑制，难以实现原地静浮. 尽管永磁电动悬浮汽车在其余5个自由度的控制研究方面已有进展，但横向运动控制的稳定性研究仍处于空白阶段，因此本文选择横向自由度展开控制研究具有重要意义.
为实现永磁电动悬浮汽车横向运动精准控制并提高横向抗扰能力，本文首先提出基于同侧磁轮对称偏转的横向运行模式，建立横向单自由度非线性动力学模型；然后，针对系统非线性特点设计非线模型预测控制器（NMPC）实现系统横向轨迹跟踪控制，并针对模型不精确与存在外部扰动的问题引入扩张状态观测器（ESO），针对ESO对传感器测量噪声较敏感的问题，引入扩展卡尔曼滤波器（EKF）对传感器原始数据滤波处理，建立NMPC-ESO-EKF控制器；最后，通过MATLAB/Simulink与Simscape联合仿真及原理样机实验验证运行模式与控制器的有效性.
1. 横向运行原理及动力学建模
1.1 车辆横向运行原理
本文中永磁电动悬浮汽车的车载永磁轮采用径向环形Halbach阵列结构，结构示意如图1所示. 图中：Ro、Rw分别为永磁轮的外半径和内半径，h为悬浮间隙，s为导体板厚度， vx为汽车速度. 该永磁轮在导体板上方以转速n高速旋转时，其产生的磁场与导体板在有限空间内产生相对运动，导体板中磁通量发生变化，产生感应电动势，从而形成较强的涡流. 涡流产生的感应磁场与永磁轮源磁场相互作用，使永磁轮产生垂直向上的悬浮力Fl和阻碍其自身旋转的径向磁阻力Fr. Bird等[19-21]已经建立了该悬浮力与磁阻力的二维、三维稳态解析模型，悬浮力用来克服自身重力，而磁阻力可以类比于传统汽车轮胎与路面的摩擦力，推动永磁电动悬浮汽车前行，实现系统悬浮驱动一体化.
将车载永磁轮沿水平方向可控地偏转一定角度α，如图2所示，永磁轮的径向磁阻力便被分解到x、y方向上，表现为横向力和纵向驱动力，从而实现横向力的主动补偿，提高系统横向抗扰能力. 驱动力Fd、横向力Fg与径向磁阻力Fr之间的关系为
图 2 常规永磁轮与偏转永磁轮示意Figure 2. Conventional permanent magnet wheel and deflecting permanent magnet wheel{Fd=Frcosα,Fg=Frsinα. (1) 图3为磁轮可偏转的永磁电动悬浮汽车原理样机结构示意，主要是由4个无刷直流电机带动4个磁轮在导体板上方高速旋转，由4个舵机带动磁轮沿水平方向偏转，通过控制4个磁轮的转速与偏转角度调节整车的悬浮力、驱动力与横向力，实现车辆浮沉、前进后退、左右横移、偏航、点头与侧滚的全维运动，从而调整车辆的位置与姿态，本文针对其中的横向运动控制展开研究.
永磁电动悬浮汽车横向运行模式原理如图4所示. 假设永磁电动悬浮汽车为刚体且结构对称，其重心在几何中心处. nj（j=1,2,3,4）为j号磁轮的转速，αj（j=2,3）为2、3号磁轮偏转角度，Flj（j=1,2,3,4）表示j号磁轮的悬浮力，Frj（j=1,2,3,4）表示j号磁轮的磁阻力. 选取逆时针和向前运动为正方向，图中前侧1、2号磁轮逆时针旋转，后侧3、4号磁轮顺时针旋转，各磁轮转速大小始终保持一致，则各磁轮磁阻力大小相同，即Fr1Fr2=−Fr3Fr4=Fr0. 将左侧2、3号磁轮对称偏转相同角度，即α2=−α3=α0, 可以将其磁阻力分解产生横向力使车辆横移，此时，其他方向的合力与合力矩均为0. 保持各磁轮转速不变，通过控制2、3号磁轮偏转角度的大小与方向便可控制车辆的横向运动，采用该横向运动模式可以将系统简化为一个单输入单输出系统.
1.2 车辆横向动力学建模
准确的永磁电动悬浮汽车横向动力学模型是设计控制器的基础. 为避免其他自由度变化时给磁轮电磁力特性带来影响，限制其他自由度的运动，以简化磁悬浮汽车横向动力学模型. 永磁轮横向运动时会始终产生阻碍其运行的横向磁阻力，其大小与横向平移速度$ \dot {\boldsymbol{x}} $有关，用横向阻尼力将其简化代替，假设车辆横向跑动时横向阻尼系数c处处相等，同时忽略磁轮偏转引起的横向动量变化，选取向左运动为正方向，则车辆受到的总横向力Fgx为
Fgx=Fr2sinα2+Fr3sinα3−c˙x. (2) 根据上文可知，2、3号磁轮转速与偏转角度大小相同，则式（1）简化为
Fgx=2Fr0sinα0−c˙x. (3) 根据控制器设计的需要，本文利用牛顿定律得到永磁电动悬浮汽车横向自由度的非线性动力学模型，如式（4）所示.
¨x=2Fr0msinα0−cm˙x, (4) 式中：$ \ddot {\boldsymbol{x}} $为永磁电动悬浮汽车的横向加速度，m为永磁电动悬浮汽车的总质量.
2. NMPC-ESO-EKF控制器设计
NMPC控制的基本思想是利用当前时刻系统的状态及约束条件，对未来一段时间内的状态、输入变量进行预测，并求解出1组最优的控制输入序列，但只选取该序列中的第1组结果作用于系统中，如此循环直到系统达到期望状态，主要包括预测模型、滚动优化和反馈校正三部分[22]. 为解决控制中模型不精确和存在外部干扰的问题，引入ESO对系统扰动量进行实时观测与反馈补偿，采用EKF对传感器测量数据进行滤波降噪，提高观测精度，降低环境噪声对ESO的影响.
2.1 NMPC控制器设计
2.1.1 预测模型
预测模型可根据当前时刻的状态量和控制量迭代地预测出系统未来的状态. 针对永磁电动悬浮汽车非线性的特点，结合式（4），本文采用非线性离散状态空间方程作为系统的预测模型.
定义系统的状态量为X[k]=[x1 x2]T，其中，x1＝x，x2＝$ \dot {\boldsymbol{x}} $；系统的输入量为U[k]=u， 其中，u=α；系统的输出量为Y[k]. 则永磁电动悬浮汽车的非线性离散状态空间模型为
{X[k+1]=f(X[k],U[k]),Y[k]=X[k], (5) 式中：$ f({\boldsymbol{X}}_{[k]},U_{[k]}) $为非线性离散映射函数，如式（6）所示.
f(X[k],U[k])=[x1+x2Tcx2+(2Fr0msinu−cmx2)Tc], (6) 式中：Tc模型预测离散步长.
2.1.2 滚动优化
为保证永磁电动悬浮汽车能更快速与稳定地达到横向参考期望值，本文采用二次型函数设计代价函数，使其最小化来求解最优控制序列U*，在预测区间Np和控制区间Nc内，定义代价函数J如式（7）所示.
J=Np∑i=1‖Xpre[k+i]−Xref[k+i]‖2Q+Nc−1∑i=0‖U[k+i]‖2R, (7) 式中：Q为运行状态权重矩阵，R为输入权重，Q∈R2 × 2，R∈R1 × 1；$ \|{\text{•}}\|_{s}^{2} $为欧几里得范数，S∈{Q，R}；Xpre[k]为预测状态量；Xref[k]为期望状态量.
式（7）中，第1项表示系统的跟踪误差，使其在预测区间Np内的预测状态量Xpre[k + 1]、Xpre[k + 2]、…、Xpre[k + Np]尽可能接近期望状态量Xref[k + 1]、Xref[k + 2]、 …、Xref[k + Np]，第2项为输入控制量以反映系统能量损耗情况. 另一方面，考虑控制过程中磁轮偏转角度受机械限位的限制，还需要对控制输入量进行幅值约束，如式（8）所示.
Umin⩽ (8) 式中：Umin和Umax分别为控制量的最小值和最大值.
根据预测模型、代价函数和系统约束，将NMPC控制问题转化有限时域内的非线性规划求解问题（NLP）[23]，如式（9）所示.
\begin{split} & \mathrm{min}\text{ }J,\\ & {\mathrm{s.t.}}{{\boldsymbol{X}}}_{\text{pre}[k]}={{\boldsymbol{X}}}_{[k]},\\ &\quad {{\boldsymbol{X}}}_{\text{pre}[k + i + 1]}=f({{\boldsymbol{X}}}_{\text{pre}[k + i]},{U}_{[k + i]}),\quad i\in [0,{N}_{\text{p}}-1],\\ &\quad {U}_{\mathrm{min}}\leqslant{U}_{[k + i]}\leqslant{U}_{\mathrm{max}}. \end{split} (9) 本文采用序列二次规划法（SQP）求解上述NLP问题，将求解得到的最优控制序列U*中第1项元素作为永磁电动悬浮汽车的控制输入U[k]，使其状态更新到下一时刻，并将更新后的状态作为新的初始条件，继续求解下一时刻的优化问题，上述过程循环进行.
2.1.3 反馈校正
MPC和NMPC中的反馈校正模块通过利用预测模型输出的预测值Xpre[k]与状态量测量滤波值（后验状态估计值）$ \hat {\boldsymbol{X}}_{[k]} $之间的误差e[k]，对预测值Xpre[k]进行实时反馈修正，得到预测输出校正值$ {\boldsymbol{X}}_{{\mathrm{c}}[k]} $，此值用作新的滚动优化的输入，以修正控制器输出U[k]，从而应对外部干扰和模型不精确导致的轨迹跟踪性能下降. 这类似于ESO对输入控制量的观扰反馈补偿. 但ESO的抗扰效果取决于带宽，过小的带宽可能降低系统的抗扰性能，而过大的带宽则可能对传感器测量噪声过于敏感，因此，需要牺牲部分抗扰效果，选取适当的带宽[24]. 为实现更好抗扰效果，采用上述反馈校正模块与ESO协同作用，对控制量进行辅助补偿，反馈校正模块如式（10）所示.
\left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{e}}_{[k]}} = {{\hat {\boldsymbol{X}}}_{[k]}} - {{\boldsymbol{X}}_{{\text{pre}}[k]}}}, \\ {{{\boldsymbol{X}}_{{\mathrm{c}}[k]}} = {{\boldsymbol{X}}_{{\text{pre}}[k]}} + {{\boldsymbol{K}}_{\mathrm{c}}}{{\boldsymbol{e}}_{[k]}}}, \end{array}} \right. (10) 式中：Kc为反馈校正系数.
2.2 扩张状态观测器设计
ESO可以将永磁电动悬浮汽车模型不精确部分、内部参数变化和时变干扰视为总扰动，将其进行估计并提前作用在控制器的输出部分进行反馈补偿，提高系统的鲁棒性[25]，ESO观测器设计如式（11）所示.
\left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{e}}_{\textit{z}}} = {{\boldsymbol{z}}_1} - {{\hat {\boldsymbol{x}}}_1}}, \\ {{{\dot {\boldsymbol{z}}}_1} = {{\boldsymbol{z}}_2} - {\beta _1}{{\boldsymbol{e}}_{\textit{z}}}}, \\ {{{\dot {\boldsymbol{z}}}_2} = {{\boldsymbol{z}}_3} - {\beta _2}f_{\mathrm{al}}({{\boldsymbol{e}}_{\textit{z}}},{\gamma _1},\delta ) + b{u_{\text{t}}}}, \\ {{{\dot {\boldsymbol{z}}}_3} = - {\beta _3}f_{\mathrm{al}}({{\boldsymbol{e}}_{\textit{z}}},{\gamma _2},\delta )} , \end{array}} \right. (11) f_{\mathrm{al}}({{\boldsymbol{e}}}_{{\textit{z}}},\gamma ,\delta )=\left\{\begin{array}{l} |{{\boldsymbol{e}}}_{{\textit{z}}}{|}^{\gamma }\text{sgn}({{\boldsymbol{e}}}_{{\textit{z}}}),\quad |{{\boldsymbol{e}}}_{{\textit{z}}}| > \delta, \\ {{\boldsymbol{e}}}_{{\textit{z}}},\quad \left|{{\boldsymbol{e}}}_{{\textit{z}}}\right|\leqslant\delta, \end{array}\right. (12) 式中：z1、z2、z3分别为永磁电动悬浮汽车的横向位移、速度与总扰动的状态观测值，ez为横向位移的观测值与传感器测量滤波值之间的误差，$ {\hat {\boldsymbol{x}}_1} $为传感器测量并滤波后的位移反馈值，β1、β2、β3为恒正增益系数，ut为作用到系统上的控制量，b为控制增益. $ f_{\mathrm{al}}({{\boldsymbol{e}}_{\textit{z}}},\gamma ,\delta ) $为非线性函数，γ与$ \delta $为恒正参数，γ的取值与fal函数的非线性形状有关[17].
为方便计算，将式（11）离散化为
\left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{e}}_{{\mathrm{z}}[k]}} = {{\boldsymbol{z}}_{1[k]}} - {{\hat {\boldsymbol{x}}}_{1[k]}}}, \\ {{{\boldsymbol{z}}_{1[k + 1]}} = {{\boldsymbol{z}}_{1[k]}} + {T_{\mathrm{z}}}\left( {{{\boldsymbol{z}}_{2[k]}} - {\beta _1}{{\boldsymbol{e}}_{{\mathrm{z}}[k]}}} \right)}, \\ {{{\boldsymbol{z}}_{2[k + 1]}} = {{\boldsymbol{z}}_{2[k]}} + {T_{\mathrm{z}}}\left( {{{\boldsymbol{z}}_{3[k]}} - {\beta _2}f_{\mathrm{al}}\left( {{{\boldsymbol{e}}_{{\mathrm{z}}[k]}},{\gamma _1},\delta } \right) + b{u_{{\mathrm{t}}[k]}}} \right)}, \\ {{{\boldsymbol{z}}_{3[k + 1]}} = {{\boldsymbol{z}}_{3[k]}} - {T_{\mathrm{z}}}{\beta _3}f_{\mathrm{al}}\left( {{{\boldsymbol{e}}_{{\mathrm{z}}[k]}},{\gamma _2},\delta } \right)}, \end{array}} \right. (13) 式中：Tz为ESO离散步长.
2.3 扩展卡尔曼滤波器设计
测量传感器本身存在一定误差，为准确获取永磁电动悬浮汽车系统状态量，降低测量噪声对控制器的影响，提高控制系统的表现，本文采用最优化的、递归的、数字处理的扩展卡尔曼滤波算法[26]对噪声信号进行滤波降噪，用其提供的状态估计值作为控制器的反馈输入.
结合式（5），同时为系统加入过程噪声$ {\boldsymbol{w}} = [ {w_1}\;\;\; {w_2} ]^{\text{T}} $，其符合正态分布，期望为0，协方差矩阵为Le，即$ {\boldsymbol{w}}\sim {\boldsymbol{N}}(0,{{\boldsymbol{L}}_{\mathrm{e}}}) $，得到带有过程噪声的非线性状态空间方程（式（14））.
{{\boldsymbol{X}}_{[k + 1]}} = g({{\boldsymbol{X}}_{[k]}},{U_{[k]}},{{\boldsymbol{w}}_{[k]}}), (14) 式中：
\begin{split} & g({{\boldsymbol{X}}_{[k]}},{U_{[k]}},{{\boldsymbol{w}}_{[k]}}) = \\ &\quad \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{x}}_1} + {{\boldsymbol{x}}_2}{T_{\mathrm{s}}} + {w_{1[k]}}} \\ {{{\boldsymbol{x}}_2} + \left( {\dfrac{{2{{\boldsymbol{F}}_{{\mathrm{r}}0}}}}{m}\sin\; u - \dfrac{c}{m}{{\boldsymbol{x}}_2}} \right){T_{\mathrm{s}}} + {w_{2[k]}}} \end{array}} \right], \end{split} (15) Ts为EKF采样时间.
使用传感器测量系统的位移x1与速度x2，并考虑符合正态分布、期望为0和协方差矩阵为Me的测量误差$ {\boldsymbol{\lambda}} = {[{\lambda _1}\;\;\;{\lambda _2}]^{\text{T}}},{\boldsymbol{\lambda}} \sim {\boldsymbol{N}}(0,{{\boldsymbol{M}}_{\text{e}}}) $，定义测量方程为
{{\boldsymbol{Z}}_{[k]}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{x}}_1} + {\lambda _{1[k]}}} \\ {{{\boldsymbol{x}}_2} + {\lambda _{2[k]}}} \end{array}} \right] \text{，} (16) 式中：Z[k]为状态观测量，λ1[k]和λ2[k]分别为位移和速度的测量误差.
首先，求时刻k + 1的g(·)对X和w偏导的雅可比矩阵A和W，分别如式（17）、（18）所示.
{{\boldsymbol{A}}_{[k + 1]}} = \frac{{\partial g}}{{\partial {\boldsymbol{X}}}}({\hat {\boldsymbol{X}}_{[k]}},{U_{[k]}},0) = \left[ {\begin{array}{*{20}{c}} 1&{{T_{\mathrm{s}}}} \\ 0&{1 - \dfrac{c}{m}{T_{\mathrm{s}}}} \end{array}} \right], (17) {{\boldsymbol{W}}_{[k + 1]}} = \frac{{\partial g}}{{\partial {\boldsymbol{w}}}}({\hat {\boldsymbol{X}}_{[k]}},{U_{[k]}},0) = \left[ {\begin{array}{*{20}{c}} 1&0 \\ 0&1 \end{array}} \right]. (18) 将式（15）、 （17）、式（18）代入时间更新方程（式（19）、（20）），求时刻k + 1的先验状态估计值$ \bar {\hat {\boldsymbol{X}}}_{[k + 1]} $和先验状态估计误差协方差矩阵$ \bar {\boldsymbol{P}}_{[k + 1]} $.
{\bar {\hat {\boldsymbol{X}}}_{[k + 1]}} = g({\hat {\boldsymbol{X}}_{[k]}},{U_{[k]}},0) (19) {\bar {\boldsymbol{P}}_{[k + {\text{l}}]}} = {{\boldsymbol{A}}_{[k + {\text{l}}]}}{{\boldsymbol{P}}_{[k + {\text{l}}]}}{\boldsymbol{A}}_{[k]}^{\text{T}} + {{\boldsymbol{W}}_{[k + {\text{l}}]}}{{\boldsymbol{L}}_{\text{e}}}{\boldsymbol{W}}_{[k + {\text{l}}]}^{\text{T}} (20) 式中：P[k + 1]为时刻k + 1后验状态估计误差协方差矩阵，由时刻k计算得到.
然后，再求时刻k + 1的h(·)对X和$ {\boldsymbol{\lambda}} $偏导的雅可比矩阵H和$ {\boldsymbol{\varLambda}} $.
{{\boldsymbol{H}}_{[k + 1]}} = \frac{{\partial h}}{{\partial {\boldsymbol{X}}}}({\bar {\hat {\boldsymbol{X}}}_{[k + 1]}},0) = \left[ {\begin{array}{*{20}{c}} 1&0 \\ 0&1 \end{array}} \right], (21) {{\boldsymbol{\varLambda}} _{[k + 1]}} = \frac{{\partial h}}{{\partial {\boldsymbol{\lambda}} }}({\bar {\hat {\boldsymbol{X}}}_{[k + 1]}},0) = \left[ {\begin{array}{*{20}{c}} 1&0 \\ 0&1 \end{array}} \right]. (22) 将式（21）和式（22）计算的结果代入扩展卡尔曼的测量更新方程（式（23）～（25）），求时刻k + 1的卡尔曼增益K[k + 1]、状态量测量滤波值（后验状态估计值）$ \hat {\boldsymbol{X}}_{[k + 1]} $，并更新后验状态估计误差协方差矩阵P[k + 2]. 上述过程进行迭代循环运行.
{{\boldsymbol{K}}_{[k + 1]}} = \frac{{{{\bar {\boldsymbol{P}}}_{[k + 1]}}{\boldsymbol{H}}_{[k + 1]}^{\text{T}}}}{{{{\boldsymbol{H}}_{[k + 1]}}{{\bar {\boldsymbol{P}}}_{[k + 1]}}{\boldsymbol{H}}_{[k + 1]}^{\text{T}} + {{\boldsymbol{V}}_{[k + 1]}}{{\boldsymbol{M}}_{\text{e}}}{\boldsymbol{\varLambda}} _{[k + 1]}^{\text{T}}}}, (23) {\hat {\boldsymbol{X}}_{[k + 1]}} = {\bar {\hat {\boldsymbol{X}}}_{[k + 1]}} + {{\boldsymbol{K}}_{[k + 1]}}\left({{\boldsymbol{Z}}_{[k + 1]}} - h({\bar {\hat {\boldsymbol{X}}}_{[k + 1]}},0)\right), (24) {{\boldsymbol{P}}_{[k + 2]}} = ({\boldsymbol{I}} - {{\boldsymbol{K}}_{[k + 1]}}{{\boldsymbol{H}}_{[k + 1]}}){\bar {\boldsymbol{P}}_{[k + 1]}}. (25) 综上，本文设计的NMPC-ESO-EKF控制器结构如图5所示.
2.4 收敛性证明
由于本文对u进行了小幅值约束，$ \sin\; u \approx u $，可对式（6）线性化处理来简化收敛性证明，线性化后的离散型状态空间方程如式（26）所示.
\left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{X}}_{[k + 1]}} = {{\boldsymbol{AX}}_{[k]}} + {{\boldsymbol{BU}}_{[k]}}}, \\ {{{\boldsymbol{Y}}_{[k]}} = {{\boldsymbol{CX}}_{[k]}}} , \end{array}} \right. (26) 式中：$ {\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 1&{{T_{\text{s}}}} \\ 0&{1 - \dfrac{c}{m}{T_{\text{s}}}} \end{array}} \right] $，$ {\boldsymbol{B}} = {\left[ {\begin{array}{*{20}{c}} 0&{\dfrac{{2{F_r}}}{m}{T_s}} \end{array}} \right]^T} $，$ {\boldsymbol{C}} = [\begin{array}{*{20}{c}} 1&0 \end{array}] $.
定义期望目标为Xref[k]的动态方程为
{{\boldsymbol{X}}_{{\mathrm{ref}}[k + 1]}} = {{\boldsymbol{A}}_D}{{\boldsymbol{X}}_{{\mathrm{ref}}[k]}}, (27) 式中：AD为目标状态转移矩阵.
组合式（26）和式（27），定义新的状态空间方程为
{{\boldsymbol{X}}_{\mathrm{a}[k + 1]}} = {{\boldsymbol{A}}_{\mathrm{a}}}{{\boldsymbol{X}}_{\mathrm{a}[k]}} + {{\boldsymbol{B}}_{\mathrm{a}}}{{\boldsymbol{U}}_{[k]}}, (28) 式中：$ {{\boldsymbol{X}}_{\mathrm{a}[k]}} = {[\begin{array}{*{20}{c}} {{{\boldsymbol{X}}_{[k]}}}&{{{\boldsymbol{X}}_{{\mathrm{ref}}}}_{[k]}} \end{array}]^{\mathrm{T}}} $，$ {{\boldsymbol{A}}_{\mathrm{a}}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{{{\boldsymbol{A}}_{\mathrm{D}}}} \end{array}} \right] $，$ {{\boldsymbol{B}}_{\mathrm{a}}} = {[\begin{array}{*{20}{c}} {\boldsymbol{B}}&{\boldsymbol{0}} \end{array}]^{\mathrm{T}}} $.
定义状态误差为
{{\boldsymbol{e}}_{\mathrm{d}[k]}} = {{\boldsymbol{X}}_{[k]}} - {{\boldsymbol{X}}_{{\mathrm{ref}}[k]}} = {{\boldsymbol{C}}_{\mathrm{a}}}{{\boldsymbol{X}}_{\mathrm{a}[k]}}, (29) 式中：$ {{\boldsymbol{C}}_{\mathrm{a}}} = [\begin{array}{*{20}{c}} {{{\boldsymbol{I}}_{2 \times 2}}}&{ - {{\boldsymbol{I}}_{2 \times 2}}} \end{array}] $.
将式（29）带入式（7），并令$ C_{\text{a}}^{\text{T}}Q{C_{\text{a}}} = {Q_{\text{a}}} $，则式（7）可进一步表示为
J = \sum\limits_{i = 1}^{{N_p}} \| {X_{{\text{a}}[k + i]}}\| _{ {Q_{\text{a}}}}^2 + \sum\limits_{i = 0}^{{N_c} - 1} \| {U_{[k + i]}}\| _{ R}^2. (30) 通过上述定义将轨迹跟踪控制问题转化为标准调节问题，可得到控制输入U[k]解析解[21]为
{{\boldsymbol{U}}_{[k]}} = - \left[ {\begin{array}{*{20}{c}} 1&{\mathbf{0}}& \cdots &{\mathbf{0}} \end{array}} \right]{({{\boldsymbol{\varGamma}} ^{\mathrm{T}}}{\boldsymbol{\varOmega \varGamma}} + {\boldsymbol{\varPsi}} )^{ - 1}}({\boldsymbol{\varGamma \varOmega \varPhi}} ){{\boldsymbol{X}}_{[k]}}, (31) 式中：$ {\boldsymbol{\varGamma}} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{B}}_{\mathrm{a}}}}&0& \cdots &0 \\ {{{\boldsymbol{A}}_{\mathrm{a}}}{{\boldsymbol{B}}_{\mathrm{a}}}}&{{{\boldsymbol{B}}_{\mathrm{a}}}}& \cdots &0 \\ \vdots & \vdots & & \vdots \\ {{\boldsymbol{A}}_{\mathrm{a}}^{{N_{\mathrm{p}}} - 1}{\boldsymbol{B}}}&{{\boldsymbol{A}}_{\mathrm{a}}^{{N_{\mathrm{p}}} - 2}{{\boldsymbol{B}}_{\mathrm{a}}}}& \cdots &{{{\boldsymbol{B}}_{\mathrm{a}}}} \end{array}} \right]_{2{N_{\mathrm{p}}} \times {N_{\mathrm{p}}}}} $，$ {\boldsymbol{\varOmega}} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{Q}}_{\mathrm{a}}}}& {} &{}& {} \\ {} &{{{\boldsymbol{Q}}_{\mathrm{a}}}}& {}& {} \\ {} &{}& {\ddots}& {} \\ {}& {} && {}{{{\boldsymbol{Q}}_{\mathrm{a}}}} \end{array}} \right]_{2{N_{\mathrm{p}}} \times 2{N_{\mathrm{p}}}}} $，${\boldsymbol{\varPsi}} = {\left[ {\begin{array}{*{20}{c}} R& {}& {} &{} \\ { } & {{\boldsymbol{R}}}& {} &{} \\ {}& {}& {} \ddots & {} \\ {} && {} {} &R \end{array}} \right]_{{N_{\mathrm{p}}} \times {N_{\mathrm{p}}}}}$，$ {\boldsymbol{\varPhi}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{{{\boldsymbol{A}}^2}}& \cdots &{{{\boldsymbol{A}}^{{N_{\mathrm{p}}}}}} \end{array}} \right]_{2{N_{\mathrm{p}}} \times 2}^T $.
将状态量X[k]的估计值$\hat {\boldsymbol{X}}_{[k]}$反馈输入给控制器，令${K_{{\mathrm{mpc}}}} = \left[ {\begin{array}{*{20}{c}} 1&{\mathbf{0}}& \cdots &{\mathbf{0}} \end{array}} \right]{({{\boldsymbol{\varGamma}} ^{\mathrm{T}}}{\boldsymbol{\varOmega \varGamma}} + {\boldsymbol{\varPsi}} )^{ - 1}}({\boldsymbol{\varGamma \varOmega \varPhi}} )$，则最终的控制量U[k]表示为
{U_{[k]}} = - {{\boldsymbol{K}}_{{\mathrm{mpc}}}}{\hat {\boldsymbol{X}}_{[k]}}. (32) 根据式（26）定义系统标称模型为
{{\boldsymbol{X}}_{[k + 1]}} = {{\boldsymbol{AX}}_{[k]}} + {{\boldsymbol{BU}}_{[k]}} + {{\boldsymbol{D}}_{[k]}}, (33) 式中：$ {{\boldsymbol{D}}_{[k]}} = {[ \begin{array}{*{20}{c}} 0&{{d_{[k]}}} \end{array} ]^{\mathrm{T}}} $为系统时变有界扰动值.
假设系统扰动为0，不考虑ESO，根据式（32）～（33）得到系统理想轨迹为
{{\boldsymbol{X}}_{{\mathrm{p}}[k + 1]}} = {{\boldsymbol{AX}}_{{\mathrm{p}}[k]}} + {\boldsymbol{B}}( - {{\boldsymbol{K}}_{{\mathrm{mpc}}}}{\hat {\boldsymbol{X}}_{{\mathrm{p}}[k]}}), (34) 式中：Xp[k]和$ \hat {\boldsymbol{X}}_{{\mathrm{p}}[k]} $分别为理想轨迹的系统状态量真实值和EKF估计值.
当存在时变扰动时，引入ESO来补偿控制输入，根据式（13）得到补偿后的控制量Ut[k]为
\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{{\mathrm{t}}[k]}} = {{\boldsymbol{U}}_{{\mathrm{r}}[k]}} - {{\boldsymbol{U}}_{{\mathrm{eso}}[k]}}} ,\\ {\boldsymbol{U}}_{\mathrm{eso}[k]} = \dfrac{{{T_{\mathrm{z}}}}}{b}{{\boldsymbol{z}}_{3[k]}} = \dfrac{1}{b}{{\hat d}_{[k]}} , \end{array}} \right. (35) 式中：令$ \hat d_{[k]} = T_{\mathrm{z}} {\boldsymbol{z}}_{3[k]} $，用以保证与d[k]量纲统一；Ur[k]为实际轨迹补偿前的控制输入；Ueso[k]为ESO补偿量.
联立式（32）～（33）和式（35），得到系统的实际轨迹如式（36）所示.
{{\boldsymbol{X}}_{{\mathrm{r}}[k + 1]}} = {{\boldsymbol{AX}}_{{\mathrm{r}}[k]}} + {\boldsymbol{B}}\left( - {{\boldsymbol{K}}_{{\mathrm{mpc}}}}{\hat {\boldsymbol{X}}_{{\mathrm{r}}[k]}} - \frac{1}{b}{\hat d_{[k]}}\right) + {{\boldsymbol{D}}_{[k]}}, (36) 式中：Xr[k]和$ \hat {\boldsymbol{X}}_{{\mathrm{r}} [k]} $分别为实际轨迹的系统状态量真实值和EKF估计值.
定义系统的理想轨迹与实际轨迹误差epr[k]为
{{\boldsymbol{e}}_{{\mathrm{pr}}[k + 1]}} = {{\boldsymbol{X}}_{\mathrm{p}[k + 1]}} - {{\boldsymbol{X}}_{\mathrm{r}[k + 1]}}. (37) 将式（34）和式（36）带入式（37），可得
{{\boldsymbol{e}}_{{\mathrm{pr}}[k + 1]}} = {{\boldsymbol{Ae}}_{{\mathrm{pr}}[k]}} - {{\boldsymbol{BK}}_{{\text{mpc}}}}({\hat {\boldsymbol{X}}_{{\mathrm{p}}[k]}} - {\hat {\boldsymbol{X}}_{{\mathrm{r}}[k]}}) + \frac{1}{b}{\boldsymbol{B}}{\hat d_{[k]}} - {{\boldsymbol{D}}_{[k]}}. (38) 分别定义理想和实际轨迹的EKF估计值为
\left\{ {\begin{array}{*{20}{c}} {{{\hat {\boldsymbol{X}}}_{{\mathrm{p}}[k]}} = {{\boldsymbol{X}}_{{\mathrm{p}}[k]}} + {{\boldsymbol{e}}_{{\mathrm{p}}[k]}}}, \\ {{{\hat {\boldsymbol{X}}}_{{\mathrm{r}}[k]}} = {{\boldsymbol{X}}_{{\mathrm{r}}[k]}} + {{\boldsymbol{e}}_{{\mathrm{r}}[k]}}}, \end{array}} \right. (39) 式中：ep[k]和er[k]分别为理想和实际轨迹的EKF估计误差值.
则式（38）可以进一步表示为
\begin{split} & {{\boldsymbol{e}}}_{{\mathrm{pr}}[k + 1]}= ({\boldsymbol{A}}-{{\boldsymbol{BK}}}_{\text{mpc}}){{\boldsymbol{e}}}_{{\mathrm{pr}}[k]}-{{\boldsymbol{BK}}}_{\text{mpc}}({{\boldsymbol{e}}}_{{\mathrm{p}}[k]}-{{\boldsymbol{e}}}_{{\mathrm{r}}[k]})\\ &\quad + \frac{1}{b}({\boldsymbol{B}}{\widehat{d}}_{[k]}-b{{\boldsymbol{D}}}_{[k]}). \end{split} (40) 文献[27]对EFK估计误差序列进行收敛分析，并得到期望的收敛速率，ep[k]和er[k]在一定条件下收敛，则$ {\boldsymbol{e}}_{{\mathrm{p}}[k]} - {\boldsymbol{e}}_{{\mathrm{r}}[k]} $有界. 离散ESO观测器的收敛性已被证明[28]，根据兼顾带宽扩大和噪声抑制的ESO参数整定方法[29]来选取适当参数，可使$ \hat d[k] $收敛于D[k]中的d[k].
当b适当取值时，$ {\boldsymbol{B}}\hat d_{[k]} - b{\boldsymbol{D}}_{[k]} $有界，由于$ {\boldsymbol{A}} - {\boldsymbol{BK}}_{\mathrm{mpc}} $是赫尔维兹矩阵[30]，且ep[k]-er[k]有界，则误差epr[k]关于原点渐进稳定. 因此，系统实际轨迹收敛于理想轨迹，而理想轨迹已经被证明在NMPC-EKF共同作用下系统收敛[31].
3. 仿真验证
为验证本文设计的永磁电动悬浮汽车横向运行方式和NMPC-ESO-KEF控制算法的有效性，并考虑永磁电动悬浮汽车的运行成本与安全性，首先在MATLAB中搭建Simulink-Simscape多体动力学联合仿真平台进行验证，永磁电动悬浮汽车的模型参数如表1所示. 本文设计3组仿真实验，分别为横向定常数轨迹跟踪控制、外部扰动条件下的横向抗扰控制和内外部扰动条件下的横向抗扰控制.
表 1 永磁电动悬浮汽车模型参数Table 1. Model parameters of PMEDS vehicle编号 参数 数值 整车
参数整车质量m/kg
悬浮间隙h/mm
横向阻尼系数c/（N·s·m−1）18.1
10
17.2磁轮
参数极对数P
外径Ro/mm
内径Ri/mm
宽度d/mm
磁化角q/°
磁轮转速n0/rpm
磁轮磁阻力Fr0/N4
50
32.5
35
90
2000
17.33磁轮偏转角度范围α/° −20～20 3.1 横向定常数轨迹跟踪性能分析
设定仿真初始状态为$ {{\boldsymbol{X}}_0} = {[0\;\;0]^{\mathrm{T}}} $，在第2 s处加入0.5 m的阶跃跟踪信号，采用PID、MPC与NMPC控制方法进行对比验证，该工况没有加入测量噪声扰动，所以不考虑EKF滤波效果，各控制器下的系统输出和控制输入响应如图6所示.
仿真结果表明，4种控制方法均能将系统稳定控制在0.5 m处，证明本文提出的系统横向运行模式的有效性，各控制器系统响应速度较为接近，PID超调量最大，MPC也存在一定超调，但远小于PID；NMPC与NMPC-ESO将系统超调控制到最小，两者跟踪曲线基本一致，证明ESO对系统模型不精确引起的内部扰动观测估计的值较小，系统动力学建模较为准确，整个控制过程偏转角度不超过20°，验证了输入控制量约束的有效性.
3.2 横向抗扰性能分析
3.2.1 外部扰动条件下的横向抗扰分析
为验证改进NMPC的抗外部扰动能力，在系统静浮状态下加入幅值为2 N、频率为0.1 Hz的方波信号（第5 s时加入）来模拟永磁电动悬浮汽车受到的外部横向阵风扰动，同时额外施加高斯白噪声来使扰动更符合真实情况，外部扰动如图7所示. 设定静浮初始状态为$ {{\boldsymbol{X}}_0} = {[0\;\;0]^{\mathrm{T}}} $，系统响应曲线如图8所示.
仿真结果表明，在外部扰动下，NMPC-ESO控制的系统横向位移变化幅值明显小于NMPC、MPC和PID控制，收敛速度也最快，因此可以证明NMPC-ESO的横向抗扰能力最好，MPC横向抗扰能力最差. ESO的扰动观测值如图9所示，加入ESO后对NMPC的抗扰性能有显著提升，能较为准确地对外部扰动进行实时估计，然后对输入偏转角度进行反馈补偿，从而提高系统抗扰性能.
3.2.2 内外部扰动条件下的横向抗扰分析
为验证NMPC-ESO-EKF在内外部干扰条件下的抗扰能力，在系统位移输出中加入高斯噪声来模拟永磁电动悬浮汽车传感器测量噪声，将其作为内部扰动，在第2 s加入0.5 m的阶跃参考信号，在第12 s施加2N的横向力来模拟外部扰动，同时与NMPC、NMPC-EKF、NMPC-ESO进行对比仿真，分析ESO与EKF在控制过程中的作用，系统响应结果如图10和表2所示.
表 2 系统响应结果对比Table 2. Comparison of system response results控制器 0～12 s 12 s后 平均误差/mm 性能提升/% 平均误差/mm 性能提升/% NMPC 71.31 34.41 NMPC-ESO 319.54 −348.10 183.07 −432.03 NMPC-EKF 57.99 18.70 6.12 82.20 NMPC-ESO-EKF 57.93 18.76 3.52 89.77 仿真结果表明，除NMPC-ESO以外的3种控制器上升时间基本一致，在NMPC控制下，系统可以稳定的控制在0.54 m附近，存在一定误差，而且波动较大；在NMPC-EKF控制下，系统可以稳定收敛在期望目标值，位移输出与输入偏转角度的波动幅值都较小，EKF的滤波效果如图11所示，滤波值与真实值误差较小，使系统有效消除了内部扰动；在NMPC-ESO控制下，系统不能稳定收敛，因为ESO对环境噪声比较敏感，会把传感器的测量噪声误认为扰动进行反馈补偿引起系统失稳；在NMPC-ESO-EKF控制下的系统响应表现最好，EKF滤波保证了ESO观扰的有效性，实现最佳控制效果，在12 s后同时施加内外扰动条件下，将系统横向轨迹跟踪性能提升89.77%.
4. 实验验证
4.1 实验平台
为进一步验证本文提出的永磁电动悬浮汽车横向运行模式及其改进NMPC控制器的有效性，利用已有的1∶50原理样机进行实验验证，并搭建横向跑动实验装置保证实验过程的安全性，实验平台如图12所示.
原理样机主要由4个磁轮、无刷电机、FOC （field oriented control）矢量驱动器和舵机组成，磁轮转速通过FOC矢量控制，舵机偏转角度通过PID控制，以上控制以及整车横向轨迹跟踪控制器运算均通过dSPACE以及匹配的Autodesks上位机实现，横向位移通过激光传感器实时反馈得到，并通过24 V直流电源集中供电. 横向跑动实验装置主要通过横向导轨滑块解除横向自由度的限制，底部铺设铜板.
实验参数与表1一致，依然采用PID、MPC和NMPC控制进行对比实验验证，共设计3组实验工况，分别为小距离定常数轨迹跟踪控制、方波轨迹跟踪控制与大距离定常数轨迹跟踪下的抗扰控制.
4.2 小距离定常数轨迹跟踪性能分析
设定实验初始条件为$ {{\boldsymbol{X}}_0} = {[0.2\;\;\;0]^{\mathrm{T}}} $，在第2 s施加0.2 m至0.7 m的阶跃参考信号，为避免传感器测量噪声影响各个控制器的控制性能，均采用EKF对传感器反馈数据进行在线滤波处理，EKF滤波效果如图13所示，传感器自身原始数据噪声较小，但采样频率只有10 Hz，EKF滤波性能较优.
实验结果如图14所示，在施加阶跃信号后，PID-EKF响应最快，但存在192.20 mm的超调量，调节时间最长达到9 s，稳定后也存在小幅值的持续震荡；其余3种控制器性能差距不大，其中，NMPC-ESO-EKF超调为2.12 mm，超调量最小，调节时间最短，缩减到4.7 s，性能相较于PID-EKF分别提升98.90%和47.78%，但其输入偏转角度曲线存在强烈波动现象，这是铜板表面不平整导致的，各处的悬浮间隙不一致，对磁轮电磁力产生影响，导致系统模型参数不确定，ESO将其观测并持续以偏转角度的形式补偿，使系统表现出最好跟随性能，整个控制过程偏转角度不超过20°，验证了输入控制量约束的有效性.
图 14 小距离定常数轨迹跟踪控制系统响应Figure 14. System response of short-distance constant trajectory tracking control4.3 方波信号轨迹跟踪性能分析
设定实验初始条件为$ {{\boldsymbol{X}}_0} = {[0.9\;\;\;0]^{\mathrm{T}}} $，在第5 s时施加幅值为0.25 m、周期为60 s的方波信号，实验结果如图15和表3所示. 4种控制器均能有效跟踪方波信号，使系统具备了灵活左右横移的基本功能，适用于横向变道、爬坡、停车入位等场景，大大提高了永磁电动悬浮汽车的机动性能和交通效率. 其中，PID-EKF在方波信号跳变时出现较大超调和震荡，在实际工程应用中存在安全性问题，平均跟踪误差为57.60 mm，跟踪性能表现最差；MPC-EKF、NMPC-EKF、NMPC-ESO-EKF的轨迹跟踪曲线均较为平稳，其中，NMPC-ESO-EKF的平均超调量与平均跟踪误差最小，相较于PID-EKF性能分别提升93.77%和36.13%，控制表现最佳.
图 15 方波信号轨迹跟踪控制系统响应Figure 15. System response of square wave signal trajectory tracking control表 3 方波信号轨迹跟踪控制系统响应结果对比Table 3. Comparison of system response results of square wave signal trajectory tracking controlmm 控制器 平均跟踪误差 平均超调量 PID-EKF 57.60 164.96 MPC-EKF 38.48 20.48 NMPC-EKF 36.82 15.69 NMPC-ESO-EKF 36.79 10.27 4.4 大距离定常数轨迹跟踪下的抗扰性能分析
为进一步探索系统更大距离横向轨迹跟踪性能和抗扰性，设定实验初始条件为$ {{\boldsymbol{X}}_0} = {[0.7\;\;\;0]^{\mathrm{T}}} $，在2 s时施加0.7 m到1.7 m的阶跃参考信号，并在第22 s时施加扰动，扰动形式为图4中1号和4号磁轮对称偏转10°，即α1=−α4=10°，产生水平向左的固有电磁力来模拟系统扰动，实验结果如图16和表4所示. 在定常数轨迹跟踪方面，该工况与小距离轨迹跟踪控制的结果基本一致，NMPC-ESO-EKF整体性能最好，该工况下各控制器会产生更大的超调量. 施加扰动后，各控制器均能使系统回到稳定目标值，稳定后输入偏转角度在−10°附近波动，与1、4号磁轮对称偏转10°产生的电磁扰动力相互抵消，实现抗扰. MPC-EKF的位移波动幅值最大，恢复时间最长；NMPC-EKF的位移波动幅值和恢复时间略小于PID-EKF；NMPC-ESO-EKF相较于PID-EKF波动幅值减小34.51%，恢复时间缩短42.08%，抗扰性能大幅提升，表现最佳.
图 16 大距离定常数轨迹跟踪控制系统响应Figure 16. System response of long-distance constant trajectory tracking control表 4 大距离定常数轨迹跟踪控制系统响应结果对比Table 4. Comparison of system response results of long-distance constant trajectory tracking control控制器 位移波动幅值/mm 恢复时间/s PID-EKF 40.89 6.82 MPC-EKF 45.39 10.46 NMPC-EKF 35.73 6.80 NMPC-ESO-EKF 26.78 3.95 5. 结 论
本文率先针对永磁电动悬浮汽车横向运动控制展开研究，提出同侧磁轮对称偏转的横向运行模式，设计了基于系统横向非线性动力学模型和约束条件的NMPC-ESO-EKF横向控制器，通过仿真分析和实验验证得到以下结论：
1） 验证了通过对称偏转磁轮来补偿系统横向力的横向运行模式的有效性.
2） 引入EKF滤波器消除了传感器测量噪声，解决了ESO对输入噪声的敏感性问题.
3） 在定常数和方波信号轨迹跟踪控制下，NMPC-ESO-EKF超调量和平均跟踪误差最小，轨迹跟踪性能最好.
4） NMPC-ESO-EKF的抗干扰性能更好，施加扰动后横向位移波动幅值最小，系统恢复收敛速度大幅上升.
未来将聚焦在永磁电动悬浮汽车系统多自由度耦合控制研究等工作中.
-
图 1 年发生频率-落石尺寸关系
Figure 1. Relationships between annual frequency of occurrence and rockfall size
图 2 100 a重现期落石尺寸频次与累积概率分布
Figure 2. Frequency and cumulative probability distribution of rockfall size during 100-year return period
图 7 不同vmin的E(μ)-μ0和ξ-μ0曲线
Figure 7. E(μ)-μ0 and ξ-μ0 curves with different vmin values
年份 v>3.0 m3 2.0 m3<v≤3.0 m3 1.5 m3<v≤2.0 m3 1.0 m3<v≤1.5 m3 0.5 m3<v≤1.0 m3 v≤0.5 m3 合计 1961 年 1 1 0 10 12 1962 年 4 1 2 2 17 26 1963 年 4 2 1 4 30 41 1964 年 3 4 2 3 5 54 71 1965 年 1 2 2 45 50 1966 年 5 2 3 12 99 121 1967 年 1 1 1 5 8 1969 年 3 1 5 26 35 1970 年 1 1 5 16 51 74 1971 年 1 1 2 3 3 30 40 1972 年 3 2 3 3 2 39 52 合计 23 15 10 24 52 406 530 分占比例% 4.34 2.83 1.89 4.53 9.81 76.6 100.00 
下载: 导出CSV
表 2 GPD参数σ、ξ估计
Table 2. Estimated results of GPD parameters ξ and σ
条件 ξ σ 90% 置信区间 [0.40,0.60] [0.86,1.02] 95% 置信区间 [0.38,0.63] [0.84,1.04] 均值 0.50 0.94
下载: 导出CSV
表 3 各重现期对应落石尺寸均值及其置信区间
Table 3. Mean rockfall sizes and their confidence intervals during different return periods
m3 重现期/a 均值 90% 置信区间 95% 置信区间 vmax=4.0 m3 vmax=6.0 m3 vmax=4.0 m3 vmax=6.0 m3 vmax=4.0 m3 vmax=6.0 m3 50 4.93 13.42 [4.28，5.89] [11.30，14.92] [4.13, 6.12] [11.03, 15.26] 100 5.09 15.91 [4.39，6.12] [13.34，17.86] [4.23, 6.38] [12.89，18.72] 1000 5.41 26.75 [4.60，6.69] [21.23，31.83] [4.42, 7.01] [20.28，33.38] 10000 5.58 42.99 [4.70，7.02] [31.95，53.43] [4.51, 7.39] [30.21，56.86]
下载: 导出CSV
表 4 各设计基准期落石尺寸标准值
Table 4. Standard values of rockfall size in different design reference periods
m3 设计基准期/a 设防
水准T/a 1/T/a− 1 1−1/T N 年不超越
概率
((1− 1/T)N)N 年超越
概率
(1−(1−1/T)N)与 T 对应的
落石尺寸标准值vmax=4.0 m3 vmax=6.0 m3 50 Ⅰ 50 0.0200 0.9800 0.36 0.64 4.93 13.4 Ⅱ 475 0.0021 0.9979 0.90 0.10 5.33 22.8 Ⅲ 2475 0.0004 0.9996 0.98 0.02 5.49 32.4 100 Ⅰ 100 0.0100 0.9900 0.37 0.63 5.09 15.9 Ⅱ 900 0.0011 0.9989 0.89 0.11 5.40 26.2 Ⅲ 5000 0.0002 0.9998 0.98 0.02 5.54 37.4
下载: 导出CSV
-
[1] MARCHELLI M, COLTRINARI G, ALFARO DEGAN G, et al. Towards a procedure to manage safety on construction sites of rockfall protective measures[J]. Safety Science, 2023, 168: 106307.1-106307.8. [2] 中华人民共和国住房和城乡建设部. 铁路工程结构可靠性设计统一标准:GB 50216—2019[S]. 北京:中国计划出版社,2019. [3] 中华人民共和国交通运输部. 公路工程结构可靠性统一标准:JTG 2120—2020[S]. 北京:人民交通出版社,2020. [4] 国家铁路局. 铁路隧道设计规范:TB 10003—2016[S]. 北京:中国铁道出版社, 2017. [5] 中华人民共和国交通运输部. 公路隧道设计规范 第一册 土建工程:JTG 3370.1—2018[S]. 北京:人民交通出版社, 2019. [6] 中华人民共和国住房和城乡建设部. 建筑结构可靠性设计统一标准:GB 50068—2018[S]. 北京:中国建筑工业出版社. [7] 铁道第二勘测设计院. 铁路工程设计技术手册. 隧道[M]. 北京:中国铁道出版社,1999. [8] 曾廉. 明洞顶设计荷载的研究[J]. 铁路标准设计通讯,1974(7): 3-17. [9] ILLEDITSCH M, PREH A. Determination of meaningful block sizes for rockfall modelling[J]. Natural Hazards, 2024, 120(6): 5685-5710. doi: 10.1007/s11069-024-06432-4 [10] STELZER G, BICHLER A. ONR 24810–A comprehensive guideline for building better rockfall protection structures[C]//64th Annual Highway Geology Symposium. New Hampshire: Highway Geology Symposium, 2013. [11] MARCHELLI M, DE BIAGI V, PEILA D. Reliability-based design of rockfall passive systems height[J]. International Journal of Rock Mechanics and Mining Sciences, 2021, 139: 104664.1-104664.9. [12] MARCHELLI M, DE BIAGI V, PEILA D. Reliability-based design of protection net fences: influence of rockfall uncertainties through a statistical analysis[J]. Geosciences, 2020, 10(8): 280.1-280.24. [13] DE BIAGI V, MARCHELLI M, PEILA D. Reliability analysis and partial safety factors approach for rockfall protection structures[J]. Engineering Structures, 2020, 213: 110553.1-110553.1. [14] Austrian Standards Institute. Technical protection against rockfall-terms and definitions, effects of actions, design, monitoring and maintenance: ONR 24810[S]. Vienna: [s.n.], 2017. [15] MELZNER S, ROSSI M, GUZZETTI F. Impact of mapping strategies on rockfall frequency-size distributions[J]. Engineering Geology, 2020, 272: 105639.1-105639.11. [16] LAIMER H J. Determination of rockfall design blocks in Upper Triassic limestones and Dolomites (Dachstein Formation, Northern Calcareous Alps)[J]. Bulletin of Engineering Geology and the Environment, 2020, 79(3): 1581-1590. doi: 10.1007/s10064-019-01640-w [17] DE BIAGI V, LIA NAPOLI M, BARBERO M, et al. Estimation of the return period of rockfall blocks according to their size[J]. Natural Hazards and Earth System Sciences, 2017, 17(1): 103-113. doi: 10.5194/nhess-17-103-2017 [18] CURCEAC S, ATKINSON P M, MILNE A, et al. An evaluation of automated GPD threshold selection methods for hydrological extremes across different scales[J]. Journal of Hydrology, 2020, 585: 124845.1-124845.15. [19] STORMY A. MATLAB[M]. Fifth edition. Boston: [s.n.], 2019. [20] COLES S. An Introduction to Statistical Modeling of Extreme Values[M]. London: Springer, 2001. [21] 史道济. 实用极值统计方法[M]. 天津:天津科学技术出版社,2006. [22] 仇学艳,王超,秦崇仁. 阈值法在河海工程设计中的应用[J]. 水利学报,2001(8): 32-37. doi: 10.3321/j.issn:0559-9350.2001.08.006QIU Xueyan, WANG Chao, QIN Chongren. The application of threshold procedure in design of river and coastal engineering[J]. Shuili Xuebao, 2001(8): 32-37. doi: 10.3321/j.issn:0559-9350.2001.08.006 [23] PHILIP J, DAVID R, JENNY W, et al. Uncertainties in return from extreme value analysis of peaks over threshold using the generalised Pareto distribution[J]. Ocean Engineering, 2021, 220: 107725.1-107725.17. [24] THOMAS R K, ILARIA P. Use of peak over threshold data for flood frequency estimation: An application at the UK national scale.[J]. Jounal of Hydrology, 2023, 626(5): 130235.1-130235.31. [25] 中华人民共和国住房和城乡建设部. 建筑抗震设计规范GB 50011—2010(2016版)[S]. 北京:中国建筑工业出版社,2016. -
下载:
下载:
点击查看大图
计量
- 文章访问数: 160
- HTML全文浏览量: 34
- PDF下载量: 44
- 被引次数: 0