GROUP-SMA ALGORITHM BASED JOINT ESTIMATION OF TRAIN PARAMETER AND STATE

The braking rate and train arresting operation is important in the train braking performance. It is difficult to obtain the states of the train on time because of the measurement noise and a long calculation time. A type of Group Stochastic M-algorithm (GSMA) based on Rao-Blackwellization Particle Filter (RBPF) algorithm and Stochastic M-algorithm (SMA) is proposed in this paper. Compared with RBPF, GSMA based estimation precisions for the train braking rate and the control accelerations were improved by 78% and 62%, respectively. The calculation time of the GSMA was decreased by 70% compared with SMA.


INTRODUCTION
With the high-speed development of rail transportation, the safety and efficiency requirements during the transport process are increasing.Urban train parking precision affects the alignment of the train door with the platform screen door (PSD).The train braking system is a crucial factor affecting the parking accuracy.The braking rate is affecting the performance of the train braking system.The long-time continuous operation and the variation of the loading conditions and operation system will lead to the inevitable change of the braking rate.The displacement and velocity measurement accuracy is affecting the braking accuracy and it is affected by the performance of the velocity sensors.The aforementioned precision problems can be converted to the joint estimation of unknown parameter and states by modelling the train braking process.In the train braking model, the velocity and displacement are states in the state space model, and the braking rate is the unknown parameter to be estimated.
The commonly used parameter identification methods are point estimation algorithms and filtering algorithms [1].The traditional point estimation algorithms contain the least squares [2], maximum likelihood estimation method [3] and EM algorithm (Expectation Maximization) [4].When the point estimation algorithms are used for the parameter estimation, it requires that the objective function is continuous and derivable.The estimated parameter value is searched by the gradient information, so it easily falls into the local minimum.The parameter identification method based on neural network [5] has the ability to approximate non-linear functions with arbitrary precision, but there are still several problems such as the determination of the network structure, sample data selection and network training algorithms.Genetic algorithm [6] has defects of premature or slow convergence.The Kalman filter and extended Kalman filter (EKF) have been applied successfully in the state estimation [7].The Interacting Multiple Model (IMM) algorithm [8] is based on Kalman Filter (KF) and it has been used for the state estimation of the hybrid system.Knowledge of the transition probability between different models must be achieved in advance, which is difficult in practical engineering.
The Particle Filter (PF) uses the non-parametric Monte Carlo method to achieve the recursive Bayesian filter.This algorithm is suitable for any non-linear systems that can be modelled by state space equation and its estimation precision is approximate to the optimal estimation.N. Metropolis and S. Ulam in [9] proposed a particle filter for the first time.Development of the particle filter is introduced in [10] and its application in sensor diagnosis and accident detection in [11,12].There exist two major defects of the method.Firstly, the particle diversity decreases after several steps and all the particles tend to be of the same value.The Auxiliary particle filter (APF) and Gaussian particle filter (GPF) were discussed to solve this problem in [13].
The second problem is that the particle filter needs a large number of samples to approximate the posterior density of the system state.Rao-Blackwellization technology in [14] reduces the number of particles by reducing the dimension of the state vector.In RBPF algorithm, the state space is divided into linear substate space estimated by Kalman filter algorithm and non-linear sub-state space estimated by particle filter algorithm [15][16][17][18].The RBPF has applications in the hybrid Gaussian [19][20][21], fixed parameters estimation [22], hidden Markov models (HMMs) [20][21], Dirichlet process models and Dynamic Bayesian Networks (DBN) [23].There is no transition function for the unknown parameter in the state space equations.After several steps of the iterative calculation, all values of the parameter particles will tend to the same if the RBPF algorithm is used for the joint train braking rate estimation and the train states.SMA [24] is more efficient than RBPF for the state estimation of the hybrid system.SMA originally grows from the stochastic Malgorithm [25], based on a random sampling.In the SMA algorithm, the unknown parameter is discrete in several values and each discrete value corresponds to a system mode.When each parameter particle is predicted, the particle is transferred to all the possible discrete values from the unknown parameter, leading to an exponential growth of the particles number.
The GSMA is proposed to reduce the computation complexity, keeping the estimation precision.The second part of this paper introduces the principle of the estimation algorithms for the hybrid system.The GSMA is proposed in the third part.After the train braking model is specified in the fourth part, the RBPF algorithm, the SMA algorithm and the GSMA are employed for the joint estimation of train braking rate and states in the fifth part.Advantages and limitations of the GSMA algorithm are given in section six.

