Two-Stage Optimization Method of Train Energy-Efficient Operation Based on Dynamic Programming
-
摘要: 针对地铁列车多站间节能运行优化问题,提出将列车节能驾驶优化过程和时刻表优化过程结合的两阶段优化方法,分别求解两优化过程的全局最优解,从而获取列车在多站间运行的最优操纵策略. 首先考虑节能和节时两个目标,构建列车节能驾驶多目标优化模型,结合动态规划多阶段寻优思路,建立一系列包含多个过程指标及约束的子阶段求解模型,逆序求解后获取列车站间运行最优操纵策略的Pareto前沿;其次建立时刻表优化模型,基于动态规划方法,调用各站间Pareto前沿,搜索站间运行时间最优分配方案;最后以北京地铁亦庄线为例,验证两阶段优化方法的有效性和高效性. 试验结果表明,与最速操纵策略相比,经过两个阶段优化后的列车牵引能耗分别降低了53.87%和54.69%,两阶段优化过程分别用时258.90 s和0.08 s.Abstract: Focusing on the problem of train energy-efficient operation between multi-sections in urban rail transit, a two-stage optimization method was proposed by integrating the processes of energy-efficient driving optimization and timetable optimization. To obtain the optimum train driving strategy between multi-sections, each process was solved with global optimum solutions, respectively. First, to realize energy saving and time saving, a multi-objective energy-efficient driving model was constructed. Utilizing the multistage-based dynamic programming searching approach, a series of sub-stage models that contain multiple objects and constrains were constructed. The Pareto front of the optimum driving strategy was generated by inverse order method. Then, a timetable optimization model was constructed, in which the Pareto front of sections was applied, and the optimum running time allocation of multi-sections was searched by dynamic programming approach. A case study of Yizhuang urban rail line in Beijing was conducted to verify the effectiveness and efficiency of the two-stage optimization method. Compared with the flat-out running strategy, the optimization of two stages resulted in 53.87% and 54.69% energy saving improvement respectively; the calculation time of two process was 258.90 s and 0.08 s respectively.
-
Key words:
- subways /
- energy efficiency /
- optimization /
- train operation /
- dynamic programming
-
随着中速磁悬浮列车的广泛商业化应用[1],作为其关键核心技术的电磁悬浮控制技术已经成熟稳定[2]. 商业化运营过程中随着规模和复杂性的迅速增加,对中速磁悬浮列车悬浮系统的安全性和可靠性的需求不断增长. 在过去几十年中,对磁浮列车悬浮系统的性能监测尤其是异常检测的研究引起了极大的兴趣[3-4].
磁浮列车悬浮系统的异常指在列车实际运行中,悬浮系统出现工作状态与期望状态不匹配,但没有导致故障或失效的情况. 悬浮系统的异常具有出现时刻未知,持续时间短,通过控制器调节可以恢复稳定的特点. 对悬浮系统异常的精确检测可以提前获取系统的运行情况,有利于提前安排预防性维护计划,降低工程事故率[5].
针对磁悬浮列车的异常检测研究通常可以分为针对轨道的检测和针对悬浮系统的检测. 针对磁浮列车轨道的异常检测问题,Zhang等[6]提出基于小波与分形相结合的高速磁悬浮列车长轨道磁异常检测方法,以发现高速磁悬浮列车长定子铁芯的局部短路故障;Deng等[7]对高温超导磁悬浮列车运行时永磁导轨的不规则性进行测量和表征,并设计了一种不规则性检测设备;罗茹丹[8]提出了“2+1”的检测方法实现对长定子行波主漏磁场异常情况的动态检测;林国斌等[9]提出一种基于车轨状态监测的悬浮冗余控制系统,用于在控制列车稳定悬浮的同时对车轨状态进行检测,并通过对悬浮间隙和电磁铁振动情况的分类和学习确定状态类型,具有实时性好、容错性高等优点;杨杰等[10]提出一种基于悬挂式磁悬浮轨道交通系统的轨道维护设备,在对轨道状态进行取样和检测后对轨道异常实现诊断. 文献[11]收集了振动条件下高温超导钉扎系统的水平和垂直加速度数据.
针对磁浮列车悬浮系统的异常检测问题,朱跃欧等[12]提出一种磁浮列车悬浮控制器的异常预警方法与检测系统,通过分析悬浮间隙值所在区间,实现了车辆异常状况的在线检测与预警;王平等[13]提出了一种基于改进典型相关分析的中低速悬浮系统异常检测方法;通过运营线数据验证了该方法能获得更高的检测率;王平等[14]提出了一种基于超球体高斯分布的悬浮系统异常检测方法,该方法对系统异常检测的结果较为精确. 以新型悬挂式永磁磁浮列车研究对象,马政[15]设计并实现了一套基于物联网的永磁磁浮列车远程监控系统,实现了对列车悬浮状态的异常监控;Wang等[16]引入了支持向量机方法,利用悬浮传感器采集的实时信号来检测磁悬浮列车悬浮系统的异常状态;Ma等[17]将磁悬浮列车悬浮控制系统收集的数据用于训练机器学习模型,并根据数据确定磁悬浮列车的实际运行状态,以达到实时监控的目的.
然而,上述方法存在2个问题:1) 大部分方法是基于单变量的异常检测,没有充分利用系统的运行数据. 磁浮列车悬浮系统的异常检测大多都是只基于悬浮间隙,而实际运行中电流、加速度数据同样包含系统的异常信息. 因此,如何应用列车悬浮系统运行时产生的多变量数据进行异常检测是一个关键问题. 2) 需要有异常的先验信息,而实际系统的异常先验信息是不完全可知的,这无疑会降低异常检测的准确率;此外,面向实际工程系统的异常检测模块或系统在设计时需考虑异常漏检率和异常误检率的折中,这需要从数理的角度进行分析.
基于上述分析,本文将基于参数化残差的异常检测方法引入悬浮系统的异常检测问题,以最低异常漏检率为设计目标,提出了磁悬浮列车悬浮系统基于最低漏检率的残差生成器设计方法及异常检测方法. 结合磁悬浮列车悬浮系统的运行数据,验证了该残差生成器能够在满足期望误报率的情况下,实现异常漏检率的最小化.
1. 残差生成器与异常检测
1.1 参数化残差生成器
系统出现异常时的状态空间模型[18]为
{x(k+1)=Ax(k)+Bu(k)+Bff(k)+w(k),y(k)=Cx(k)+Du(k)+Dff(k)+v(k), (1) 式中:x(k)为时刻k的系统状态;u(k)为时刻k的系统输入;y(k)为时刻k的系统输出;A、B、C、D为一定维数的系统参数矩阵;Bf和Df为异常参数矩阵;f(k)为异常向量;w(k)为系统的过程噪声;v(k)为系统的输出噪声.
则状态检测量z(k)为
{z(k)=Ψs⊥[us(k)ys(k)],us(k)=[uT(k−s)uT(k−s+1)⋯uT(k)]T,ys(k)=[yT(k−s)yT(k−sH)⋯yT(k)]T, (2) 式中:s为正整数,表示截断数据的长度;us(k)、ys(k)为长度为s+1的输入、输出数据矩阵;Ψs为基于数据驱动的系统信息矩阵[19];Ψs⊥为Ψs的正交补矩阵.
当系统出现异常时的状态检测量为
z(k)=Ψs⊥[us(k)ys(k)]=[0Hw,sWk,s+Vk,s+Hf,sFs], (3) 式中:Hw,s为过程噪声的Hankel矩阵;Wk,s为数据驱动的过程噪声矩阵;Vk,s为数据驱动的过程噪声矩阵;Hf,s为异常向量f(k)的Hankel矩阵;Fs为基于数据驱动的异常矩阵.
从式(3)可以看出,当系统出现异常时,z(k)不仅包含噪声信息,同时也包含了异常信息[20].
在基于模型的异常检测或故障检测领域,残差是通过比较实际系统的输出y(k)与估计输出ˆy(k)得到的,即定义r(k)=y(k)−ˆy(k). 残差数据中包含了系统扰动和系统异常的全部信息,残差产生器设计的主要目标是设法让异常检测系统对异常敏感,同时对其他信号具有鲁棒性.
结合Youla参数化和基于输入输出数据的状态检测量,残差生成器参数化形式为
r(k)=gTz(k)=gTΨs⊥[us(k)ys(k)], (4) 式中:g为参数向量,是残差生成器和异常检测的核心[21].
出于实际工程运营需要,通常在系统运行前就需要完成g与Ψs⊥的设计,并且在一段时间内都不能改变其参数.
1.2 基于残差的异常检测
由于残差数据中包含了系统扰动和系统异常的全部信息,因此,在系统运行时,可以通过对残差进行分析、处理,以及基于某种准则或阈值,来进一步掌握系统的运行状态信息,判断系统的运行状态.
当系统运行时,残差生成器实时获得残差数据r(k),并将残差数据送入残差评估单元,计算残差评估函数J(r)=Jk(r). 在设置评估函数阈值Jth后,根据式(5)判断系统是否发生异常
{Jk(r)>Jth,异常,Jk(r)⩽Jth,健康. (5) 2. 基于最低漏检率的残差生成器设计
2.1 基于均值和方差的置信集
传统的异常检测方案通常采用基于假设检验的方法进行残差评估,即假设残差的分布知识是准确已知的,或者可以从历史数据中估计出来[22]. 然而,对于磁浮列车悬浮系统而言,使用历史过程数据获得高度可信的概率分布估计比较繁琐,经验估计与真实概率分布的偏差可能会导致异常检测结果不可靠[23].
尽管难以获得悬浮系统的高度可信概率分布,但在无异常时的z(k)可以视为属于同一个集合. 考虑到悬浮系统运行数据的统计特性,可以依据数据的均值和协方差矩阵对其进行建模.
设无异常情况下,z(k)满足
Sh={Sz∈ζz|E=ˉzh,V=Σh}, (6) 式中:Sh为无异常数据的置信集;ζz为所有z(k)的分布集合;ˉzh和Σh为无异常时状态检测变量的均值和方差,ˉzh和Σh是通过对健康运行数据进行计算得来的,通常可以选择刚投入运行时无异常的运行数据,也可以是某一次无异常运行时的数据;E为z(k)的均值;V为z(k)的方差.
在悬浮系统运行初期,由于缺乏异常数据,难以分析出异常的数据分布特性. 因此,通常会依据正常运行时的数据特点预先设计一个异常检测系统,以实现对大部分潜在异常的检测. 假设系统运行时会出现M个典型异常,且发生第i个异常时的异常信号设为fi(k),i=1,2,⋯,M. 为不失一般性,假设每个异常发生的概率相同. 采用均值和协方差矩阵对异常进行建模. 设异常情况时的z(k)满足
Sfi={Sz∈ζz|E=ˉzfi,V=Σfi}, (7) 式中:Sfi为异常数据的置信集;ˉzfi为zfi(k)的均值,zfi(k)为第i个异常时的状态检测变量;Sfi为zfi(k)的方差.
当有可用的异常数据时,可以直接使用该数据建立置信集. 若无异常数据或异常数据少时,则可以自定义M个异常情况下的均值和协方差,进而构建置信集,当有新可用的异常数据时,记录数据并离线更新置信集.
2.2 基于最低漏检率的异常检测
针对第 i个异常,基于异常检测的残差信号为
ri(k)=gTiz(k), (8) 式中:gi为第i个异常时的参数向量.
因此,第i个异常时的残差评估函数J(r)与阈值Jth为
{J(r)=(ri(k)−ri,h)2,ri,h=gTiˉzh,Jth=1, (9) 式中:ri,h为第i个无异常场景的残差.
定义异常误检率PF=Pr,定义异常漏检率为{P_{\rm{M}}} = \Pr \left\{ {J({\boldsymbol{r}}) \leqslant 1\left| {{{{\boldsymbol{f}}}_i}(k) \ne {\boldsymbol{0}}} \right.} \right\}. 理想情况下,异常检测系统应该对所有潜在的异常检测且不存在异常的误报,即{P_{\rm{F}}} = {P_{\rm{M}}} = 0,在实际中却难以实现. 当设计的{P_{\rm{M}}}越小时,往往会导致{P_{\rm{F}}}较高,反之亦然. 因此,异常检测系统通常会采取{P_{\rm{F}}}和{P_{\rm{M}}}的折中.
本文研究的悬浮系统的最低异常漏检率可被表述为{P_{\rm{F}}}固定时的最小化{P_{\rm{M}}},即考虑当{P_{\rm{F}}}固定且可接受时,如何使得{P_{\rm{M}}}最小化[24]. 根据式(8)可知,{{{\boldsymbol{g}}}_i}是{{\boldsymbol{r}}_i}(k)生成的关键, {{{\boldsymbol{g}}}_i} 的取值会影响异常检测的{P_{\rm{F}}}和{P_{\rm{M}}}. 因此,基于最低漏检率的残差生成器的设计最终转化为求得参数向量{{{\boldsymbol{g}}}_i}({{{\boldsymbol{g}}}_i} \ne {\boldsymbol{0}}),使得
\begin{gathered} \min {\alpha _i}, \\ {\rm{ s.t.}}\left\{ \begin{gathered} \sup \Pr \left\{ {{{\left( {{{\boldsymbol{g}}}_i^{\rm{T}}\left( {{{\boldsymbol{z}}}(k) - {{\bar {\boldsymbol{z}}}_{\rm{h}}}} \right)} \right)}^2} \leqslant 1} \right\} \leqslant {\alpha _i}, \\ \sup \Pr \left\{ {{{\left( {{{\boldsymbol{g}}}_i^{\rm{T}}\left( {{{\boldsymbol{z}}}(k) - {{\bar {\boldsymbol{z}}}_{\rm{h}}}} \right)} \right)}^2} > 1} \right\} \leqslant {\beta _i}, \\ \end{gathered} \right. \end{gathered} (10) 式中: {\alpha _i} \in (0,1) , {\beta _i} \in (0,1) ,分别为{P_{\rm{M}}}和{P_{\rm{F}}}的上界, {\beta _i} 在设计时应预先给出[25].
进一步,异常检测式(10)可以被重构为
\begin{gathered} \max a({{{\boldsymbol{g}}}_i}), \\ {\rm{ s.t.}}\left\{ \begin{gathered} {k_{\beta_i}}\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}\varSigma _{\rm{h}} {{{\boldsymbol{g}}}_i}} \leqslant 1, \\ a({{{\boldsymbol{g}}}_i}) = \frac{{\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}} {{\tilde {\boldsymbol{z}}}_{{\boldsymbol{f}}_i}} {\tilde {\boldsymbol{z}}}_{{\boldsymbol{f}}_i}^{\rm{T}} {{\boldsymbol{g}}_i}} - {k_{\beta_i}}\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}\varSigma_{\rm{h}} {{{\boldsymbol{g}}}_i}} }}{{\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}\varSigma_{{{\boldsymbol{f}}_i}} {{{\boldsymbol{g}}}_i}} }} > 0 , \\ \end{gathered} \right. \\ \end{gathered} (11) 式中:{k_{\beta_i}} = \sqrt {1/{\beta _i}};{\tilde{\boldsymbol{z}}}_{{\boldsymbol{f}}_i} = {\bar {\boldsymbol{z}}_{{\boldsymbol{f}}_i}} - {\bar {\boldsymbol{z}}_{\rm{h}}}.
当获得最优{{{\boldsymbol{g}}}_i}后,只需将{{{\boldsymbol{g}}}_i}乘以一个非零常数,即可满足{k_{\beta_i}}\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}{\varSigma _{\rm{h}}}{{{\boldsymbol{g}}}_i}} \leqslant 1. 因此,先分析a({{{\boldsymbol{g}}}_i}) > 0时a({{{\boldsymbol{g}}}_i})的最大值.
令
\left\{ \begin{gathered} b({{{\boldsymbol{g}}}_i}) = \frac{{\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}} {{\boldsymbol{z}}_{{\boldsymbol{f}}_i}} {{\boldsymbol{z}}_{{\boldsymbol{f}}_i}^{\rm{T}}} {{{\boldsymbol{g}}}_i}} }}{{\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}{{\boldsymbol{\varSigma }} _{{\boldsymbol{f}}_i}}{{{\boldsymbol{g}}}_i}} }}, \\ c({{{\boldsymbol{{\boldsymbol{g}}}}}_i}) = \frac{{{k_{\beta_i}}\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}{{\boldsymbol{\varSigma }}_{\rm{h}}}{{{\boldsymbol{g}}}_i}} }}{{\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}{{\boldsymbol{\varSigma }} _{{\boldsymbol{f}}_i}}{{{\boldsymbol{g}}}_i}} }}, \\ \end{gathered} \right. (12) 式中:{{\boldsymbol{\varSigma }} _{{\boldsymbol{f}}_i}} \ne 0;b({{{\boldsymbol{g}}}_i}) > c({{{\boldsymbol{g}}}_i}) > 0.
则d({{{\boldsymbol{g}}}_i}) = c({{{\boldsymbol{g}}}_i})/b({{{\boldsymbol{g}}}_i}) < 1,因此,可以得到式(13)的关系.
\left\{ \begin{gathered} a({{{\boldsymbol{g}}}_i}) = b({{{\boldsymbol{g}}}_i}) - c({{{\boldsymbol{g}}}_i}) = (1 - d({{{\boldsymbol{g}}}_i}))b({{{\boldsymbol{g}}}_i}), \\ e({{{\boldsymbol{g}}}_i}) = {b^2}({{{\boldsymbol{g}}}_i}) - {c^2}({{{\boldsymbol{g}}}_i}) = \frac{{1 + d({{{\boldsymbol{g}}}_i})}}{{1 - d({{{\boldsymbol{g}}}_i})}}{a^2}({{{\boldsymbol{g}}}_i}). \\ \end{gathered} \right. (13) 对e({{{\boldsymbol{g}}}_i})求关于a({{{\boldsymbol{g}}}_i})的导数,得
\frac{{\partial e({{{\boldsymbol{g}}}_i})}}{{\partial a({{{\boldsymbol{g}}}_i})}} = 2a({{{\boldsymbol{g}}}_i})\frac{{1 + d({{{\boldsymbol{g}}}_i})}}{{1 - d({{{\boldsymbol{g}}}_i})}} > 0. (14) 因此,e({{\boldsymbol{g}}_i})关于a({{{\boldsymbol{g}}}_i})是单调增加的,所以式(11)可以变为求{{{\boldsymbol{g}}}_i} = \arg \mathop {\max }\limits_{{{\boldsymbol{g}}_i} \ne 0} e({{{\boldsymbol{g}}}_i})的问题[26],即
\mathop {\max }\limits_{{{\boldsymbol{g}}_i} \ne 0} e({{{\boldsymbol{g}}}_i}) = \mathop {\max }\limits_{{{\boldsymbol{g}}_i} \ne 0} \frac{{\sqrt {{{\boldsymbol{g}}}_i^{\rm{T }} ({{\tilde {\boldsymbol{z}}}_{{{\boldsymbol{f}}_i}}} { \tilde{\boldsymbol{z}}^{\rm{T}}_{{{\boldsymbol{f}}_i}}} - k_{\beta_i}^2{{{\boldsymbol{\varSigma }}} _{\rm{h}}}) {{{\boldsymbol{g}}}_i}} }}{{\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}{{{\boldsymbol{\varSigma }}} _{{{\boldsymbol{f}}_i}}}{{{\boldsymbol{g}}}_i}} }}. (15) 通过上述分析,基于最低漏检率的残差生成器的设计问题最终转换为求解式(15)的广义特征值和特征向量[27],具体求解步骤如算法1所示.
1) 算法 1:最优参数向量的求解
步骤 1 对\sqrt {{{{\boldsymbol{\varSigma }}} _{{{\boldsymbol{f}}_i}}}}做 SVD 分解,\sqrt {{{{\boldsymbol{\varSigma }}} _{{\boldsymbol{f}}_i}}} = {{{{\boldsymbol{U}}}}_i}\left[ { {{{{{\boldsymbol{S}}}}_i}}\;\; 0 } \right] {{{\boldsymbol{V}}}}_i^{\rm{T}},其中,\lambda _{{\rm{max}}} 为最大奇异值;{\boldsymbol{U}}_i、{\boldsymbol{S}}_i、{\boldsymbol{V}}_i 均为SVD分解后得出的矩阵;{\boldsymbol{v}}_i 为{\boldsymbol{V}}_i 中的元素向量.
步骤 2 通过求解方程组
\left\{ \begin{gathered} {{\boldsymbol{v}}}_i^{\rm{T}}({\lambda _{m,i}}{{{\boldsymbol{I}}}} - {{{\boldsymbol{S}}}}_i^{ - 1}{{{\boldsymbol{U}}}}_i^{\rm{T}}{{{\varXi}} _i}{{{{\boldsymbol{U}}}}_i}{{{\boldsymbol{S}}}}_i^{ - 1}){{{\boldsymbol{v}}}_i} = 0, \\ {\varXi _i} = {{\tilde {\boldsymbol{z}}}_{{{\boldsymbol{f}}_i}}} { {\boldsymbol{z}}^{\rm{T}}_{{{\boldsymbol{f}}_i}}} - k_{{\beta _i}}^2{{\boldsymbol{\varSigma}} _{\rm{h}}} , \\ {\lambda _{m,i}} = {\lambda _{\max }}\left\{ {{{{\boldsymbol{S}}}}_i^{ - 1}{{{\boldsymbol{U}}}}_i^{\rm{T}}{\varXi _i}{{{{\boldsymbol{U}}}}_i}{{{\boldsymbol{S}}}}_i^{ - 1}} \right\}, \\ {{\boldsymbol{v}}}_i^{\rm{T}}{{{\boldsymbol{v}}}_i} = 1, \\ \end{gathered} \right. 得出{{{\boldsymbol{v}}}_i}.
步骤 3 依据 {{{\boldsymbol{g}}}_i} = {{{{\boldsymbol{U}}}}_i}{{{\boldsymbol{S}}}}_i^{ - 1}{{{\boldsymbol{v}}}_i}计算 {{{\boldsymbol{g}}}_i} .
步骤 4 依据 {{\hat {\boldsymbol{g}}}_i} = \dfrac{{{{\boldsymbol{g}}}_i}}{{{k_{{\beta _i}}}\sqrt {{{\boldsymbol{g}}}_i^{\rm{T}}{{{{\boldsymbol{\varSigma}}}} _{\rm{h}}}{{{\boldsymbol{g}}}_i}} }}计算出满足式(10)的最优参数向量{{\hat {\boldsymbol{g}}}_i}.
在得出最优参数向量后,第i个异常的参数化残差生成器可以通过式(16)构建.
{{{\boldsymbol{r}}}_i}(k) = {{\hat {\boldsymbol{g}}}_i}^{\rm{T}}{{{{\boldsymbol{\varPsi}} }}_s}^ \bot \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{u}}}_s}(k)} \\ {{{{\boldsymbol{y}}}_s}(k)} \end{array}} \right]. (16) 令{{\hat {\boldsymbol{r}}}_{i,{\rm{h}}}} = {{\hat {\boldsymbol{g}}}_i}^{\rm{T}}{\bar {\boldsymbol{z}}_{\rm{h}}},异常检测器的残差评估函数和异常检测阈值为
\left\{ \begin{gathered} {\boldsymbol{J}}({{{\boldsymbol{r}}}_i}) = {\left( {{{{\boldsymbol{r}}}_i}\left( k \right) - {{{\hat {\boldsymbol{r}}}}_{i,h}}} \right)^2}, \\ {J_{{\rm{th}}}} = 1. \\ \end{gathered} \right. (17) 完整的异常检测见算法 2
2) 算法 2:磁悬浮列车悬浮系统异常检测算法
① 离线辨识
步骤 1 采集无异常情况下的输入输出数据,基于子空间辨识法离线辨识出矩阵{{{{\boldsymbol{\varPsi}} }}_{\boldsymbol{s}}}^ \bot.
步骤 2 利用无异常的输入输出数据,计算得到{\bar {\boldsymbol{z}}_{\rm{h}}}和{{\boldsymbol{\varSigma}} _{\rm{h}}},构建健康状况下的置信集{S_{\rm{h}}}.
步骤 3 建立M个异常的置信集{S_{{{\boldsymbol{f}}_i}}} . 当有M个异常时的输入输出数据,计算每个异常的 {\bar {\boldsymbol{z}}_{{{\boldsymbol{f}}_i}}} 和{{{{\boldsymbol{\varSigma}} }}_{{\boldsymbol{f}}_i}},建立 {S _{{{\boldsymbol{f}}_i}}} . 如果无异常数据,预设每个异常的{\bar {\boldsymbol{z}}_{{{\boldsymbol{f}}_i}}}和{{{{\boldsymbol{\varSigma}}}} _{{{\boldsymbol{f}}}_i}},依此建立 {S _{{{\boldsymbol{f}}_i}}} .
步骤 4 确定每个异常下的{\beta _i},根据算法1计算出{{\hat {\boldsymbol{g}}}_i}和{{\hat {\boldsymbol{r}}}_{i,h}}.
② 在线异常检测
步骤 5 在线收集时间间隔为\left[ {k-s,k} \right]的输入输出数据,并构造数据矩阵\left[ {{{{\boldsymbol{u}}}_s}(k)}\quad{{{{\boldsymbol{y}}}_s}(k)} \right]^{\rm{T}}.
步骤 6 依据式(16)计算{{{\boldsymbol{r}}}_i}\left( k \right),依据式(17)执行异常检测并输出异常检测的结果.
步骤 7 令k = k + 1,重复步骤5和步骤6,直到中止时刻.
3. 异常检测实验验证
本节结合长沙磁浮快线某辆磁浮列车的悬浮系统的健康数据和异常数据,验证所提出的基于最低漏检率的残差生成器在悬浮系统异常检测中的有效性.
3.1 悬浮系统的3类异常
以长沙磁浮快线的运营磁浮列车为例,如图1所示,作为典型的常导电磁型磁浮列车[28],通常选用间隙作为其悬浮系统异常检测的检测变量,通过比较间隙值和异常阈值的大小来判断是否出现异常. 磁浮列车悬浮系统的工程异常阈值设置为标准间隙值 ±4 mm. 该线路的标准悬浮工作间隙为9 mm,工程异常阈值分别设为5 mm和13 mm.
磁浮列车悬浮系统经常会遭受到3类异常:间隙突变异常、砸轨异常和加速度传感器异常. 间隙突变异常的现象是指在运行中间隙有明显的波动,但没有超过工程异常阈值,此时加速度、电流也对应有一定的波动,出现两个尖峰,但电压几乎没有波动. 该异常通常是由于磁浮列车经过轨道接缝引起的,如图2所示.
砸轨异常的现象是间隙有剧烈波动,在一段时间内间隙值超过工程异常阈值并且需要经历5 s左右的恢复期,此时加速度、电流、电压均在一段时间内出现剧烈波动,如图3所示.
加速度计异常的现象为间隙没有超过工程异常阈值,也没有较为明显的波动,但加速度的波动十分明显且变为负值,电压也出现两个极小值尖峰. 由于实际加速度计存在着1g的偏置,当加速度计输出信号为负时,表明此刻的加速度已经低于 - 10\;{\rm{m/{s^2}}},如图4所示.
3.2 异常检测验证与分析
以磁浮列车悬浮系统的第13个悬浮点的运行数据作为检测对象,如图5所示.
列车正线行驶时存在着站内悬浮和站间行驶两个运行工况,考虑到上述3类异常均出现在站间行驶运行工况,因此,以速度为划分基准,分段提取站间行驶工况的运行数据. 本节选择一个无异常站间行驶数据段、3个包含间隙异常的数据段、2个包含砸轨异常的运行数据段和包含3个加速度传感器异常的数据段作为检测数据.
选择悬浮系统的电流数据作为输入数据,间隙数据和加速度数据作为输出数据,如图6所示. 依据正常运行时的输入输出数据建立{S_{\rm{h}}}并离线辨识出{{{{\boldsymbol{\varPsi}} }}_s}^ \bot.
依据3种异常的输入输出数据建立异常置信集 {S _{{{\boldsymbol{f}}_i}}} . 出于工程运营需求,设期望的故障误报率{P_{\rm{F}}} = 5{\text{%}} ,执行算法1可得
{{\hat {\boldsymbol{g}}}_i} =(0.257\;5, 0.257\;5,\cdots, 0.257\;5)_{1\times9}. (18) 基于最低漏检率的残差生成器为
{{{\boldsymbol{r}}}_i}(k) = 0.257\;5 {{{{\boldsymbol{\varPsi}} }}_s}^ \bot \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{u}}}_s}(k)} \\ {{{{\boldsymbol{y}}}_s}(k)} \end{array}} \right]. (19) 依据式(17)算出9个检测数据段的J({{\boldsymbol{r}}_i})并送入检测单元. 图7为无异常数据段的检测结果. 从图中可以看出,当系统无异常时,J({{\boldsymbol{r}}_i})小于0.4,均在异常检测阈值{J_{{\rm{th}}}}之下,不存在异常误检.
含间隙突变异常的3个数据段的检测结果如图8所示. 在每个数据段,大部分数据都低于异常检测阈值{J_{{\rm{th}}}},但是在出现异常的时刻,存在着J({{{\boldsymbol{r}}}_i})大于{J_{{\rm{th}}}}的情况,个数分别为8、6和5,说明该方案能够实现对间隙异常的检测. 此外,间隙突变异常的J({{{\boldsymbol{r}}}_i})的峰值略高于检测阈值{J_{{\rm{th}}}},一般在2以下,且不存在第二个大于{J_{{\rm{th}}}}的次尖峰. 表1为间隙突变异常的检测位置与异常实际发生位置的对比,从表中可以看出,间隙突变异常1的检测时刻落后实际异常发生时刻 0.2 s,其他两个数据段可以及时检测出异常.
表 1 间隙突变异常的检测位置与异常实际发生位置Table 1. Detection position and actual occurrence position of gap mutation anomalies异常 检测位置/s 实际位置/s 间隙异常 1 5.9 5.7 间隙异常 2 568.7 568.5 间隙异常 3 21.1 21.1 含砸轨异常的两个数据段的检测结果如图9所示. 从图中可以看出,该方法能够实现对砸轨异常的检测. 图9(a)中有20个样本的J({{{\boldsymbol{r}}}_i})>{J_{{\rm{th}}}},图9(b)中有20个样本的J({{{\boldsymbol{r}}}_i})>{J_{{\rm{th}}}},此外,砸轨异常位置的J({{{\boldsymbol{r}}}_i})存在着两个高于{J_{{\rm{th}}}}的峰值. 表2为砸轨异常的检测位置与异常实际发生位置对比. 从表2可以看出,砸轨异常的检测时刻落后于实际异常的发生时刻,落后约0.2 s. 结合图3可以发现,造成这一现象的主要原因是砸轨异常发生时刻数据变化不明显,间隙和加速度都有一个缓慢减小的阶段,没有大范围变化.
表 2 砸轨异常的检测位置与异常实际发生位置Table 2. Detection position and actual occurrence position of rail smashing anomalies异常 检测位置/s 实际位置/s 砸轨异常 1 7.1 7.0 砸轨异常 2 20.2 20.0 含加速度计异常的3个数据段的检测结果如图10所示. 在每个数据段,大部分数据都低于异常检测阈值{J_{{\rm{th}}}},但是在出现加速度计异常的时刻,存在着J({{{\boldsymbol{r}}}_i})>{J_{{\rm{th}}}}的情况,个数分别为5、7和4,且J({{{\boldsymbol{r}}}_i})的值较大,说明该方案能够实现对加速度计异常的检测. 表3为加速度计异常的检测位置与异常实际发生位置对比. 从表3中可以看出,一旦出现加速度计异常,本文所提出的方法可以立即检测出该异常,及时性较好.
表 3 加速度计异常的检测位置与异常实际发生位置Table 3. Detection position and actual occurrence position of accelerometer anomalies异常 检测位置/s 实际位置/s 加速度计异常 1 52.2 52.2 加速度计异常 2 52.1 52.1 加速度计异常 3 40.5 40.5 经过上述分析,基于残差生成器的异常检测方法可以用于磁悬浮列车悬浮系统的异常检测. 基于最低漏检率的异常检测方法能够实现磁悬浮列车悬浮系统3类典型异常的完全检测,不存在漏检的情况且不存在误检.
4. 结 论
本文针对电磁悬浮系统的异常检测问题开展了研究,以基于参数化残差的异常检测为核心,依据数据的均值和方差建立置信集,设计了基于最低异常漏检率的残差生成器,依据生成的残差数据在线进行异常检测. 最后,以长沙磁浮快线的运行数据为例,对悬浮系统常见的3类异常进行检测,实验结果表明,在给定异常误报率为5%时,本文所提出的方法能够实现对3类典型异常的全部检测,没有异常的漏检和误检. 检测时间对比结果进一步说明,本文提出的方法检测实时性较好.
-
表 1 符号定义
Table 1. Symbol definition
符号 定义 符号 定义 x 列车运行位置 P 列车最大牵引力 v 列车运行速度 Q 列车最大制动力 t 列车运行时间 \Delta v 速度离散间隔 M 列车总质量 \overline V \left( {{x_k}} \right) 位置{x_k}处的列车最大运行速度 \alpha 列车回转质量系数 e\left( {{s_{k,i}},{u_k}} \right) {s_{k,i}} 处施加 {u_k}后的列车牵引能耗 r 列车运行阻力 t\left( {{s_{k,i}},{u_k}} \right) {s_{k,i}} 处施加{u_k} 后的列车运行时间 a,b,c 戴维斯方程系数 F\left( {{s_{k,i}},{u_k}} \right) {s_{k,i}} 处施加{u_k} 后的列车状态转移方程 {T_{{\rm{run }},n}} 第n站间的运行时间 N 全线站间个数 g 重力加速度 {T_{{\rm{dwell }},n}} 第 n 站的停站时间 \theta 线路坡度 {T_{\min,n}} 第n站间的最小运行时间 X 站间距离 {E_{{\rm{total}}}} 列车全线总牵引能耗 {V_{\rm{limit}}} 线路限速 {T_{{\rm{total}}}} 列车全线总运行时间 K 站间子阶段划分个数 \Delta t 时间离散间隔 U 列车运行操纵策略 E_n^*\left( {{T_{{\rm{run }},n}}} \right) 第n站间条件{T_{{\rm{run }},n}}下的最小牵引能耗 {U^{\rm{*}}} 列车运行 Pareto 最优操纵策略 {A_n}、{B_n}、{C_n} 第n站间 Pareto 前沿函数拟合参数 {s_{k,i}} 节能驾驶优化第k子阶段第 i开始状态 L 全线运行时间分配策略 U\left( {{s_{k,i}}} \right) {s_{k,i}} 处的操纵子策略 {L^{\rm{*}}} 全线运行时间最优分配策略 U_w^*\left( {{s_{k,i}}} \right) {s_{k,i}} 处的第w个 Pareto 最优子策略 E\left( L \right) 时间分配策略L的牵引能耗 W\left( {{s_{k,i}}} \right) {s_{k,i}} 处的 Pareto 最优子策略个数 {D_{n,i}} 时刻表优化第n子阶段第 i开始状态 E\left( U \right) 操纵策略 U的牵引能耗 F\left( {{D_{n,i}},{T_{{\rm{run }},n}}} \right) {D_{n,i}} 运行{T_{{\rm{run }},n}}时间后的状态转移方程 T\left( U \right) 操纵策略 U的运行时间 L\left( {{D_{n,i}}} \right) {D_{n,i}}处的时间分配子策略 {u_k} 第k子阶段的控制工况 {L^*}\left( {{D_{n,i}}} \right) {D_{n,i}} 处的最优时间分配子策略 {x_k} 第k子阶段的开始位置 q\left( {{u_k}} \right) 工况{u_k}下的列车制动力 p\left( {{u_k}} \right) 工况 uk下的列车牵引力 表 2 不同控制工况下的牵引制动力输出
Table 2. Tractive and braking force under different control commands
控制工况 输出 最大牵引 p\left( {{u_k}} \right) = P,\;\;q\left( {{u_k}} \right) = 0, 巡航 p\left( {{u_k}} \right) = r,\;\;q\left( {{u_k}} \right) = 0, 惰行 p\left( {{u_k}} \right) = 0,\;\;q\left( {{u_k}} \right) = 0, 最大制动 p\left( {{u_k}} \right) = 0,\;\;q\left( {{u_k}} \right) = Q. 表 3 最优驾驶Pareto曲线拟合结果
Table 3. Fitting results of each Pareto front
区间 站间距离/m {A_n} {B_n} {C_n} 最优操纵策略解数/个 求解时间/s 亦庄站—次渠站 1 334 2.68 × 1014 6.83 3.28 116 8.92 次渠站—次渠南 1 286 2.88 × 1013 6.54 12.13 113 8.19 次渠南—经海路 2 086 2.49 × 1014 6.49 21.97 214 23.24 经海路—同济南 2 265 2.57 × 1011 4.93 8.33 247 44.83 同济南—荣昌东 2 338 2.75 × 1010 4.38 7.00 256 35.69 荣昌东—荣京东 1 354 4.90 × 109 4.43 3.85 106 10.20 荣京东—万源街 1 280 3.55 × 1010 4.94 3.92 109 7.01 万源街—文化园 1 538 1.84 × 1011 5.15 6.44 141 11.11 文化园—亦庄桥 993 3.90 × 1012 6.33 5.22 98 4.25 亦庄桥—旧宫站 1 982 3.56 × 107 3.07 2.22 167 16.56 旧宫站—小红门 2 366 3.79 × 109 3.96 1.95 165 47.60 小红门—肖村站 1 275 2.20 × 1011 5.36 6.23 126 4.71 肖村站—宋家庄 2 631 1.72 × 1014 6.06 10.28 195 36.59 总计 258.90 表 4 列车节能运行两阶段优化结果
Table 4. Results of two-stage optimization method
区间 优化第 1 阶段 优化第 2 阶段 最速牵引能耗/(kW•h) 时刻表运行时间/s 牵引能耗/
(kW•h)节能率/% 优化后运行时间/s 牵引能耗/
(kW•h)节能率/% 节能效果
对比/%亦庄站—次渠站 18.60 105 7.59 59.21 112 6.05 67.46 8.25 次渠站—次渠南 29.42 102 14.23 51.65 99 14.68 50.11 −1.54 次渠南—经海路 38.07 140 24.92 34.55 136 25.53 32.95 −1.60 经海路—同济南 22.79 150 13.17 42.20 149 13.33 41.49 −0.71 同济南—荣昌东 28.51 164 12.36 56.63 159 13.14 53.90 −2.73 荣昌东—荣京东 22.63 104 9.50 58.03 111 8.08 64.29 6.26 荣京东—万源街 21.59 103 7.98 63.02 106 7.45 65.51 2.49 万源街—文化园 24.22 114 11.26 53.52 119 10.30 57.46 3.94 文化园—亦庄桥 21.35 90 6.86 67.86 86 7.41 65.28 −2.57 亦庄桥—旧宫站 25.95 135 12.68 51.15 147 10.27 60.41 9.26 旧宫站—小红门 21.64 157 9.54 55.90 161 8.82 59.23 3.33 小红门—肖村站 23.35 108 9.07 61.17 105 9.53 59.19 −1.98 肖村站—宋家庄 32.07 190 12.93 59.67 172 15.13 52.82 −6.85 总计 330.19 1662 152.05 53.87 1662 149.76 54.69 0.82 -
中国轨道交通协会. 城市轨道交通2019年度统计和分析报告[EB/OL]. (2020-05-18)[2020-07-22]. http://www.camet.org.cn/tjxx/5133. GONZALEZ-GIL, PALACIN R, BATTY P, et al. A systems approach to reduce urban rail energy consumption[J]. Energy Conversion and Management, 2014, 80: 509-524. doi: 10.1016/j.enconman.2014.01.060 ICHIKAWA K. Application of optimization theory for bounded state variable problems to the operation of train[J]. Bulletin of JSME, 1968, 11(47): 857-865. doi: 10.1299/jsme1958.11.857 HOWLETT P G, PUDNEY P J, VU X. Local energy minimization in optimal train control[J]. Automatica, 2009, 45(11): 2692-2698. doi: 10.1016/j.automatica.2009.07.028 LIU R, GOLOVITCHER I M. Energy-efficient operation of rail vehicles[J]. Transportation Research Part A:Policy and Practice, 2003, 37(10): 917-932. doi: 10.1016/j.tra.2003.07.001 KHMELNITSKY E. On an optimal control problem of train operation[J]. IEEE Transactions on Automatic Control, 2000, 45(7): 1257-1266. doi: 10.1109/9.867018 CHANG C S, SIM S S. Optimising train movements through coast control using genetic algorithms[J]. IEE Proceedings-Electric Power Applications, 1997, 144(1): 65-73. doi: 10.1049/ip-epa:19970797 卢启衡,冯晓云,王青元. 基于遗传算法的追踪列车节能优化[J]. 西南交通大学学报,2012,47(2): 265-270.LU Qiheng, FENG Xiaoyun, WANG Qingyuan. Energy-saving optimal control of following trains based on genetic algorithm[J]. Journal of Southwest Jiaotong University, 2012, 47(2): 265-270. 李诚,王小敏. 基于粒子群优化的ATO控制策略[J]. 铁道学报,2017,39(3): 53-58. doi: 10.3969/j.issn.1001-8360.2017.03.010LI Cheng, WANG Xiaomin. An ATO control strategy based on particle swarm optimization[J]. Journal of the China Railway Society, 2017, 39(3): 53-58. doi: 10.3969/j.issn.1001-8360.2017.03.010 余进,何正友,钱清泉. 基于微粒群算法的多目标列车运行过程优化[J]. 西南交通大学学报,2010,45(1): 70-75. doi: 10.3969/j.issn.0258-2724.2010.01.012YU Jin, HE Zhengyou, QIAN Qingquan. Multi-objective train operation optimization based on particle swarm algorithm[J]. Journal of Southwest Jiaotong University, 2010, 45(1): 70-75. doi: 10.3969/j.issn.0258-2724.2010.01.012 MIYATAKE M, KO H. Optimization of train speed profile for minimum energy consumption[J]. IEEJ Transactions on Electrical and Electronic Engineering, 2010, 5(3): 263-269. doi: 10.1002/tee.20528 LU Shaofeng, HILLMANSEN S, HO T K, et al. Single-train trajectory optimization[J]. IEEE Transactions on Intelligent Transportation Systems, 2013, 14(2): 743-750. doi: 10.1109/TITS.2012.2234118 刘炜,王栋,李群湛,等. 基于时间逼近搜索算法的城轨列车运行节能优化研究[J]. 西南交通大学学报,2016,51(5): 918-924. doi: 10.3969/j.issn.0258-2724.2016.05.014LIU Wei, WANG Dong, LI Qunzhan, et al. A novel time-approaching search algorithm for energy-saving optimization of urban rail train[J]. Journal of Southwest Jiaotong University, 2016, 51(5): 918-924. doi: 10.3969/j.issn.0258-2724.2016.05.014 唐海川,朱金陵,王青元,等. 一种可在线调整的列车正点运行节能操纵控制算法[J]. 中国铁道科学,2013,34(4): 89-94. doi: 10.3969/j.issn.1001-4632.2013.04.15TANG Haichuan, ZHU Jinling, WANG Qingyuan, et al. An on-line adjustable control algorithm for on-time and energy saving operations of trains[J]. China Railway Science, 2013, 34(4): 89-94. doi: 10.3969/j.issn.1001-4632.2013.04.15 ALBRECHT T, OETTICH S. A new integrated approach to dynamic schedule synchronization and energy-saving train control[C]//Computers in Railways. Southampton: WIT Press, 2002: 847-856. 丁勇,刘海东,栢赟,等. 地铁列车节能运行的两阶段优化模型算法研究[J]. 交通运输系统工程与信息,2011,11(1): 96-101. doi: 10.3969/j.issn.1009-6744.2011.01.016DING Yong, LIU Haidong, BAI Yun, et al. A two-level optimization model and algorithm for energy-efficient urban train operation[J]. Journal of Transportation Systems Engineering and Information Technology, 2011, 11(1): 96-101. doi: 10.3969/j.issn.1009-6744.2011.01.016 SU Shuai, LI Xiang, TANG Tao, et al. A subway train timetable optimization approach based on energy-efficient operation strategy[J]. IEEE Transactions on Intelligent Transportation Systems, 2013, 14(2): 883-893. doi: 10.1109/TITS.2013.2244885 张惠茹,贾利民,王莉,等. 面向列车节能控制的时刻表优化[J]. 铁道学报,2019,41(2): 8-15.ZHANG Huiru, JIA Limin, WANG Li, et al. Study of timetable optimization based on train energy saving control[J]. Journal of the China Railway Society, 2019, 41(2): 8-15. SICRE C, CUCALAAP, FERNANDEZ A, et al. A method to optimise train energy consumption combining manual energy efficient driving and scheduling[J]. WIT Transactions on the Built Environment, 2010, 114: 549-560. CUCALA A P, FERNANDEZ A, SICRE C, et al. Fuzzy optimal schedule of high speed train operation to minimize energy consumption with uncertain delays and driver's behavioral response[J]. Engineering Applications of Artificial Intelligence, 2012, 25(8): 1548-1557. doi: 10.1016/j.engappai.2012.02.006 黄友能,宫少丰,曹源,等. 基于粒子群算法的城轨列车节能驾驶优化模型[J]. 交通运输工程学报,2016,16(2): 118-124. doi: 10.3969/j.issn.1671-1637.2016.02.014HUANG Youneng, GONG Shaofeng, CAO Yuan, et al. Optimization model of energy-efficient driving for train in urban rail transit based on particle swarm algorithm[J]. Journal of Traffic and Transportation Engineering, 2016, 16(2): 118-124. doi: 10.3969/j.issn.1671-1637.2016.02.014 陈荣武,诸昌钤,刘莉. 基于CBTC的城市轨道交通列车能耗算法及仿真[J]. 计算机应用研究,2011,28(6): 2126-2129. doi: 10.3969/j.issn.1001-3695.2011.06.034CHEN Rongwu, ZHU Changqian, LIU Li. CBTC based urban rail transit train energy consumption algorithm and simulation[J]. Application Research of Computers, 2011, 28(6): 2126-2129. doi: 10.3969/j.issn.1001-3695.2011.06.034 杨欣. 面向节能的城市轨道交通列车运行图优化研究[D]. 北京: 北京交通大学, 2016. HUANG Youneng, YANG Chen, GONG Shaofeng. Energy optimization for train operation based on an improved ant colony optimization methodology[J]. Energies, 2016, 9(8): 626. doi: 10.3390/en9080626 期刊类型引用(5)
1. 赵超,付斌,林立. 基于改进樽海鞘算法的含电动汽车微电网经济优化调度. 控制理论与应用. 2025(01): 167-180 . 百度学术
2. 马龙飞,张宝群,王立永,曾佳妮,焦然,宫成. 基于数字孪生混合储能的电动汽车参与微电网负荷功率波动平抑研究. 汽车工程. 2024(06): 1045-1053+1084 . 百度学术
3. 罗世刚,刘正英,李威武,许青,妥建军. 新能源接入交直流混合微电网的设备调度方法研究. 自动化仪表. 2024(06): 111-115 . 百度学术
4. 张懿璞,金雨哲,柯吉,茹锋,张成朔. 考虑电动汽车负荷的高速公路服务区离网微电网容量配置. 长安大学学报(自然科学版). 2024(05): 126-140 . 百度学术
5. 何玉朝,郭枫,钱辉敏,姚胜红,张能. 微电网源网荷储分布鲁棒优化调度方法研究. 电子设计工程. 2023(17): 123-127 . 百度学术
其他类型引用(4)
-