2.793

2018影响因子

(CJCR)

  • 中文核心
  • EI
  • 中国科技核心
  • Scopus
  • CSCD
  • 英国科学文摘

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于混合滤波最大期望算法的高速列车建模

王呈 陈晶 荀径 李开成

王呈, 陈晶, 荀径, 李开成. 基于混合滤波最大期望算法的高速列车建模. 自动化学报, 2019, 45(12): 2260−2267. doi: 10.16383/j.aas.c190193
引用本文: 王呈, 陈晶, 荀径, 李开成. 基于混合滤波最大期望算法的高速列车建模. 自动化学报, 2019, 45(12): 2260−2267. doi: 10.16383/j.aas.c190193
Wang Cheng, Chen Jing, Xun Jing, Li Kai-Cheng. Hybrid filter based expectation maximization algorithm for high-speed train modeling. Acta Automatica Sinica, 2019, 45(12): 2260−2267. doi: 10.16383/j.aas.c190193
Citation: Wang Cheng, Chen Jing, Xun Jing, Li Kai-Cheng. Hybrid filter based expectation maximization algorithm for high-speed train modeling. Acta Automatica Sinica, 2019, 45(12): 2260−2267. doi: 10.16383/j.aas.c190193

基于混合滤波最大期望算法的高速列车建模


DOI: 10.16383/j.aas.c190193
详细信息
    作者简介:

    江南大学物联网工程学院副教授. 2014年获得北京交通大学交通信息工程及控制专业博士学位. 主要研究方向为先进列车控制技术, 非线性系统建模与控制, 机器学习与数据挖掘. 本文通信作者. E-mail: wangc@jiangnan.edu.cn

    江南大学理学院副教授. 2013年获得江南大学控制理论与控制工程博士学位. 主要研究方向为系统辨识和过程控制. E-mail: chenjing1981929@126.com

    北京交通大学轨道交通控制与安全国家重点实验室副教授. 2012年获北京交通大学交通信息工程与控制博士学位. 主要研究方向为先进列车控制方法, 轨道交通优化问题, 交通流理论, 元胞自动机和加强学习. E-mail: jxun@bjtu.edu.cn

    北京交通大学轨道交通运行控制系统国家工程研究中心研究员. 主要研究方向为轨道交通列车运行控制, 智能控制技术及应用. E-mail: kchli@bjtu.edu.cn

  • 基金项目:  国家自然科学基金(61603156, 61973137), 高速铁路基础研究联合基金(U1734210), 北京交通大学教育基金会基金(9907006519)资助

Hybrid Filter Based Expectation Maximization Algorithm for High-speed Train Modeling