PROBLEM FORMULATION
Train braking system main function is to achieve consistent braking performance by a brake control-ler.The driver completes the vehicle control using the vehicle traction control system and the brake control system, but the driver cannot directly manipulate the train power actuators.The brake model is a vehicle operating model including train braking control system [26][27][28].The target acceleration has the fixed linear relationship with the drive input so it is chosen as the input of the braking controller and represented by a t Tar ^h .The brake controller will track a t Tar ^h by feedback regulation.The dynamic regulation process is described as Equation ( 1): where a t t^h denotes the control acceleration and it is the output of the controller, a t t o ^h means the derivation of a t t^h after value t and f represents the braking rate which is defined as the ratio of the actual braking power with the expected braking power.The nominal value of f is 1, meaning the actual braking power is ideally equal to the expected braking power.Response time of the power actuator x represents the time needed for the braking power to reach the desired value.Transmission delay time from the generation of the target acceleration to the execution of the power actuator is T, so a t T  ing that x D v a az T = t 7 A and y D v az = 7 A, the state space equations of the model can be described as: Here x represents the state vector of the system, y denotes the observation vector of the system and m is a parameter related to the transmission delay time T. The system coefficient matrices are A, b, E and C and X represents the process Gaussian noise.The train braking model was demonstrated and verified with the practical operation data of YIZHUANG subway line in Beijing in literature [28].As can be seen in the state space model, there are parameters m , x , f to be estimated.Static parameters m and x are estimated through offline methods.

PRESENTATION OF GSMA ALGORITHM
The hybrid system, particle filter algorithm, SMA and GSMA algorithms are presented further in the text.

Hybrid system
The state of the hybrid system , x k k i ^h is composed of parameter k i and the state xk .Each system model corresponds to a set of state space equations containing state transition equations and observation equations:  [26].Here value m is an introduced coefficient calculated from T with the Pade^h function in Matlab Software and a s z ^h is considered as the value derived from a s Tar ^h [26].
In this formula T represents the transposition for the rows with columns.
i + represent the value of the discrete system parameter at discrete moment k 1 + .The input value of the system at moment k is uk .The probability transition matrix of discrete parameters is p k . The system coefficient matrices are Process noise and measurement noise at moment k are k ~ and vk and both are the independent white Gaussian noise.Parameter k 1 i + is calculated from the probability transition matrix (5).The continuous state transition equation and observation equation are ( 6) and ( 7), respectively.The state estimation is in calculating the parameter and state posterior probability density , p x Y : where Y :k 1 describes the system measurements from moment 1 to moment k in steps.

Particle Filter algorithm
(1) Filter is sequential to estimate parameter i and state x of a system given a set of observation variables.Variable z is used to represent the parameters and states variables set , x i " ,.The posterior probability density of z for , , , , is estimated by the given observation process values yi , which presents the value of y at moment i, , , , , ) Particle filters use a large random number of sample points sampled from the sampling function in the state space to approximate the posterior probability distribution of the estimated variables.(3) The procedures of particle filters are as follows [20]: Here N is the number of the sample particles.

2) For
, , , normalize the importance weights by ) Compute an estimation of the effective number of particles Neff t by:

