# High-Speed Train Emergency Brake Modeling and Online Identification of Time-Varying Parameters

Authors: Yongze Jin; Guo Xie; Pang Chen; Xinhong Hei; Wenjiang Ji; Jinwei Zhao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

By analyzing the mechanism of pure air emergency brake for high-speed train, the discrete emergency brake model is established. Aiming at the problem that time-varying hidden parameters cannot be observed directly, the sliding window-based expectation maximization is proposed, and the unobserved time-varying brake parameters are identified. Firstly, the position and size of the sliding window are selected; then, the sliding window-based expectation maximization is used for brake parameter identification; finally, combined with the gradient optimization, the optimal identifications of emergency brake parameters are obtained. The simulation results show that the brake parameters can be identified quickly and accurately by the proposed method. Under uniform noise, the identification errors of friction coefficient and braking ratio are ±0.0068 and ±0.0349, respectively, and the maximum relative errors between the identifications and true values are 2.4807% and 1.3154%, respectively, which can meet the actual requirements of the brake system. The effectiveness and practicability of the proposed model and method are verified.


---

## Body

## 1. Introduction

As an important part of the rail transportation system [1–3], high-speed train plays an important role in national economic development, cultural exchange, and urban construction [4–6]. As a crucial part to ensure the safe operation of the train, the braking performance of the train braking system has received great attention. Emergency braking, the last safety barrier of train braking system, should be paid more attention. Accurate dynamic modeling is the basis of precise braking. Through the real-time accurate estimation of the performance parameters of braking model, the real-time performance state of the train can be obtained, which greatly ensures the driving safety and improves the driving efficiency [7–9]. Therefore, it is of great significance to study on high-speed train emergency brake modeling and online identification of time-varying parameters.At present, important progress has been made in the research of braking modeling for high-speed trains. Considering resistance friction, aerodynamic resistance, and inherent nonlinear braking characteristics of trains, a single coordinate multiparticle braking model was established by Song et al., which reflects the resistance and transient effects on the trains [10]. Integrating the train system dynamics and braking dynamics, the train multiparticle braking dynamics model was established by Song et al., which has strong practicability in reflecting train linear/nonlinear parameter uncertainties and possible braking faults [11]. In order to effectively utilize the regenerative energy of braking train and reduce the actual energy consumption during running, a cooperative train braking control model was put forward by Su et al., which can minimize the actual energy consumption effectively [12]. In order to solve the problem of tracking and controlling the position and speed of high-speed train connected by couplers, an elastic collision dynamics model reflecting the nonlinear traction/braking and driving faults of asymmetric vehicles was proposed by Song and Song [13]. On the basis of researching the braking principle of EMU train, the Hammerstein model of EMU braking system was put forward by Guo and Xie, according to the transmission of braking instructions and the generation of braking power [14].In the parameter estimation of train modeling, deep research has been conducted. For the dynamic model of high-speed train, a generalized weighted parameter estimation was proposed by Wang et al., and the controller was designed according to the estimation error [15]. Zhang et al. established an alternative model based on neural network for the complex structure and numerous parameters of trains [16]. By analyzing the composition and characteristics of the braking system of urban rail train and its interface with the driver, the braking model suitable for controller design and the identification of model parameters was proposed by Yu and Chen [17]. Zhong et al. constructed the discrete state space system of single particle model of high-speed train. Aiming at the non-Gaussian characteristics of train noise, the nonlinearity of the train model, and the incomplete measurability of train state, a maximum likelihood estimation suitable for the nonlinear model of high-speed train was proposed, and the constant performance parameters of the train were estimated [18]. By analyzing the longitudinal dynamic characteristics of high-speed train, combining particle filter and Bayesian posterior probability, based on expectation maximization criterion, an online estimation method for time-varying parameters of the longitudinal dynamic model of high-speed train was proposed [19].The above research has made great progress in the establishment of train model and parameter identification. However, the previous research mainly focused on the constant model parameters of the train system, only estimating the system state or constant parameters and taking offline estimation as the main method. However, due to the difference of train characteristics, the complexity of running environment, and mechanical wear, there must be uncertain and time-varying parameters in the dynamic model of train. Moreover, only Gaussian noise is used to simulate the interference of random factors such as pedestrians, road conditions, and weather, which inevitably leads to the problem of information loss, leading to the inaccuracy of the train model [20–22]. In view of the above problems, the emergency braking model of high-speed train is established by analyzing the mechanism of train emergency braking. Aiming at the time-varying hidden variable parameters which are difficult to observe in the model, an online identification of hidden variable parameters based on sliding window and maximum expectation is proposed, which realizes the real-time acquisition of dynamic characteristics and running state of high-speed train.


## 2. Modeling of Emergency Braking for High-Speed Trains