More Information
  • Fund Project:  Supported by National Natural Science Foundation of China (61603156, 61973137), the Joint Fund for Basic Research of High Speed Railway (U1734210), and Beijing Jiaotong University Education Foundation (9907006519)
  • 摘要: 针对高速列车非线性单质点模型的特殊结构及含有隐含变量问题, 提出一种基于混合滤波的最大期望辨识方法. 借助递阶辨识理论, 将高铁列车状态空间模型分解为线性子系统模型和非线性子系统模型. 进而, 分别利用卡尔曼滤波和粒子滤波对速度和位移状态进行联合估计. 最后, 使用最大期望方法辨识高铁列车子系统模型参数, 解决了隐含变量辨识问题. 和传统方法相比, 本文所提出方法计算量小, 且具有较高的辨识精度. 仿真对比实验结果验证了该方法的有效性.
  • 图  1  参数误差$\tau$$k$变化曲线

    Fig.  1  Parameter estimation error $\tau$ versus $k$

    图  2  位移估计变化曲线($+$: 估计位移; $-$: 真实位移)

    Fig.  2  Displacement estimation curve

    图  3  速度估计变化曲线($+$: 估计速度; $-$: 真实速度)

    Fig.  3  Velocity estimation curve

    表  1  模型参数的估计(混合滤波方法)

    Table  1  Parameters and their estimates (Hybrid filter)

    $k$$\bar{d}$$\bar{d}a$$\bar{d}b$$\bar{d}c$$\delta\ ({\text{%}})$
    50.009200120.004912630.000039910.000000980.33063
    80.009200150.004923160.000039920.000000980.37553
    100.009203140.004943210.000039910.000000970.59953
    200.009202870.004942950.000039870.000000970.58897
    300.009199720.004838220.000038460.000000990.74688
    真值0.009200000.004900000.000040000.00000100
    下载: 导出CSV

    表  2  模型参数的估计(拓展的卡尔曼滤波方法)

    Table  2  Parameters and their estimates (Extended Kalman filter)

    $k$$\bar{d}$$\bar{d}a$$\bar{d}b$$\bar{d}c$$\delta ({\text{%}})$
    50.009216320.005482360.000031240.000000865.52150
    80.009215940.005444690.000032620.000000895.19838
    100.009216210.005453210.000031730.000000915.22873
    200.009218450.005461510.000032920.000000915.39818
    300.009217060.005450760.000031890.000000915.31865
    真值0.009200000.004900000.000040000.00000100
    下载: 导出CSV

    表  3  模型参数的估计(粒子滤波方法)

    Table  3  Parameters and their estimates (Particle filter)

    $k$$\bar{d}$$\bar{d}a$$\bar{d}b$$\bar{d}c$$\delta\ ({\text{%}})$
    50.009176520.005021360.000038650.000000961.14182
    80.009191380.004987540.000039730.000000980.82781
    100.009191410.004987630.000040160.000000980.85850
    200.009191360.005016270.000040980.000000961.02914
    300.009191390.004900320.000040740.000000970.94999
    真值0.009200000.004900000.000040000.00000100
    下载: 导出CSV
  • [1] 衷路生, 颜真, 杨辉, 齐叶鹏, 张坤鹏, 樊晓平. 数据驱动的高速列车子空间预测控制方法. 铁道学报, 2013, 35(4): 77−83 doi:  10.3969/j.issn.1001-8360.2013.04.012

    1 Zhong Lu-Sheng, Yan Zhen, Yang Hui, Zhang Kun-Peng, Fan Xiao-Ping. Predictive control of high-speed train based on data driven subspace approach. Journal of the China Railway Society, 2013, 35(4): 77−83 doi:  10.3969/j.issn.1001-8360.2013.04.012
    [2] 2 Dong H R, lin X, Yao X M, Bai W Q, Ning B. Composite disturbance-observer-based control and H control for high speed trains with actuator faults. Asian Journal of Control, 2018, 20(2): 735−745 doi:  10.1002/asjc.1590
    [3] 王呈, 唐涛, 罗仁士. 基于UIO方法的列车自动驾驶信号故障检测研究. 铁道学报, 2013, 35(6): 48−52 doi:  10.3969/j.issn.1001-8360.2013.06.008

    3 Wang Cheng, Tang Tao, Luo Ren-Shi. ATO fault detection based on UIO method. Journal of the China Railway Society, 2013, 35(6): 48−52 doi:  10.3969/j.issn.1001-8360.2013.06.008
    [4] 4 Zhuan X, Xia X. Cruise control scheduling of heavy haul trains. IEEE Transactions on Control Systems Technology, 2006, 33(1): 757−766
    [5] 5 Xia X, Zhang J. Modeling and control of heavy haul trains. Control Systems IEEE, 2011, 31(4): 18−31 doi:  10.1109/MCS.2011.941403
    [6] 6 Chou M, Xia X, Kayser C. Modelling and model validation of heavy-haul trains equipped with electronically controlled pneumatic brake systems. Control Engineering Practice, 2007, 15(4): 501−509 doi:  10.1016/j.conengprac.2006.09.006
    [7] 7 Chen D W, Gao C H. Soft computing methods applied to train station parking in urban rail transit. Applied Soft Computing, 2012, 12(2): 759−767 doi:  10.1016/j.asoc.2011.10.016
    [8] 辛斌, 白永强, 陈杰. 基于偏差消除最小二乘估计和Durbin方法的两阶段ARMAX参数辨识. 自动化学报, 2012, 38(3): 491−496

    8 Xin Bin, Bai Yong-Qiang, Chen Jie. Two-stage ARMAX parameter identification based on bias-eliminated least squares estimation and Durbin′s method. Acta Automatica Sinica, 2012, 38(3): 491−496
    [9] 9 Chen J, Jiang B. Modified stochastic gradient parameter estimation algorithms for a nonlinear two-variable difference system. International Journal of Control, Automation, and Systems, 2016, 14(6): 1493−1500 doi:  10.1007/s12555-015-0185-x
    [10] 10 Wang C, Li K C. Aitken-based stochastic gradient algorithm for ARX models with time delay. Circuits Systems and Signal Processing, 2019, 38(6): 2863−2876 doi:  10.1007/s00034-018-0998-y
    [11] 11 Li J H, Zheng W X, Gu J P, Hua L. Parameter estimation algorithms for Hammerstein output error systems using Levenberg-Marquardt optimization method with varying interval measurements. Journal of the Franklin Institute, 2017, 354(1): 316−331 doi:  10.1016/j.jfranklin.2016.10.002
    [12] 12 Xia H F, Yang Y Q, Ding F, Alsaedi A, Hayat T. Maximum likelihood-based recursive least-squares estimation for multivariable systems using the data filtering technique. International Journal of Systems Science, 2019, 50(6): 1121−1135 doi:  10.1080/00207721.2019.1590664
    [13] 13 Yang X Q, Yin S. Robust global identification and output estimation for LPV dual-rate systems subjected to random output time-delays. IEEE Transactions on Industrial Informatics, 2017, 13(6): 2876−2885 doi:  10.1109/TII.2017.2702754
    [14] 14 Chen J, Huang B, Ding F, Gu Y. Variational Bayesian approach for ARX systems with missing observations and varying time-delays. Automatica, 2018, 94: 194−204 doi:  10.1016/j.automatica.2018.04.003
    [15] 15 Xiong W L, Yang X Q, Ke L, Xu B G. EM algorithm-based identification of a class of nonlinear Wiener systems with missing output data. Nonlinear Dynamics, 2015, 80(1): 329−339
    [16] 16 Zhao Y J, Fatehi A, Huang B. Robust estimation of ARX models with time-varying time delays using variational Bayesian approach. IEEE Transactions on Cybernetics, 2018, 48(2): 532−542 doi:  10.1109/TCYB.2016.2646059
    [17] 衷路生, 李兵, 龚锦红, 张永贤, 祝振敏. 高速列车非线性模型的极大似然辨识. 自动化学报, 2014, 40(12): 2950−2958

    17 Zhong Lu-Sheng, Li Bing, Gong Jin-Hong, Zhang Yong-Xian, Zhu Zhen-Min. Maximum likelihood identification of nonlinear models for high-speed train. Acta Automatica Sinica, 2014, 40(12): 2950−2958
    [18] Simon D. Optimal State Estimation: Kalman, H, and Nonlinear Approaches. John Wiley & Sons, Inc., New Jersey, 2006.
    [19] 19 Simon D, Shmaliy Y S. Unified forms for Kalman and finite impulse response filtering and smoothing. Automatica, 2013, 49(6): 1892−1899 doi:  10.1016/j.automatica.2013.02.026
    [20] 20 Chen J, Li J, Liu Y J. Gradient iterative algorithm for dual-rate nonlinear systems based on a novel particle filter. Journal of the Franklin Institute, 2017, 354(11): 4425−4437 doi:  10.1016/j.jfranklin.2017.04.003
    [21] 顾晓清, 倪彤光, 张聪, 戴臣超, 王洪元. 结构辨识和参数优化协同学习的概率TSK模糊系统. 自动化学报, 2019. DOI 10.16383/j.aas.c180298

    Gu Xiao-Qing, Ni Tong-Guang, Zhang Cong, Dai Chen-Chao, Wang Hong-Yuan. Probabilistic TSK fuzzy system in the simulataneous learning of structure identification and parameter optimization. Acta Automatica Sinica, 2019. DOI 10.16383/j.aas.c180298
  • [1] 谢国, 金永泽, 黑新宏, 姬文江, 高士根, 高桥圣, 望月宽. 列车动力学模型时变环境参数自适应辨识[J]. 自动化学报, 2019, 45(12): 2268-2280. doi: 10.16383/j.aas.c190215
    [2] 陆耿虹, 冯冬芹. 基于粒子滤波的工业控制网络态势感知建模[J]. 自动化学报, 2018, 44(8): 1405-1412. doi: 10.16383/j.aas.2017.c160830
    [3] 张继文, 梁桐, 张淑平. 实验小鼠运动参数的模板匹配及粒子滤波提取方法[J]. 自动化学报, 2018, 44(1): 25-34. doi: 10.16383/j.aas.2018.c160573
    [4] 江涛, 钱富才, 杨恒占, 胡绍林. 具有双重不确定性系统的联合滤波算法[J]. 自动化学报, 2016, 42(4): 535-544. doi: 10.16383/j.aas.2016.c150486
    [5] 田梦楚, 薄煜明, 陈志敏, 吴盘龙, 赵高鹏. 萤火虫算法智能优化粒子滤波[J]. 自动化学报, 2016, 42(1): 89-97. doi: 10.16383/j.aas.2016.c150221
    [6] 赵光琼, 陈绍刚, 付奎, 唐忠樑, 贺威. 基于变尺度变换减少Sigma点的粒子滤波算法研究[J]. 自动化学报, 2015, 41(7): 1350-1355. doi: 10.16383/j.aas.2015.c140833
    [7] 杜党波, 张伟, 胡昌华, 周志杰, 司小胜, 张建勋. 含缺失数据的小波-卡尔曼滤波故障预测方法[J]. 自动化学报, 2014, 40(10): 2115-2125. doi: 10.3724/SP.J.1004.2014.02115
    [8] 宋宇, 李庆玲, 康轶非, 闫德立. 平方根容积Rao-Blackwillised粒子滤波SLAM算法[J]. 自动化学报, 2014, 40(2): 357-367. doi: 10.3724/SP.J.1004.2014.00357
    [9] 周振威, 方海涛. 在不准确方差下带随机系数矩阵的卡尔曼滤波稳定性[J]. 自动化学报, 2013, 39(1): 43-52. doi: 10.3724/SP.J.1004.2013.00043
    [10] 欧阳成, 姬红兵, 郭志强. 改进的多模型粒子PHD和CPHD滤波算法[J]. 自动化学报, 2012, 38(3): 341-348. doi: 10.3724/SP.J.1004.2012.00341
    [11] 甘敏, 彭辉, 黄云志, 董学平. 自组织状态空间模型参数初始分布搜索算法[J]. 自动化学报, 2012, 38(9): 1538-1543. doi: 10.3724/SP.J.1004.2012.01538
    [12] 王帅, 杨文, 侍洪波. 带丢包一致性滤波算法研究[J]. 自动化学报, 2010, 36(12): 1689-1696. doi: 10.3724/SP.J.1004.2010.01689
    [13] 宋宇, 孙富春, 李庆玲. 移动机器人的改进无迹粒子滤波蒙特卡罗定位算法[J]. 自动化学报, 2010, 36(6): 851-857. doi: 10.3724/SP.J.1004.2010.00851
    [14] 赵玲玲, 马培军, 苏小红. 一种快速准蒙特卡罗粒子滤波算法[J]. 自动化学报, 2010, 36(9): 1351-1356. doi: 10.3724/SP.J.1004.2010.01351
    [15] 胡昌华, 张琪, 乔玉坤. 强跟踪粒子滤波算法及其在故障预报中的应用[J]. 自动化学报, 2008, 34(12): 1522-1528. doi: 10.3724/SP.J.1004.2008.01522
    [16] 段琢华, 蔡自兴, 于金霞. 不完备多模型混合系统故障诊断的粒子滤波算法[J]. 自动化学报, 2008, 34(5): 581-587. doi: 10.3724/SP.J.1004.2008.00581
    [17] 安德玺, 梁彦, 周东华. 一种基于滤波参数在线辨识的鲁棒自适应滤波器[J]. 自动化学报, 2004, 30(4): 560-566.
    [18] 丁锋, 陈通文, 萧德云. 一般双率随机系统状态空间模型及其辨识[J]. 自动化学报, 2004, 30(5): 652-663.
    [19] 丁锋, 杨家本. 大系统的递阶辨识[J]. 自动化学报, 1999, 25(5): 647-654.
    [20] 南志远, 王瑞申. 奇异摄动型卡尔曼滤波算法及其在互联电力系统负荷频率控制中的应用[J]. 自动化学报, 1990, 16(5): 385-392.
  • 加载中
