# Optimal Operation of High-Speed Trains Using Hybrid Model Predictive Control

Authors: Yingze Yang; Zheng Xu; Weirong Liu; Heng Li; Rui Zhang; Zhiwu Huang
Journal: Journal of Advanced Transportation
 (2018)

---

## Abstract

The high-speed train operation process is highly nonlinear and has multiple constraints and objectives, which lead to a requirement for the automatic train operation (ATO) system. In this paper, a hybrid model predictive control (MPC) framework is proposed for the controller design of the ATO system. Firstly, a piecewise linear system with state and input constraints is constructed through piecewise linearization of the high-speed train’s nonlinear dynamics. Secondly, the piecewise linear system is transformed into a mixed logical dynamical (MLD) system by introducing the auxiliary binary variables. For the transformed MLD system, a hybrid MPC controller is designed to realize the precise control under hard constraints. To reduce the online computation complexity, the explicit control law is computed offline by employing the mixed-integer linear programming (MILP) technique. Simulation results validate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

With rapid development of high-speed railway, the operation safety, punctuality, and energy consumption of train have received more and more attentions. Nowadays, in China, the high-speed railway generally adopts distributed traction power to improve speed and traction efficiency, where locomotives and wagons are marshalled together to form a high-speed train [1]. For distributed traction power trains, there are some critical issues that need to be solved, such as vehicle traction and brake force allocation, operation safety, punctuality, and energy consumption [2].To achieve the safety and multiobjective optimization, automatic train control system is developed for high-speed railway, and it is one of the key on-board equipment with high safety assurance [3]. The automatic train operation (ATO) system is the essential component that plays a key role in the train operation process [4]. The ATO system controls the train traction/brake force to follow the reference speed. Moreover, ATO system can reduce the energy consumption and improve riding comfort. Now, it is still a challenge to develop an efficient modeling and control method for ATO system [5].Currently, the researches of automatic train control are mainly based on single-point model and multipoint model [6]. For single-point model, the whole high-speed train is simplified as a single rigid mass point, which ignores the coupling characteristics among the vehicles and the allocation of the distributed traction/brake force. The dynamic properties of single-point model are simple and easy to design. However, high-speed train adopts distributed traction power mode; thus the single-point model is obviously not able to describe the coupling characteristics and hydrodynamic dispersion characteristics of the train. Thus multipoint model is favoured to study the automatic train control strategy.Yang and Sun [7, 8] analyzed the coupling characteristics among vehicles of high-speed train and established train multipoint dynamical model. For both push-pull driving and distributed driving train types, the hybrid H2/H∞ automatic train controller was designed and was compared by H2 controller and H∞ controller, respectively. Lin et al. [9] proposed a control strategy based on improved sliding mode control to overcome the effect of the disturbance in the process of train operation, which had good robustness and effectively restrained the high frequency chattering phenomenon. Thus the work reduced the impact of the frequent switching of control input and increased riding comfort. To simplify the controller design, Song et al. [10] utilized the geometric topology to reduce multiple positions to one position for multipoint model. In their work, to address the nonlinear saturation constraints of traction and brake forces, the robust adaptive controller was designed with low computing load, which had robustness for train disturbance and uncertainty of model parameters.Based on the multipoint model, Zhuan and Xia [11] have achieved constructive results on the automatic train control strategy of the heavy-duty train in South Africa. Specifically, Zhuan and Xia [12] put forward three kinds of offline open-loop optimal operation strategies, whose goals are to minimize the workshop bonding force, ensure the safety of train driving, and reduce the maintenance costs of the coupler buffer. Then, Chou and Xia [13] proposed a closed-loop cruising linear quadratic regulator controller to minimize the running cost of electric air brake system in heavy-haul train, in which the control objective includes the speed tracking, coupling force, and energy consumption. In order to overcome the communication limitation, the vehicle barrier was introduced and the controller was reconstructed based on the current orbital slope. In addition, Zhuan and Xia [14] adopted the output feedback adjustment method to regulate the speed of the heavy train and verified the application conditions of their proposed method.As a summary, Scheepmaker et al. [15] presented an extensive review on energy-efficient train control and the related topic of energy-efficient train timetabling, from the first simplest model of a train running on a level track to the advanced models and algorithms of the last decade dealing with varying gradients and speed limits and including regenerative braking. And there have appeared various theories and technologies to realize the efficient and safety train operation, such as the passivity-based cruise control [16], the robust adaptive control [17], and the iterative learning and fault detection approaches [18, 19].Li et al. [20] considered the optimal guaranteed cost of cruise control for high-speed train movement. The sufficient condition for the existence of guaranteed cost cruise control law is given in terms of linear matrix inequalities. And a convex optimization problem is formulated to determine the optimal guaranteed cost that minimizes the performance upper bound of the cruise control law. Ye and Liu [21] propose a novel approach to solve the complex optimal train control problems in closed loop by introducing some simplifications. The operation sequence consists of maximum traction, speed holding, coasting, and maximum braking on each subsection of the track or a constant force is applied on each subsection. Zhao et al. [22] proposed a new cluster consensus technique to design the distributed control law, by which the trains can track the desired speeds asymptotically and the distance between the neighboring cars can be kept in the ideal range.However, with the improvement of the requirement for ATO control performance, the above methods lack the ability to address the nonlinear, multiconstraint, and multiobjective problem based on the multipoint model, which leads to their limitations in practical applications. Especially for online operation process of high-speed trains, it raises the nonlinear resistance force which makes many classical control methods based on the linear resistance force model or the single equilibrium point linearizable model difficult to be implemented. And it is necessary to consider the safe speed, the saturation characteristics of traction and braking units, and the work conditions such as maximization of bonding forces, operating punctuality, energy-efficiency, and reduction of coupling wear forces. The model predictive control has the advantages of fully considering the input and state constraints of the system and the multiobjective optimization problems in the design of the controller, which is more suitable for the design of ATO controller [23, 24].In this paper, an optimal automatic train control strategy is proposed for high-speed trains based on model predictive control. Due to the high computation complexity of nonlinear train model, the common method is to linearize the nonlinear system model at the equilibrium points [25]. However, if the nonlinear system has a wide range of operating condition, there are many deviations on entire running process for singular linearized model with one equilibrium point. Thus, the piecewise linear systems are utilized to approximate the original system, which can describe complex nonlinear systems accurately [26]. In this paper, the nonlinear running resistance is fitted by multiple line segments to construct the piecewise linear system model of train.In the traditional controller design based on piecewise linear model, the linear submodel will be determined according to the initial state of receding horizon; then the predicted control is solved based on the determined sublinear model, which is fixed in the entire receding horizon [27]. However, for each actual control step, the system state changes and the corresponding submodel switch may also occur. The fixed submodel is difficult to approximate the original nonlinear model, especially near the model switching point. Therefore, the traditional model prediction control just according to current linear submodel can result in a large prediction error.In order to solve the above problems, the logic variables are introduced to describe the switch of different linear submodels in this paper. By using this method, the piecewise linear system can be transformed into a hybrid logic dynamical system [28, 29]. For the hybrid logic dynamical system, the model switch can be described by changing the logical variables. Then the hybrid model predictive control (MPC) can be utilized to generate the control inputs based on the constructed mixed logic dynamical system.The hybrid MPC controller uses the mixed logic dynamical model to predict the future evolution of the system at each time step, and a certain performance index is optimized under operating constraints with respect to a sequence of future input moves. The first of such optimal solution is the control action applied to the plant. For the hybrid MPC controller, the optimal problem that needs to be solved is a mixed integer linear program (MILP) due to the hybrid system description. Because the mixed logic dynamical model ensures that the submodel can switch on the entire receding horizon. Thus, more accurate approximation can be obtained by hybrid MPC and the prediction error is reduced.Compared to the existing work, the proposed hybrid MPC makes a superior trade-off between the model complexity and control performance. The proposed controller can obtain the optimal input sequence without excessive simplifying assumption, while the piecewise linear model reduces the computation complexity under the sufficient precision compared to the classical MPC method. The main contributions in this paper can be summarized as follows.(1) The nonlinear multiple mass point dynamical model of high-speed trains is constructed. In the model, the complex nonlinear characteristics of running resistance are piecewise linearized and the piecewise linear description is transformed to mixed logic dynamical system.(2) Based on the MLD system model, a hybrid MPC framework is introduced and the MLD model is used as the prediction model. Based on the MLD prediction model, the control problem can be formulated as a unified mixed integer linear programming (MILP) problem, even the piecewise linear characteristics of train model.(3) For the formulated MILP problem containing logical variables and continuous variables, the branch and bound method (B&B) is utilized to transform the original MILP problem into linear programming (LP) by relaxing the logic variable constraints.The remainder of the paper is organized as follows. In Section2, we establish the piecewise linear model of the high-speed train and transform the piecewise linear model to MLD model. In Section 3, a mixed logic dynamical system is developed. In Section 4, a hybrid MPC controller is designed by analyzing the performance indices and constraints of high-speed trains. Then, the simulation results are provided to validate the effectiveness of the controller in Section 5. Finally, the major conclusion of this paper is given in Section 6.


## 2. The Model of High-Speed Train

In this section, we establish the nonlinear multiple mass-point dynamic model of high-speed trains by analyzing their dynamical characteristics. Because the nonlinear model is complex for designing the system controller, the nonlinear curve of the running resistance is piecewise linearized and the piecewise linear model of the train is constructed.


### 2.1. The Multiple Mass-Point Model of Trains