Emergency braking is one of the most critical devices in the high-speed train braking system. When there are major failures or accidents affecting train safety during train operation, it is the last safeguard to ensure the safety of high-speed train and passengers. Therefore, it is of great significance to carry out deep research on emergency braking modeling for high-speed trains. In this paper, pure air braking is regarded as the only effective braking mode for train emergency braking. By studying the braking mechanism of the brake device, the force analysis of the train during braking is carried out. Referring to theRegulations on Railway Train Traction Calculation, the discrete emergency braking model of high-speed train is established as follows:(1)st+1vt+1=st+Tvt−1−vt2×3.6vt−3.6TBPt,ηt,γt,μ+R1v+R2v,xM+wt,yt=10st+1vt+1+et,where s and v indicate the braking distance and running speed of the train, respectively, x is the train location, T is the train sampling time, M is the sum of train quality and rotation quality, wt is the random interference on train speed and displacement, which is caused by track conditions, weather conditions, and other factors, et is used to characterize measurement errors, and BPt,ηt,γt,μ represents the braking force of a train. The braking force is not only affected by the braking performance parameters of the train but also affected by the adhesion force between wheels and rails.(2)BPt,ηt,γt,μ=π4rRcd2⋅Pt⋅ηt⋅γt⋅μ⋅N⋅10−6,BPt,ηt,γt,μ≤μA⋅M⋅g,where d is brake cylinder diameter, Pt is air pressure of brake cylinder, ηt is transmission efficiency of foundation brake device, γt is the braking ratio, r is the friction radius of the brake disc, Rt is the diameter of the wheel, N is the total number of vehicle brake pads, μ and μA are friction coefficient of brake disc and adhesion coefficient between wheel and rail, respectively (the values of them vary with the braking speed v), mA is the axle load of train, and g is the gravity acceleration 9.8m/s2.R1v represents the basic running resistance of the train; it is positively correlated with the train speed v. It is usually expressed as(3)R1v=M⋅c0+c1⋅vt+c2⋅v2⋅g⋅10−3,where c0 is the rolling resistance coefficient, c1 is the other mechanical resistance coefficient, which is proportional to the train speed v, and c2 is air resistance coefficient, which is proportional to the square of train speed v2.R2v,x represents the aggregation of additional resistance of ramps, curves, and tunnels. According to the Regulations on Railway Train Traction Calculation, it can be obtained as(4)R2v,x=Rrx+Rcx+Rtx.Rrx represents the additional resistance of the ramp at x:(5)Rrx=Mgsinsin−1hl,where h and l represent the height and length of the slope, respectively. Rcx represents the additional resistance of the curve at x:(6)Rcx=0.6MgR,Lo<LC,0.0105MgαLo,Lo≥Lc,where R represents the curve radius of track, Lo and Lc are the lengths of train and curve, respectively, and α is the deflection angle.Rtx represents the additional resistance of the tunnel at x:(7)Rtx=0.00013MgL,where L represents the length of the tunnel.The train braking parameters are estimated accurately and in real time; it is of great significance to accurately grasp braking distance, ensure operation safety, and improve operation efficiency. However, in the actual operation process of trains, the range of braking parameters is determined by the train operation experience, and then the so-called maximum braking distance can be obtained. The time-varying, drift, and unobserved characteristics of braking parameters cannot be considered by the above method, which poses a great threat to the safe operation of trains. At present, the estimation of constant parameter in linear model with Gaussian noise has been a mature solution. However, as a typically time-varying parameter identification problem of nonlinear and non-Gaussian complex system, it is still lacking an effective theory and method. Therefore, this paper improves the traditional expectation maximization, and a sliding window-based expectation maximization was proposed for online identification of time-varying braking parameters of high-speed train emergency braking model.


## 3. Online Identification of Emergency Braking Model for High-Speed Trains

EM (expectation maximization) [23, 24] follows the maximum likelihood criterion; unobservable fixed parameters can be estimated by EM from incomplete datasets; the offline estimation of parameters is obtained [25]. For the identification of the train braking model, EM has some shortcomings, such as slow identification speed, low computational efficiency, and offline non-time-varying identification. Therefore, the sliding window-based expectation maximization is proposed, and the unobserved parameters are identified by this method.


### 3.1. EM Identification of Braking Model for High-Speed Train Based on Sliding Window

Sliding window is a common data processing model [26, 27]. A large number of continuous, infinite, and fast observation data of high-speed train can be processed by sliding window. The starting and ending positions of windows are set, and the current window data are selected for analysis. The starting and ending positions of the window will change with the arrival of the new data of train measurement. Old data will be removed from the window and new data will be added to the window. Therefore, the incoming data can be processed by the sliding window. And only the latest train data in the window can be processed, which accords with the actual measurement scene of the train braking process.Equation (1) shows that the train model has Markov probabilistic properties, namely,(8)xt+1∼pθxt+1|xt,yt∼pθyt|xt,where pθxt+1|xt denotes the probability density of xt+1 when the train state xt is given, pθyt|xt denotes the probability density of yt when the train state xt is given, and θ is a vector of unknown train braking parameters which may be nonlinear and time-varying, such as(9)θ=μγ.This paper divides the sliding window size based on the data record count. The measurement att1 is set as the front of the window, and the measurement at t2 is set as the back of the window. The window size is S, and the sliding length is L. Taking the mth window as the object, the identification result of braking parameter set θ in m−1th interval is initialized as θ^m−1:(10)θm=θ^m−1.The sequence of window output vectors and their likelihood functions are defined as(11)Yt1:t2=yt1,…,yt2,(12)pθYt1:t2=pθyt1,…,yt2=pθyt1∏t=t1+1t2pθyt|Yt1:t−1,where pθytYt1:t−1 is the probability density of yt when the train displacement vector sequence Yt1:t−1=yt1,…,yt−1 is given. Considering the characteristic of maximum likelihood identification, the displacement vector sequence is logarithmized on both sides of equation (12), and the logarithmic form of the sequence of output vectors is obtained.(13)LθYt1:t2=lnpθYt1:t2=lnpθyt1+∑t=t1+1t2lnpθyt|Yt1:t−1.Aiming at equation (13), combined with the Bayesian theorem and Markov characteristics of model (1), pθytYt1:t−1 is obtained as follows:(14)pθyt|Yt1:t−1=∫pθyt|xtpθxt|Yt1:t−1dxt,where pθytxt denotes the probability density of yt when the train state xt is given. Because the train model contains an unmeasurable braking parameter set θ, it is difficult to calculate pθytxt. In this paper, expectation maximization is used to indirectly optimize the logarithmic likelihood function of equation (14). The train state sequence Xt1:t2=xt1,…,xt2 is regarded as incomplete measurable data, and the joint probability density function of all data is obtained from Xt1:t2 and train output observation data Yt1:t2.(15)LθXt1:t2,Yt1:t2=lnpθXt1:t2,Yt1:t2=lnpθXt1:t2|Yt1:t2+lnpθYt1:t2.Considering the idea of minimum variance estimation, assuming that the current estimated value of train parameters isθk, the expected value Qθ,θk of LθXt1:t2,Yt1:t2 is calculated by combining θk with output observation sequence Yt1:t2.(16)Qθ,θk=∫LθXt1:t2,Yt1:t2pθkXt1:t2|Yt1:t2dXt1:t2=LθYt1:t2+∫lnpθXt1:t2|Yt1:t2pθkXt1:t2|Yt1:t2dXt1:t2,where pθkXt1:t2|Yt1:t2 is the probability density of state Xt1:t2 when the train parameter is θk and output sequence Yt1:t2 is given.


#### 3.1.1. Construction of Conditional Expectation of High-Speed Train