RBPF algorithm
RBPF algorithm is commonly used for parameter and state estimation of hybrid system.From the Bayesian theory, the state posterior probability density function , Y p x The numerical solution of Y p : k k 1 i ^h can be approximately calculated through the particle filter algorithm.Given the available parameter values and the state coefficient matrices of the system, the analytic solution of , Y p x ^h can be derived using Kalman filter algorithm.If the particle filter is used to estimate the unknown parameter, ,Y p is chosen as the sampling function to approximate the unknown posterior probability density function of the parameter.The parameter particles are obtained based on ,Y p and then the probability density function of Y p denotes the probability density function of the unknown parameter k i given the history values of the parameter particle :k According to the theory of probability, ,Y p can be simplified as [12]: ,Y p py x dx In formula ( 9), , ,Y p x is derived through the Kalman filter algorithm as [12]:  9) to (11), the sampling function can be parsed as: where and the weights of the particles are calculated from: , ,Y p y Here w k j ^h represents the weight of the j-th particle at moment k and the system hybrid state estimation is obtained from the weighted sum of particle weights and particle values.
The mean and variance of continuous state corresponding to each particle can be derived using the Kalman filtering algorithm.The posterior probability density function of xk can be expressed as a mixture of Gaussian function.The particle will show the degradation property if RBPF is used for the joint estimation of unknown parameter and system states because most of the parameter particles tend to transfer to the same discrete value after sampling step according to the sampling function, as the unknown parameter does not have the transition function.During the resampling process, the particles with large weights are replicated many times while the particles with small weights are rejected, and the replicated particles represent the samples.After several re-sampling steps, all the values of the parameter particles tend to be the same.

SMA algorithm
The SMA algorithm is derived from the M-algorithm to solve the particle degradation problem occurring in the RBPF algorithm.The unknown parameter i in RBPF is assumed as discrete in . Discrete values represent the value space and any particle for i should belong to the value space.In each moment step from 0 to k, the parameter particle needs to be updated based on the sampling function.If the parameter particle , is calculated first and then the parameter value corresponding to the largest transition probability is sampled as the parameter particle k i i ^h .In each particle trajectory, although there are M possible values to which parameter particle k i 1 i - ^h can be transferred at moment k, only one value of them is finally retained.Since the diversity of the particles is important for obtaining accurate posterior probability density, it is preferable to retain more than one value which has relatively large weight.Compared with RBPF algorithm, the SMA algorithm retains all the possible values for each particle trajectory and the particle trajectories grow like a tree.After each sampling step, one particle trajectory will be extended into M particle trajectories with the results that the total number of particles at time k become M times to that of particles at precious moment k 1 -.If the number of particles is initialized as N at moment t0 , the whole number of the particles will be M N k 1 # + at moment tk 1 + .