图(3) / 表(3)
计量
  • 文章访问数:  616
  • HTML全文浏览量:  134
  • PDF下载量:  34
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-03-19
  • 录用日期:  2019-08-15
  • 刊出日期:  2019-12-01

基于混合滤波最大期望算法的高速列车建模

doi: 10.16383/j.aas.c190193
    作者简介:

    江南大学物联网工程学院副教授. 2014年获得北京交通大学交通信息工程及控制专业博士学位. 主要研究方向为先进列车控制技术, 非线性系统建模与控制, 机器学习与数据挖掘. 本文通信作者. E-mail: wangc@jiangnan.edu.cn

    江南大学理学院副教授. 2013年获得江南大学控制理论与控制工程博士学位. 主要研究方向为系统辨识和过程控制. E-mail: chenjing1981929@126.com

    北京交通大学轨道交通控制与安全国家重点实验室副教授. 2012年获北京交通大学交通信息工程与控制博士学位. 主要研究方向为先进列车控制方法, 轨道交通优化问题, 交通流理论, 元胞自动机和加强学习. E-mail: jxun@bjtu.edu.cn

    北京交通大学轨道交通运行控制系统国家工程研究中心研究员. 主要研究方向为轨道交通列车运行控制, 智能控制技术及应用. E-mail: kchli@bjtu.edu.cn

基金项目:  国家自然科学基金(61603156, 61973137), 高速铁路基础研究联合基金(U1734210), 北京交通大学教育基金会基金(9907006519)资助

摘要: 针对高速列车非线性单质点模型的特殊结构及含有隐含变量问题, 提出一种基于混合滤波的最大期望辨识方法. 借助递阶辨识理论, 将高铁列车状态空间模型分解为线性子系统模型和非线性子系统模型. 进而, 分别利用卡尔曼滤波和粒子滤波对速度和位移状态进行联合估计. 最后, 使用最大期望方法辨识高铁列车子系统模型参数, 解决了隐含变量辨识问题. 和传统方法相比, 本文所提出方法计算量小, 且具有较高的辨识精度. 仿真对比实验结果验证了该方法的有效性.