Based on the above deduction, combined with Bayesian theorem and Markov characteristics of braking model (1),LθXt1:t2,Yt1:t2 is expressed as follows:(17)LθXt1:t2,Yt1:t2=lnpθYt1:t2|Xt1:t2+lnpθXt1:t2=lnpθxt1+∑t=t1t2−1lnpθxt+1|xt+∑t=t1t2lnpθyt|xt.By substituting equation (17) into equation (16), the conditional expectation Qθ,θk can be obtained:(18)Qθ,θk=I1+I2+I3,(19)I1=∫lnpθxt1pθkxt1Yt1:t2dxt1,(20)I2=∑t=t1t2−1∬lnpθxt+1xtpθkxt+1,xtYt1:t2dxtdxt+1,(21)I3=∑t=t1t2∫lnpθytxtpθkxtYt1:t2dxt,where I1, I2, and I3 are related to the conditional smoothing probability density pθkxt1Yt1:t2, pθkxt+1,xtYt1:t2, and pθkxtYt1:t2, respectively. However, the train state contains unmeasurable variables, which make it difficult to analytically calculate by integration. Using the idea of particle filter and particle smoothing for reference, the integral calculation is solved by the particle method, pθxtYt is estimated by particle filter, and pθkxtYt1:t2 and pθkxt+1,xtYt1:t2 are estimated by particle smoother; ωti denotes the weight of particles after particle filtering, and ωtt2i denotes the weight of particles after smoothing.(22)I1=∑i=1Mωt1t2ilnpθx˜t1i,(23)I3=∑t=t1t2∑i=1Mωtt2ilnpθytx˜ti,where I1 and I3 can be solved by equations (22) and (23). For I2, however, equation (20) involves the joint conditional smoothing density pθkxt+1,xtYt1:t2 of two continuous time state variables xt and xt+1. Therefore, two groups of particles are selected to discretize the integral variables xt and xt+1 of I2 so that I2 also can be expressed as the sum of particles:(24)I2=∑t=t1t2−1∑i=1M∑j=1Mωtt2ijlnpθx˜t+1jx˜ti,where(25)ωtt2ij=ωtiωt+1t2jpθkx˜t+1jx˜ti∑l=1Mωtlpθkx˜t+1jx˜tl.So far, the three components ofI1, I2, and I3 of conditional expectation Qθ,θk are solved:(26)Qθ,θk=I1+I2+I3.


#### 3.1.2. Maximization of Conditional Expectations of High-Speed Train

By analyzing the train braking model, it can be seen that the relationship between conditional mathematical expectationQ^θ,θk and train braking parameters is nonlinear and nonconvex, and it is difficult to obtain its closed solution. Therefore, the gradient descent optimization [28] is used to find θ which maximizes Q^θ,θk, where θ represents a set of unknown braking parameters, and the partial derivative of Q^θ,θk with respect to parameter θ is as follows:(27)∂∂θQ^θ,θk=∂I^1∂θ+∂I^2∂θ+∂I^3∂θ,where(28)∂I^1∂θ=∑i=1Mωt1t2i∂lnpθx˜t1i∂θ,∂I^2∂θ=∑t=t1t2−1∑i=1M∑j=1Mωtt2ij∂lnpθx˜t+1jx˜ti∂θ,∂I^3∂θ=∑t=t1t2∑i=1Mωtt2i∂lnpθytx˜ti∂θ.After calculating the partial derivative ofQ^θ,θk with respect to parameter θ, the gradient method is used to optimize the parameters iteratively. The iteration and updating process are as follows:(29)αj+1=αj+βjpj,pj=Hjgj,gj=∂∂θQ^θ,θk,where Hj is a positive definite matrix to determine the gradient search direction of pj and βj represents a scalar step to ensure Q^αj+βjpj,θk>Q^αj,θk. Through the iteration process above, when the increment of Q^α,θk reaches the set range, the process of iteration optimization is cut off. The parameter estimation θk+1 makes Qθk+1,θk>Qθk,θk and ensures Lθk+1>Lθk so that θk+1 is a better estimation than θk.When the identification effect of braking parameters satisfiesθk+1−θk≤ɛ, the algorithm terminates and regards θk+1 as the final identification result θ^m of the braking parameters set of the mth window; otherwise, the algorithm continues to iterate:(30)θ^m=θk+1.


### 3.2. EM Identification for Sliding Window of Train Braking Model

Based on the analysis of the previous sections, an online sliding window-based expectation maximization identification of the high-speed train braking model is obtained. The flowchart of the algorithm is shown in Figure1. The specific steps of the algorithm are summarized as follows:(i)
Step 1. The observation data Y of the braking model of high-speed train are collected and recorded. The braking parameter θ identified is initialized. The size of sliding window is S, the sliding length of window is L, and the serial number of window is m=1. The train observation data Y1:S are identified by EM.(ii)
Step 2. Construct train conditional expectation (E): implement particle filter and particle smoothing to calculate I1,I2,I3, and the train expectation is Q^θ,θk=I1+I2+I3.(iii)
Step 3. Maximize train expectations Q^θ,θk (M): the parameter is estimated as θk+1=argmaxθQθ,θk. When the identification accuracy satisfies θk+1−θk≤ɛ, the iteration stops and the braking parameter identification result θ^m=θk+1 of the mth window is obtained. Otherwise, it returns to Step 2 to continue the iteration optimization.(iv)
Step 4. The window moves L length backwards: let m=m+1 and take θ^m−1 as the initial value of braking parameters. The existing window data are analyzed. Step 2 and Step 3 are iteratively executed. The result of braking parameter identification θ^m is obtained and then transferred to Step 5.(v)
Step 5. If the length of the remaining data is longer than the window size, turn to Step 4. Otherwise, the remaining data will be used as identification data. Step 2 and Step 3 are executed to obtain the identification result θ^m of the braking parameters of the last window.Figure 1
The flowchart of time-varying parameter identification.


## 3.1. EM Identification of Braking Model for High-Speed Train Based on Sliding Window