Proposed GSMA algorithm
The unknown parameter is usually static or changes slowly in the industry domain.During the estimation process, when each particle is sampled, it is enough to retain the possible values closer to the parameter value instead of retaining all the possible values.
For instance, assuming that one parameter is discrete as ., ., ., ., . 1 0 0 8 0 5 0 3 0 1 " , .If the current estimated parameter value is 0.8, only the particle value 1.0 or 0.5 will be sampled for the next step, because the parameter is assumed to change slowly.So the estimated parameter will be first discrete and then all discrete values will be divided into different groups in steps as follows.
During each sampling step, parameter particles can transfer only within the group.For example, if discrete space of i is assumed to contain five values as , , , , , , then at least two different groups could be generated: , , i i i ^^ĥ h h " , .The parameter particle value 1 i ^h at moment k 1 -, can transfer to a value in the first group, not in the second one.
The member of the group should include the similar value in the discrete space.
Each discrete value is allowed to exist in different groups.The particle should have the transferring possibility between different groups, such as , , i i i ^^ĥ h h " , .If the value of the parameter is appointed as 1 i ^h initially and since 2 i ^h exists both in these two groups, the sampling scope for 2 i ^h is possibly in the second group.
The group, in which the average value of the member is the closest to the parameter particle, will be selected.The average values of the groups should cover all the discrete values in the state space, except the two boundary values in the group.The number of members in a group is kept as small as possible, but at least three.If the state space of estimated param- , the largest number of the groups is M 2 -.The groups could be specified as , , i i i i i i ^^ĥ h h " , ,g, , , The particle filter procedure is described in Figure 5.The discrete space of estimated values is assumed as , , , , 1 2 3 4 5 " , and the space is grouped into three groups shown as , , 1 2 3 " ,, , , 2 3 4 " ,, , , 3 4 5 " ,.The initial value of the parameter is set as 3 and the number of the initial particles is described as N0 .Figure 5 shows the increased procedure from moment 0 to moment t n + for one particle.For any t, re-sampling scheme will be executed if the number of the particles reaches the boundary maximum number Q and N particles with higher weights will be kept for the next step.
The calculation steps of the GSMA algorithm are shown as Figure 8, where S denotes the particles, N0 is set as 1 and N represents the number of resampling particles; if the space is set as , , , , , , i 1 2 3 !" ,.
Step 1: Parameter particles initialization: From In the particle filter algorithm, p k i ^h is sampled from the initial defined particle values of parameter i .The weights of the particle are set as /N 1 0 .
Step 2: States particles initialization: , where the particles for x i 0 ^h and the particles for P i 0 ^h are sam-pled on the normal distribution with the mean value x0 t and variance p0 , at moment 0.
Step 3: For , , k 1 2 f = , repeat the following steps: For , , , i N 1 2 f = , the state particle x k k i 1 -^h of the system and the mean square error matrix where is estimation mean value of the parameter particle values k i i ^h given by (5).
is the estimation mean value of y k i ^h given by (7) and it is calculated by: R k i ^h denotes the estimation standard deviation of the parameter and it is derived by: The weights for the parameter and state particles are calculated and normalized: S t The mean values of the parameters and states are calculated by:

Figure -The filter process of GSMA algorithm
where Nk means the number of the particles at moment k.Kalman gain K k i ^h is calculated by: The state particle x k i ^h of the system and the mean square error matrix P k i ^h are updated by observation value yk from (7): The mean values of the states are calculated by: The group, in which the parameter particle exists must be determined before the sampling action and the particle is extended within the group.

SIMULATION RESULTS
The RBPF algorithm, SMA algorithm and GSMA algorithm are used for the joint estimation of the train braking rate and the states of the urban train.The state space equation of the simulation system is described in (4) .The braking rate f in (4) represents the system parameter and displacement D, velocity v, control acceleration a t .The filter value of the target acceleration az denotes the estimated status of the system.
The initial condition of the system parameter and states are defined as follows: (1) It is assumed that the braking rate of the urban railway changes as: .
where initialization value f0 is set as 1.0.The scope of the braking rate f is , 0 1 6 @ and f, x and m are used to calculate the "real" state value x D v a az A by ( 4).(2) For the displacement D, velocity v, control acceleration a t and intermediate variable az in the state space, the initial values are 0 m, 20 m/s, (-1) m/s 2 and (-2) m/ s 2 .
These values are not calculated from (23).The values of parameter particles will be assigned from them.
The initial numbers of parameter and state particles are both represented by N0 and they are both set as 50.The boundary maximum upper limit particle number Q and the re-sampling particle number N are set as 500 and 50, respectively.
The initial numbers of parameter and state particles are both represented by N0 and they are both set as 50; Q and N are set as 500 and 50, respectively.
The initial numbers of parameter and state particles are both represented by N0 and they are both set as 250.
During the estimation calculating procedure, the discrete time integration is 1.0 second.The values of the observations y are generated based on the system state space equations by the simulation tool MATLAB.The hardware is PC ThinkPad R400, with Intel T6570, 2.1GHz and 2GB DDR3 memory.
In (19) and (22), k i t represents braking rate f and A. Thus f, D, v, a t and az can be up- dated and calculated by the observations generated in (23) and (4).The estimation errors of the parameters and states, calculated by ( 23) and ( 4) are shown from Figure 6 to Figure 10.The corresponding estimated mean errors are specified in Table 1.In (4), y is observable output and y D v a 0 z T = 7 A Displacement D of the train is measured by displacement sensor and velocity v is calculated from the displacement at the moment.The az is derived from the target acceleration and can be calculated from the system input a t Tar ^h .From the system parameter and status shown in Figure 7, Figure 8, Figure 10 and in Table 1, the error curves of urban train displacement D, velocity v and az are almost equal because these three variables are observable.Figure 6 and Figure 9 show that the estima-tion precisions of the train braking rate f and the control acceleration a t using GSMA and SMA are similar.The estimation accuracy of GSMA and SMA are much higher than that of RBPF algorithm.The computation times of the GSMA, SMA and the RBPF algorithm in 150 recursive steps are shown in Table 2.The real-time performance of the GSMA is better than that of SMA and weaker than that of RBPF.Although the RBPF has the best real-time estimation, its estimation error for the braking rate and control acceleration is relatively larger, which will bring about the unexpected hidden danger during the train operation process.Compared with RBPF, GSMA has the better estimation precision the GSMA is most suitable for the joint estimation of the train braking rate and the train operation states.