English Abstract

  • 随着高速列车进一步提速, 高铁已成为当今国人出行首选. 高铁列控系统是保障列车安全、精确、节能和舒适运行的中枢大脑. 为进一步提升安全性与可靠性, 满足对各类恶劣环境的适应性, 需要对高速列车进行控制和故障检测与诊断[1-3]. 而控制器的设计以及故障的检测与诊断离不开高速列车模型. 因此对高速列车模型进行建模并获得精确模型是开展控制与故障检测研究的基础.

    列车模型可分为多质点模型和单质点模型. 多质点模型最早用于刻画重载列车制动动态. 原因在于重载列车的主风管在排气减压时, 各车厢制动性能无法完全一致, 使得某些载货车厢会产生前冲动作, 因此建模时必须考虑车钩内力及其耦合作用[4-6]. 针对模型参数辨识和滤波估计, 目前, 已出现了许多方法, 如: 最小二乘算法[7-8]、随机梯度算法[9-10]、极大似然辨识算法等[11-12]. 其中, 极大似然算法采用概率的手段来对模型进行辨识, 其核心思想是: 假设有若干参数, 每个参数都对应固定样本出现的概率, 其中最大概率对应参数即为估计的真实值. 极大似然方法存在所有样本需可测的前提, 当部分样本不可测时, 该方法将无能为力.

    最大期望(Expectation maximization, EM)方法作为极大似然方法的一个拓展, 可以处理样本中存在隐含变量的问题. EM方法分为两步, 即E步和M步. 通过迭代计算, 在E步, 估计系统未知变量(隐含变量); 在M步, 辨识系统参数, 两步交替进行[13-14]. 如, Xiong等利用EM方法对具有损失数据的非线性系统提出了辨识方法[15], 在E步用辅助模型方法辨识系统未知输出, 在M步辨识系统参数. Zhao等针对具有马尔科夫时延的ARX模型, 设计EM辨识算法, 在E步辨识出系统的未知时延, 在M步利用辨识出的时延估计系统参数[16]. 最近, 衷路生等借助粒子滤波思想, 利用EM算法对具有非线性特性的高速列车模型进行参数辨识研究, 开创了从概率角度研究高速列车建模的新思路[17], 本文在文献[17]的基础上, 进一步针对非线性高速列车模型特征, 提出了一种基于混合滤波的EM辨识算法.

    卡尔曼滤波估计器能对线性状态空间模型进行优化估计, 当模型强非线性时, 粒子滤波可以取代卡尔曼滤波, 但计算量较大[18-21]. 递阶辨识方法能很好地解决结构复杂的耦合系统辨识问题, 从而减少算法计算量. 其思想是通过模型分解简化辨识问题复杂性, 从而减小辨识计算量. 本文借助递阶辨识思想, 将高速列车模型分为两个子模型. 其一是关于位移的线性状态空间子模型, 其二是关于速度的非线性状态空间子模型. 针对线性模型利用卡尔曼滤波方法辨识位移变量, 而对非线性速度状态空间模型, 利用粒子滤波估计速度变量. 最后借助EM方法辨识系统参数. 本文的主要贡献可以归纳如下:

    1)单质点高速列车模型具有其特殊结构, 借助递阶辨识思想和EM算法, 建立了针对线性与非线性耦合模型的混合滤波辨识理论框架.

    2)提出的混合滤波方法, 解决了扩展卡尔曼滤波处理非线性辨识精度较低, 以及粒子滤波计算量大的问题, 提高了模型的辨识效率.

    3)提出的混合滤波理论框架及方法, 为精确列车自动驾驶控制和面向节能的优化驾驶控制提供模型基础, 并可以进一步拓展到多工况条件下的高速列车多模型辨识和重载列车多质点模型辨识研究.

    本文结构组织如下: 第1节重点介绍高速列车状态空间模型和EM算法框架, 第2节提出基于混合滤波的EM辨识方法, 第3节通过仿真例子验证提出方法的有效性, 最后在第4节给出结论和未来研究方向.

    • 考虑如下高速列车单质点模型:

      $$ \frac{{\rm d}s}{{\rm d}t} = v \hspace{60pt}$$ (1)
      $$ \frac{{\rm d}v}{{\rm d}t} = \varepsilon(f(v)-w(v)) $$ (2)

      其中, $ s $表示列车位移, $ v $表示运行速度, $ \varepsilon = {0.0098}/ $$({1+d}) $($ d $为列车回转质量系数)表示列车加速度系数, $ f(v) $$ w(v) $分别为列车的单位动力和单位运行阻力. 将高铁连续模型转换为离散模型

      $$ \left[\begin{array}{l} s(t)\\ v(t) \end{array}\right] = \left[\begin{array}{l} s(t-1)+Tv(t-1)\\ v(t-1)+\dfrac{0.0098T}{1+d}\varphi(t) \end{array}\right] \hspace{9pt}$$ (3)
      $$ \varphi(t) = u(t)-(a+bv(t-1)+cv^2(t-1))\\ $$ (4)

      其中, $ u(t) = f(v(t-1)) $可测, 不失一般性, 高铁离散模型可进一步表示为如下状态空间模型

      $$ \begin{split} \left[\begin{array}{l} s(t)\\v(t)\end{array}\right] = \;& \left[\begin{array}{l} s(t-1)+Tv(t-1)\\v(t-1)+\dfrac{0.0098T}{1+d}\varphi(t)\end{array}\right]+\\ & \left[\begin{array}{l} w_1(t)\\w_2(t)\end{array}\right] \end{split}$$ (5)
      $$ y_1(t) = s(t)+e_1(t) \hspace{85pt}$$ (6)
      $$ y_2(t) = v(t)+e_2(t) \hspace{85pt}$$ (7)

      其中, $ w_1(t), w_2(t) $为系统的过程噪声, 满足均值为零方差为$ \delta_1^2 $$ \delta_2^2 $的高斯分布, $ e_1(t) $$ e_2(t) $分别是系统的位移和速度测量噪声, 满足均值为零, 方差为$ \sigma_1^2 $$ \sigma_2^2 $的高斯分布, 且$ \delta_i^2 $$ \sigma_i^2,\; i = 1,2 $已知. 本文目标是利用可测输出$ y_1(t) $, $ y_2(t) $以及输入$ u(t) $, 辨识系统真实状态$ v(t) $, $ s(t) $以及未知参数.

      目前, 状态空间模型辨识已相对成熟, 其中应用最为广泛的是将状态空间模型转化为输入输出模型, 进而利用自回归模型辨识方法来辨识系统参数. 然而, 本文所述系统状态空间模型是非线性的, 很难转化为输入输出模型, 因而常规方法难以适用高速列车模型建模. 定义

      $$ X(t) = [s(t),v(t)]^{\rm T} $$ (8)

      由于量测噪声存在, 式(8)真实模型状态$ s(t) $$ v(t) $不可测, 即辨识过程中存在隐含变量, 而最大期望(EM)方法在处理具有隐含变量的辨识问题中具有独特优势. 本文即利用EM方法来对高铁模型进行参数辨识. 定义参数向量为

      $$ {{\Theta}} = [a,b,c,d]^{\rm T} $$ (9)

      设共采集$ n $组数据, 则系统的可测变量$ C_{\rm obs} $

      $$\begin{split} C_{\rm obs} =\;& \{y_1(1),\cdots,y_1(n),y_2(1),\cdots,y_2(n),\\ &u(1),\cdots,u(n)\} \end{split} $$ (10)

      系统的隐含变量$ C_{\rm mis} $

      $$ C_{\rm mis} = \{v(1),\cdots,v(n),s(1),\cdots,s(n)\} $$ (11)

      $ {\rm EM} $算法分为$ {\rm E} $步和$ {\rm M} $步, 步骤如下:

      1. 初始化: 给定一个初始的系统参数$ {{\Theta}}_{{{0}}} $.

      2. E步: 根据辨识出的参数$ {{\Theta}}_{{k}} $, 计算下列Q函数

      $$ Q({{\Theta}}|{{\Theta}}_k) = E{_{C_{{\rm mis}}|C_{\rm obs},{{\Theta}}_k}{\{{\rm ln} p(C_{\rm obs},C_{\rm mis}|{{\Theta}}})\} }$$

      3. M步: 极大化Q函数, 求得使Q函数最大的参数$ {{\Theta}}_{k+1} $

      $$ {{\Theta}}_{k+1} = {\rm arg} {\rm max} \;Q({{\Theta}}|{{\Theta}}_k) $$
    • EM算法通过极大化下列概率密度函数来求解未知参数

      $$ p(C_{\rm obs}|{{\Theta}}) $$

      将上式转换为

      $$ p(C_{\rm obs}|{{\Theta}}) = \int p(C_{\rm obs}, C_{\rm mis}|{{\Theta}}){\rm d}C_{\rm mis} $$

      等式两边取对数

      $$ {\rm ln} p(C_{\rm obs}|{{\Theta}}) = {\rm{ ln}} \int p(C_{\rm obs}, C_{\rm mis}|{{\Theta}}){\rm d}C_{\rm mis} $$

      引入一个新的密度函数$ q(C_{\rm mis}|C_{\rm obs},{{\Theta}}) $

      $$ \begin{split} &{\rm ln} p(C_{\rm obs}|{{\Theta}}) = \\ &\qquad{\rm ln} \int q(C_{\rm mis}|C_{\rm obs},{{\Theta}})\frac{p(C_{\rm obs}, C_{\rm mis}|{{\Theta}})}{q(C_{\rm mis}|C_{\rm obs},{{\Theta}})}{\rm d}C_{\rm mis} \end{split} $$ (12)

      根据Jensen不等式, 可得

      $$ \begin{split} & {\rm ln} p(C_{\rm obs}|{{\Theta}})\ge \\ &\qquad \int q(C_{\rm mis}|C_{\rm obs},{{\Theta}}){\rm ln}\frac{p(C_{\rm obs}, C_{\rm mis}|{{\Theta}})}{q(C_{\rm mis}|C_{\rm obs},{{\Theta}})}{\rm d}C_{\rm mis} \end{split}$$ (13)

      等式成立的条件是

      $$ q(C_{\rm mis}|C_{\rm obs},{{\Theta}}) = p(C_{\rm obs}, C_{\rm mis}|{{\Theta}}) $$

      其中, 后验概率

      $$ \begin{split}& p(C_{\rm obs},C_{\rm mis}|{{\Theta}}) = \\ &\qquad\frac{p(C_{\rm obs}|C_{\rm mis},{{\Theta}})p(C_{\rm mis}|{{\Theta}})}{\int p(C_{\rm obs}|C_{\rm mis},{{\Theta}})p(C_{\rm mis}|{{\Theta}}){\rm d}C_{\rm mis}} \end{split} $$ (14)

      然而, 上述后验概率中$ p(C_{\rm obs}|C_{\rm mis},{{\Theta}}) $含有未知变量$ C_{\rm mis} $, 导致其无法计算. 重新将$ Q $函数表示为

      $$ \begin{split} Q({{\Theta}}&|{{\Theta}}_k) = E{{_{C_{{\rm mis}}|C_{\rm obs},{{\Theta}}_k}}\{{\rm ln} p(C_{\rm obs},C_{\rm mis}|{{\Theta}}})\}=\\ & \int p(C_{\rm mis}|C_{\rm obs},{{\Theta}}_k){\rm ln} p(C_{\rm obs},C_{\rm mis}|{{\Theta}}){\rm d}C_{\rm mis}=\\ & \int p(C_{\rm mis}|C_{\rm obs},{{\Theta}}_k)[{\rm ln} p(C_{\rm mis}|C_{\rm obs},{{\Theta}})+\\ & {\rm ln} p(C_{\rm obs}|{{\Theta}})]{\rm d}C_{\rm mis}=\\ & \int p(C_{\rm mis}|C_{\rm obs},{{\Theta}}_k){\rm ln}p(C_{\rm mis}|C_{\rm obs},{{\Theta}}){\rm d}C_{\rm mis}+\\ &{\rm ln}C=\\ & \sum\limits_{t = 1}^n\int p(x(t)|,x(t-1),{{\Theta}}_k)\times \\ & {\rm ln} p(x(t)|x(t-1),{{\Theta}}){\rm d}C_{\rm mis}+{\rm ln}C \end{split}$$

      显然, 由于含有隐含变量$ x(t-1) $, 使得上式积分无法计算. 假设在第$ k-1 $次的参数已经获得, 则可以通过该参数估计和可测$ y,u $来估计所有隐含变量.

    • 将高速列车模型重写为

      $$ \begin{split} \left[\begin{array}{l} s(t)\\v(t) \end{array}\right] =\;& \left[\begin{array}{l} s(t-1)+Tv(t-1)\\v(t-1)+\dfrac{0.0098T}{1+\hat{d}_{k-1}}{\varphi}(t) \end{array}\right] +\\ & \left[\begin{array}{l} w_1(t)\\w_2(t) \end{array}\right] \end{split}$$ (15)

      由于式(15)是线性项和非线性项的组合, 直接采用卡尔曼滤波方法无法估计隐含变量. 若使用文献[14]的方法, 即粒子滤波算法, 则由于隐含变量是二维的, 会导致计算量增大. 此外, 由于位移模型为线性模型, 若使用粒子滤波会影响计算精度. 本文借助递阶辨识思想, 将上述状态空间模型分解为两个子系统. 线性子系统为

      $$ s(t) = s(t-1)+Tv(t-1)+w_1(t) $$ (16)
      $$ y_1(t) = s(t)+e_1(t) \hspace{67pt}$$ (17)

      非线性子系统为

      $$ v(t) = v(t-1)+\frac{0.0098T}{1+d_{k-1}}\varphi(t)+w_2(t) $$ (18)
      $$ y_2(t) = v(t)+e_2(t) \hspace{84pt}$$ (19)

      假设在第$ k $次, 参数估计为$ {{\Theta}}_{k-1} = [a_{k-1},b_{k-1},$$c_{k-1},d_{k-1}] ^{\rm T} $. 对线性子系统, 利用卡尔曼滤波方法估计最优隐含变量$ s(t) $, 设为$ \bar{s}(t) $, 并用非线性子系统估计出的最优$ \bar{v}(t-1) $来代替$ v(t-1) $. 得到关于隐含变量$ s(t) $估计的卡尔曼滤波算法如下:

      $$ \bar{s}(t) = \hat{s}(t)+K(t)\left[y_1(t)-\hat{s}(t)\right] $$ (20)
      $$ \hat{s}(t) = \bar{s}(t-1)+T\bar{v}(t-1) \hspace{20pt}$$ (21)
      $$ K(t) = P^{-}(t)[P^{-1}(t)+\sigma^2_1]^{-1} \hspace{10pt} $$ (22)
      $$ P^{-}(t) = P^{+}(t-1)+\delta^2_1 \hspace{35pt}$$ (23)
      $$ P^{+}(t) = [I-K(t)]P^{-}(t) \hspace{32pt} $$ (24)

      由于关于$ v(t) $的状态模型是非线性的, 所以卡尔曼滤波方法无法直接应用. 本文主要通过粒子滤波来估计状态$ v(t) $, 与文献[17]相比, 由于速度模型是一维的, 因而计算量会大大减小. 根据粒子滤波理论, 在第$ k $次, 第$ t $时刻, 状态$ v(t) $可由如下等式近似

      $$ \bar{v}(t) = \beta_i\hat{v}_i(t),\; i = 1,\cdots,n $$ (25)

      其中, $ \hat{v}_i(t) $是第$ i $个粒子的预测值, 并由以下状态方程产生

      $$ v(t) = v(t-1)+\frac{0.0098T}{1+d}\varphi(t)+w_2(t) $$

      $ \beta_i $是粒子对应的权值, 由输出方程产生

      $$ y_2(t) = v(t)+e_2(t) $$

      假设第$ i $个粒子的预测值为$ \hat{v}_i(t) $, 则

      $$ \begin{split} \hat{v}_i(t) =\;& \hat{v}_i(t-1)+\frac{0.0098T}{1+d_{k-1}}\times[u(t)-(a_{k-1}+\\ &{b_{k-1}\hat{v}_i(t-1)+c_{k-1}\hat{v}^2_i(t-1))]} \end{split} $$ (26)

      由于$ e_2(t) $满足均值为零, 方差为$ \sigma_2^2 $的高斯分布, 则根据粒子滤波理论可得

      $$ \begin{split}\hat{\beta}_i =\;& p(y_2(t)|\hat{v}_i(t))=\\& \frac{1}{\sqrt{2{\text{π}} }\sigma_2}\exp\left(\frac{[y_2(t)-\hat{v}_i(t)]^2}{2{\delta_2^2}}\right) \end{split} $$ (27)

      对上述$ \hat{\beta}_i $归一化, 得

      $$ {\beta}_i = \frac{\hat{\beta}_i}{\displaystyle\sum\limits_{j = 1}^{n}\hat{\beta}_i} $$ (28)

      则隐含变量的卡尔曼/粒子滤波, 即混合滤波(Hybrid filter)联合估计算法总结如下:

      1. 假设在第$ k $次, 已获得参数的估计为$ {{\Theta}}_{k-1} $, 给定初始$ \bar{s}(0) = 0 $, $ P^{+}(0) = 10 $.

      2. $ n $个初始粒子为$ v_i(0) = 0.2randn\,(n,1) $, $ i = $$ 1,\cdots,n $, $randn (n, 1)$为生成正态分布随机数的函数, 并为每一个粒子赋初始权值为$ 1/n $.

      3. 根据式(25)计算$ \bar{v}(0) $. 初始取$ t = 1 $.

      卡尔曼滤波:

      4. 根据式(21)计算$ \hat{s}(t) $.

      5. 由式(23)计算$ P^{-}(t) $.

      6. 根据式(22)计算$ K(t) $.

      7. 由式(24)和式(20)分别计算$ P^{+}(t) $$ \bar{s}(t) $.

      粒子滤波:

      8. 根据式(26)计算$ \hat{v}_i(t) $, $ i = 1,\cdots,n $.

      9. 由式(27)计算$ \hat{\beta}_i $, $ i = 1,\cdots,n $.

      10. 根据式(28)计算归一化权值$ \bar{\beta}_i $, $ i = $$ 1,\cdots,n $.

      11. 由式(25)计算$ \bar{v}(t) $.

      12. 使$ t = t+1 $, 如果$ t<n $则转步骤4; 否则结束循环.

      注 1. 与文献[14]相比, 本文将复杂的二维高速列车单质点模型分解为相对简单的一维线性位移模型和非线性速度模型, 并分别对每个一维模型单独进行辨识, 因此本文所提出的方法计算量较小.

    • $ Q $函数可简化为

      $$ \begin{split}&{E{{_{C_{{\rm mis}}|C_{\rm obs},{{\Theta}}_{k-1}}}}\{{\rm ln} p(C_{\rm obs},C_{\rm mis}|{{\Theta}}})\}=\\ &\qquad \sum\limits_{t = 1}^n\int p(v(t)|v(t-1),{{\Theta}}_{k-1})\times\\ &\qquad {\rm ln} p(v(t)|v(t-1),{{\Theta}}){\rm d}C_{\rm mis}+{\rm ln}C=\\ &\qquad -\sum\limits_{t = 1}^n\frac{1}{2}{\rm ln} (2{\text{π}} \delta^2_2)-\sum\limits_{t = 1}^n\frac{1}{2\delta^2_2}\times\\ &\qquad \int p(v(t)|v(t-1),{{\Theta}}_{k-1})\times (v(t)-\psi(t))^2{\rm d}v(t)\end{split} $$ (29)

      由概率论中的二阶矩定义可得

      $$ \begin{split} & \int p(v(t)|v(t-1),{{\Theta}}_{k-1})\times (v(t)-\psi(t))^2{\rm d}v(t)=\\ &\qquad \delta^2_2+\psi^2_{k-1}(t)-2\psi^2_{k-1}(t)\psi(t)+\psi^2(t)=\\ &\qquad \delta^2_2+[\psi(t)-\psi_{k-1}(t)]^2 \\[-12pt] \end{split}$$ (30)

      其中

      $$ \begin{split} {\psi _{k - 1}}(t) =\;& {{\bar v}_{k - 1}}(t - 1) + \frac{{0.0098T}}{{1 + d}} \times \\ &{\left[u(t) - (a + b{{\bar v}_{k - 1}}(t - 1) + c\bar v_{k - 1}^2(t - 1))\right]} \end{split} $$

      将式(30)代入式(29), 并分别对$ a $, $ b $, $ c $, $ d $求偏导可得参数$ a $, $ b $, $ c $, $ d $的估计. 为了进一步解决由于参数$ d $在分母上, 导致求偏导时计算量大的问题, 定义

      $$ \phi_k(t) = \bar{v}_{k-1}(t)-\bar{v}_{k-1}(t-1) $$ (31)
      $$ \frac{0.0098T}{1+d} = \bar{d} \hspace{70pt}$$ (32)
      $$ \begin{split} \Phi_k(t) =\;& [u(t),-1,-\bar{v}_k(t-1),\\ &-\bar{v}^2_k(t-1)]^{\rm T} \end{split}$$ (33)
      $$ \vartheta = [\bar{d},\bar{d}a,\bar{d}b,\bar{d}c]^{\rm T} \hspace{46pt}$$ (34)

      则式(29)可简化为

      $$ \begin{split} &{E{_{C_{{\rm mis}}|C_{\rm obs},{{\Theta}}_{k-1}}}}\{{\rm ln} p(C_{\rm obs},C_{\rm mis}|{{\Theta}})\}=\\ &\qquad \sum\limits_{t = 1}^n-\frac{1}{2}{\rm ln} (2{\text{π}}\delta^2_2)-\sum\limits_{t = 1}^n\frac{1}{2\delta^2_2}\times\\ &\qquad [\delta^2_2+(\phi_k(t)-\Phi^{\rm T}_k(t)\vartheta)^2] \end{split}$$ (35)

      极大化式(35), 得到第$ k $次的参数估计

      $$ \vartheta_k = \left[{\sum\limits_{t = 1}^n}\Phi_k(t)\Phi^{\rm T}_k(t)\right]^{-1}\left[{\sum\limits_{t = 1}^n}\Phi_k(t)\phi_k(t)\right] $$ (36)

      一旦计算得到$ \vartheta_k $, 则可根据式(32)和式(34)获得参数估计$ \{a_k,b_k,c_k,d_k\} $.

      基于混合滤波EM算法的参数和隐含变量计算步骤如下:

      1. 初始化$ u(-j) = 0,v(-j) = 0,s(-j) = 0 $, $ j = $$ 0,1,2,\cdots,n-1 $, 并给定一个小的正常数$ \varepsilon $.

      2. 定义$ {{\Theta}}_0 = [0,0,0,0]^{\rm T} $, 并令$ k = 1 $.

      3. 分别收集输入和输出数据$ y_1(t),y_2(t), u(t),$$ t = 1,2, \cdots, n $.

      4. 利用卡尔曼滤波和粒子滤波分别计算隐含变量$ \bar{s}(t),\bar{v}(t) $.

      5. 根据式(31)计算$ \phi_k(t) $以及式(33)构建$ \Phi_k(t) $.

      6. 根据式(36)更新参数$ \vartheta_k $, 由$ \vartheta_k $计算获得$ {{\Theta}}_k $.

      7. 比较$ {{\Theta}}_k $$ {{\Theta}}_{k-1} $: 如果$ \|{{\Theta}}_k-{{\Theta}}_{k-1}\| \le\varepsilon $, 则得到最终的参数估计$ {{\Theta}}_k $; 否则, 令$ k = k+1 $并返回步骤4.

      注 2. 针对线性位移模型利用卡尔曼滤波估计未知位移, 非线性速度模型采用粒子滤波估计未知速度. 由于卡尔曼滤波是线性状态空间模型最优估计器, 且计算量远小于粒子滤波. 因此本文所采用的混合滤波方法具有计算量小、精度高的优点.

    • 考虑CRH3高速列车动车组(四动四拖)进行仿真实验, 则高速列车运行控制模型可以写为

      $$ \begin{split} \left[\begin{array}{l} s(t)\\ v(t) \end{array}\right] & = \left[\begin{array}{l} s(t-1)+v(t-1)\\ v(t-1)+\frac{0.0098}{1+0.06}\varphi(t) \end{array}\right]+\\ & \quad \; \left[\begin{array}{l} w_1(t)\\ w_2(t) \end{array}\right]\\ y_1(t) &= s(t)+e_1(t)\\ y_2(t) &= v(t)+e_2(t)\\ \varphi(t) &= u(t)- [0.53+0.0039v(t-1)+\\& \;\;\;\;\;\; 0.000114v^2(t-1)] \end{split}$$

      将系统分解为如下两个子模型, 关于位移的线性子模型为

      $$ s(t) = s(t-1)+v(t-1)+w_1(t) $$ (37)
      $$ y_1(t) = s(t)+e_1(t) \hspace{60pt}$$ (38)

      关于速度的非线性子模型为

      $$ \begin{split} v(t)-\;&v(t-1) = 0.0092u(t)-0.0049-\\ &0.00004v(t-1)-0.000001v(t-1)+w_2(t) \end{split}$$ (39)
      $$ y_2(t) = v(t)+e_2(t) \hspace{110pt}$$ (40)

      定义

      $$ \begin{split} \;\phi(t) &= v(t)-v(t-1)\\ \Phi(t) &= [u(t),-1,-v(t-1),-v^2(t-1)]^{\rm T}\\ { \vartheta} &= [\bar{d},\bar{d}a,\bar{d}b,\bar{d}c]^{\rm T}=\\ &\;\;\;\;\;[0.0092,0.0049,0.00004,0.000001]^{\rm T}.\end{split} $$

      所有噪声均采用零均值方差为$ 0.10^2 $的白噪声序列.

      首先采用本文所提出的混合滤波最大期望方法辨识模型参数、真实位移和速度. 参数估计误差$ \tau: = \|{ \vartheta}_k-{ \vartheta}\|/\|{ \vartheta}\| $随迭代次数$ k $变化的曲线如图1表1所示. 由图1表1可知, 参数估计值在第5次就已经逼近真实值, 且估计误差从第5次迭代往后基本控制在2 %以内, 与文献[14]相比, 本文的参数收敛速度和精度更有优势. 为避免对分母求导, 本文采用简化的参数模型, 故辨识参数并非原始参数. 为此, 仍需通过辨识参数$ \vartheta_k $计算初始参数$ { \Theta}_k $, 由表1可知

      图  1  参数误差$\tau$$k$变化曲线

      Figure 1.  Parameter estimation error $\tau$ versus $k$

      表 1  模型参数的估计(混合滤波方法)

      Table 1.  Parameters and their estimates (Hybrid filter)

      $k$$\bar{d}$$\bar{d}a$$\bar{d}b$$\bar{d}c$$\delta\ ({\text{%}})$
      50.009200120.004912630.000039910.000000980.33063
      80.009200150.004923160.000039920.000000980.37553
      100.009203140.004943210.000039910.000000970.59953
      200.009202870.004942950.000039870.000000970.58897
      300.009199720.004838220.000038460.000000990.74688
      真值0.009200000.004900000.000040000.00000100
      $$ \begin{split}{ \vartheta}_{30} =\;& [0.00919972, 0.00483822, 0.00003846, \\ &0.00000099]^{\rm T} \end{split} $$

      根据式(32)和式(34)计算可得

      $$\begin{split} {{\Theta}}_{30} =& [a_{30}, b_{30}, c_{30}, d_{30}]^{\rm T}=\\ & [0.06524981,0.52590948,0.00418056,\\ & 0.00010761]^{\rm T} \end{split}$$

      取采样点(t = 50:100), 则位移真值和估计值变化如图2所示, 速度的估计值和真值变化如图3所示. 由图2图3可知, 通过卡尔曼滤波和粒子滤波并结合递阶辨识原理所得到的位移和速度估计能很好地跟踪真实的位移和速度.

      图  2  位移估计变化曲线($+$: 估计位移; $-$: 真实位移)

      Figure 2.  Displacement estimation curve

      图  3  速度估计变化曲线($+$: 估计速度; $-$: 真实速度)

      Figure 3.  Velocity estimation curve

      进一步通过扩展卡尔曼滤波(Extended Kalman filter)最大期望方法以及粒子滤波(Particle filter)最大期望方法对CRH3高速列车动车组进行建模辨识, 参数估计值及相应的误差分别如图1表2表3所示. 由图1表1$\sim $3可知, 基于扩展的卡尔曼滤波方法辨识精度最低, 粒子滤波方法辨识效果几乎和本文所提出的混合滤波方法效果相当, 但其计算量远远超过本文所提出方法. 因此, 本文所提出的混合滤波最大期望方法具有辨识精度高、计算量小的特点.

      表 2  模型参数的估计(拓展的卡尔曼滤波方法)

      Table 2.  Parameters and their estimates (Extended Kalman filter)

      $k$$\bar{d}$$\bar{d}a$$\bar{d}b$$\bar{d}c$$\delta ({\text{%}})$
      50.009216320.005482360.000031240.000000865.52150
      80.009215940.005444690.000032620.000000895.19838
      100.009216210.005453210.000031730.000000915.22873
      200.009218450.005461510.000032920.000000915.39818
      300.009217060.005450760.000031890.000000915.31865
      真值0.009200000.004900000.000040000.00000100

      表 3  模型参数的估计(粒子滤波方法)

      Table 3.  Parameters and their estimates (Particle filter)

      $k$$\bar{d}$$\bar{d}a$$\bar{d}b$$\bar{d}c$$\delta\ ({\text{%}})$
      50.009176520.005021360.000038650.000000961.14182
      80.009191380.004987540.000039730.000000980.82781
      100.009191410.004987630.000040160.000000980.85850
      200.009191360.005016270.000040980.000000961.02914
      300.009191390.004900320.000040740.000000970.94999
      真值0.009200000.004900000.000040000.00000100
    • 针对单质点高速列车模型线性与非线性耦合的特点, 基于递阶辨识原理, 将高速列车状态空间模型分解为线性位移子模型和非线性速度子模型. 利用最大期望方法处理模型中的隐含变量, 并分别利用卡尔曼滤波器估计位移, 用粒子滤波器估计速度, 提出基于混合滤波的最大期望辨识方法. 由于综合应用两类滤波器的优点, 提出的方法具有计算量小和辨识精度高的优势.

      本文提出的基于混合滤波的最大期望理论框架及方法, 为实现自动列车驾驶精确控制和面向节能的优化驾驶控制提供模型, 并可以作为进一步拓展到多工况条件下的高速列车多模型辨识和重载列车多质点模型辨识的研究基础.

参考文献 (21)

目录

    /

    返回文章
    返回