Sliding window is a common data processing model [26, 27]. A large number of continuous, infinite, and fast observation data of high-speed train can be processed by sliding window. The starting and ending positions of windows are set, and the current window data are selected for analysis. The starting and ending positions of the window will change with the arrival of the new data of train measurement. Old data will be removed from the window and new data will be added to the window. Therefore, the incoming data can be processed by the sliding window. And only the latest train data in the window can be processed, which accords with the actual measurement scene of the train braking process.Equation (1) shows that the train model has Markov probabilistic properties, namely,(8)xt+1∼pθxt+1|xt,yt∼pθyt|xt,where pθxt+1|xt denotes the probability density of xt+1 when the train state xt is given, pθyt|xt denotes the probability density of yt when the train state xt is given, and θ is a vector of unknown train braking parameters which may be nonlinear and time-varying, such as(9)θ=μγ.This paper divides the sliding window size based on the data record count. The measurement att1 is set as the front of the window, and the measurement at t2 is set as the back of the window. The window size is S, and the sliding length is L. Taking the mth window as the object, the identification result of braking parameter set θ in m−1th interval is initialized as θ^m−1:(10)θm=θ^m−1.The sequence of window output vectors and their likelihood functions are defined as(11)Yt1:t2=yt1,…,yt2,(12)pθYt1:t2=pθyt1,…,yt2=pθyt1∏t=t1+1t2pθyt|Yt1:t−1,where pθytYt1:t−1 is the probability density of yt when the train displacement vector sequence Yt1:t−1=yt1,…,yt−1 is given. Considering the characteristic of maximum likelihood identification, the displacement vector sequence is logarithmized on both sides of equation (12), and the logarithmic form of the sequence of output vectors is obtained.(13)LθYt1:t2=lnpθYt1:t2=lnpθyt1+∑t=t1+1t2lnpθyt|Yt1:t−1.Aiming at equation (13), combined with the Bayesian theorem and Markov characteristics of model (1), pθytYt1:t−1 is obtained as follows:(14)pθyt|Yt1:t−1=∫pθyt|xtpθxt|Yt1:t−1dxt,where pθytxt denotes the probability density of yt when the train state xt is given. Because the train model contains an unmeasurable braking parameter set θ, it is difficult to calculate pθytxt. In this paper, expectation maximization is used to indirectly optimize the logarithmic likelihood function of equation (14). The train state sequence Xt1:t2=xt1,…,xt2 is regarded as incomplete measurable data, and the joint probability density function of all data is obtained from Xt1:t2 and train output observation data Yt1:t2.(15)LθXt1:t2,Yt1:t2=lnpθXt1:t2,Yt1:t2=lnpθXt1:t2|Yt1:t2+lnpθYt1:t2.Considering the idea of minimum variance estimation, assuming that the current estimated value of train parameters isθk, the expected value Qθ,θk of LθXt1:t2,Yt1:t2 is calculated by combining θk with output observation sequence Yt1:t2.(16)Qθ,θk=∫LθXt1:t2,Yt1:t2pθkXt1:t2|Yt1:t2dXt1:t2=LθYt1:t2+∫lnpθXt1:t2|Yt1:t2pθkXt1:t2|Yt1:t2dXt1:t2,where pθkXt1:t2|Yt1:t2 is the probability density of state Xt1:t2 when the train parameter is θk and output sequence Yt1:t2 is given.


### 3.1.1. Construction of Conditional Expectation of High-Speed Train

Based on the above deduction, combined with Bayesian theorem and Markov characteristics of braking model (1),LθXt1:t2,Yt1:t2 is expressed as follows:(17)LθXt1:t2,Yt1:t2=lnpθYt1:t2|Xt1:t2+lnpθXt1:t2=lnpθxt1+∑t=t1t2−1lnpθxt+1|xt+∑t=t1t2lnpθyt|xt.By substituting equation (17) into equation (16), the conditional expectation Qθ,θk can be obtained:(18)Qθ,θk=I1+I2+I3,(19)I1=∫lnpθxt1pθkxt1Yt1:t2dxt1,(20)I2=∑t=t1t2−1∬lnpθxt+1xtpθkxt+1,xtYt1:t2dxtdxt+1,(21)I3=∑t=t1t2∫lnpθytxtpθkxtYt1:t2dxt,where I1, I2, and I3 are related to the conditional smoothing probability density pθkxt1Yt1:t2, pθkxt+1,xtYt1:t2, and pθkxtYt1:t2, respectively. However, the train state contains unmeasurable variables, which make it difficult to analytically calculate by integration. Using the idea of particle filter and particle smoothing for reference, the integral calculation is solved by the particle method, pθxtYt is estimated by particle filter, and pθkxtYt1:t2 and pθkxt+1,xtYt1:t2 are estimated by particle smoother; ωti denotes the weight of particles after particle filtering, and ωtt2i denotes the weight of particles after smoothing.(22)I1=∑i=1Mωt1t2ilnpθx˜t1i,(23)I3=∑t=t1t2∑i=1Mωtt2ilnpθytx˜ti,where I1 and I3 can be solved by equations (22) and (23). For I2, however, equation (20) involves the joint conditional smoothing density pθkxt+1,xtYt1:t2 of two continuous time state variables xt and xt+1. Therefore, two groups of particles are selected to discretize the integral variables xt and xt+1 of I2 so that I2 also can be expressed as the sum of particles:(24)I2=∑t=t1t2−1∑i=1M∑j=1Mωtt2ijlnpθx˜t+1jx˜ti,where(25)ωtt2ij=ωtiωt+1t2jpθkx˜t+1jx˜ti∑l=1Mωtlpθkx˜t+1jx˜tl.So far, the three components ofI1, I2, and I3 of conditional expectation Qθ,θk are solved:(26)Qθ,θk=I1+I2+I3.


### 3.1.2. Maximization of Conditional Expectations of High-Speed Train

