A robust recursive speed tracking estimation for subway trains under uncertain bounded disturbance

Abstract

This article investigates the speed tracking estimation problem for the subway system under uncertain bounded disturbance. Firstly, the nonlinear running model of the train with dynamic descriptions is formulated, which contains traction/braking force, basic resistance, gradient, and uncertain disturbance under different conditions. By modeling the linearization error as uncertain parameters, a convex optimization problem of Kalman-like nonlinear estimator is established. Then, a robust speed tracking method is designed based on linear matrix inequalities and bounded recursive optimization idea, and the sufficient conditions to ensure the bounded stability of the nonlinear estimator are derived. Therefore, an effective solution to the stability problem of the designed nonlinear estimator with unknown disturbance statistical characteristics is given. Finally, the actual operational data of Ningbo Metro Line 7 is employed to verify the effectiveness of the proposed method.

Introduction

With the rapid and all-encompassing development of smart cities in China, the public’s increasing demand for transportation capacity has brought unprecedented challenges to urban transportation systems, and the subway system has become the key to solving the dilemma of urban public transportation with its many superiorities. ¹ In particular, the automatic train operation (ATO) system enables fine management and optimal scheduling of train operations by utilizing advanced control strategies, precise sensor technologies, and efficient communication technologies. ² Most importantly, the safe operation of the train is ensured through comprehensive monitoring and rapid fault diagnosis.

Related work

The trajectory tracking problem is a crucial issue in train operation monitoring, especially speed profile tracking. In fact, the automatic train speed regulation module is a key component of the ATO system. Based on the off-line speed profile, the speed tracking controller can be designed to realize the operation optimization of the train running time, running speed, stopping position, and departure interval. ³ In order to design stable and effective controllers, various train speed profile tracking control methods have been developed, such as proportional-integral-derivative (PID) control, ⁴ fuzzy control, ⁵ adaptive control, ⁶ model predictive control (MPC). ⁷ Though PID control is widely used for train tracking control, its performance degradation is inevitable in dynamic nonlinear train motion models with uncertain parameters. ⁸ Then, a comprehensive optimization method was proposed in Ref. ⁹ to track the speed profile of the train. Particularly, the nonlinear layered neural network was employed for speed profile tracking in Ref. ¹⁰ By combing neural network and PID control methodologies, a model-free adaptive controller was designed in Ref. ¹¹ for achieving precise speed profile tracking. Consider various braking conditions and parameter uncertainties, ¹² constructed an adaptive neuro-fuzzy inference system model to capture the characteristics of running train. Moreover, a MPC-based profile tracking controller was proposed in Ref. ¹³ to deal with uncertain parameters and varying speed, and Liu et al. ¹⁴ developed a predictive fuzzy proportional-integral-derivative control method for suspended permanent magnetic maglev train.

To improve the performance of the speed controller, the observer or estimator can provide accurate speed tracking information for controller designing. ¹⁵ Since the actual modeling of train operation system is nonlinear, the extended Kalman filter (EKF) ¹⁶ was usually employed for nonlinear target tracking problems. ¹⁷ Since EKF utilizes Taylor series for linearizing operation, while the linearization errors and the uncertainties are ignored in this process. Though the second order EKF ¹⁸ reduces the linearization errors, the Hessian matrix is difficult to be calculated and computationally intensive. Then, the unscented Kalman filter (UKF) ¹⁹ and cubature Kalman filter (CKF) ²⁰ were proposed by sampling sigma points to approximate the probability distribution of the system. In Ref., ²¹ an adaptive maximum correntropy UKF tracking method was developed based on multi-layer isogradient sound speed profile. However, the UKF and CKF performed a sensitive to nonpositive definite covariance matrices. Moreover Hu et al., ²² adopted a statical linear regression method to linearize the nonlinear systems. Unfortunately, the above filtering methods require the system noises to obey Gaussian distribution, and the noise covariances need to be known in advance.

In fact, the uncertain disturbance during train operation often does not obey any random distribution, ²³ and the statistical information is difficult to obtain accurately, but it is always bounded. Under this case, various $H_{\infty }$ filter^24,25 methods for dealing with energy-bounded noises have been developed. Note that the energy-bounded noises will gradually decay to 0 in the $H_{\infty }$ framework over time, but the noises in the subway systems are always persistent and do not fade away. Moreover, set-member filter ²⁶ can also deal with the bounded disturbance. Based on this method, a proportional integral estimator was proposed in Ref., ²⁷ and an iterative proportional-integral interval estimation method was developed in Ref. ²⁸ to deal with bounded noises for linear systems. Meanwhile Shen et al., ²⁹ constructed improved proportional-integral framework for linear system with unknown disturbances, and an iterative interval estimator was designed to obtain tighter estimation intervals. However, the above-mentioned works are all within the linear framework. Though a set membership estimation method was proposed in Ref. ³⁰ for nonlinear dynamic systems with bounded noises, the noises are still confined to specified ellipsoidal sets. Moreover, an approach utilizing bounded recursive optimization was devised in Ref. ³¹ to tackle both nonlinearities and unknown bounded noises, while the stability criteria associated with this method are rather conservative and fall short of being comprehensive.

Motivation and contributions

Based on the above analysis, it remains a challenge to study the stable tracking algorithm for nonlinear systems with more general unknown but bounded disturbance. Therefore, this paper shall study the speed tracking estimation issue for metro trains under uncertain bounded disturbance. The main contributions of this study are summarized as follows:

• (i) A robust recursive optimization tracking estimation method is proposed to deal with uncertain bounded disturbance, and the sufficient conditions for ensuring the bounded stability of the proposed nonlinear estimator are derived by modeling the linearization error as the state-dependent matrix with uncertain parameters.

• (ii) Based on the actual operational data that was collected from the YLZ-YFL segment of Ningbo Metro Line 7, the tracking performance comparison of different disturbances and different nonlinear estimators are presented to show the effectiveness and the advantages of the proposed algorithm.

Notations: ${\mathbb{R}}^{n\times q}$ represents the set of real matrices with $n\times q$ dimension, and $I$ stands for the identity matrix with corresponding dimension. $A^{\top }$ represents the transpose of matrix $A$ , and $A<0$ denotes the negative definite matrix. $\parallel ·{\parallel }_2$ is the Euclidean norm, $\mathrm{diag}\{·\}$ denotes the block diagonal matrix, and “*” stands for the symmetric term of a symmetric matrix.

Dynamics modeling and problem formulation

Motion equation and parameter

When analyzing the motion of a train along a track, the train is usually considered a mass point. ³² Since the modeling of train kinematics equations usually only considers the forces acting in the forward direction of the train along the track, the longitudinal motion equation of the train can be modeled by Newton’s laws and described as the following time-varying model:

$$\overline{\mathrm{M}}\left(t\right)\stackrel{·}{s}\left(t\right)=\mathrm{F}\left(t\right)-\mathrm{R}\left(t\right)$$

(1)

where $t$ is the running time, $s\left(t\right)$ denotes the speed of the train, $\overline{\mathrm{M}}\left(t\right)$ denotes the rotational mass, $\mathrm{F}\left(t\right)$ denotes the traction/braking force that is generated by the corresponding traction/braking system, while $\mathrm{R}\left(t\right)$ denotes the resistance that is affected by the characteristics of the train and environmental conditions.

Notice that, the moment of inertia will cause a mismatch between the mass of the train and its static mass during train operation. In this case, the train mass shall be amended as rotational mass in the dynamic calculation. Then, the computation of the rotational mass is given by:

$$\overline{\mathrm{M}}\left(t\right)=(1+\gamma )\mathrm{M}\left(t\right)$$

(2)

where $\gamma$ is the known rotational mass coefficient. $\mathrm{M}\left(t\right)$ denotes the total static mass of the train, which includes the vehicle mass and the passenger mass.

In practice, the motor drives the train through traction/braking control, where the maximum output traction and braking force of the motor varies with the running speed of the train. Therefore, the maximum traction/braking force can be determined by the following equations at a specific operating speed:

$$\{\begin{array}{c}{\mathrm{F}}_{\max }^t\left(t\right)=f_{\max }^t\left(s\left(t\right)\right)\hfill \\ {\mathrm{F}}_{\max }^b\left(t\right)=f_{\max }^b\left(s\left(t\right)\right)\hfill \end{array}$$

(3)

where ${\mathrm{F}}_{\max }^t\left(t\right)$ represents the maximum traction force, and $f_{\max }^t\left(s\left(t\right)\right)$ is a function that describes the relationship between the train speed and the maximum traction force; ${\mathrm{F}}_{\max }^b\left(t\right)$ represents the maximum braking force, and $f_{\max }^b\left(s\left(t\right)\right)$ is a function for describing the relationship between the train speed and the maximum braking force.

It has been verified that the ATO system determines the percentage output of the force signal. ³³ In this case, the percentage of maximum traction/braking force output is employed for train controlling, and then the traction/braking force $\mathrm{F}\left(t\right)$ can be calculated by:

$$\begin{array}{l}\mathrm{F}\left(t\right)=u\left(t\right){\mathrm{F}}_{\max }\left(t\right)\\ =\{\begin{array}{cc}u\left(k\right)f_{\max }^t\left(s\left(t\right)\right),& {\mathrm{\Theta }}_1\left(t\right)=1,{\mathrm{\Theta }}_2\left(t\right)=0;\\ u\left(k\right)f_{\max }^b\left(s\left(t\right)\right),& {\mathrm{\Theta }}_1\left(t\right)=0,{\mathrm{\Theta }}_2\left(t\right)=1;\\ 0,& {\mathrm{\Theta }}_1\left(t\right)=0,{\mathrm{\Theta }}_2\left(t\right)=0.\end{array}\end{array}$$

(4)

where $u\left(t\right)$ is the control signal, ${\mathrm{F}}_{\max }\left(t\right)$ denotes the corresponding maximum traction/braking force, ${\mathrm{\Theta }}_1\left(t\right)$ and ${\mathrm{\Theta }}_2\left(t\right)$ are command indicators. Since $u\left(t\right)$ is a percentage of force output, and then $\left|u\left(t\right)\right|\le 1$ holds. It should be noted that traction and braking forces are in opposite directions, and $0<u\left(t\right)\le 1$ implies traction force action, while $-1\le u\left(t\right)<0$ implies braking force action. In other words, ${\mathrm{\Theta }}_1\left(t\right)=1,{\mathrm{\Theta }}_2\left(t\right)=0$ indicates traction command, ${\mathrm{\Theta }}_1\left(t\right)=0,{\mathrm{\Theta }}_2\left(t\right)=1$ indicates braking command, and ${\mathrm{\Theta }}_1\left(t\right)=0,{\mathrm{\Theta }}_2\left(t\right)=0$ indicates no control command. Particularly, $\sum _{i=1}^2{\mathrm{\Theta }}_i\left(t\right)\le 1$ , that is, traction and braking commands cannot be executed at the same time.

On the other hand, resistance is also an important component of train operation, which including basic resistance, gradient resistance and uncertain disturbances. Then, the resistance $\mathrm{R}\left(t\right)$ can be modeled by ³⁴ :

$$\mathrm{R}\left(t\right)={\mathrm{R}}_b\left(t\right)+{\mathrm{R}}_g\left(t\right)+\varpi \left(t\right)$$

(5)

where ${\mathrm{R}}_b\left(t\right)$ denotes the basic resistance, ${\mathrm{R}}_g\left(t\right)$ denotes the gradient resistance, and $\varpi \left(t\right)$ is the bounded uncertain disturbance.

In fact, basic resistance is mainly caused by mechanical and aerodynamic resistance, which is usually molded by the Davis formula ² :

$${\mathrm{R}}_b\left(t\right)=\alpha \left(t\right)s^2\left(t\right)+\beta \left(t\right)s\left(t\right)+\zeta \left(t\right)$$

(6)

where $\alpha \left(t\right)$ , $\beta \left(t\right)$ and $\zeta \left(t\right)$ are empirical parameters.

Meanwhile, gradient resistance is usually arisen from track gradients, which is determined by:

$${\mathrm{R}}_g\left(t\right)=M\left(t\right)g\sin \left(\arctan \left(\delta \left(t\right)\right)\right)$$

(7)

where $g$ is gravity, $\delta \left(t\right)$ denotes the slope gradient of the subway tarck.

According to the requirements of track construction standards, the slope gradient is small enough, and then $\sin \left(\arctan \left(\delta \left(t\right)\right)\right)\approx \delta \left(t\right)$ holds by using Taylor series expansion. In terms of resistance, the direction of basic resistance is always opposite to the train operation, the direction of the gravity resistance is determined by the sign of the gradient, while the direction of the perturbation is indeterminate. Notice that the uncertainty in the direction of gravity resistance and perturbation force implies that these two forces may favor train operation.

Model foundation

It is noted above that trains are also affected by traction force, braking Force and resistance. Then, from (1)–(7), the longitudinal motion model is modeled as

$$\begin{array}{c}(1+\gamma )\mathrm{M}\left(t\right)\stackrel{·}{s}\left(t\right)\hfill \\ =u\left(t\right)\mathcal{F}\left(t\right)\mathcal{E}\left(t\right)-\mathcal{H}\left(t\right)\mathcal{X}\left(t\right)\hfill \\ -\mathrm{M}\left(t\right)g\delta \left(t\right)-\varpi \left(t\right)\hfill \end{array}$$

(8)

Where

$$\{\begin{array}{c}\mathcal{F}\left(t\right)\stackrel{\Delta }{=}\left[f_{\max }^t\left(s\left(t\right)\right)f_{\max }^b\left(s\left(t\right)\right)\right]\hfill \\ \multicolumn{1}{c}{\mathcal{E}\left(t\right)\stackrel{\Delta }{=}{\left[{\mathrm{\Theta }}_1\left(t\right){\mathrm{\Theta }}_2\left(t\right)\right]}^{\top }}\\ \multicolumn{1}{c}{\begin{array}{c}\mathcal{H}\left(t\right)\stackrel{\Delta }{=}\left[\alpha \left(t\right)\beta \left(t\right)\zeta \left(t\right)\right]\hfill \\ \multicolumn{1}{c}{\mathcal{X}\left(t\right)\stackrel{\Delta }{=}{\left[1s\left(t\right)s^2\left(t\right)\right]}^{\top }}\end{array}}\end{array}$$

(9)

It follows from (8) that

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{s}\left(t\right)=\frac{1}{1+\gamma }(u\left(t\right){\mathcal{F}}_m\left(t\right)\mathcal{E}\left(k\right)}\\ \multicolumn{1}{c}{-{\mathcal{H}}_m\left(t\right)\mathcal{X}\left(k\right)-g\delta \left(t\right)-{\varpi }_m\left(t\right))}\end{array}$$

(10)

where ${\mathcal{F}}_m\left(t\right)=\mathcal{F}\left(t\right)/\mathrm{M}\left(t\right)$ is the maximum traction/braking force per unit mass, ${\mathcal{H}}_m\left(t\right)=\mathcal{H}\left(t\right)/\mathrm{M}\left(t\right)$ is the corresponding basic resistance parameters per unit mass, and ${\varpi }_m\left(t\right)$ is the corresponding disturbance.

To facilitate numerical calculation and analysis in the computer, the continuous system (10) is discretized as

$$\begin{array}{c}s(k+1)=s\left(k\right)+\frac{T}{1+\gamma }({\mathcal{F}}_m\left(k\right)\mathcal{E}\left(k\right)u\left(k\right)\hfill \\ -{\mathcal{H}}_m\left(k\right)\mathcal{X}\left(k\right)-g\delta \left(k\right))+\mathrm{Tw}\left(k\right)\hfill \end{array}$$

(11)

where $T$ is the sampling period, $k$ is the running step, and $w\left(k\right)={\varpi }_m\left(k\right)/(1+\gamma )$ .

Meanwhile, according to the relationship between position and velocity, the position model can be modeled as

$$p(k+1)=p\left(k\right)+\mathrm{Ts}\left(k\right)+\frac{1}{2}{T}^{2}w\left(k\right)$$

(12)

Then, the motion model of the train can be described by

$$\chi (k+1)=\mathcal{A}\left(k\right)\chi \left(k\right)+f(\chi \left(k\right),u\left(k\right))+\mathcal{B}\left(k\right)w\left(k\right)$$

(13)

Where

$$\{\begin{array}{c}\begin{array}{c}\mathcal{A}\left(k\right)\stackrel{\Delta }{=}\left[\begin{array}{cc}1\hfill & T\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]\hfill \\ \mathcal{B}\left(k\right)\stackrel{\Delta }{=}{[0.5T^2,T]}^{\top }\hfill \\ \chi \left(k\right)\stackrel{\Delta }{=}{[p\left(k\right),s\left(k\right)]}^{\top }\hfill \\ f(\chi \left(k\right),u\left(k\right))\stackrel{\Delta }{=}[0\frac{T}{1+\gamma }({\mathcal{F}}_m\left(k\right)\mathcal{E}\left(k\right)u\left(k\right)\hfill \\ {-{\mathcal{H}}_m\left(k\right)\mathcal{X}\left(k\right)-g\delta \left(k\right)]}^{\top }\hfill \end{array}\hfill \end{array}$$

(14)

Moreover, the sensor can be used to monitor the motion of train, where the measurement model of the sensor can be described by

$$\mathcal{Y}\left(k\right)=\mathcal{C}\left(k\right)\chi \left(k\right)+v\left(k\right)$$

(15)

where $\mathcal{C}\left(k\right)=\left[10\right]$ , and measurement noise $v\left(k\right)$ usually does not follow the gaussian distribution in a real train operating environment, while it always bounded, that is,

$$v^2\left(k\right)\le {\sigma }_v$$

(16)

where ${\sigma }_v$ is an unknown scalar.

Problem formulation

Based on the measurements $\{\mathcal{Y}\left(1\right),\mathcal{Y}\left(2\right),\dots ,\mathcal{Y}\left(k\right)\}$ , an estimator to obtain the train’s speed is designed as

$$\{\begin{array}{c}\begin{array}{c}{\hat{\chi }}^p\left(k\right)=\mathcal{A}(k-1)\hat{\chi }(k-1)\hfill \\ +f(\hat{\chi }(k-1),u(k-1))\hfill \\ \hat{\chi }\left(k\right)={\hat{\chi }}^p\left(k\right)+\mathcal{K}\left(k\right)[\mathcal{Y}\left(k\right)-\mathcal{C}\left(k\right){\hat{\chi }}^p\left(k\right)]\hfill \end{array}\hfill \end{array}$$

(17)

where ${\hat{\chi }}^p\left(k\right)$ denotes one-step prediction, $\hat{\chi }\left(k\right)$ indicates the tracking state of the train, and $\mathcal{K}\left(k\right)$ is the gain matrix to be determined.

Consequently, the problems to be addressed in this paper can be summarized in the following aspects:

1. To obtain the estimator gain $\mathcal{K}\left(k\right)$ such that the upper bound of the square error (i.e. ${(\chi \left(k\right)-\hat{\chi }\left(k\right))}^{\top }(\chi \left(k\right)-\hat{\chi }\left(k\right))$ ) of the designed estimator (17) is minimal at step $k$ .

2. To construct stability conditions such that the square error of the designed estimator (17) is asymptotically bounded.

Main results

In this section, a nonlinear estimator is designed for tracking the speed of the train. First of all, a useful lemma is introduced as follows before deriving the main results:

• Lemma 1. ³⁵ Given appropriate dimension matrices ${\mathrm{\Sigma }}_1$ , ${\mathrm{\Sigma }}_2$ and ${\mathrm{\Sigma }}_3$ , where ${\mathrm{\Sigma }}_1^{\top }{\mathrm{\Sigma }}_1$ then

$$\begin{array}{c}{\mathrm{\Sigma }}_1+{\mathrm{\Sigma }}_3\mathrm{P}{\mathrm{\Sigma }}_2+{\mathrm{\Sigma }}_2^{\top }{\mathrm{P}}^{\top }{\mathrm{\Sigma }}_3^{\top }<0\hfill \end{array}$$

holds for all matrices $\mathrm{P}$ satisfying ${\mathrm{P}}^{\top }\mathrm{P}\le \mathrm{I}$ if and only if for some $ϵ>0$

$$\begin{array}{c}{\mathrm{\Sigma }}_1+{ϵ}^{-1}{\mathrm{\Sigma }}_3{\mathrm{\Sigma }}_3^{\top }+ϵ{\mathrm{\Sigma }}_2^{\top }{\mathrm{\Sigma }}_2<0\hfill \end{array}$$

□

Define $\tilde{\chi }\left(k\right)\stackrel{\Delta }{=}\chi \left(k\right)-\hat{\chi }\left(k\right)$ is tracking error, and it follows from (13) and (17) that

$$\begin{array}{c}\tilde{\chi }\left(k\right)=\mathcal{A}(k-1)\chi (k-1)+f(\chi (k-1),u(k-1))\hfill \\ +\mathcal{B}(k-1)w(k-1)-\mathcal{A}(k-1)\hat{\chi }(k-1)\hfill \\ -f(\hat{\chi }(k-1),u(k-1))\hfill \\ -\mathcal{K}\left(k\right)[\mathcal{Y}\left(k\right)-\mathcal{C}\left(k\right){\hat{\chi }}^p\left(k\right)]\hfill \end{array}$$

(18)

By using the first-order Taylor series to expand the nonlinear function $f(\chi \left(k\right),u\left(k\right))$ near $\hat{\chi }\left(k\right)$ , and modeling the linearization errors by state-dependent matrix with uncertain parameter, ³⁶ one has

$$\begin{array}{c}f(\chi (k-1),u(k-1))\hfill \\ =f(\hat{\chi }(k-1),u(k-1))+{\mathcal{A}}_f(k-1)\hfill \\ \times \tilde{\chi }(k-1)+\mathcal{P}(k-1)\mathcal{O}(k-1)\tilde{\chi }(k-1)\hfill \end{array}$$

(19)

where $\mathcal{P}(k-1)$ is the state-dependent scaling matrix and the determination criterion was presented in Ref., ³⁶ $\mathcal{O}(k-1)$ is the uncertain parameter that satisfying $\parallel \mathcal{O}(k-1){\parallel }_2\le 1$ and

$${{\mathcal{A}}_f(k-1)=\frac{\partial f(\chi ,u)}{\partial \chi }|}_{\chi =\hat{\chi }(k-1)}$$

(20)

It follows from (18)–(20) that

$$\begin{array}{c}\tilde{\chi }\left(k\right)={\mathcal{K}}_{\mathcal{C}}\left(k\right)[\mathcal{A}(k-1)+{\mathcal{A}}_f(k-1)+\mathcal{P}(k-1)\hfill \\ \times \mathcal{O}(k-1)]\tilde{\chi }(k-1)+{\mathcal{K}}_{\mathcal{C}}\left(k\right)\hfill \\ \times \mathcal{B}(k-1)w(k-1)-\mathcal{K}\left(k\right)v\left(k\right)\hfill \end{array}$$

(21)

where ${\mathcal{K}}_{\mathcal{C}}\left(k\right)\stackrel{\Delta }{=}I-\mathcal{K}\left(k\right)\mathcal{C}\left(k\right)$ .

Whereafter, the estimator gain $\mathcal{K}\left(k\right)$ in (17) shall be determined in following Theorem.

• Theorem 1. For the dynamic system (13)–(15) of a moving train, a convex optimization problem in terms of linear matrix inequalities is proposed to determine estimator gain $\mathcal{K}\left(k\right)$ as follows:

$$\begin{array}{c}\underset{\begin{array}{c}\mathfrak{P}\left(k\right)>0,\mathfrak{Q}\left(k\right)>0\hfill \\ \multicolumn{1}{c}{\mathfrak{G}\left(k\right),\mathcal{K}\left(k\right),{ϵ}_1\left(k\right),{ϵ}_2\left(k\right)}\end{array}}{\min }\mathrm{Tr}\{\mathfrak{P}\left(k\right)\}+\mathrm{Tr}\{\mathfrak{Q}\left(k\right)\}\hfill \\ \multicolumn{1}{c}{\mathrm{s}.\mathrm{t}.\{\begin{array}{c}\left[\begin{array}{ccc}-{ϵ}_1\left(k\right)I\hfill & {ϵ}_1\left(k\right){\mathcal{I}}_{\mathcal{O}}\left(k\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{*}& {\mathcal{K}}_{\mathcal{B}}^{\mathcal{A}}\left(k\right)\hfill & {\mathcal{P}}_{\mathcal{I}}\left(k\right)\hfill \\ \multicolumn{1}{c}{*}& *\hfill & -{ϵ}_1\left(k\right)I\hfill \end{array}\right]<0\hfill \\ \multicolumn{1}{c}{\left[\begin{array}{cccc}-{ϵ}_2\left(k\right)I\hfill & 0\hfill & {ϵ}_2\left(k\right)I\hfill & 0\hfill \\ \multicolumn{1}{c}{*}& -I\hfill & {\mathcal{K}}_{\mathcal{C}}^{\mathcal{A}}\left(k\right)\hfill & {\mathcal{K}}_{\mathcal{C}}\left(k\right)\mathcal{P}(k-1)\hfill \\ \multicolumn{1}{c}{*}& *\hfill & -I\hfill & 0\hfill \\ \multicolumn{1}{c}{*}& *\hfill & *\hfill & -{ϵ}_2\left(k\right)I\hfill \end{array}\right]<0}\\ \multicolumn{1}{c}{\left[\begin{array}{cc}-I\hfill & {\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)\hfill \\ \multicolumn{1}{c}{*}& -{\xi }_{\mathcal{K}}I\hfill \end{array}\right]\le 0}\end{array}}\end{array}$$

(22)

where

$$\{\begin{array}{c}\multicolumn{1}{c}{\begin{array}{c}{\mathcal{K}}_{\mathcal{B}}^{\mathcal{A}}\left(k\right)\stackrel{\Delta }{=}\left[\begin{array}{ccc}-I\hfill & {\mathcal{K}}_{\mathcal{C}}^{\mathcal{A}}\left(k\right)\hfill & {\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)\hfill \\ \multicolumn{1}{c}{*}& -\mathfrak{P}\left(k\right)\hfill & -\mathfrak{G}\left(k\right)\hfill \\ \multicolumn{1}{c}{*}& *\hfill & -\mathfrak{Q}\left(k\right)\hfill \end{array}\right]\hfill \\ \multicolumn{1}{c}{{\mathcal{K}}_{\mathcal{C}}^{\mathcal{A}}\left(k\right)\stackrel{\Delta }{=}{\mathcal{K}}_{\mathcal{C}}\left(k\right)[\mathcal{A}(k-1)+{\mathcal{A}}_f(k-1)]}\\ \multicolumn{1}{c}{{\mathcal{P}}_{\mathcal{I}}\left(k\right)\stackrel{\Delta }{=}{\left[\begin{array}{ccc}{\mathcal{P}}^{\top }(k-1){\mathcal{K}}_{\mathcal{C}}^{\top }\left(k\right)\hfill & 0\hfill & 0\hfill \end{array}\right]}^{\top }}\\ \multicolumn{1}{c}{{\mathcal{I}}_{\mathcal{O}}\left(k\right)\stackrel{\Delta }{=}\left[\begin{array}{ccc}0\hfill & I\hfill & 0\hfill \end{array}\right]}\end{array}}\end{array}$$

(23)

In this case, the square error of the estimate $\hat{\chi }\left(k\right)$ shall be asymptotically bounded, that is,

$$\begin{array}{c}\underset{k\to \infty }{\lim }{\tilde{\chi }}^{\top }\left(k\right)\tilde{\chi }\left(k\right)<\xi \hfill \end{array}$$

(24)

where $\xi$ is an unknown scalar sufficiently small.

Proof. In general, the tracking error (21) of the train can be rewritten as

$$\begin{array}{c}\tilde{\chi }\left(k\right)={\mathcal{A}}_{\mathcal{K}}\left(k\right)\tilde{\chi }(k-1)+{\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)w_v\left(k\right)\hfill \end{array}$$

(25)

where

$$\{\begin{array}{c}\begin{array}{c}{\mathcal{A}}_{\mathcal{K}}\left(k\right)\stackrel{\Delta }{=}{\mathcal{K}}_{\mathcal{C}}^{\mathcal{A}}\left(k\right)+{\mathcal{K}}_{\mathcal{C}}\left(k\right)\mathcal{P}(k-1)\mathcal{O}(k-1)\hfill \\ {\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)\stackrel{\Delta }{=}[{\mathcal{K}}_{\mathcal{C}}\left(k\right)\mathcal{B}(k-1),-\mathcal{K}\left(k\right)]\hfill \\ w_v\left(k\right)\stackrel{\Delta }{=}{[w^{\top }(k-1),v^{\top }\left(k\right)]}^{\top }\hfill \end{array}\hfill \end{array}$$

(26)

Then, introducing an index performance ³⁷ to construct the upper bound of the square error of the designed estimator (17):

$$\begin{array}{c}\mathfrak{J}\left(k\right)\stackrel{\Delta }{=}{\tilde{\chi }}^{\top }\left(k\right)\tilde{\chi }\left(k\right)-{\tilde{\chi }}^{\top }(k-1)\mathfrak{P}\left(k\right)\tilde{\chi }(k-1)\hfill \\ -2{\tilde{\chi }}^{\top }(k-1)\mathfrak{G}\left(k\right)w_v\left(k\right)-w_v^{\top }\left(k\right)\mathfrak{Q}\left(k\right)w_v\left(k\right)\hfill \end{array}$$

(27)

where $\mathfrak{P}\left(k\right)>0$ , $\mathfrak{Q}\left(k\right)>0$ and $\mathfrak{G}\left(k\right)$ introduced optimization matrices to construct an upper bound of the square error of the estimator (17).

Let ${\tilde{\chi }}_w\left(k\right)\stackrel{\Delta }{=}[{\tilde{\chi }}^{\top }(k-1),w_v^{\top }\left(k\right){]}^{\top }$ , and substituting (25) into (27) yields

$$\begin{array}{c}\multicolumn{1}{c}{\mathfrak{J}\left(k\right)={\tilde{\chi }}_w^{\top }\left(k\right)\underset{\mathfrak{X}\left(k\right)}{\underbrace{\left[\begin{array}{cc}{\mathfrak{X}}_1\left(k\right)\hfill & {\mathfrak{X}}_2\left(k\right)\hfill \\ \multicolumn{1}{c}{*}& {\mathfrak{X}}_3\left(k\right)\hfill \end{array}\right]}}{\tilde{\chi }}_w\left(k\right)}\end{array}$$

(28)

where

$$\{\begin{array}{c}\begin{array}{c}{\mathfrak{X}}_1\left(k\right)={\mathcal{A}}_{\mathcal{K}}^{\top }\left(k\right){\mathcal{A}}_{\mathcal{K}}\left(k\right)-\mathfrak{P}\left(k\right)\hfill \\ {\mathfrak{X}}_2\left(k\right)={\mathcal{A}}_{\mathcal{K}}^{\top }\left(k\right){\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)-\mathfrak{G}\left(k\right)\hfill \\ {\mathfrak{X}}_3\left(k\right)={\left({\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)\right)}^{\top }{\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)-\mathfrak{Q}\left(k\right)\hfill \end{array}\hfill \end{array}$$

(29)

Meanwhile, it can be deduced from $\mathfrak{J}\left(k\right)<0$ that

$$\begin{array}{c}{\tilde{\chi }}^{\top }\left(k\right)\tilde{\chi }\left(k\right)<{\tilde{\chi }}^{\top }(k-1)\mathfrak{P}\left(k\right)\tilde{\chi }(k-1)\hfill \\ +2{\tilde{\chi }}^{\top }(k-1)\mathfrak{G}\left(k\right)w_v\left(k\right)+w_v^{\top }\left(k\right)\mathfrak{Q}\left(k\right)w_v\left(k\right)\hfill \end{array}$$

(30)

Obviously, the right term of (30) can be viewed as an upper bounded of ${\tilde{\chi }}^{\top }\left(k\right)\tilde{\chi }\left(k\right)$ . By using Schur complement lemma, ³⁵ $\mathfrak{J}\left(k\right)<0$ holds is equivalent to following inequality holds:

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{ccc}-I\hfill & {\mathcal{A}}_{\mathcal{K}}(k-1))\hfill & {\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)\hfill \\ \multicolumn{1}{c}{*}& -\mathfrak{P}\left(k\right)\hfill & -\mathfrak{G}\left(k\right)\hfill \\ \multicolumn{1}{c}{*}& *\hfill & -\mathfrak{Q}\left(k\right)\hfill \end{array}\right]}\\ \multicolumn{1}{c}{={\mathcal{K}}_{\mathcal{B}}^{\mathcal{A}}\left(k\right)+{\mathcal{P}}_{\mathcal{I}}\left(k\right)\mathcal{O}(k-1){\mathcal{I}}_{\mathcal{O}}\left(k\right)}\\ \multicolumn{1}{c}{+{\mathcal{I}}_{\mathcal{O}}^{\top }\left(k\right){\mathcal{O}}^{\top }(k-1){\mathcal{P}}_{\mathcal{I}}^{\top }\left(k\right)<0}\end{array}$$

(31)

It is derived from Lemma 1 that the first inequality in (20) holds if and only if (31) holds.

Notice that

$$\begin{array}{c}{\tilde{\chi }}^{\top }(k-1)\mathfrak{P}\left(k\right)\tilde{\chi }(k-1)\hfill \\ \multicolumn{1}{c}{+2{\tilde{\chi }}^{\top }(k-1)\mathfrak{G}\left(k\right)w_v\left(k\right)+w_v^{\top }\left(k\right)\mathfrak{Q}\left(k\right)w_v\left(k\right)}\\ \multicolumn{1}{c}{\le {\lambda }_{\max }\left({\tilde{\chi }}_w\left(k\right){\tilde{\chi }}_w^{\top }\left(k\right)\right)(\mathrm{Tr}\{\mathfrak{P}\left(k\right)\}+\mathrm{Tr}\{\mathfrak{Q}\left(k\right)\})}\end{array}$$

(32)

Based on (30) and (32), “ $min(\mathrm{Tr}\{\mathfrak{P}\left(k\right)\}+\mathrm{Tr}\{\mathfrak{Q}\left(k\right)\})$ ” can be selected as the optimization objective to determine $\mathcal{K}\left(k\right)$ .

Moreover, the stability of the tracking error should be considered. It follows from (25) and the similar derivation in Ref., ³⁸ ${\tilde{\chi }}^{\top }\left(k\right)\tilde{\chi }\left(k\right)$ is asymptotically bounded when the following inequalities hold:

$$\{\begin{array}{c}\begin{array}{c}{\Vert {\mathcal{A}}_{\mathcal{K}}\left(k\right)\Vert }_2<1\hfill \\ \multicolumn{1}{c}{{\Vert {\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)\Vert }_2\le {\xi }_{\mathcal{K}}}\end{array}\hfill \end{array}$$

(33)

where ${\xi }_{\mathcal{K}}$ is a scalar, ${\mathcal{A}}_{\mathcal{K}}\left(k\right)$ and ${\mathcal{K}}_{\mathcal{B}}^{\mathcal{C}}\left(k\right)$ have been defined in (26). Similarly, by using Schur complement lemma and Lemma 1 (33) shall be hold if and only if the second and the third inequalities in (22) hold. Since $\varpi \left(k\right)$ and $v\left(k\right)$ are bounded noises, the $w_v\left(k\right)$ in (25) is also bounded. Therefore, the convex optimization problem (22) is established to determine the stability estimator gain $\mathcal{K}\left(k\right)$ .

Simulation example

The Ningbo Metro Line 7 is used in this section to validate the tracking performance of the designed nonlinear estimator (17).

In detail, this line consists of 25 stations and has a total length of approximately 40 km. Meanwhile, the train has an unladen weight of 212,904 kg and a fully loaded weight of 339,864 kg. Since different operating sections have different passenger loads, the train mass data simulated during off-peak hours are shown in Table 1. By using the train’s mass, rotational mass, maximum traction force, and maximum braking force, the maximum and minimum accelerations of the train are 16,848 and −14,139.36 km/ ${\mathrm{h}}^2$ , respectively. Moreover, the safety speed, reference speed, and slope gradient of the motion train throughout the whole line are presented in Figure 1.

Table 1.

Train mass of off-peak hours.

Section no.	Running section	Passenger mass (kg)	Train mass (kg)
1	YLZ-XYJ	10,340	223,240
2	XYJ-DQH	12,180	225,080
3	DQH-CXL	13,970	226,870
4	CXL-FQD	15,830	228,730
5	FQD-BZD	17,640	230,540
6	BZD-SML	19,450	232,350
7	SML-BMC	21,130	234,030
8	BMC-MAD	22,790	235,690
9	MAD-XTL	24,420	237,320
10	XTL-FMG	25,910	238,810
11	FMG-TYG	27,260	240,160
12	TYG-SGL	28,470	241,370
13	SGL-WTD	29,540	242,440
14	WTD-DJY	28,300	241,200
15	DJY-HCB	27,050	239,950
16	HCB-KQN	25,700	238,600
17	KQN-NCL	24,250	237,150
18	NCL-BZX	22,700	235,600
19	BZX-JLD	21,150	234,050
20	JLD-ZHD	19,500	232,400
21	ZHD-JHL	17,750	230,650
22	JHL-GSZ	16,000	228,900
23	GSZ-MHD	14,250	227,150
24	MHD-YFL	12,500	225,400

Station name abbreviation: YLZ: YunLong Station; XYJ: XiaoYangJiang Station, etc.

Figure 1.

The safety speed, reference speed and gradient of the Ningbo line 7.

The basic resistance parameter in (10) is set as ${\mathcal{H}}_m=[0.01,0.03,0.08]$ , and the state-dependent matrix in (19) is taken as $\mathcal{P}\left(k\right)=\mathrm{diag}\{0.1,0.1,0.1\}$ . Moreover, ${\mathcal{A}}_f(k-1)$ can be calculated by (20) that

$${\mathcal{A}}_f(k-1)={\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \mathcal{A}\left(x^{*}\left(k\right)\right)\hfill \end{array}\right]}_{x^{*}\left(k\right)=\hat{\chi }(k-1)}$$

(34)

where $\mathcal{A}\left(x^{*}\left(k\right)\right)\stackrel{\Delta }{=}\frac{T}{(1+\gamma )}({\stackrel{·}{\mathcal{F}}}_m\left(k\right)\mathcal{E}\left(k\right)u\left(k\right)-2\alpha \left(k\right)x^{*}-\beta \left(k\right))$ , and $\gamma =0.1$ . The uncertain disturbance of unit mass is $w\left(k\right)=0.03{\sigma }_w\left(k\right)-0.01$ and $v\left(k\right)=0.003{\sigma }_v-0.001$ , where ${\sigma }_w\left(k\right)$ and ${\sigma }_v$ are uniformly distributed random variables over [0, 1].

By implementing the tracking algorithm that proposed in Section 3, the actual position trajectory of running section YLZ-XYJ, the tracking trajectory of the proposed method, EKF, UKF and CKF methods are shown in Figure 2, which shows that all nonlinear tracking methods can track the train’s motion trajectory well. Figure 3 shows the safety speed, true speed of the train, and the speed tracking effect of the proposed method, EKF, UKF and CKF methods are also plotted. It can be seen from this figure that the proposed method has a better tracking performance.

Figure 2.

The position tracking trajectory of proposed method and different nonlinear estimation methods for section YLZ-XYJ.

Figure 3.

The speed tracking trajectory of proposed method and different nonlinear estimation methods for section YLZ-XYJ.

To clearly display the tracking precious of the proposed method, the root mean square error (RMSE) is employed to assess the tracking performance. Meanwhile, the classical nonlinear estimation methods (such as: EKF, ¹⁶ UKF, ¹⁹ CKF ²⁰ ) are presented to show the advantages of the proposed methods. Then, the RMSEs of these nonlinear estimation methods with 100 Monte Carlo simulations are given in Figure 4. It can be visualized from this figure that the proposed method in this paper has better performance in both position and velocity tracking. In particular, a numerical table of RMSEs of different methods is presented in Table 2, which shows the proposed method has a lower tracking error intuitively.

Figure 4.

The comparison of tracking performance under bounded disturbances: (a) The position error of the proposed method and EKF, UKF and CKF methods; (b) The speed error of the proposed method and EKF, UKF and CKF methods.

Table 2.

The RMSEs of the position and speed tracking errors of proposed and EKF, UKF and CKF methods, under bounded disturbances.

RMSE	Proposed method	EKF 16	UKF 19	CKF 20
Position	0.0146	0.0875	0.1385	0.1385
Speed	0.0535	0.0863	0.0665	0.0665

Since the proposed method is mainly dealing with bounded disturbances, the EKF, UKF, and CKF requires a priori statistical information of the disturbance. Therefore, assuming that the motion and measurement disturbance of a running train obey Gaussian distribution, and statistical information are $\mathcal{N}(0,0.3)$ and $\mathcal{N}(0.03)$ , respectively. Under this case, the position and speed tracking errors of different tracking methods are given in Figure 5 and Table 3. It is seen from the figure that the tracking performance of the position is similar, and EKF has a better performance in speed tracking. Meanwhile, the RMSEs in Table 3 more intuitively verify the above results. This is mainly because the design of conventional nonlinear filters is highly dependent on precisely known system noises and thus is more applicable to Gaussian noise with known covariance. Note that uncertain disturbances in train operation often do not obey a Gaussian distribution, and their statistical properties are difficult to obtain precisely but are bounded. The method proposed in this paper is not only applicable to Gaussian disturbances but also to bounded disturbances and thus has better applicability for practical applications.

Figure 5.

The comparison of tracking performance under Gaussian disturbances: (a) The position error of the proposed method and EKF, UKF and CKF methods; (b) The speed error of the proposed method and EKF, UKF and CKF methods.

Table 3.

The RMSEs of the position and speed tracking errors of proposed and EKF, UKF and CKF methods, under Gaussian disturbances.

RMSE	Proposed method	EKF 16	UKF 19	CKF 20
Position	0.0402	0.0399	0.0434	0.0434
Speed	0.1109	0.0699	0.1096	0.1096

Conclusion

In this article, a robust recursive optimization estimation method was proposed for speed tracking of subway system. Specifically, the nonlinear model with dynamic descriptions was formulated to describe the actual train running. Then, a Kalman-like nonlinear estimator was designed by modeling the linearization error as state-dependent matrix with uncertain parameters. Particularly, a convex optimization problem was established to determine estimator gain, and the bounded stability conditions of the designed nonlinear estimator were also derived. Finally, the effectiveness of the proposed method was verified by the actual operational data of Ningbo Metro Line 7.

Noted that the key to the proposed robust estimation method relies on solving the convex optimization problem, the computational complexity is mainly determined by the dimension and the introduced unknown parameters of the linear matrix inequalities. In this case, the proposed algorithm is more computationally intensive in high-dimensional systems. Therefore, how to design robust estimation algorithms with low computational complexity is one of our future works. Moreover, the limitation of network resources makes it difficult to transmit large amounts of information concurrently, and uncertainties such as delay, packet loss, and timing errors are inevitably incurred in the transmission process; how to design effective networked estimation algorithms will also be our future work.