Figure1 presents the multiple mass-point model structure for the simplified electric multiple power units. The total number of locomotives and wagons is denoted by n. These vehicles are connected by the vehicle hook couplers. The mechanical properties of the couplers are usually described as the “elastic-damping” system. Thus, the hook coupler force between vehicle i and vehicle i+1 is as follows:(1)fiin=kiΔLi+diΔL˙i,where ΔLi=xi-xi+1-hi+1-hi-L0 is the offset value of the relative equilibrium position of the hook coupler between vehicle i and vehicle i+1, L0 is the natural length of the coupler and hi is the half length of vehicle i, and xi is the displacement of vehicle i. Since the length of each vehicle is essentially the same and constant (denoted by h) and the hook coupler has the same type, the balance position between two vehicles can be denoted by Le=2h+L0. Thus ΔLi=xi-xi+1-Le. And the hook coupler force (1) can be transferred to (2)fiin=kxi-xi+1-Le+dx˙i-x˙i+1.Figure 1
The multiple mass-point model structure of the electric multiple units.The multiple mass-point dynamic equation of electric multiple power units can be described as(3)m1v˙1=u1-f1in-f1r,miv˙i=ui+fi-1in-fiin-fir,i=2,3,…,n-1mnv˙n=un+fn-1in-fnr,x˙i=vi,i=1,2,…,n,where ui is the control input of vehicle i. When the train is in traction state or cruise state, ui>0 is traction force. When the train is in brake state, ui<0 is the brake force. When the train is in coasting state, ui=0, implying that there is no traction force and brake force. For locomotive vehicles, they can provide both traction force and brake force, while for wagon vehicles, they can just provide brake force. Assume that fir denotes the running resistance of vehicle i; that is, (4)fir=fib+fia=mic0+micvvi+micavi2+mig·sin⁡θi+0.004·miDi,where fib and fia are the basic resistance and appended resistance, respectively, and mi and vi are the mass and velocity of vehicle i, θi is the ramp angle, and Di is the curvature of rail.


### 2.2. The Multiple Mass-Point Piecewise Linearized Model of Trains

According to the dynamic analysis of high-speed trains, the basic resistance of the train is composed of mechanical resistance and air resistance. When the train is running at low speed, for instance, below 44.44 m/s (160 km/h), the basic resistance is mainly mechanical resistance, and air resistance can be ignored. However, when the train runs at high speed, air resistance will be dominant, and it is related to the square item of the train speed with nonlinear characteristics. As the speed of high-speed trains becomes higher, this nonlinear relationship is becoming stronger. In order to facilitate the controller design of nonlinear model, we consider linearization of the nonlinear model of the train.For nonlinear systems, the linear model is usually obtained by using the equilibrium linearization method. However, for high-speed trains with varying speeds, it is not appropriate to choose one equilibrium point. Piecewise linear systems approximate the original system by using multiple linear subsystems, which can approximate most complex nonlinear systems accurately. Thus, the least square method is adopted to fit the train running resistance linearly. By using the method, the piecewise linear system model of the train is established.The train running parameters are selected from Chinese electric multiple units high-speed train CRH3, which are shown in Table1. Basic resistance after the piecewise linearization is presented in Figure 2, where the train speed range is from 0 m/s to 100 m/s (0–360 km/h). Through lots of experiments, we set va=35m/s and vb=70m/s as the boundaries of piecewise linearization. The original nonlinear equations are divided into three linear equations, and the piecewise linear parameters of basic resistance are shown in Table 2.Table 1
The basic resistance parameters of Chinese electric multiple units high-speed train CRH3.
Parameter Value Unit mi 47.5 t c0 7.75×10-3 Nkg-1 cv 2.28×10-4 Ns(mkg)-1 ca 1.66×10-5 Ns2(m2kg)-1Table 2
The piecewise linear parameters of basic resistance.
j α j  ( s - 1 ) β j  ( N k g - 1 ) Speed (m⋅s−1) 1 0.0008 0.0044 [ 0,35 ] 2 0.0020 −0.0363 [ 35,70 ] 3 0.0031 −0.1100 [ 70,100 ]Figure 2
The piecewise linearization of train basic resistance.It can be seen from Figure2 that the running resistance of the train unit after linearization is already a good approximation to the original nonlinear curve. In fact, using more sublinear model is better for the approximation effect of the original nonlinear model, but it will increase the calculation of the controller. In order to further verify the approximate effect of the linearized train running resistance, the error of the approximate linear curve and the original curve is drawn, which is shown in Figure 3. From Figure 3, the error between the piecewise linearization of running basic resistance and the practical basic resistance is less than 3.3 N/kg, which is just 0.55% compared to the biggest traction force of single vehicle, 600 N/kg.Figure 3
The error between the piecewise linearization of train basic resistance and the original true basic resistance.Thus, the linearization of train basic running resistance can be approximated as(5)ω0=αjv+βj,j=1,2,3.In fact, due to the force constraints of hook couplers, there is small difference among the operative velocities of vehicles. Thus, the velocity of vehicle 1 can be taken as the piecewise point of piecewise linear model. Taking it into the train multiple mass-points model, the piecewise linear dynamic model of train can be described as(6)x˙t=A1cxt+Bcut+f1cvi≤vax˙t=A2cxt+Bcut+f2cva<vi≤vbx˙t=A3cxt+Bcut+f3cvi>vb,where the state variable x=[v1,v2,…,vn,x1,x2,…,xn]T, control input u=[u1,u2,…,un]T, system matrix Ajc∈R2n×2n, Bc∈R2n×n, and fjc∈R2n×1 are as follows:(7)
A
j
c
=
A
11
j
A
12
j
I
n
×
n
0
n
×
n
,
B
c
=
B
11
0
n
×
n
,
f
j
c
=
f
j
1
,
f
j
2
,
…
,
f
j
n
,
0,0
,
…
,
0
T
,
A
11
j
=
-
d
m
1
-
α
j
d
m
1
0
⋯
0
0
0
d
m
2
-
2
d
m
2
-
α
j
d
m
2
⋯
0
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
d
m
n
-
1
-
2
d
m
n
-
1
-
α
j
d
m
n
-
1
0
0
0
⋯
0
d
m
n
-
d
m
n
-
α
j
,
A
12
j
=
-
k
m
1
k
m
1
0
⋯
0
0
0
k
m
2
-
2
k
m
2
k
m
2
⋯
0
0
0
⋮
⋮
⋱
⋮
⋮
⋮
0
0
0
⋯
k
m
n
-
1
-
2
k
m
n
-
1
k
m
n
-
1
0
0
0
⋯
0
k
m
n
-
k
m
n
,and B11=diag1/m1,1/m2,…,1/mn, and fji=k1Le/mi-βj-9.98sin⁡θi-0.004Di. In practical controller design, Le can be eliminated by state transformation, so the equilibrium position is set to zero for each subsystem.Then, the zero-order holder is used to discretize the continuous time model. By setting the sampling period asTs, the discrete-time piecewise linear model can be obtained as(8)xk+1=Ajdxk+Bjduk+fjd,j=1,2,3, where Ajd=eAjcTs, Bjd=∫0TseAjcτdτBc, and fjd=∫0TseAjcτdτfjc.


## 2.1. The Multiple Mass-Point Model of Trains

Figure1 presents the multiple mass-point model structure for the simplified electric multiple power units. The total number of locomotives and wagons is denoted by n. These vehicles are connected by the vehicle hook couplers. The mechanical properties of the couplers are usually described as the “elastic-damping” system. Thus, the hook coupler force between vehicle i and vehicle i+1 is as follows:(1)fiin=kiΔLi+diΔL˙i,where ΔLi=xi-xi+1-hi+1-hi-L0 is the offset value of the relative equilibrium position of the hook coupler between vehicle i and vehicle i+1, L0 is the natural length of the coupler and hi is the half length of vehicle i, and xi is the displacement of vehicle i. Since the length of each vehicle is essentially the same and constant (denoted by h) and the hook coupler has the same type, the balance position between two vehicles can be denoted by Le=2h+L0. Thus ΔLi=xi-xi+1-Le. And the hook coupler force (1) can be transferred to (2)fiin=kxi-xi+1-Le+dx˙i-x˙i+1.Figure 1
The multiple mass-point model structure of the electric multiple units.The multiple mass-point dynamic equation of electric multiple power units can be described as(3)m1v˙1=u1-f1in-f1r,miv˙i=ui+fi-1in-fiin-fir,i=2,3,…,n-1mnv˙n=un+fn-1in-fnr,x˙i=vi,i=1,2,…,n,where ui is the control input of vehicle i. When the train is in traction state or cruise state, ui>0 is traction force. When the train is in brake state, ui<0 is the brake force. When the train is in coasting state, ui=0, implying that there is no traction force and brake force. For locomotive vehicles, they can provide both traction force and brake force, while for wagon vehicles, they can just provide brake force. Assume that fir denotes the running resistance of vehicle i; that is, (4)fir=fib+fia=mic0+micvvi+micavi2+mig·sin⁡θi+0.004·miDi,where fib and fia are the basic resistance and appended resistance, respectively, and mi and vi are the mass and velocity of vehicle i, θi is the ramp angle, and Di is the curvature of rail.


## 2.2. The Multiple Mass-Point Piecewise Linearized Model of Trains