By analyzing the train braking model, it can be seen that the relationship between conditional mathematical expectationQ^θ,θk and train braking parameters is nonlinear and nonconvex, and it is difficult to obtain its closed solution. Therefore, the gradient descent optimization [28] is used to find θ which maximizes Q^θ,θk, where θ represents a set of unknown braking parameters, and the partial derivative of Q^θ,θk with respect to parameter θ is as follows:(27)∂∂θQ^θ,θk=∂I^1∂θ+∂I^2∂θ+∂I^3∂θ,where(28)∂I^1∂θ=∑i=1Mωt1t2i∂lnpθx˜t1i∂θ,∂I^2∂θ=∑t=t1t2−1∑i=1M∑j=1Mωtt2ij∂lnpθx˜t+1jx˜ti∂θ,∂I^3∂θ=∑t=t1t2∑i=1Mωtt2i∂lnpθytx˜ti∂θ.After calculating the partial derivative ofQ^θ,θk with respect to parameter θ, the gradient method is used to optimize the parameters iteratively. The iteration and updating process are as follows:(29)αj+1=αj+βjpj,pj=Hjgj,gj=∂∂θQ^θ,θk,where Hj is a positive definite matrix to determine the gradient search direction of pj and βj represents a scalar step to ensure Q^αj+βjpj,θk>Q^αj,θk. Through the iteration process above, when the increment of Q^α,θk reaches the set range, the process of iteration optimization is cut off. The parameter estimation θk+1 makes Qθk+1,θk>Qθk,θk and ensures Lθk+1>Lθk so that θk+1 is a better estimation than θk.When the identification effect of braking parameters satisfiesθk+1−θk≤ɛ, the algorithm terminates and regards θk+1 as the final identification result θ^m of the braking parameters set of the mth window; otherwise, the algorithm continues to iterate:(30)θ^m=θk+1.


## 3.1.1. Construction of Conditional Expectation of High-Speed Train

Based on the above deduction, combined with Bayesian theorem and Markov characteristics of braking model (1),LθXt1:t2,Yt1:t2 is expressed as follows:(17)LθXt1:t2,Yt1:t2=lnpθYt1:t2|Xt1:t2+lnpθXt1:t2=lnpθxt1+∑t=t1t2−1lnpθxt+1|xt+∑t=t1t2lnpθyt|xt.By substituting equation (17) into equation (16), the conditional expectation Qθ,θk can be obtained:(18)Qθ,θk=I1+I2+I3,(19)I1=∫lnpθxt1pθkxt1Yt1:t2dxt1,(20)I2=∑t=t1t2−1∬lnpθxt+1xtpθkxt+1,xtYt1:t2dxtdxt+1,(21)I3=∑t=t1t2∫lnpθytxtpθkxtYt1:t2dxt,where I1, I2, and I3 are related to the conditional smoothing probability density pθkxt1Yt1:t2, pθkxt+1,xtYt1:t2, and pθkxtYt1:t2, respectively. However, the train state contains unmeasurable variables, which make it difficult to analytically calculate by integration. Using the idea of particle filter and particle smoothing for reference, the integral calculation is solved by the particle method, pθxtYt is estimated by particle filter, and pθkxtYt1:t2 and pθkxt+1,xtYt1:t2 are estimated by particle smoother; ωti denotes the weight of particles after particle filtering, and ωtt2i denotes the weight of particles after smoothing.(22)I1=∑i=1Mωt1t2ilnpθx˜t1i,(23)I3=∑t=t1t2∑i=1Mωtt2ilnpθytx˜ti,where I1 and I3 can be solved by equations (22) and (23). For I2, however, equation (20) involves the joint conditional smoothing density pθkxt+1,xtYt1:t2 of two continuous time state variables xt and xt+1. Therefore, two groups of particles are selected to discretize the integral variables xt and xt+1 of I2 so that I2 also can be expressed as the sum of particles:(24)I2=∑t=t1t2−1∑i=1M∑j=1Mωtt2ijlnpθx˜t+1jx˜ti,where(25)ωtt2ij=ωtiωt+1t2jpθkx˜t+1jx˜ti∑l=1Mωtlpθkx˜t+1jx˜tl.So far, the three components ofI1, I2, and I3 of conditional expectation Qθ,θk are solved:(26)Qθ,θk=I1+I2+I3.


## 3.1.2. Maximization of Conditional Expectations of High-Speed Train

By analyzing the train braking model, it can be seen that the relationship between conditional mathematical expectationQ^θ,θk and train braking parameters is nonlinear and nonconvex, and it is difficult to obtain its closed solution. Therefore, the gradient descent optimization [28] is used to find θ which maximizes Q^θ,θk, where θ represents a set of unknown braking parameters, and the partial derivative of Q^θ,θk with respect to parameter θ is as follows:(27)∂∂θQ^θ,θk=∂I^1∂θ+∂I^2∂θ+∂I^3∂θ,where(28)∂I^1∂θ=∑i=1Mωt1t2i∂lnpθx˜t1i∂θ,∂I^2∂θ=∑t=t1t2−1∑i=1M∑j=1Mωtt2ij∂lnpθx˜t+1jx˜ti∂θ,∂I^3∂θ=∑t=t1t2∑i=1Mωtt2i∂lnpθytx˜ti∂θ.After calculating the partial derivative ofQ^θ,θk with respect to parameter θ, the gradient method is used to optimize the parameters iteratively. The iteration and updating process are as follows:(29)αj+1=αj+βjpj,pj=Hjgj,gj=∂∂θQ^θ,θk,where Hj is a positive definite matrix to determine the gradient search direction of pj and βj represents a scalar step to ensure Q^αj+βjpj,θk>Q^αj,θk. Through the iteration process above, when the increment of Q^α,θk reaches the set range, the process of iteration optimization is cut off. The parameter estimation θk+1 makes Qθk+1,θk>Qθk,θk and ensures Lθk+1>Lθk so that θk+1 is a better estimation than θk.When the identification effect of braking parameters satisfiesθk+1−θk≤ɛ, the algorithm terminates and regards θk+1 as the final identification result θ^m of the braking parameters set of the mth window; otherwise, the algorithm continues to iterate:(30)θ^m=θk+1.


## 3.2. EM Identification for Sliding Window of Train Braking Model

