Hybrid Filter Based Expectation Maximization Algorithm for High-speed Train Modeling
-
摘要: 针对高速列车非线性单质点模型的特殊结构及含有隐含变量问题, 提出一种基于混合滤波的最大期望辨识方法. 借助递阶辨识理论, 将高铁列车状态空间模型分解为线性子系统模型和非线性子系统模型. 进而, 分别利用卡尔曼滤波和粒子滤波对速度和位移状态进行联合估计. 最后, 使用最大期望方法辨识高铁列车子系统模型参数, 解决了隐含变量辨识问题. 和传统方法相比, 本文所提出方法计算量小, 且具有较高的辨识精度. 仿真对比实验结果验证了该方法的有效性.Abstract: For the special high-speed train model structure with hidden variables in the form of the single mass-point, a hybrid filter based expectation maximization (EM) algorithm is proposed. By employing the hierarchical identification theory, the high-speed train state-space model is decomposed into a linear subsystem and a nonlinear subsystem. Furthermore, the Kalman filter and the particle filter are provided to estimate the velocity and displacement, respectively. Finally, the parameters of subsystems are identified by using the EM algorithm. Compared to the classical methods, the proposed algorithm can produce high accuracy estimation with less computational effort. The simulation results verify the effectiveness of the algorithm.
-
随着高速列车进一步提速, 高铁已成为当今国人出行首选. 高铁列控系统是保障列车安全、精确、节能和舒适运行的中枢大脑. 为进一步提升安全性与可靠性, 满足对各类恶劣环境的适应性, 需要对高速列车进行控制和故障检测与诊断[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节给出结论和未来研究方向.
1. 问题描述
考虑如下高速列车单质点模型:
$ \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) $
2. 基于混合滤波的EM算法
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 $ 来估计所有隐含变量.2.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+\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]相比, 本文将复杂的二维高速列车单质点模型分解为相对简单的一维线性位移模型和非线性速度模型, 并分别对每个一维模型单独进行辨识, 因此本文所提出的方法计算量较小.
2.2 EM算法
$ 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. 针对线性位移模型利用卡尔曼滤波估计未知位移, 非线性速度模型采用粒子滤波估计未知速度. 由于卡尔曼滤波是线性状态空间模型最优估计器, 且计算量远小于粒子滤波. 因此本文所采用的混合滤波方法具有计算量小、精度高的优点.
3. 仿真例子
考虑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 模型参数的估计(混合滤波方法)Table 1 Parameters and their estimates (Hybrid filter)$k$ $\bar{d}$ $\bar{d}a$ $\bar{d}b$ $\bar{d}c$ $\delta\ ({\text{%}})$ 5 0.00920012 0.00491263 0.00003991 0.00000098 0.33063 8 0.00920015 0.00492316 0.00003992 0.00000098 0.37553 10 0.00920314 0.00494321 0.00003991 0.00000097 0.59953 20 0.00920287 0.00494295 0.00003987 0.00000097 0.58897 30 0.00919972 0.00483822 0.00003846 0.00000099 0.74688 真值 0.00920000 0.00490000 0.00004000 0.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可知, 通过卡尔曼滤波和粒子滤波并结合递阶辨识原理所得到的位移和速度估计能很好地跟踪真实的位移和速度.
进一步通过扩展卡尔曼滤波(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{%}})$ 5 0.00921632 0.00548236 0.00003124 0.00000086 5.52150 8 0.00921594 0.00544469 0.00003262 0.00000089 5.19838 10 0.00921621 0.00545321 0.00003173 0.00000091 5.22873 20 0.00921845 0.00546151 0.00003292 0.00000091 5.39818 30 0.00921706 0.00545076 0.00003189 0.00000091 5.31865 真值 0.00920000 0.00490000 0.00004000 0.00000100 表 3 模型参数的估计(粒子滤波方法)Table 3 Parameters and their estimates (Particle filter)$k$ $\bar{d}$ $\bar{d}a$ $\bar{d}b$ $\bar{d}c$ $\delta\ ({\text{%}})$ 5 0.00917652 0.00502136 0.00003865 0.00000096 1.14182 8 0.00919138 0.00498754 0.00003973 0.00000098 0.82781 10 0.00919141 0.00498763 0.00004016 0.00000098 0.85850 20 0.00919136 0.00501627 0.00004098 0.00000096 1.02914 30 0.00919139 0.00490032 0.00004074 0.00000097 0.94999 真值 0.00920000 0.00490000 0.00004000 0.00000100 4. 总结
针对单质点高速列车模型线性与非线性耦合的特点, 基于递阶辨识原理, 将高速列车状态空间模型分解为线性位移子模型和非线性速度子模型. 利用最大期望方法处理模型中的隐含变量, 并分别利用卡尔曼滤波器估计位移, 用粒子滤波器估计速度, 提出基于混合滤波的最大期望辨识方法. 由于综合应用两类滤波器的优点, 提出的方法具有计算量小和辨识精度高的优势.
本文提出的基于混合滤波的最大期望理论框架及方法, 为实现自动列车驾驶精确控制和面向节能的优化驾驶控制提供模型, 并可以作为进一步拓展到多工况条件下的高速列车多模型辨识和重载列车多质点模型辨识的研究基础.
-
表 1 模型参数的估计(混合滤波方法)
Table 1 Parameters and their estimates (Hybrid filter)
$k$ $\bar{d}$ $\bar{d}a$ $\bar{d}b$ $\bar{d}c$ $\delta\ ({\text{%}})$ 5 0.00920012 0.00491263 0.00003991 0.00000098 0.33063 8 0.00920015 0.00492316 0.00003992 0.00000098 0.37553 10 0.00920314 0.00494321 0.00003991 0.00000097 0.59953 20 0.00920287 0.00494295 0.00003987 0.00000097 0.58897 30 0.00919972 0.00483822 0.00003846 0.00000099 0.74688 真值 0.00920000 0.00490000 0.00004000 0.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{%}})$ 5 0.00921632 0.00548236 0.00003124 0.00000086 5.52150 8 0.00921594 0.00544469 0.00003262 0.00000089 5.19838 10 0.00921621 0.00545321 0.00003173 0.00000091 5.22873 20 0.00921845 0.00546151 0.00003292 0.00000091 5.39818 30 0.00921706 0.00545076 0.00003189 0.00000091 5.31865 真值 0.00920000 0.00490000 0.00004000 0.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{%}})$ 5 0.00917652 0.00502136 0.00003865 0.00000096 1.14182 8 0.00919138 0.00498754 0.00003973 0.00000098 0.82781 10 0.00919141 0.00498763 0.00004016 0.00000098 0.85850 20 0.00919136 0.00501627 0.00004098 0.00000096 1.02914 30 0.00919139 0.00490032 0.00004074 0.00000097 0.94999 真值 0.00920000 0.00490000 0.00004000 0.00000100
下载: 导出CSV
-
[1] 衷路生, 颜真, 杨辉, 齐叶鹏, 张坤鹏, 樊晓平. 数据驱动的高速列车子空间预测控制方法. 铁道学报, 2013, 35(4): 77−83 doi: 10.3969/j.issn.1001-8360.2013.04.0121 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.0083 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−4968 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−295817 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.c180298Gu 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 期刊类型引用(3)
1. 丁锋，徐玲，张霄. 大规模系统高效辨识方法研究. 青岛科技大学学报(自然科学版). 2025(01): 1-16 . 
百度学术2. 陈晶，朱全民. 有理模型辨识的两类新方法—–混合迭代与柔性最小二乘法. 控制与决策. 2022(01): 58-66 . 百度学术3. 刘晓宇，荀径，高士根，阴佳腾. 高速列车精确停车的鲁棒自触发预测控制. 自动化学报. 2022(01): 171-181 . 本站查看其他类型引用(5)
-
下载:
计量
- 文章访问数: 2586
- HTML全文浏览量: 446
- PDF下载量: 190
- 被引次数: 8