According to the dynamic analysis of high-speed trains, the basic resistance of the train is composed of mechanical resistance and air resistance. When the train is running at low speed, for instance, below 44.44 m/s (160 km/h), the basic resistance is mainly mechanical resistance, and air resistance can be ignored. However, when the train runs at high speed, air resistance will be dominant, and it is related to the square item of the train speed with nonlinear characteristics. As the speed of high-speed trains becomes higher, this nonlinear relationship is becoming stronger. In order to facilitate the controller design of nonlinear model, we consider linearization of the nonlinear model of the train.For nonlinear systems, the linear model is usually obtained by using the equilibrium linearization method. However, for high-speed trains with varying speeds, it is not appropriate to choose one equilibrium point. Piecewise linear systems approximate the original system by using multiple linear subsystems, which can approximate most complex nonlinear systems accurately. Thus, the least square method is adopted to fit the train running resistance linearly. By using the method, the piecewise linear system model of the train is established.The train running parameters are selected from Chinese electric multiple units high-speed train CRH3, which are shown in Table1. Basic resistance after the piecewise linearization is presented in Figure 2, where the train speed range is from 0 m/s to 100 m/s (0–360 km/h). Through lots of experiments, we set va=35m/s and vb=70m/s as the boundaries of piecewise linearization. The original nonlinear equations are divided into three linear equations, and the piecewise linear parameters of basic resistance are shown in Table 2.Table 1
The basic resistance parameters of Chinese electric multiple units high-speed train CRH3.
Parameter Value Unit mi 47.5 t c0 7.75×10-3 Nkg-1 cv 2.28×10-4 Ns(mkg)-1 ca 1.66×10-5 Ns2(m2kg)-1Table 2
The piecewise linear parameters of basic resistance.
j α j  ( s - 1 ) β j  ( N k g - 1 ) Speed (m⋅s−1) 1 0.0008 0.0044 [ 0,35 ] 2 0.0020 −0.0363 [ 35,70 ] 3 0.0031 −0.1100 [ 70,100 ]Figure 2
The piecewise linearization of train basic resistance.It can be seen from Figure2 that the running resistance of the train unit after linearization is already a good approximation to the original nonlinear curve. In fact, using more sublinear model is better for the approximation effect of the original nonlinear model, but it will increase the calculation of the controller. In order to further verify the approximate effect of the linearized train running resistance, the error of the approximate linear curve and the original curve is drawn, which is shown in Figure 3. From Figure 3, the error between the piecewise linearization of running basic resistance and the practical basic resistance is less than 3.3 N/kg, which is just 0.55% compared to the biggest traction force of single vehicle, 600 N/kg.Figure 3
The error between the piecewise linearization of train basic resistance and the original true basic resistance.Thus, the linearization of train basic running resistance can be approximated as(5)ω0=αjv+βj,j=1,2,3.In fact, due to the force constraints of hook couplers, there is small difference among the operative velocities of vehicles. Thus, the velocity of vehicle 1 can be taken as the piecewise point of piecewise linear model. Taking it into the train multiple mass-points model, the piecewise linear dynamic model of train can be described as(6)x˙t=A1cxt+Bcut+f1cvi≤vax˙t=A2cxt+Bcut+f2cva<vi≤vbx˙t=A3cxt+Bcut+f3cvi>vb,where the state variable x=[v1,v2,…,vn,x1,x2,…,xn]T, control input u=[u1,u2,…,un]T, system matrix Ajc∈R2n×2n, Bc∈R2n×n, and fjc∈R2n×1 are as follows:(7)
A
j
c
=
A
11
j
A
12
j
I
n
×
n
0
n
×
n
,
B
c
=
B
11
0
n
×
n
,
f
j
c
=
f
j
1
,
f
j
2
,
…
,
f
j
n
,
0,0
,
…
,
0
T
,
A
11
j
=
-
d
m
1
-
α
j
d
m
1
0
⋯
0
0
0
d
m
2
-
2
d
m
2
-
α
j
d
m
2
⋯
0
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
d
m
n
-
1
-
2
d
m
n
-
1
-
α
j
d
m
n
-
1
0
0
0
⋯
0
d
m
n
-
d
m
n
-
α
j
,
A
12
j
=
-
k
m
1
k
m
1
0
⋯
0
0
0
k
m
2
-
2
k
m
2
k
m
2
⋯
0
0
0
⋮
⋮
⋱
⋮
⋮
⋮
0
0
0
⋯
k
m
n
-
1
-
2
k
m
n
-
1
k
m
n
-
1
0
0
0
⋯
0
k
m
n
-
k
m
n
,and B11=diag1/m1,1/m2,…,1/mn, and fji=k1Le/mi-βj-9.98sin⁡θi-0.004Di. In practical controller design, Le can be eliminated by state transformation, so the equilibrium position is set to zero for each subsystem.Then, the zero-order holder is used to discretize the continuous time model. By setting the sampling period asTs, the discrete-time piecewise linear model can be obtained as(8)xk+1=Ajdxk+Bjduk+fjd,j=1,2,3, where Ajd=eAjcTs, Bjd=∫0TseAjcτdτBc, and fjd=∫0TseAjcτdτfjc.


## 3. Hybrid MPC Based ATO Design

A discrete-time piecewise linear model of high-speed trains has been established, which consists of multiple piecewise linear submodels. For the traditional predictive controller design based on the piecewise linear model, the linear submodel based on the current state is determined at the beginning of the prediction and then the predictive control is operated according to the determined linear submodel. However, the predictive submodel has been fixed and does not change throughout the prediction horizon, which cannot reflect the system nonlinear characteristics by only using a fixed submodel. Actually, the predictive system state may change, and the corresponding submodel may also switch during the solution process. Using a fixed submodel is not enough to express the nonlinearity of the original model, especially near the switching point between submodels. Therefore, the traditional predictive control based on piecewise linear model may produce large prediction errors.In order to address the above issues, a mixed logic dynamical system is introduced in this paper, which is proposed by Bemporad and Morari in the context of MPC technology [28]. By introducing logical variables to represent different linear submodels, the piecewise linear system can be transformed into a mixed logic dynamical system, which is a global dynamical system covering the whole nonlinear dynamics. Based on the model predictive control of the mixed logic dynamical system, the switching of the model can be expressed by the change of the logic quantity in the model. The model can be switchable throughout the predicted horizon, which ensures more accurate approximation of the original nonlinear model and reduces the prediction error. According to the propositional transformation of the mixed logic dynamic system and mixed integer linear inequality conversion rules, we transform the train piecewise linear model into a mixed integer dynamic model by introducing the binary variables [28].Two binary variablesδ1k and δ2k are given by(9)vi≤va⟺δ1k=1vi≤vb⟺δ2k=1.Based onvi∈vmin,vmax, (9) can be written as follows: (10)δ1kva-vmin≥vik-vminδ1kva-vmax≤va-vik-εδ2kvb-vmin≥vik-vminδ2kvb-vmax≤vb-vik-ε,where ε is a small positive number and generally takes 10-6. It is used to transform a strict equality into a nonstrict inequality.By definingΦi(k)=Aidx(k)+Bidu(k)+fid, the relationship between δ1(k), δ2(k), and x(k+1) can be shown in Table 3.Table 3
The relationship betweenδ1(k), δ2(k), and x(k+1).
δ 1 ( k ) δ 2 ( k ) x ( k + 1 ) 1 1 Φ1(k) 0 1 Φ2(k) 0 0 Φ3(k)Then, the piecewise linear system (8) can be described as(11)xk+1=δ1kδ2kΦ1k+1-δ1kδ2kΦ2k+1-δ1k1-δ2kΦ3k.We introduce a new variableδ3(k)=δ1(k)δ2(k), which satisfies(12)-δ1k+δ3k≤0-δ2k+δ3k≤0δ1k+δ2k-δ3k≤1.In addition, the auxiliary mixed logical vectors and variables are introduced, that is,z1x(k)=δ1(k)x(k), z2x(k)=δ2(k)x(k), z3x(k)=δ3(k)x(k), z1u(k)=δ1(k)u(k), z2u(k)=δ2(k)u(k), and z3u(k)=δ3(k)u(k), which can be expressed as(13)xminδik≤zixk≤xmaxδikxk-xmax1-δik≤zixk≤xk-xmin1-δikuminδik≤ziuk≤umaxδikuk-umax1-δik≤ziuk≤uk-umin1-δik,for i=1,2,3, where x∈xmin,xmax and u∈umin,umax. The maximum tract input umax is a nonincreasing function of speed, which is approximated by a piecewise linear, quadratic, and/or hyperbolic function of speed [15]. Thus it can be described as a group of hyperbolic or parabolic formulas and each formula approximates the actual traction force for a certain speed interval [30], such as(14)umaxv=c1,j+c2,jv+c3,jv2,v∈vj,vj+1 or(15)umaxv=ch,jv,v∈vj,vj+1,forj=1,2,…,M-1,where the piecewise point vj and vj+1 and the parameters c1j, c2j, c3j, and chj can be determined by train operation experiments.Then (8) can be transformed into the following mixed logical dynamical (MLD) system:(16)xk+1=Axk+B1uk+B2δk+B3zk+B4E2δk+E3zk≤E1uk+E4xk+E5,where zk=z1x(k)Tz2x(k)Tz3x(k)Tz1u(k)z2u(k)z3u(k)T, δ(k)=δ1(k)δ2(k)δ3(k)T, A=A3d, B1=B3d, B4=f3d, B2=[-f3d,f2d-f3d,f3d-f2d+f1d], and B3=[-A3d,A2d-A3d,A1d-A2d+A3d,-B3d,B2d-B3d,B1d-B2d+B3d].Specifically, the resulting MLD model (14) is obtained by using the HYSDEL compiler [31]. Matrices E1,…,E5 include constraints (16). They can be obtained and analyzed in Matlab using, for example, the Hybrid Toolbox [32].


## 4. Controller Design

Based on the MLD system model, a hybrid MPC framework is introduced in this section. Within this, the MLD model is used as the prediction model. Then the control problem is formulated as a constrained finite time optimal control problem and converted into mixed integer linear programming (MILP) problem, for which efficient solvers are available.


### 4.1. Formulation of Optimal Control Problem