Based on the analysis of the previous sections, an online sliding window-based expectation maximization identification of the high-speed train braking model is obtained. The flowchart of the algorithm is shown in Figure1. The specific steps of the algorithm are summarized as follows:(i)
Step 1. The observation data Y of the braking model of high-speed train are collected and recorded. The braking parameter θ identified is initialized. The size of sliding window is S, the sliding length of window is L, and the serial number of window is m=1. The train observation data Y1:S are identified by EM.(ii)
Step 2. Construct train conditional expectation (E): implement particle filter and particle smoothing to calculate I1,I2,I3, and the train expectation is Q^θ,θk=I1+I2+I3.(iii)
Step 3. Maximize train expectations Q^θ,θk (M): the parameter is estimated as θk+1=argmaxθQθ,θk. When the identification accuracy satisfies θk+1−θk≤ɛ, the iteration stops and the braking parameter identification result θ^m=θk+1 of the mth window is obtained. Otherwise, it returns to Step 2 to continue the iteration optimization.(iv)
Step 4. The window moves L length backwards: let m=m+1 and take θ^m−1 as the initial value of braking parameters. The existing window data are analyzed. Step 2 and Step 3 are iteratively executed. The result of braking parameter identification θ^m is obtained and then transferred to Step 5.(v)
Step 5. If the length of the remaining data is longer than the window size, turn to Step 4. Otherwise, the remaining data will be used as identification data. Step 2 and Step 3 are executed to obtain the identification result θ^m of the braking parameters of the last window.Figure 1
The flowchart of time-varying parameter identification.


## 4. Simulation Results and Analysis

In order to verify the effectiveness of the identification method, the CRH3 four-motor and four-trailer high-speed train is selected as the research object, and the simulation experiment is carried out.


### 4.1. Explanation of Experimental Parameters

In the train braking simulation experiment, the main parameters of the braking device are shown in Table1.Table 1
Main braking parameters of CRH3.
Braking parameterValueTotal weight of train (t)536Maximum operating speed (km/h)350Continuous operating speed (km/h)300Brake cylinder diameter (mm)203Air pressure of brake cylinder (kPa)410Transmission efficiency0.85Braking ratio2.55Friction coefficient0.28Brake disc friction radius (mm)297.6Wheel rolling radius (mm)460Table1 shows that the friction coefficient of the brake disc and the braking ratio of the train are constant parameters, but in the actual braking process of the train, the friction coefficient of the brake disc will change nonlinearly with the change of speed, and the braking ratio of the train will fluctuate slightly with the train running. Therefore, the unobserved friction coefficient of brake disc and the braking ratio of train will be identified in this paper.According to TSI (Technical Specification for Interoperability), the brake adhesion coefficient is shown in Table2.Table 2
Braking adhesion coefficient.
Speed (km/h)0–200200–240240–300300–350Adhesion coefficient0.150.130.110.1Assuming that the train state equation is disturbed by non-Gaussian noise, taking uniform noise as an example, the observation equation is disturbed by Gaussian noise.(31)wt=w1tw2tw1t∼U−0.1,0.1w2t∼U−0.1,0.1,et∼N0,0.1.In the train model identification, the number of iterations is 50, the number of particles is 20, and the initial values of friction coefficient and braking ratio are set to random numbers within 100% of the real value. Based on the above simulation environment, 200 Monte Carlo experiments are carried out to identify the braking parameters of high-speed trains in real time.


### 4.2. Identification of Friction Coefficient of Brake Disc

In the actual braking process of train, the friction coefficient of brake disc will change nonlinearly with the change of speed. Figure2 shows the real-time identification curve of friction coefficient of brake disc. The red solid line is the real change curve of the train friction coefficient varying with running speed. The black, blue, and purple solid lines are the friction coefficient identification curves with initial values of 0.40, 0.35, and 0.03, respectively. The black, blue, and purple dashed lines are the friction coefficient identification curves with initial values of 0.15, 0.20, and 0.25, respectively.Figure 2
Real-time identification of friction coefficient by online EM.It can be seen that in the three deceleration stages of350∼320km/h, 320∼210km/h, and 150∼0km/h, the proposed method can accurately identify the friction coefficient of the brake disc in real time. In the 210∼150km/h deceleration stage, because the friction coefficient is too small and there is inflection point with higher change rate, the identification accuracy decreases slightly, but it can accurately predict the changing trend of the friction coefficient of the brake disc, which provides a reference for the identification of braking parameters and the monitoring of the braking performance of trains.Figures3 and 4 show the identification error and relative error curves of braking ratio, respectively. It can be seen from the figures that the identification error and relative error under different initial values gradually decrease with the increase of time, and the fluctuations of the two curves gradually decrease and remain almost unchanged in a fixed range after 320km/h.Figure 3
Real-time identification error of friction coefficient by online EM.Figure 4
Relative error of friction coefficient by online EM.In order to further prove the effectiveness of this method proposed in this paper, this method is compared with extended Kalman filter (EKF). EKF is a typical extension of linear Kalman filter in theory and application. Its theoretical core is to linearize the nonlinear function, that is, take the first-order expression of Taylor expansion to realize the estimation of the optimal linear problem. The innovation of its application in state estimation is the accurate calculation of the gain matrix. In this process, system noise and measurement noise are integrated into the update of the gain matrix so that the algorithm has better estimation accuracy and anti-interference ability. Figures5–7 show the real-time identification results, estimation errors, and relative errors of the friction coefficient, respectively, based on EKF under different initial conditions.Figure 5
Real-time identification of friction coefficient by EKF.Figure 6
Real-time identification error of friction coefficient by EKF.Figure 7
Relative error of friction coefficient by EKF.As can be seen from Figures5–7, compared with EKF, the identification results obtained by the proposed method in this paper can converge to the true value of the friction coefficient rapidly. In the process of convergence, the fluctuation of identification error is smaller.Table3 shows the comparison of the estimation results of friction coefficient between the method in this paper and EKF under different initial values.Table 3
Identification results of friction coefficient.
Online EMEKFInitial valueIdentification errorRelative error (%)Identification errorRelative error (%)0.15±0.00742.6730±0.01734.94440.20±0.00762.8845±0.01233.52010.25±0.00492.1075±0.01825.22550.30±0.00572.1446±0.01875.33590.35±0.00642.1186±0.01795.12630.40±0.00862.9561±0.01905.4146Mean±0.00682.4807±0.01724.9278From Table3, it can be seen that the average values of the estimation error and relative error of friction coefficient obtained by the method in this paper are ±0. 0068 and 2.4807%, respectively, and the average values of the estimation error and relative error of friction coefficient obtained by EKF are ±0.0172 and 4.9278%, respectively. The method proposed in this paper has higher estimation accuracy and better real-time performance, and the estimation results meet the actual needs of the brake system.Figure8 shows the identification curve of train braking speed by online EM and the calculation curve of braking displacement. It takes 94.2 seconds from normal operation to stop, and the braking distance is 4647.16 meters, which meets the national standard.Figure 8
Braking speed and displacement curve under the friction coefficient.