DISCUSSION AND CONCLUSION
In order to improve the parking precision of the urban train, the GSMA is proposed, based on the RBPF algorithm and SMA.Adequate estimation for the dynamic braking rate, the displacement and velocity measurement accuracy are the key factors affecting the parking precision.Compared with the RBPF and SMA, GSMA reduces the estimation error of braking rate greatly so that the parking precision of the urban train is improved.In addition, GSMA keeps the realtime performance by decreasing the number of the calculating particles during the estimation process.
In general, the braking rate of urban train changes occasionally and slowly so that the longer parameter estimation time of GSMA, compared with RBPF, is tolerable.
In the urban railway braking domain the benchmark is the accuracy requirement, so that the proposed GSMA takes precedence on the accuracy requirements.If the range of the unknown parameter is wide and the number of the discrete values for the discrete parameter is big, the computation time advantage of the GSMA will become more apparent, because the computation complexity of the proposed algorithm is only associated with the number of discrete values in one group.Whether the unknown parameter is static or dynamic, GSMA can obtain accurate estimation result if the range of the parameter is known in advance.If there is little knowledge of the unknown parameter, it is difficult to use the GSMA.
In the future work, the GSMA will be improved for estimation of the unknown parameter with little a priori knowledge by combining it with an algorithm that can determine the range of the parameter.

Figure 2 -ure 1
Figure 2 -Laplace transformation blocks of urban train braking model

Figure 3 Figure 3
Figure 3(a) -Approximation transfer block of urban train braking model

/ 5 ) 6 )
If the effective number of particles Neff t is less than the given threshold Nthr , then resample in steps: (a) Draw N particles from the current particle set with probabilities proportional to their weights.Replace the current particle set with the new one.(b) For , , go back to step 2).
state particles and parameter particles of the system at moment k.Here, state mean and variance of state x k j ^h , respectively for the given state x k j 1 -^h .The variance of process noise wk in (6) is Q~.From Equations (

Figure - 2 Figure 10 -
Figure -Comparison of the estimation error of velocity with different algorithms 8 v By inverse Laplace transformation, values fromFigure 3(b) are transformed to Figure 4.
D t ^h velocity v t ^h , control acceleration a t t^h and intermediate variables a t z ^h are selected as the state variables.
border of the integral is the value space of xk .
. Zheng, J. Han, W. Kong, L. Wang: Group-SMA Algorithm Based Joint Estimation of Train Parameter and State W . The static parameters x and m are set as .

Table 1 -
Comparison of the parameter and status estimation using RBPF algorithm, SMA and GSMA

Table 2 -
Comparison of computing time using RBPF algorithm, SMA and GSMA