Automatic train control needs to track the desired velocity profile quickly such that the timeliness of the interstation travel can be ensured. With the development of high-speed railway, energy-saving driving has become an important goal of automatic train control. As the running distance of the train between stations is fixed, the energy consumption directly related to the traction and braking force. Thus, the control input is used to express energy consumption through this paper. In addition, the in-train force should be kept in the desired range to ensure the safety of the train.Therefore, the optimal problem can be established as follows, which considers the speed tracking error, control input, and the in-train force.(17)JuNk=∑i=0N-1εk+i∣kQ1+uk+i∣kR1+xink+i∣kH1,where k is the current sampling time; uN(k)=[u(k),…,u(k+N-1)] is the sequence of N control inputs. ε((k+i)∣k)=v((k+i)∣k)-vr(k+i)∈Rn×1 represents the difference between actual speed v and reference speed vr, u((k+i)∣k)=u1(k+i),u2(k+i),…,un(k+i)T represents the control input, and xin((k+i)∣k)=(xin1((k+i)∣k,…,xinn-1(((k+i)∣k)) represents the relative displacements between neighboring cars. 1 means 1-norm and Q, R, and H are the weight matrixes, respectively.In order to ensure the safe and efficient operation of the train, some constraints must be satisfied as follows.(1) Constraint on Control Input. The tractive and braking forces are bounded due to the nature physical characteristics of the traction motor. Thus, the control input is set as -FmaxB≤uik≤FmaxT.(2) Constraint on Maximum Operational Velocity. The maximum speed of the train depends first on the maximum speed allowed by the train itself. Second, when the line conditions and operating conditions change, the maximum allowable speed of the train may change.(3) Constraint on In-Train Force. The in-train force should be in the range -fmaxin,fmaxin to keep the operational safety of the trains.(4) Constraint on Coupler Deformation. The coupler deformation should be in the range -ΔLmax,ΔLmax to keep the operational safety of the trains.Then, the constraints considered can be formulated by(18)-FmaxB≤uik≤FmaxT0≤vik≤vmax-fmaxin≤fi,i+1in≤fmaxin-ΔLmax≤ΔLi≤ΔLmax.


### 4.2. Hybrid MPC Controller Design

Based on (17) and (18), a constrained optimal control problem can be constructed as follows:(19)minJuNk=∑i=0N-1εk+i∣kQ1+uk+i∣kR1+xink+i∣kH1s.t.xk+1=Axk+B1uk+B2δk+B3zk+B4E2δk+E3zk≤E1uk+E4xk+E5-FmaxB≤uik≤FmaxT0≤vik≤vmax-fmaxin≤fi,i+1in≤fmaxin-ΔLmax≤ΔLi≤ΔLmax.In order to facilitate the solution of the optimal problem, (19) can be transformed to the linear objective function (20) by introducing the auxiliary variables wε(k+i∣k)=ε((k+i)∣k)Q1, wu(k+i∣k)=u((k+i)∣k)R1, and wxin(k+i∣k)=xin((k+i)∣k)H1.(20)JuNk=wεk+N∣k+∑i=0N-1wεk+i∣k+wuk+i∣k+wxink+i∣kwith some linear inequalities:(21)wεk+i∣k≥±εk+i∣kwuk+i∣k≥±uk+i∣kwxink+i∣k≥±xink+i∣k.In addition, the state update equation (22) can be obtained by the MLD model. (22)xk+i=Aixk+∑j=0i-1AjB1uk+i-j-1+B2δk+i-j-1+B3zk+i-j-1E2δk+E3zk≤E1uk+E4xk+E5.Then, based on (20), (21), and (22), the optimal control problem can be converted to a mixed integer linear programming (MILP) problem as follows:(23a)JV,θ=minVCTV(23b)st.GV≤W+Sθk,where θ(k)=[x(k)T,vr(k),…,vr(k+Np-1)]T is commonly represented as the parameter vector and V=u′Np(k),δ~(k)T,z~(k)T,w~(k)TT is the optimization vector. Let (24)δ~k=δkδk+1⋮δk+Np-1,z~k=zkzk+1⋮zk+Np-1,w~k=wεk,wu′k,…,wεk+Np-1,wu′k+Np-1T;CT=0,…,0,1,…,1. The matrices G, W, and S are obtained according to (16), (18), and (21).At each sampling stepk, the MILP problem is solved based on the current parameter vector θ(k), and an optimal sequence of future control input uNp′(k) can be obtained. At time k+1, the parameter vector is updated to θ(k+1), and a new MILP problem is solved to get the new optimal control input. This procedure is repeated at each sampling step, and the problems can be solved in real time. In this paper, only Nc=Np is considered.


### 4.3. Solution to MILP Problem

The problem of automatic train control is established as a constrained finite time optimal control problem and further transformed into mixed integer linear programming (MILP) problem by analyzing the operation target and constraints. Since the constructed MILP problem contains logical variables and continuous variables, this section uses the branch and bound method (B&B) to solve this problem so that, in each control cycle, the optimal traction of each vehicle or braking force output can be obtained.The basic idea of branch and bound method is to relax the logic variable constraints in the optimization problem.δi(k)∈0,1 in the MLD model of formula (16) can be relaxed as δi(k)∈0,1, which transforms the original MILP problem into linear programming (LP). Then, the traditional linear programming method is used to solve the optimal subproblem. Algorithm 1 is presented to solve the MILP problem and some notations are defined and presented in Table 4.Table 4
Definition of some notations.
Notation Definition P The MILP problem J0∗ The minimum upper bound of optimal performance of train operation V∗ The optimal solution vector of optimal performance index for train operation Π The object table of recording all LP subproblems that need to be solved J(V-) The optimal performance of the current relaxed LP problem V- The optimal solution vector of the relaxed LP subproblem Qnum Number of LP subproblems that have been solved M a x n u m  L P ( j ) [ P 0 ] , L P ( j ) [ P 1 ] The maximum number of LP problemsAlgorithm 1:Solution of MILP problem.
(1) initialize J0∗=+∞, Qnum=0, Π= Null; (2) for each i∈[1,n]do (3)  δi(k)∈0,1←δi(k)∈0,1; (4)  obtain LP-subproblem LPi; (5)  Π←Π∪LPi (6) end for (7) if there are no LPi∈Π that need to be solved then (8)  obtain optimal performance index J0∗; (9)  obtain optimal solution vector V∗; (10)  obtain optimal control input uk←V∗; (11) else (12)  whileQnum<Maxnumdo (13)   for each j∈[n,1]do (14)    extracting LPj from Π; (15)    ifLPj has feasible solution then (16)     obtain J(V-) and V-; (17)    else (18)     ifLPj satisfies δi(k)∈0,1 and J(V-)< J0∗then (19)      J0∗←J(V-) and V∗←V-; (20)      Qnum=Qnum+1; (21)     end if (22)     ifJ(V-)>J0∗then (23)      Qnum=Qnum+1; (24)     end if (25)     ifLPj dose not satisfy δi(k)∈0,1 and J(V-)<J0∗ then (26)      choose δi(k); (27)      LPjP0←δi(k)=1; (28)      LPjP1←δi(k)=0; (29)      Π←ΠLPj∪LPjP0∪LPjP1; (30)     else (31)      Qnum=Qnum+1; (32)     end if (33)     Qnum=Qnum+1; (34)    end if (35)   end for (36)  end while (37) end if


## 4.1. Formulation of Optimal Control Problem

Automatic train control needs to track the desired velocity profile quickly such that the timeliness of the interstation travel can be ensured. With the development of high-speed railway, energy-saving driving has become an important goal of automatic train control. As the running distance of the train between stations is fixed, the energy consumption directly related to the traction and braking force. Thus, the control input is used to express energy consumption through this paper. In addition, the in-train force should be kept in the desired range to ensure the safety of the train.Therefore, the optimal problem can be established as follows, which considers the speed tracking error, control input, and the in-train force.(17)JuNk=∑i=0N-1εk+i∣kQ1+uk+i∣kR1+xink+i∣kH1,where k is the current sampling time; uN(k)=[u(k),…,u(k+N-1)] is the sequence of N control inputs. ε((k+i)∣k)=v((k+i)∣k)-vr(k+i)∈Rn×1 represents the difference between actual speed v and reference speed vr, u((k+i)∣k)=u1(k+i),u2(k+i),…,un(k+i)T represents the control input, and xin((k+i)∣k)=(xin1((k+i)∣k,…,xinn-1(((k+i)∣k)) represents the relative displacements between neighboring cars. 1 means 1-norm and Q, R, and H are the weight matrixes, respectively.In order to ensure the safe and efficient operation of the train, some constraints must be satisfied as follows.(1) Constraint on Control Input. The tractive and braking forces are bounded due to the nature physical characteristics of the traction motor. Thus, the control input is set as -FmaxB≤uik≤FmaxT.(2) Constraint on Maximum Operational Velocity. The maximum speed of the train depends first on the maximum speed allowed by the train itself. Second, when the line conditions and operating conditions change, the maximum allowable speed of the train may change.(3) Constraint on In-Train Force. The in-train force should be in the range -fmaxin,fmaxin to keep the operational safety of the trains.(4) Constraint on Coupler Deformation. The coupler deformation should be in the range -ΔLmax,ΔLmax to keep the operational safety of the trains.Then, the constraints considered can be formulated by(18)-FmaxB≤uik≤FmaxT0≤vik≤vmax-fmaxin≤fi,i+1in≤fmaxin-ΔLmax≤ΔLi≤ΔLmax.


## 4.2. Hybrid MPC Controller Design

Based on (17) and (18), a constrained optimal control problem can be constructed as follows:(19)minJuNk=∑i=0N-1εk+i∣kQ1+uk+i∣kR1+xink+i∣kH1s.t.xk+1=Axk+B1uk+B2δk+B3zk+B4E2δk+E3zk≤E1uk+E4xk+E5-FmaxB≤uik≤FmaxT0≤vik≤vmax-fmaxin≤fi,i+1in≤fmaxin-ΔLmax≤ΔLi≤ΔLmax.In order to facilitate the solution of the optimal problem, (19) can be transformed to the linear objective function (20) by introducing the auxiliary variables wε(k+i∣k)=ε((k+i)∣k)Q1, wu(k+i∣k)=u((k+i)∣k)R1, and wxin(k+i∣k)=xin((k+i)∣k)H1.(20)JuNk=wεk+N∣k+∑i=0N-1wεk+i∣k+wuk+i∣k+wxink+i∣kwith some linear inequalities:(21)wεk+i∣k≥±εk+i∣kwuk+i∣k≥±uk+i∣kwxink+i∣k≥±xink+i∣k.In addition, the state update equation (22) can be obtained by the MLD model. (22)xk+i=Aixk+∑j=0i-1AjB1uk+i-j-1+B2δk+i-j-1+B3zk+i-j-1E2δk+E3zk≤E1uk+E4xk+E5.Then, based on (20), (21), and (22), the optimal control problem can be converted to a mixed integer linear programming (MILP) problem as follows:(23a)JV,θ=minVCTV(23b)st.GV≤W+Sθk,where θ(k)=[x(k)T,vr(k),…,vr(k+Np-1)]T is commonly represented as the parameter vector and V=u′Np(k),δ~(k)T,z~(k)T,w~(k)TT is the optimization vector. Let (24)δ~k=δkδk+1⋮δk+Np-1,z~k=zkzk+1⋮zk+Np-1,w~k=wεk,wu′k,…,wεk+Np-1,wu′k+Np-1T;CT=0,…,0,1,…,1. The matrices G, W, and S are obtained according to (16), (18), and (21).At each sampling stepk, the MILP problem is solved based on the current parameter vector θ(k), and an optimal sequence of future control input uNp′(k) can be obtained. At time k+1, the parameter vector is updated to θ(k+1), and a new MILP problem is solved to get the new optimal control input. This procedure is repeated at each sampling step, and the problems can be solved in real time. In this paper, only Nc=Np is considered.


## 4.3. Solution to MILP Problem

The problem of automatic train control is established as a constrained finite time optimal control problem and further transformed into mixed integer linear programming (MILP) problem by analyzing the operation target and constraints. Since the constructed MILP problem contains logical variables and continuous variables, this section uses the branch and bound method (B&B) to solve this problem so that, in each control cycle, the optimal traction of each vehicle or braking force output can be obtained.The basic idea of branch and bound method is to relax the logic variable constraints in the optimization problem.δi(k)∈0,1 in the MLD model of formula (16) can be relaxed as δi(k)∈0,1, which transforms the original MILP problem into linear programming (LP). Then, the traditional linear programming method is used to solve the optimal subproblem. Algorithm 1 is presented to solve the MILP problem and some notations are defined and presented in Table 4.Table 4
Definition of some notations.
Notation Definition P The MILP problem J0∗ The minimum upper bound of optimal performance of train operation V∗ The optimal solution vector of optimal performance index for train operation Π The object table of recording all LP subproblems that need to be solved J(V-) The optimal performance of the current relaxed LP problem V- The optimal solution vector of the relaxed LP subproblem Qnum Number of LP subproblems that have been solved M a x n u m  L P ( j ) [ P 0 ] , L P ( j ) [ P 1 ] The maximum number of LP problemsAlgorithm 1:Solution of MILP problem.
(1) initialize J0∗=+∞, Qnum=0, Π= Null; (2) for each i∈[1,n]do (3)  δi(k)∈0,1←δi(k)∈0,1; (4)  obtain LP-subproblem LPi; (5)  Π←Π∪LPi (6) end for (7) if there are no LPi∈Π that need to be solved then (8)  obtain optimal performance index J0∗; (9)  obtain optimal solution vector V∗; (10)  obtain optimal control input uk←V∗; (11) else (12)  whileQnum<Maxnumdo (13)   for each j∈[n,1]do (14)    extracting LPj from Π; (15)    ifLPj has feasible solution then (16)     obtain J(V-) and V-; (17)    else (18)     ifLPj satisfies δi(k)∈0,1 and J(V-)< J0∗then (19)      J0∗←J(V-) and V∗←V-; (20)      Qnum=Qnum+1; (21)     end if (22)     ifJ(V-)>J0∗then (23)      Qnum=Qnum+1; (24)     end if (25)     ifLPj dose not satisfy δi(k)∈0,1 and J(V-)<J0∗ then (26)      choose δi(k); (27)      LPjP0←δi(k)=1; (28)      LPjP1←δi(k)=0; (29)      Π←ΠLPj∪LPjP0∪LPjP1; (30)     else (31)      Qnum=Qnum+1; (32)     end if (33)     Qnum=Qnum+1; (34)    end if (35)   end for (36)  end while (37) end if


## 5. Performance Evaluation

In this section, we will evaluate the performance of the proposed hybrid MPC ATO algorithm using numerical results. We will first introduce the parameter settings of the ATO system. Specifically, the parameters of the train model are from the CRH-3 high-speed train in China, which are given in Table5. Therein, FmaxT and FmaxB are the maximum tractive force and braking force, respectively.Table 5
Parameters of the CRH-3 high-speed train.
Parameter Value Unit mi 47.5 t c0 7.75×10-3 Nkg-1 cv 2.28×10-4 Ns(mkg)-1 ca 1.66×10-5 Ns2(m2kg)-1 FmaxT 300 kN FmaxB 300 kN k 1 × 1 0 4 kNm-1 d 1000 kNsm-1 θ 0 − D 0 − Δ L m a x 0.2 m f m a x i n 2000 kN


### 5.1. Simulation Scenario

In order to evaluate the performance of the controller, the desired velocity curve including accelerating, decelerating, velocity step increase, velocity step decrease, and constant-velocity stages is characterized by the following equation:(25)vrt=0.4t0≤t≤10040100≤t≤40075400≤t≤90050900≤t≤1150-t+12001150≤t≤1200.In (25), the acceleration is 0.4m/s, and the deceleration is 0.5m/s, both of which are within the threshold of acceleration and velocity. In the operation of trains, the basic resistance is the main resistance for the trains. Thus, we only consider the basic resistance in this paper.The designed MPC based ATO controller considers the train velocity tracking performance, energy consumption, and intertrain force. With the train parameters, we provide two simulation scenarios to evaluate the performance of the proposed algorithm under different operating conditions. Specifically, we evaluate the performance of the controller with different train types: all-locomotive train and partial-locomotive train.


#### 5.1.1. Scenario A

In this scenario, we evaluate the performance of the ATO controller with the all-locomotive train. To simplify the simulation, we choose the number of trains asn=4, and the train is comprised of all locomotives. This implies that each locomotive has the tracking motor and braking unit and thus can provide both the tracking force and the braking force. The parameters of the train are shown in Table 5, where mi is the mass of locomotive i. Since the numbers of passengers in each locomotive may be different, we choose m1=48.5t,m2=52.5t,m3=51.0t,m4=47.5t. For the MPC controller, the prediction horizon Np and the control horizon Nc also affect the controller performance. Without loss of generality, we choose Np=Nc=3. The weighting matrix gain is chosen as(26)Q=diag4,…,4︸n,R=diag1,…,1︸n,H=diag1,…,1︸n-1.


#### 5.1.2. Scenario B

In this scenario, we evaluate the effectiveness of the ATO controller in partial-locomotive trains. To simplify the simulation, we choose the number of trains asn=4, and the train is comprised of 2 locomotives and 2 wagons. The locomotive can provide both the tracking force and the braking force, while the wagons can only supply the brake forces. In order to guarantee the total force the same as Scenario A, the maximum tracking force of each locomotive is set as FmaxT=600kN, and the other train parameters are the same as in Table 5.


### 5.2. Simulation Results

We will provide simulation results of the ATO controller in different operating conditions. The simulations are conducted in Matlab/SIMULINK. We fist consider the simulations in Scenario A, which is based on all-locomotive train. Then we will further consider the simulation results in Scenario B, which is based on partial-locomotive trains. And in Scenario C, we focus on more realistic condition where the track force limit is dynamic.


#### 5.2.1. Scenario A

The velocity curve of each locomotive is shown in Figure4, where x-axis represents the time andy-axis represents the velocity. From the figure, we can see that the velocities of four locomotives almost remain the same during the operation time. Specifically, in the accelerating and decelerating stages, each locomotive can operate correspondingly according to the given acceleration and deceleration. At t=400 s, we notice that the reference velocity has the step increase. From the zoomed-in figure, we can find that the locomotives move with the almost same velocity, where the maximum velocity error is 0.003m/s. From the zoomed-in figure from t=470 s to t=490 s, the velocities of locomotives reach the steady value at about t=490 s, where the settling time is 90 s and the steady error is 0.5m/s. In the steady state, the velocities of locomotives remain fixed. Based on the simulation results, we can find that the train can track the time-varying reference with small settling time and steady error, which verifies the effectiveness of the algorithm.Figure 4
The velocity curve of each locomotive in Scenario A.The tracking and braking forces of locomotives are shown in Figure5. In the accelerating stage, the resistance of the train increases with the velocity. In order to guarantee the constant acceleration, the tracking forces of locomotives increase gradually. In the deceleration stage, the resistance of the train decelerates with the velocity. In order to guarantee the constant deceleration, the braking forces of locomotives increase gradually. Compared with the acceleration or deceleration stages, the tracking forces of locomotives are small and fixed during the constant velocity stage. This is because the resistance of train keeps having a relatively small value during the constant velocity stage. At t=400 s, the reference has a step increase. In order to accelerate to the desired velocity, the tracking forces of locomotives increase to the maximum rapidly and maintain the maximum tracking force for about 72 s. At t=900 s, the braking forces increase to the maximum to decrease the velocity and maintain about 32 s to decrease to the desired velocity. From the figure, we can find that the control outputs of locomotives are almost the same and thus guarantee the reasonable allocation of braking force and tracking force. When the reference velocity curve changes, each locomotive can regulate the tracking force and braking force based on the difference between the actual velocity and desired velocity. The forces for locomotives are almost the same and also satisfy the constraints, which verify the effectiveness of the proposed strategy.Figure 5
The force output of each locomotive in Scenario A.The coupler force among each locomotive is shown in Figure6. When the coupler force is larger than zero, it means that coupler is pulled. When the coupler force is smaller than zero, the coupler is pushed. From the figure, we can see that the coupler forces among locomotives are always within the safety range. Moreover, the coupler among L1 and L2 is mainly pulled, while the coupler among L3 and L4 is mainly pushed. This is because when we choose the masses of locomotives, the masses of L3 and L4 are relatively smaller than L1 and L2. In order to guarantee the velocity consensus of locomotives, the total forces should be the same and the coupler force can balance the total force of each locomotive. It can be found that, during the constant-velocity stage, constant-accelerating stage, and constant-decelerating stage, the coupler force among locomotives is small. At t=200 s and t=400 s, the coupler force changes dramatically, with the largest coupler pulling force 317 kN and the largest coupler pushing force 241 kN. From the figure, we can see that the ATO controller can guarantee that the coupler force is within the safety range, which verifies the effectiveness.Figure 6
The coupler force among locomotives in Scenario A.The operation safety of high-speed trains can be evaluated by the position difference (coupler deformation) between neighboring locomotives. If the deformation is larger than the maximum0.2m or less than the minimum -0.2m, the coupler may have been broken so the safety cannot be guaranteed. The operation safety of high-speed trains can be further illustrated from Figure 7, which shows that the coupler deformation is always with the threshold [-0.2m,0.2m]. This is because the in-train force guarantees that the position difference between neighboring locomotives cannot be too large or too small.Figure 7
The coupler deformation between two neighbouring locomotives in Scenario A.


#### 5.2.2. Scenario B

The simulation in Scenario A verifies the effectiveness of the hybrid model predictive controller for all-locomotive case. The metrics include velocity tracking, energy consumption, and coupler force.In Scenario B, we verified the effectiveness of the ATO controller with 2 locomotives and 2 wagons. The velocity curves are shown in Figure8. With the same controller parameters, at t=400 s, where the reference velocity has a step change, the convergence time of the controller from the change at t=400 s to the steady state at t=508 s implies that the convergence time increases by 18 s compared with the all-locomotive case. The steady error is 1.2 m/s, which increases by 0.7m/s compared with the all-locomotive case. Moreover, the velocities at the steady stage have some oscillations around the steady value.Figure 8
The velocity curve of each vehicle in Scenario B.The force output of each vehicle is shown in Figure9, where the left subfigure shows the force output of the locomotive, and the right subfigure shows the force output of the wagon. From the figure, we can see that, during the acceleration and constant-velocity stages, only the locomotive has the tracking force output, and the tracking force of wagons is zero. In the deceleration stage, both the locomotives and wagons supply the braking forces.Figure 9
The forces of locomotives and wagons in Scenario B.Figure10 shows the coupler force among vehicles. Compared with the all-locomotive trains, the coupler forces between L1 and W1 and between L2 and W2 are kept at large values. This is because locomotive L1 supplies the pulling force and locomotive L2 supplies the pushing force, while wagons W1 and W2 cannot supply the tracking force. The maximum coupler force is 578 kN, which is within the safety threshold. From the above simulations, we can find that the proposed hybrid MPC ATO controller has a good control effect.Figure 10
The coupler force among vehicles in Scenario B.The operation safety of high-speed trains can be illustrated from Figure11, which depicts the spring deformation between neighboring locomotives/wagons. We can see that the largest deformation is 0.14 m, and the smallest deformation is -0.15 m, which is with the safety threshold [-0.2m,0.2m]. As illustrated in Figure 7, this is because the in-train force guarantees that the position difference between neighboring locomotives cannot be too large or too small. Simulation results show that the proposed ATO controller can guarantee the safe operation of high-speed trains.Figure 11
The coupler deformation between neighboring locomotives/wagons in Scenario B.


### 5.3. Further Discussions

In this subsection, we further discuss the performance of the proposed ATO controller in terms of practicality, computation efficiency, superiority, and robustness. For simplicity, we only consider Scenario A, where the train is comprised of four locomotives. Scenario B can also be analyzed similarly. The train parameters are set the same as in Table5, unless otherwise specified.


#### 5.3.1. Practicality

In the previous simulation setup, we assume that the maximum traction forceFmaxT and the maximum braking force FmaxB are constants. Despite simplifying the implementation, the constant maximum traction/braking force are not realistic as the operation is subject to the power limit of the train. This implies that the maximum traction and braking forces evolve with the operation speed; that is, the maximum traction and braking forces for the CRH-3 train are represented as in [33].(27)FmaxT=300-0.284v,0≤v≤119.7km/h266×119.7v,119.7<v≤300km/hFmaxB=300-0.281v,0≤v≤106.7km/h266×119.7v,106.7<v≤300km/h.The force output of each locomotive with power limit is shown in Figure12, where we find that the output force of each locomotive evolves with time. It is interesting to note that at t=400 s the maximum traction force decreases at the transition stage, where the speed of the train increases from 40 m/s to 75 m/s. Due to the constant-power operation, the maximum traction force will decrease with the increase of velocity. At t=900 s, it is shown that the maximum braking force increases when the velocity decreases from 75 m/s to 50 m/s.Figure 12
The force output of each locomotive with power limit.


#### 5.3.2. Computational Time

In this subsection, we consider the computation time of the proposed method with different prediction horizons.Within different control horizons, the nonlinear MPC and the proposed piecewise MPC are both implemented on the original nonlinear system. The simulations are performed under Windows 7 operating system with Intel Core i5-2400 CPU, 4 GB RAM, on a desktop computer. Both the nonlinear MPC and the proposed piecewise MPC provide comparable dynamic performance. Their computational times are summarized in Table6. From Table 6, it can be seen that the computational time of each proposed MPC step is less than the nonlinear MPC on all control horizon lengths. It is more obvious when the horizon length is extended. The comparative results present the computation efficiency of the proposed method. To obtain the average computational time of each MPC step, the repeated numbers of MPC steps are 1200.Table 6
Computational time (s) of each MPC step with and without linearization.
MPC with linearization MPC without linearization Min. Max. Average Min. Max. Average Nc=3 1.6 e - 3 3.8 e - 3 2.5 e - 3 2.1 e - 1 3.6 e - 1 2.4 e - 1 Nc=5 1.8 e - 3 6.4 e - 3 3.5 e - 3 2.5 e - 1 8.3 e - 1 4.7 e - 1 Nc=7 1.9 e - 3 7.3 e - 3 4.6 e - 3 3.9 e - 1 9.8 e - 1 6.5 e - 1


#### 5.3.3. Superiority

In this subsection, we verify the superiority of the proposed method through comparison with the classical PID method.The velocity difference between each locomotive and the reference velocity value is plotted in the following figures. Figures13(a) and 13(b) show the simulation results under PID and the proposed MPC methods, respectively. From Figure 13(a), the velocities of four locomotives have severe chattering at time 400 s and 900 s when the reference velocity varies, which can further cause drastic changes of the coupler and endanger safe operation of high-speed trains. The proposed method in this paper presents good control performance as shown in Figure 13(b). From this figure, we can see that the velocities of four locomotives do not change frequently and can converge to the desired velocity quickly. This control performance is conducive to the safe operation of the train.Figure 13
Comparison of MPC method with PID method in velocity error.
(a)
PID method (b)
MPC method


#### 5.3.4. Stochastic Perturbations

In this subsection, we consider the performance of the proposed method with external stochastic perturbations. The stochastic perturbations are added in the system dynamics (6). The maximum amplitude of the perturbation is assumed to be 1. Considering the constant term βj=0.0044, the perturbation is large enough to evaluate the robustness of the proposed method. The external stochastic perturbations are produced by using the Matlab function randn(n), which can create ann-by-n matrix of normally distributed random numbers with mean equal to zero and standard deviation equal to 1. The simulation result of the system with the stochastic perturbations is shown in Figure 14. From this figure, we can see that the practical velocities of locomotives can track the reference values well, though the stochastic perturbations exist. The simulation results show the good robustness of the proposed method.Figure 14
The velocity of each vehicle under MPC method with external stochastic perturbations.


## 5.1. Simulation Scenario

In order to evaluate the performance of the controller, the desired velocity curve including accelerating, decelerating, velocity step increase, velocity step decrease, and constant-velocity stages is characterized by the following equation:(25)vrt=0.4t0≤t≤10040100≤t≤40075400≤t≤90050900≤t≤1150-t+12001150≤t≤1200.In (25), the acceleration is 0.4m/s, and the deceleration is 0.5m/s, both of which are within the threshold of acceleration and velocity. In the operation of trains, the basic resistance is the main resistance for the trains. Thus, we only consider the basic resistance in this paper.The designed MPC based ATO controller considers the train velocity tracking performance, energy consumption, and intertrain force. With the train parameters, we provide two simulation scenarios to evaluate the performance of the proposed algorithm under different operating conditions. Specifically, we evaluate the performance of the controller with different train types: all-locomotive train and partial-locomotive train.


### 5.1.1. Scenario A

In this scenario, we evaluate the performance of the ATO controller with the all-locomotive train. To simplify the simulation, we choose the number of trains asn=4, and the train is comprised of all locomotives. This implies that each locomotive has the tracking motor and braking unit and thus can provide both the tracking force and the braking force. The parameters of the train are shown in Table 5, where mi is the mass of locomotive i. Since the numbers of passengers in each locomotive may be different, we choose m1=48.5t,m2=52.5t,m3=51.0t,m4=47.5t. For the MPC controller, the prediction horizon Np and the control horizon Nc also affect the controller performance. Without loss of generality, we choose Np=Nc=3. The weighting matrix gain is chosen as(26)Q=diag4,…,4︸n,R=diag1,…,1︸n,H=diag1,…,1︸n-1.


### 5.1.2. Scenario B

In this scenario, we evaluate the effectiveness of the ATO controller in partial-locomotive trains. To simplify the simulation, we choose the number of trains asn=4, and the train is comprised of 2 locomotives and 2 wagons. The locomotive can provide both the tracking force and the braking force, while the wagons can only supply the brake forces. In order to guarantee the total force the same as Scenario A, the maximum tracking force of each locomotive is set as FmaxT=600kN, and the other train parameters are the same as in Table 5.


## 5.1.1. Scenario A

In this scenario, we evaluate the performance of the ATO controller with the all-locomotive train. To simplify the simulation, we choose the number of trains asn=4, and the train is comprised of all locomotives. This implies that each locomotive has the tracking motor and braking unit and thus can provide both the tracking force and the braking force. The parameters of the train are shown in Table 5, where mi is the mass of locomotive i. Since the numbers of passengers in each locomotive may be different, we choose m1=48.5t,m2=52.5t,m3=51.0t,m4=47.5t. For the MPC controller, the prediction horizon Np and the control horizon Nc also affect the controller performance. Without loss of generality, we choose Np=Nc=3. The weighting matrix gain is chosen as(26)Q=diag4,…,4︸n,R=diag1,…,1︸n,H=diag1,…,1︸n-1.


## 5.1.2. Scenario B

In this scenario, we evaluate the effectiveness of the ATO controller in partial-locomotive trains. To simplify the simulation, we choose the number of trains asn=4, and the train is comprised of 2 locomotives and 2 wagons. The locomotive can provide both the tracking force and the braking force, while the wagons can only supply the brake forces. In order to guarantee the total force the same as Scenario A, the maximum tracking force of each locomotive is set as FmaxT=600kN, and the other train parameters are the same as in Table 5.


## 5.2. Simulation Results

We will provide simulation results of the ATO controller in different operating conditions. The simulations are conducted in Matlab/SIMULINK. We fist consider the simulations in Scenario A, which is based on all-locomotive train. Then we will further consider the simulation results in Scenario B, which is based on partial-locomotive trains. And in Scenario C, we focus on more realistic condition where the track force limit is dynamic.


### 5.2.1. Scenario A

The velocity curve of each locomotive is shown in Figure4, where x-axis represents the time andy-axis represents the velocity. From the figure, we can see that the velocities of four locomotives almost remain the same during the operation time. Specifically, in the accelerating and decelerating stages, each locomotive can operate correspondingly according to the given acceleration and deceleration. At t=400 s, we notice that the reference velocity has the step increase. From the zoomed-in figure, we can find that the locomotives move with the almost same velocity, where the maximum velocity error is 0.003m/s. From the zoomed-in figure from t=470 s to t=490 s, the velocities of locomotives reach the steady value at about t=490 s, where the settling time is 90 s and the steady error is 0.5m/s. In the steady state, the velocities of locomotives remain fixed. Based on the simulation results, we can find that the train can track the time-varying reference with small settling time and steady error, which verifies the effectiveness of the algorithm.Figure 4
The velocity curve of each locomotive in Scenario A.The tracking and braking forces of locomotives are shown in Figure5. In the accelerating stage, the resistance of the train increases with the velocity. In order to guarantee the constant acceleration, the tracking forces of locomotives increase gradually. In the deceleration stage, the resistance of the train decelerates with the velocity. In order to guarantee the constant deceleration, the braking forces of locomotives increase gradually. Compared with the acceleration or deceleration stages, the tracking forces of locomotives are small and fixed during the constant velocity stage. This is because the resistance of train keeps having a relatively small value during the constant velocity stage. At t=400 s, the reference has a step increase. In order to accelerate to the desired velocity, the tracking forces of locomotives increase to the maximum rapidly and maintain the maximum tracking force for about 72 s. At t=900 s, the braking forces increase to the maximum to decrease the velocity and maintain about 32 s to decrease to the desired velocity. From the figure, we can find that the control outputs of locomotives are almost the same and thus guarantee the reasonable allocation of braking force and tracking force. When the reference velocity curve changes, each locomotive can regulate the tracking force and braking force based on the difference between the actual velocity and desired velocity. The forces for locomotives are almost the same and also satisfy the constraints, which verify the effectiveness of the proposed strategy.Figure 5
The force output of each locomotive in Scenario A.The coupler force among each locomotive is shown in Figure6. When the coupler force is larger than zero, it means that coupler is pulled. When the coupler force is smaller than zero, the coupler is pushed. From the figure, we can see that the coupler forces among locomotives are always within the safety range. Moreover, the coupler among L1 and L2 is mainly pulled, while the coupler among L3 and L4 is mainly pushed. This is because when we choose the masses of locomotives, the masses of L3 and L4 are relatively smaller than L1 and L2. In order to guarantee the velocity consensus of locomotives, the total forces should be the same and the coupler force can balance the total force of each locomotive. It can be found that, during the constant-velocity stage, constant-accelerating stage, and constant-decelerating stage, the coupler force among locomotives is small. At t=200 s and t=400 s, the coupler force changes dramatically, with the largest coupler pulling force 317 kN and the largest coupler pushing force 241 kN. From the figure, we can see that the ATO controller can guarantee that the coupler force is within the safety range, which verifies the effectiveness.Figure 6
The coupler force among locomotives in Scenario A.The operation safety of high-speed trains can be evaluated by the position difference (coupler deformation) between neighboring locomotives. If the deformation is larger than the maximum0.2m or less than the minimum -0.2m, the coupler may have been broken so the safety cannot be guaranteed. The operation safety of high-speed trains can be further illustrated from Figure 7, which shows that the coupler deformation is always with the threshold [-0.2m,0.2m]. This is because the in-train force guarantees that the position difference between neighboring locomotives cannot be too large or too small.Figure 7
The coupler deformation between two neighbouring locomotives in Scenario A.


### 5.2.2. Scenario B

The simulation in Scenario A verifies the effectiveness of the hybrid model predictive controller for all-locomotive case. The metrics include velocity tracking, energy consumption, and coupler force.In Scenario B, we verified the effectiveness of the ATO controller with 2 locomotives and 2 wagons. The velocity curves are shown in Figure8. With the same controller parameters, at t=400 s, where the reference velocity has a step change, the convergence time of the controller from the change at t=400 s to the steady state at t=508 s implies that the convergence time increases by 18 s compared with the all-locomotive case. The steady error is 1.2 m/s, which increases by 0.7m/s compared with the all-locomotive case. Moreover, the velocities at the steady stage have some oscillations around the steady value.Figure 8
The velocity curve of each vehicle in Scenario B.The force output of each vehicle is shown in Figure9, where the left subfigure shows the force output of the locomotive, and the right subfigure shows the force output of the wagon. From the figure, we can see that, during the acceleration and constant-velocity stages, only the locomotive has the tracking force output, and the tracking force of wagons is zero. In the deceleration stage, both the locomotives and wagons supply the braking forces.Figure 9
The forces of locomotives and wagons in Scenario B.Figure10 shows the coupler force among vehicles. Compared with the all-locomotive trains, the coupler forces between L1 and W1 and between L2 and W2 are kept at large values. This is because locomotive L1 supplies the pulling force and locomotive L2 supplies the pushing force, while wagons W1 and W2 cannot supply the tracking force. The maximum coupler force is 578 kN, which is within the safety threshold. From the above simulations, we can find that the proposed hybrid MPC ATO controller has a good control effect.Figure 10
The coupler force among vehicles in Scenario B.The operation safety of high-speed trains can be illustrated from Figure11, which depicts the spring deformation between neighboring locomotives/wagons. We can see that the largest deformation is 0.14 m, and the smallest deformation is -0.15 m, which is with the safety threshold [-0.2m,0.2m]. As illustrated in Figure 7, this is because the in-train force guarantees that the position difference between neighboring locomotives cannot be too large or too small. Simulation results show that the proposed ATO controller can guarantee the safe operation of high-speed trains.Figure 11
The coupler deformation between neighboring locomotives/wagons in Scenario B.


## 5.2.1. Scenario A

The velocity curve of each locomotive is shown in Figure4, where x-axis represents the time andy-axis represents the velocity. From the figure, we can see that the velocities of four locomotives almost remain the same during the operation time. Specifically, in the accelerating and decelerating stages, each locomotive can operate correspondingly according to the given acceleration and deceleration. At t=400 s, we notice that the reference velocity has the step increase. From the zoomed-in figure, we can find that the locomotives move with the almost same velocity, where the maximum velocity error is 0.003m/s. From the zoomed-in figure from t=470 s to t=490 s, the velocities of locomotives reach the steady value at about t=490 s, where the settling time is 90 s and the steady error is 0.5m/s. In the steady state, the velocities of locomotives remain fixed. Based on the simulation results, we can find that the train can track the time-varying reference with small settling time and steady error, which verifies the effectiveness of the algorithm.Figure 4
The velocity curve of each locomotive in Scenario A.The tracking and braking forces of locomotives are shown in Figure5. In the accelerating stage, the resistance of the train increases with the velocity. In order to guarantee the constant acceleration, the tracking forces of locomotives increase gradually. In the deceleration stage, the resistance of the train decelerates with the velocity. In order to guarantee the constant deceleration, the braking forces of locomotives increase gradually. Compared with the acceleration or deceleration stages, the tracking forces of locomotives are small and fixed during the constant velocity stage. This is because the resistance of train keeps having a relatively small value during the constant velocity stage. At t=400 s, the reference has a step increase. In order to accelerate to the desired velocity, the tracking forces of locomotives increase to the maximum rapidly and maintain the maximum tracking force for about 72 s. At t=900 s, the braking forces increase to the maximum to decrease the velocity and maintain about 32 s to decrease to the desired velocity. From the figure, we can find that the control outputs of locomotives are almost the same and thus guarantee the reasonable allocation of braking force and tracking force. When the reference velocity curve changes, each locomotive can regulate the tracking force and braking force based on the difference between the actual velocity and desired velocity. The forces for locomotives are almost the same and also satisfy the constraints, which verify the effectiveness of the proposed strategy.Figure 5
The force output of each locomotive in Scenario A.The coupler force among each locomotive is shown in Figure6. When the coupler force is larger than zero, it means that coupler is pulled. When the coupler force is smaller than zero, the coupler is pushed. From the figure, we can see that the coupler forces among locomotives are always within the safety range. Moreover, the coupler among L1 and L2 is mainly pulled, while the coupler among L3 and L4 is mainly pushed. This is because when we choose the masses of locomotives, the masses of L3 and L4 are relatively smaller than L1 and L2. In order to guarantee the velocity consensus of locomotives, the total forces should be the same and the coupler force can balance the total force of each locomotive. It can be found that, during the constant-velocity stage, constant-accelerating stage, and constant-decelerating stage, the coupler force among locomotives is small. At t=200 s and t=400 s, the coupler force changes dramatically, with the largest coupler pulling force 317 kN and the largest coupler pushing force 241 kN. From the figure, we can see that the ATO controller can guarantee that the coupler force is within the safety range, which verifies the effectiveness.Figure 6
The coupler force among locomotives in Scenario A.The operation safety of high-speed trains can be evaluated by the position difference (coupler deformation) between neighboring locomotives. If the deformation is larger than the maximum0.2m or less than the minimum -0.2m, the coupler may have been broken so the safety cannot be guaranteed. The operation safety of high-speed trains can be further illustrated from Figure 7, which shows that the coupler deformation is always with the threshold [-0.2m,0.2m]. This is because the in-train force guarantees that the position difference between neighboring locomotives cannot be too large or too small.Figure 7
The coupler deformation between two neighbouring locomotives in Scenario A.


## 5.2.2. Scenario B

The simulation in Scenario A verifies the effectiveness of the hybrid model predictive controller for all-locomotive case. The metrics include velocity tracking, energy consumption, and coupler force.In Scenario B, we verified the effectiveness of the ATO controller with 2 locomotives and 2 wagons. The velocity curves are shown in Figure8. With the same controller parameters, at t=400 s, where the reference velocity has a step change, the convergence time of the controller from the change at t=400 s to the steady state at t=508 s implies that the convergence time increases by 18 s compared with the all-locomotive case. The steady error is 1.2 m/s, which increases by 0.7m/s compared with the all-locomotive case. Moreover, the velocities at the steady stage have some oscillations around the steady value.Figure 8
The velocity curve of each vehicle in Scenario B.The force output of each vehicle is shown in Figure9, where the left subfigure shows the force output of the locomotive, and the right subfigure shows the force output of the wagon. From the figure, we can see that, during the acceleration and constant-velocity stages, only the locomotive has the tracking force output, and the tracking force of wagons is zero. In the deceleration stage, both the locomotives and wagons supply the braking forces.Figure 9
The forces of locomotives and wagons in Scenario B.Figure10 shows the coupler force among vehicles. Compared with the all-locomotive trains, the coupler forces between L1 and W1 and between L2 and W2 are kept at large values. This is because locomotive L1 supplies the pulling force and locomotive L2 supplies the pushing force, while wagons W1 and W2 cannot supply the tracking force. The maximum coupler force is 578 kN, which is within the safety threshold. From the above simulations, we can find that the proposed hybrid MPC ATO controller has a good control effect.Figure 10
The coupler force among vehicles in Scenario B.The operation safety of high-speed trains can be illustrated from Figure11, which depicts the spring deformation between neighboring locomotives/wagons. We can see that the largest deformation is 0.14 m, and the smallest deformation is -0.15 m, which is with the safety threshold [-0.2m,0.2m]. As illustrated in Figure 7, this is because the in-train force guarantees that the position difference between neighboring locomotives cannot be too large or too small. Simulation results show that the proposed ATO controller can guarantee the safe operation of high-speed trains.Figure 11
The coupler deformation between neighboring locomotives/wagons in Scenario B.


## 5.3. Further Discussions

In this subsection, we further discuss the performance of the proposed ATO controller in terms of practicality, computation efficiency, superiority, and robustness. For simplicity, we only consider Scenario A, where the train is comprised of four locomotives. Scenario B can also be analyzed similarly. The train parameters are set the same as in Table5, unless otherwise specified.


### 5.3.1. Practicality

In the previous simulation setup, we assume that the maximum traction forceFmaxT and the maximum braking force FmaxB are constants. Despite simplifying the implementation, the constant maximum traction/braking force are not realistic as the operation is subject to the power limit of the train. This implies that the maximum traction and braking forces evolve with the operation speed; that is, the maximum traction and braking forces for the CRH-3 train are represented as in [33].(27)FmaxT=300-0.284v,0≤v≤119.7km/h266×119.7v,119.7<v≤300km/hFmaxB=300-0.281v,0≤v≤106.7km/h266×119.7v,106.7<v≤300km/h.The force output of each locomotive with power limit is shown in Figure12, where we find that the output force of each locomotive evolves with time. It is interesting to note that at t=400 s the maximum traction force decreases at the transition stage, where the speed of the train increases from 40 m/s to 75 m/s. Due to the constant-power operation, the maximum traction force will decrease with the increase of velocity. At t=900 s, it is shown that the maximum braking force increases when the velocity decreases from 75 m/s to 50 m/s.Figure 12
The force output of each locomotive with power limit.


### 5.3.2. Computational Time

In this subsection, we consider the computation time of the proposed method with different prediction horizons.Within different control horizons, the nonlinear MPC and the proposed piecewise MPC are both implemented on the original nonlinear system. The simulations are performed under Windows 7 operating system with Intel Core i5-2400 CPU, 4 GB RAM, on a desktop computer. Both the nonlinear MPC and the proposed piecewise MPC provide comparable dynamic performance. Their computational times are summarized in Table6. From Table 6, it can be seen that the computational time of each proposed MPC step is less than the nonlinear MPC on all control horizon lengths. It is more obvious when the horizon length is extended. The comparative results present the computation efficiency of the proposed method. To obtain the average computational time of each MPC step, the repeated numbers of MPC steps are 1200.Table 6
Computational time (s) of each MPC step with and without linearization.
MPC with linearization MPC without linearization Min. Max. Average Min. Max. Average Nc=3 1.6 e - 3 3.8 e - 3 2.5 e - 3 2.1 e - 1 3.6 e - 1 2.4 e - 1 Nc=5 1.8 e - 3 6.4 e - 3 3.5 e - 3 2.5 e - 1 8.3 e - 1 4.7 e - 1 Nc=7 1.9 e - 3 7.3 e - 3 4.6 e - 3 3.9 e - 1 9.8 e - 1 6.5 e - 1


### 5.3.3. Superiority

In this subsection, we verify the superiority of the proposed method through comparison with the classical PID method.The velocity difference between each locomotive and the reference velocity value is plotted in the following figures. Figures13(a) and 13(b) show the simulation results under PID and the proposed MPC methods, respectively. From Figure 13(a), the velocities of four locomotives have severe chattering at time 400 s and 900 s when the reference velocity varies, which can further cause drastic changes of the coupler and endanger safe operation of high-speed trains. The proposed method in this paper presents good control performance as shown in Figure 13(b). From this figure, we can see that the velocities of four locomotives do not change frequently and can converge to the desired velocity quickly. This control performance is conducive to the safe operation of the train.Figure 13
Comparison of MPC method with PID method in velocity error.
(a)
PID method (b)
MPC method


### 5.3.4. Stochastic Perturbations

In this subsection, we consider the performance of the proposed method with external stochastic perturbations. The stochastic perturbations are added in the system dynamics (6). The maximum amplitude of the perturbation is assumed to be 1. Considering the constant term βj=0.0044, the perturbation is large enough to evaluate the robustness of the proposed method. The external stochastic perturbations are produced by using the Matlab function randn(n), which can create ann-by-n matrix of normally distributed random numbers with mean equal to zero and standard deviation equal to 1. The simulation result of the system with the stochastic perturbations is shown in Figure 14. From this figure, we can see that the practical velocities of locomotives can track the reference values well, though the stochastic perturbations exist. The simulation results show the good robustness of the proposed method.Figure 14
The velocity of each vehicle under MPC method with external stochastic perturbations.


## 5.3.1. Practicality

In the previous simulation setup, we assume that the maximum traction forceFmaxT and the maximum braking force FmaxB are constants. Despite simplifying the implementation, the constant maximum traction/braking force are not realistic as the operation is subject to the power limit of the train. This implies that the maximum traction and braking forces evolve with the operation speed; that is, the maximum traction and braking forces for the CRH-3 train are represented as in [33].(27)FmaxT=300-0.284v,0≤v≤119.7km/h266×119.7v,119.7<v≤300km/hFmaxB=300-0.281v,0≤v≤106.7km/h266×119.7v,106.7<v≤300km/h.The force output of each locomotive with power limit is shown in Figure12, where we find that the output force of each locomotive evolves with time. It is interesting to note that at t=400 s the maximum traction force decreases at the transition stage, where the speed of the train increases from 40 m/s to 75 m/s. Due to the constant-power operation, the maximum traction force will decrease with the increase of velocity. At t=900 s, it is shown that the maximum braking force increases when the velocity decreases from 75 m/s to 50 m/s.Figure 12
The force output of each locomotive with power limit.


## 5.3.2. Computational Time

In this subsection, we consider the computation time of the proposed method with different prediction horizons.Within different control horizons, the nonlinear MPC and the proposed piecewise MPC are both implemented on the original nonlinear system. The simulations are performed under Windows 7 operating system with Intel Core i5-2400 CPU, 4 GB RAM, on a desktop computer. Both the nonlinear MPC and the proposed piecewise MPC provide comparable dynamic performance. Their computational times are summarized in Table6. From Table 6, it can be seen that the computational time of each proposed MPC step is less than the nonlinear MPC on all control horizon lengths. It is more obvious when the horizon length is extended. The comparative results present the computation efficiency of the proposed method. To obtain the average computational time of each MPC step, the repeated numbers of MPC steps are 1200.Table 6
Computational time (s) of each MPC step with and without linearization.
MPC with linearization MPC without linearization Min. Max. Average Min. Max. Average Nc=3 1.6 e - 3 3.8 e - 3 2.5 e - 3 2.1 e - 1 3.6 e - 1 2.4 e - 1 Nc=5 1.8 e - 3 6.4 e - 3 3.5 e - 3 2.5 e - 1 8.3 e - 1 4.7 e - 1 Nc=7 1.9 e - 3 7.3 e - 3 4.6 e - 3 3.9 e - 1 9.8 e - 1 6.5 e - 1


## 5.3.3. Superiority

In this subsection, we verify the superiority of the proposed method through comparison with the classical PID method.The velocity difference between each locomotive and the reference velocity value is plotted in the following figures. Figures13(a) and 13(b) show the simulation results under PID and the proposed MPC methods, respectively. From Figure 13(a), the velocities of four locomotives have severe chattering at time 400 s and 900 s when the reference velocity varies, which can further cause drastic changes of the coupler and endanger safe operation of high-speed trains. The proposed method in this paper presents good control performance as shown in Figure 13(b). From this figure, we can see that the velocities of four locomotives do not change frequently and can converge to the desired velocity quickly. This control performance is conducive to the safe operation of the train.Figure 13
Comparison of MPC method with PID method in velocity error.
(a)
PID method (b)
MPC method


## 5.3.4. Stochastic Perturbations

In this subsection, we consider the performance of the proposed method with external stochastic perturbations. The stochastic perturbations are added in the system dynamics (6). The maximum amplitude of the perturbation is assumed to be 1. Considering the constant term βj=0.0044, the perturbation is large enough to evaluate the robustness of the proposed method. The external stochastic perturbations are produced by using the Matlab function randn(n), which can create ann-by-n matrix of normally distributed random numbers with mean equal to zero and standard deviation equal to 1. The simulation result of the system with the stochastic perturbations is shown in Figure 14. From this figure, we can see that the practical velocities of locomotives can track the reference values well, though the stochastic perturbations exist. The simulation results show the good robustness of the proposed method.Figure 14
The velocity of each vehicle under MPC method with external stochastic perturbations.


## 6. Conclusion

In this paper, we investigated the automatic cruising control problem of high-speed trains with the aim of tracking the reference speed and minimizing the energy consumption and guaranteeing the safe operation of high-speed trains. First, we built the piecewise linearization of train resistance and developed the piecewise multiple-point model of trains. Due to the multiple constraints and performance metrics, we proposed a model predictive control method to design the ATO controller. In order to improve the accuracy of the method, we introduced the binary variable and transformed the piecewise model into a mixed logical dynamical model and guarantee that the submodel in the prediction process is still switchable. Then we analyze the operation objective and constraints and transformed the tracking control problem into a constrained finite-time optimal control problem. The optimization problem is solved by the integer linear planning method, and the optimal control solution is obtained. Simulations in the Matlab verified the effectiveness of the proposed algorithm.In this work, we just consider the nonlinear characteristics of resistance when constructing the train model and neglect the control delay. In the future work, more realistic train models will be investigated to incorporate more nonlinear characteristics and the controller robustness will be improved under time delay.



---

*Source: 7308058-2018-05-15.xml*