### 4.3. Identification of Train Braking Ratio

In the actual braking process of the train, the braking ratio of the train is not constant, and its value will fluctuate slightly with the train running. Figure9 shows the real-time identification curve of braking ratio of high-speed train. The red solid line is the real change curve of the train braking ratio. The black, blue, and purple solid lines are the braking ratio identification curves with initial values of 5, 4, and 3, respectively. The black, blue, and purple dashed lines are the braking ratio identification curves with initial values of 0, 1, and 2, respectively. It can be seen from the graph that the method proposed in this paper can accurately identify the braking ratio of trains in real time, and the accurate identification results of the braking ratio can be obtained after 320km/h.Figure 9
Real-time identification of braking ratio by online EM.Figures10 and 11 show the identification error and relative error curves of braking ratio, respectively. It can be seen from the figures that the identification error and relative error under different initial values gradually decrease with the increase of time, and the fluctuations of the two curves gradually decrease and remain almost unchanged in a fixed range after 320km/h.Figure 10
Identification error of braking ratio by online EM.Figure 11
Relative error of braking ratio by online EM.Figures12–14 show the real-time identification results, estimation errors, and relative errors of the braking ratio, respectively, based on EKF under different initial conditions.Figure 12
Real-time identification of braking ratio by EKF.Figure 13
Identification error of braking ratio by EKF.Figure 14
Relative error of braking ratio by EKF.As can be seen from Figures12–14, compared with EKF, the identification results obtained by the proposed method in this paper can converge to the true value of the braking ratio rapidly. In the process of convergence, the fluctuation of identification error is smaller.Table4 shows the comparison of the estimation results of braking ratio between the method in this paper and EKF under different initial values.Table 4
Identification results of braking ratio.
Online EMEKFInitial valueIdentification errorRelative error (%)Identification errorRelative error (%)0±0.05731.9890±0.05232.02121±0.02580.9980±0.04601.77672±0.04501.7380±0.04871.88143±0.03531.3833±0.04261.64864±0.02300.8903±0.04551.75825±0.02310.8940±0.04121.5932Mean±0.03491.3154±0.04611.7799From Table4, it can be seen that the average values of the estimation error and relative error of braking ratio obtained by the method in this paper are ±0. 0349 and 1.3154%, respectively, and the average values of the estimation error and relative error of braking ratio obtained by EKF are ±0.0461 and 1.7799%, respectively. The method proposed in this paper has higher estimation accuracy and better real-time performance, and the estimation results meet the actual needs of the brake system.Figure15 shows the identification curve of train braking speed by online EM and the calculation curve of braking displacement. It takes 94.7 seconds from normal operation to stop, and the braking distance is 4687.34 meters, which meets the national standard.Figure 15
Braking speed and displacement curve under the braking ratio.By analyzing the identification results of friction coefficient and braking ratio, it can be concluded that when the initial values of braking parameters are set in a specified range, the braking parameters can be identified quickly and accurately in real time, which verifies the effectiveness of the method.


## 4.1. Explanation of Experimental Parameters

In the train braking simulation experiment, the main parameters of the braking device are shown in Table1.Table 1
Main braking parameters of CRH3.
Braking parameterValueTotal weight of train (t)536Maximum operating speed (km/h)350Continuous operating speed (km/h)300Brake cylinder diameter (mm)203Air pressure of brake cylinder (kPa)410Transmission efficiency0.85Braking ratio2.55Friction coefficient0.28Brake disc friction radius (mm)297.6Wheel rolling radius (mm)460Table1 shows that the friction coefficient of the brake disc and the braking ratio of the train are constant parameters, but in the actual braking process of the train, the friction coefficient of the brake disc will change nonlinearly with the change of speed, and the braking ratio of the train will fluctuate slightly with the train running. Therefore, the unobserved friction coefficient of brake disc and the braking ratio of train will be identified in this paper.According to TSI (Technical Specification for Interoperability), the brake adhesion coefficient is shown in Table2.Table 2
Braking adhesion coefficient.
Speed (km/h)0–200200–240240–300300–350Adhesion coefficient0.150.130.110.1Assuming that the train state equation is disturbed by non-Gaussian noise, taking uniform noise as an example, the observation equation is disturbed by Gaussian noise.(31)wt=w1tw2tw1t∼U−0.1,0.1w2t∼U−0.1,0.1,et∼N0,0.1.In the train model identification, the number of iterations is 50, the number of particles is 20, and the initial values of friction coefficient and braking ratio are set to random numbers within 100% of the real value. Based on the above simulation environment, 200 Monte Carlo experiments are carried out to identify the braking parameters of high-speed trains in real time.


## 4.2. Identification of Friction Coefficient of Brake Disc

In the actual braking process of train, the friction coefficient of brake disc will change nonlinearly with the change of speed. Figure2 shows the real-time identification curve of friction coefficient of brake disc. The red solid line is the real change curve of the train friction coefficient varying with running speed. The black, blue, and purple solid lines are the friction coefficient identification curves with initial values of 0.40, 0.35, and 0.03, respectively. The black, blue, and purple dashed lines are the friction coefficient identification curves with initial values of 0.15, 0.20, and 0.25, respectively.Figure 2
Real-time identification of friction coefficient by online EM.It can be seen that in the three deceleration stages of350∼320km/h, 320∼210km/h, and 150∼0km/h, the proposed method can accurately identify the friction coefficient of the brake disc in real time. In the 210∼150km/h deceleration stage, because the friction coefficient is too small and there is inflection point with higher change rate, the identification accuracy decreases slightly, but it can accurately predict the changing trend of the friction coefficient of the brake disc, which provides a reference for the identification of braking parameters and the monitoring of the braking performance of trains.Figures3 and 4 show the identification error and relative error curves of braking ratio, respectively. It can be seen from the figures that the identification error and relative error under different initial values gradually decrease with the increase of time, and the fluctuations of the two curves gradually decrease and remain almost unchanged in a fixed range after 320km/h.Figure 3
Real-time identification error of friction coefficient by online EM.Figure 4
Relative error of friction coefficient by online EM.In order to further prove the effectiveness of this method proposed in this paper, this method is compared with extended Kalman filter (EKF). EKF is a typical extension of linear Kalman filter in theory and application. Its theoretical core is to linearize the nonlinear function, that is, take the first-order expression of Taylor expansion to realize the estimation of the optimal linear problem. The innovation of its application in state estimation is the accurate calculation of the gain matrix. In this process, system noise and measurement noise are integrated into the update of the gain matrix so that the algorithm has better estimation accuracy and anti-interference ability. Figures5–7 show the real-time identification results, estimation errors, and relative errors of the friction coefficient, respectively, based on EKF under different initial conditions.Figure 5
Real-time identification of friction coefficient by EKF.Figure 6
Real-time identification error of friction coefficient by EKF.Figure 7
Relative error of friction coefficient by EKF.As can be seen from Figures5–7, compared with EKF, the identification results obtained by the proposed method in this paper can converge to the true value of the friction coefficient rapidly. In the process of convergence, the fluctuation of identification error is smaller.Table3 shows the comparison of the estimation results of friction coefficient between the method in this paper and EKF under different initial values.Table 3
Identification results of friction coefficient.
Online EMEKFInitial valueIdentification errorRelative error (%)Identification errorRelative error (%)0.15±0.00742.6730±0.01734.94440.20±0.00762.8845±0.01233.52010.25±0.00492.1075±0.01825.22550.30±0.00572.1446±0.01875.33590.35±0.00642.1186±0.01795.12630.40±0.00862.9561±0.01905.4146Mean±0.00682.4807±0.01724.9278From Table3, it can be seen that the average values of the estimation error and relative error of friction coefficient obtained by the method in this paper are ±0. 0068 and 2.4807%, respectively, and the average values of the estimation error and relative error of friction coefficient obtained by EKF are ±0.0172 and 4.9278%, respectively. The method proposed in this paper has higher estimation accuracy and better real-time performance, and the estimation results meet the actual needs of the brake system.Figure8 shows the identification curve of train braking speed by online EM and the calculation curve of braking displacement. It takes 94.2 seconds from normal operation to stop, and the braking distance is 4647.16 meters, which meets the national standard.Figure 8
Braking speed and displacement curve under the friction coefficient.


## 4.3. Identification of Train Braking Ratio

In the actual braking process of the train, the braking ratio of the train is not constant, and its value will fluctuate slightly with the train running. Figure9 shows the real-time identification curve of braking ratio of high-speed train. The red solid line is the real change curve of the train braking ratio. The black, blue, and purple solid lines are the braking ratio identification curves with initial values of 5, 4, and 3, respectively. The black, blue, and purple dashed lines are the braking ratio identification curves with initial values of 0, 1, and 2, respectively. It can be seen from the graph that the method proposed in this paper can accurately identify the braking ratio of trains in real time, and the accurate identification results of the braking ratio can be obtained after 320km/h.Figure 9
Real-time identification of braking ratio by online EM.Figures10 and 11 show the identification error and relative error curves of braking ratio, respectively. It can be seen from the figures that the identification error and relative error under different initial values gradually decrease with the increase of time, and the fluctuations of the two curves gradually decrease and remain almost unchanged in a fixed range after 320km/h.Figure 10
Identification error of braking ratio by online EM.Figure 11
Relative error of braking ratio by online EM.Figures12–14 show the real-time identification results, estimation errors, and relative errors of the braking ratio, respectively, based on EKF under different initial conditions.Figure 12
Real-time identification of braking ratio by EKF.Figure 13
Identification error of braking ratio by EKF.Figure 14
Relative error of braking ratio by EKF.As can be seen from Figures12–14, compared with EKF, the identification results obtained by the proposed method in this paper can converge to the true value of the braking ratio rapidly. In the process of convergence, the fluctuation of identification error is smaller.Table4 shows the comparison of the estimation results of braking ratio between the method in this paper and EKF under different initial values.Table 4
Identification results of braking ratio.
Online EMEKFInitial valueIdentification errorRelative error (%)Identification errorRelative error (%)0±0.05731.9890±0.05232.02121±0.02580.9980±0.04601.77672±0.04501.7380±0.04871.88143±0.03531.3833±0.04261.64864±0.02300.8903±0.04551.75825±0.02310.8940±0.04121.5932Mean±0.03491.3154±0.04611.7799From Table4, it can be seen that the average values of the estimation error and relative error of braking ratio obtained by the method in this paper are ±0. 0349 and 1.3154%, respectively, and the average values of the estimation error and relative error of braking ratio obtained by EKF are ±0.0461 and 1.7799%, respectively. The method proposed in this paper has higher estimation accuracy and better real-time performance, and the estimation results meet the actual needs of the brake system.Figure15 shows the identification curve of train braking speed by online EM and the calculation curve of braking displacement. It takes 94.7 seconds from normal operation to stop, and the braking distance is 4687.34 meters, which meets the national standard.Figure 15
Braking speed and displacement curve under the braking ratio.By analyzing the identification results of friction coefficient and braking ratio, it can be concluded that when the initial values of braking parameters are set in a specified range, the braking parameters can be identified quickly and accurately in real time, which verifies the effectiveness of the method.


## 5. Conclusion

The mechanism of emergency braking of high-speed trains was analyzed, and the emergency braking model of high-speed train was established in this paper. Aiming at the characteristics of train braking parameters that cannot be observed intuitively which change with time, an improved expectation maximization identification of the high-speed train braking model based on sliding window was proposed, and the time-varying braking parameters of trains were identified online. The simulation results showed that the real-time values of braking parameters can be identified quickly and accurately. It can be concluded that the application of this method will greatly reduce the cost of daily maintenance and improve the efficiency of maintenance. In the follow-up study, the algorithm will be further improved and optimized to adjust the size of sliding window adaptively, and the speed and accuracy of parameter identification will be further improved.



---

*Source: 3872852-2020-06-15.xml*