1. Introduction
In the modern metro system, the quality of operation, including energy saving, punctuality, and comfortableness, is the core factor for operation [
1]. With the popularization of the Automatic Train Operation (ATO) system, the automation of train operation has been achieved. However, the performance of the algorithm in ATO system is not good enough, and it can be further improved [
2]. Inaccuracy and high energy consumption seriously affect the quality of train operation.
The frequent switching of operation mode is needed for the metro system because the start–stop is frequent. The ATO system needs to refer to the planned speed profile as a reference to issue control commands, such as acceleration, braking, and coasting [
3]. In practice, the PID controller is widely used to realize speed control, and the fluctuation of control force occurs when the coasting mode is needed [
4]. This problem eventually leads to larger speed errors and poor energy saving effects. To improve the tracking quality, a large number of papers are published to try to find a perfect substitute of the PID controller. Predictive controller is presented to track the speed profile. In Reference [
5], an integral sliding mode control (ISMC) is used to design distributed adaptive control strategies. In Reference [
6], a tracking controller is formulated to resolve errors caused by unmolded effects. In References [
7,
8], a robust adaptive controller is proposed to adapt unknown system parameters and nonlinear traction/braking notches. Wang et al. proposed an approach to design a distributed and fault-tolerant controller which is applied to develop braking control schemes for high-speed trains subject to traction and braking failures [
9]. In Reference [
10], to deal with input saturation and unknown disturbance in high-speed trains, a neural adaptive fault-tolerant controller is designed. In all the abovementioned studies, either the convergence cannot be proved by Lyapunov theory or the convergence speed is not good enough. In addition, large speed errors in those controllers also result in large energy consumption. Therefore, they are not suitable for metro trains that require fast convergence and high-precision tracking.
Recently, sliding mode control is becoming a research hotspot. Because it is easy to design and did not require accurate parameter identification, it has been used in actual physical systems, such as power systems and industrial robot manipulators [
11,
12,
13]. In addition, the controller with a nonlinear sliding surface has been proven to converge in a limited time [
14]. Especially for the adaptive terminal sliding mode controller (ATSMC), in addition to the advantages of fast convergence, it can also be adaptive to the internal parameter disturbance and external interference of the system [
15]. Inspired by the above articles, we introduce the additional resistance caused by slope and the basic running resistance of the train into the sliding mode controller, in the form of disturbance, and an adaptive terminal sliding mode controller for train-speed-tracking control (ATSMC) is designed in this paper. In order to further suppress the disturbance and enhance the control accuracy, which are important to energy saving, comfortableness, and punctuality, this paper introduces a disturbance observer (DOB) on the basis of ATSMC to track the speed profile; finally, ATSMC+DOB is designed, and the convergence is proven by Lyapunov theory.
As the premise of speed-tracking control, the reference speed profile needs to be generated by the speed optimization algorithm. Similar to the speed-tracking controller, the speed optimization algorithm used in engineering is not good enough either. As a widely used cruise optimization algorithm (cruising mode), it is easy to calculate, in theory. However, it is not only poor in energy-saving, but also difficult to achieve precise control in industry. Therefore, the coast optimization algorithm is adopted by a large number of studies, to reduce energy consumption. It provides a reference speed curve that is convenient for the controller to track. The energy saving of coast optimization has been proven with Pontryagin maximum principle by some studies [
16,
17,
18]. In References [
19,
20], based on coast optimization rules, the genetic algorithm is used to search the optimal coast point. Chuang et al. used an optimization method based on neural networks to calculate the optimal solution [
21]. Khmelnitsky et al. used the maximum–minimum ant system (MMAS) in the ant colony optimization algorithm (ACO) to optimize the train’s running speed trajectory [
22]. In Reference [
23], a dual heuristic programming method is used to optimize the coast time and dwell time at the same time. In Reference [
24], a kind of dual speed profile optimization is designed, to save energy. Cheng et al. proposed an improved ant colony optimization method to study train operation [
25]. Rocha et al. used the coast mode to design a n approach for real time train optimization [
26]. In References [
27,
28], multi-objective optimization is considered to study speed profile. Shang et al. proposed an online energy-saving driving strategy for metro [
29]. In Reference [
30], a switched nonlinear model predictive controller (NMPC) is designed for collaborative ecodrive control of railway vehicles. Furthermore, to use more renewable energy, a dissension-based adaptive law (DAL) is proposed, to adjust the parameters of the NMPC cost. In Reference [
31], a kind of online distributed cooperative model predictive controller is proposed, to generate a speed profile for multiple high-speed train movement. In Reference [
32], a kind of energy-efficient driving strategy based on transmission losses is proposed.
In the above research, most optimization models focus on energy saving. The lack of comfortableness and punctuality is fatal to the metro. Therefore, a multi-objective optimization model is established in this paper. In this model, in addition to energy consumption, comfortableness and punctuality are taken into consideration as penalty factors to comprehensively optimize operation quality.
In addition, for the intelligent evolutionary algorithm mentioned above, there are always two problems to be solved. Firstly, maintaining good population diversity to avoid a local optimal solution while evolving rapidly is always a problem [
33]. Secondly, the initial value of the population has a significant impact on the optimization result, so how to obtain a better initial value is worth studying. To solve these problems, the artificial bee colony (ABC) algorithm, which shows good global search ability, may be used in speed profile optimization [
34]. However, though good at global search, the ABC algorithm is poor in local search. When it is close to the optimal solution, the search efficiency drops sharply in a local region. To make up for this deficiency, a regional reinforcement artificial bee colony (RR-ABC) algorithm is proposed, to optimize the speed profile. By using chaotic mirror initialization, local search strategy, and the elite group guidance strategy, the RR-ABC algorithm maintains better population diversity and enhances search efficiency.
The specific structure is as follows.
Section 2 shows the dynamic model and some special restrictions of metro.
Section 3 establishes a speed profile optimization calculation model which considers energy consumption, comfortableness, and punctuality. In
Section 4, the RR-ABC algorithm is proposed, to optimize the speed profile. In
Section 5, to track the speed, we propose the ATSMC+DOB controller. In
Section 6, the validity of the method is verified based on the practical data of Qingdao Metro Line 6. Finally, we conclude the study in
Section 7.
2. Dynamic Characteristics of Metro Trains
There are two common dynamic models: multi-mass point model and single-particle model. For the multi-particle model, each vehicle of the train is treated as an independent mass point. As a result, the force deviation caused by train length and the power transmission between different vehicles are considered in the model. However, this model has a complicated modeling process and is also difficult to calculate [
35]. Due to the lighter weight and shorter length of the metros, this accuracy is not necessary. Compared with multi-particle model, the single-particle model regards the entire train as a mass point which has a simple modeling process. A large number of related studies have adopted the single-particle model [
36,
37].
Train single-particle dynamics model in Reference [
36] is used in this paper. The schematic diagram of the train force is shown in
Figure 1. In this model,
is the train mass, which is considered as a constant value throughout the journey. About the decomposition of the forces, denotes the traction force,
denotes the braking force,
denotes the additional resistance of the line that mainly comes from the slope, and
denotes the basic running resistance of the train, mainly from the track and air. The dynamics model can be written as follows:
where
denotes the speed, and
is the position;
is the force of traction,
is the force of braking,
is the additional resistance, and
is the basic resistance; and
is the rotary mass coefficient which is usually a given constant.
About the rotary mass coefficient, , it is a coefficient to compensate the moment of inertia. The train kinetic energy has two parts, one is the horizontal movement energy of the train, and the other is the moment of inertia of the wheels. Therefore, when only the horizontal movement is concerned, the moment of inertia needs to be compensated. This value can be obtained by dividing the converted weight of the moment of inertia by the total weight. Since this value is difficult to determine, empirical values are generally used in simulation. Therefore, is the total weight of the train after correction.
About the traction, the maximum that the train can provide is as follows:
where
is the real force, and
control factor.
About the braking, similar to traction, the calculation formula can be written as follows:
About the basic running resistance, it can be written as follows:
where
,
, and
are empirical parameters determined by experiment.
The additional resistance can be written as follows:
where
is the additional resistance of the slope,
is the additional resistance in the curve, and
is the additional resistance in the tunnel. Generally, the resistance in the curve and tunnel are not considered in the study and the slope resistance is given as:
where
is the train mass,
is the acceleration of gravity, and
is the slope value.
The external characteristic curves (max force) of train traction and braking force are shown in
Figure 2, and Equations (2), (4), and (6) are examined by the actual experiment from Qingdao Metro Line 6.
4. Speed Profile Optimization with RR-ABC
4.1. Emergency Braking Intervention (EBI) Curve and Warning Speed Limit
Resulted from some safety reasons, the metro line should consider various speed-limit protections.
Firstly, the metro line has a fixed speed limit from the line. In addition, the metro also has a security protection mechanism to ensure that emergency braking is feasible in any position that the train cannot cross the fixed speed limit; this mechanism is called emergency braking intervention (EBI). To ensure that the metro does not exceed the fixed speed limit during emergency braking under any circumstances, the EBI protection considers such an extreme case:
At first, the train is in the process of maximum acceleration. After the train issued an emergency braking command, the train still maintained the maximum acceleration during the transmission delay and traction cutoff delay, because the traction cannot be cut off instantly and transmission delay exists objectively at any moment.
Since the braking system needs some time to receive instructions and gradually generate braking force, it usually takes some time to reach the maximum force from 0. During this period, the train is considered to keep coasting.
Therefore, as shown in
Figure 7, the emergency braking process is divided into three stages:
Stage1: In this stage, the train accelerates with the maximum traction acceleration, and the acceleration time is the sum of the traction cutoff delay and the braking system action delay.
Stage2: After the train traction is cut off, the emergency braking force needs a certain time to reach the nominal value, and the train is considered to coast during this period.
Stage3: The train follows the emergency braking curve to brake.
Finally, as a result of the single-mass model, we ignore the size and length of the train. However, the ignorance will bring security risks to trains, so we design the warning speed limit based on the EBI curve.
As shown in
Figure 8, when the train moves from the low-speed-limit zone to the high-speed-limit zone, the train needs to reserve a buffer equal to the length of the train. When the train moves from the high-speed-limit zone to the low-speed-limit zone, the same buffer is required.
According to the above principle, taking Qingdao Metro Line 6 grabbing Zhuomashan station to Heluobu station as an example, these kinds of speed limit curves are obtained as shown in
Figure 9.
4.2. Energy-Saving Strategies and Coast Choices
According to actual operating conditions and theory, coast travel is an effective strategy for train energy saving. The coast strategy we designed is shown in the
Figure 10, in which two different coast intervals are set.
4.2.1. Coast Interval
Features:
- (1)
Driving from high speed limit zone to low speed limit zone
- (2)
Maintain constant speed at the next low speed limit
Strategy selection:
From the coast point, let the train maintain the coast condition until point , then we get the speed of point .
- (1)
If , then brake to ;
- (2)
If , then when , let the train drive at a constant speed .
4.2.2. Coast Interval
Features:
- (1)
Driving from high speed limit zone to low speed limit zone;
- (2)
Braking in the next low speed limit.
Strategy selection:
From the coast point, let the train maintain the coast condition until point , then we get the speed of point .
- (1)
If , then brake to ;
- (2)
If , then keep coasting to maintain the coast condition until it contacts the braking curve.
4.3. Regional Reinforcement Artificial Bee Colony (RR-ABC) Algorithm
The ABC algorithm simulates the process of the bee colony searching for the honey source with the highest honey content. It has few control parameters, and it has good global convergence. However, when the local search is approached to converge to the optimal solution, the search efficiency drops sharply [
39,
40]. With the deepening of evolution, population diversity is seriously lacking, and the search efficiency is significantly reduced at the end of evolution. To balance the global and local search capabilities and to improve efficiency of the ABC, a more efficient ABC algorithm is proposed called RR-ABC algorithm. The principles and steps of the RR-ABC algorithm are as follows:
- (1)
Model initialization: Chaos mirror initialization
Because of the randomness of the chaotic map and the sensitivity to the initial conditions, we extracted the decision space information to grow the diversity of the group. In this part, the sine formula is used to iteratively generate chaotic variables. The chaotic initialization Equation is as follows:
where
is random number between 0 and 1;
is number of iterations;
, and
;
is the j-dimensional parameter of the i-th honey source; and
,
are the top and bottom limits of the source.
Use Equation (19) to perform the chaos initialization operation to obtain the honey source , then mirror it to get the mirror honey source , then compare the fitness value of the initial solution one by one, and select the better individual as the initial population.
- (2)
Employed bee stage: Evolutionary dimension adjustment strategy
In the stage of employed bee and onlooker bee, they need to update the honey source through Equation (20); the original honey source has only one dimension for each change. Obviously, this method is extremely inefficient. To solve the problem, we use a new parameter
, to control the dimensions that need to be changed in each evolution. The new search Equation is as follows:
where
, parameter
controls the number of dimensions of each individual change, and
.
Every time an employed bee searches for a new honey source, it will estimate the potential of the source (calculate the fitness value of the objective function):
where
is the objective function value of the i-th honey source;
is the fitness value of i-th honey source. The employed bee chooses the better honey source by comparing the fitness value of
and
(Greedy choice). The larger the fitness value, the smaller the corresponding objective function value, indicating that the honey source is also better, and with a higher probability of getting better honey sources around it.
- (3)
Onlooker bees stage: Search radius adjustment strategy
As mentioned in “Model of employed bees”, the employed bee will search within a certain radius of the initial solution. In the early stage of evolution, because the honey source is far from the real solution, the use of a larger search radius can accelerate the convergence rate; when the honey source evolves to the final stage, it is closer to the real solution, and reducing the search radius can improve the search accuracy. Therefore, using a variable search radius and adjusting the size of the search radius can significantly speed up the convergence rate.
In this paper, the strategy of changing the search radius with the individual threshold, that is, as the threshold,
(the digit of times that honey source
was not selected), increases, the search radius continues to decrease. The Equation is as follows:
where
is search radius honey source
,
is the basic value of search radius,
is the threshold of the honey source
, and
is the maximum allowed threshold. The new search Equation is as follows:
Remark 1. For the traditional ABC algorithm, with the deepening of evolution, more nectar sources are continuously explored on a global scale. However, due to the unknown quality of the preliminary search results, each explored nectar needs to be searched more deeply, which seriously affects the search efficiency in the later stages of evolution, even if it can avoid the generation of local optimal solutions. In addition, for classic biomimetic algorithms, the quality of the initial value usually greatly affects the evolution process, while a simple random process can easily produce low-quality initial values. Through the flowchart about the comparison of classic method and reinforcement method (Figure 11), we can see that dynamically adjusted evolutionary dimensions and search radius are adopted. In the later stages of the search, fewer high-quality nectar sources will be paid attention, and a better efficiency can be got. 4.4. Speed Profile Optimization Model
Based on the above discussion, the speed curve optimization model is summarized as follows:
5. Adaptive Terminal Sliding Mode Controller
This section needs to design a speed tracking controller to track the speed profile optimized before. As mentioned before, compared with PID controller, which is widely used in modern industrial control, the sliding mode controller is quick to respond and completely robust to system parameter changes and external interference. However, sliding mode control will produce chattering when dealing with uncertain factors. This kind of high-frequency chattering will not only affect the accuracy of the control, but also cause the system to oscillate or become unstable.
Adaptive terminal sliding mode control not only has strong robustness, but also lets the system state converge to the ideal trajectory within a limited time. In addition, it also has a strong parameter adaptive processing function, to ensure that the system will not have unnecessary continuous switching when the parameters are uncertain. In order to design a speed-tracking controller with excellent comprehensive performance, this paper designs an adaptive terminal sliding mode controller with disturbance observer.
5.1. Dynamics Model of Speed Tracking
This part models and transforms the known fixed time delay. Here, the time delay is compensated by the sliding mode control coefficient, so that the system can be used as a non-delay system to design the sliding mode control rate in the subsequent chapters.
As mentioned in
Section 2, the ideal dynamics model is Equation (1); for convenience of expression, the model with speed
as its independent variable is rewritten as the model with time,
, as its independent variable. Furthermore, the traction force factor,
, and braking force factor,
, are replaced by the actual generated traction/braking force,
, and the model is summarized as follows:
In the actual control system, due to response delay, the model needs to be written as a first-order lag model:
where
is the control acceleration, which is produced by the controller;
is the target acceleration;
is the response time constant; and
is controlling transmission delay.
For subsequent kinetic description, we define the following Equation:
Therefore, the acceleration formed by the basic resistance of the train and the additional resistance can be expressed as follows:
Furthermore, according to Equation (30), we can get the following:
where
is the resistance acceleration.
The existing metro controllers are mostly stepless control modes, and the control system model is shown in
Figure 12.
In order to facilitate the design of subsequent controllers, pade approximation is used to rewrite the delay parameters
, in Equation (28):
where
is a constant related to delay,
.
Laplace transform of Equation (28):
In order to facilitate expression, Equation (33) is rewritten as follows:
As we can see, if the response time constant and the delay are known and kept unchanged, then , , and are all constants. In the following content ,, and are used to design the sliding mode control rate, so the delay will be compensated. The time delay is compensated to the sliding mode control rate coefficient through the abovementioned method, and the effect of delay is eliminated in this way.
5.2. Design of Sliding Mode Terminal Controller
We define the state error as follows:
where
is the train position error,
is the train speed error,
is the reference position of optimized profile, and
is the reference speed of optimized profile.
Differentiate Equation (35) and substitute it into Equation (31), and the resulting error state space Equation is as follows:
To improve the tracking accuracy, the designed sliding mode surface needs to introduce the train position error
, and speed error
, to ensure fast error convergence. The designed terminal sliding mode function is as follows:
where
; and
and
are positive odds,
.
Design the sliding mode controller below.
Differentiate the sliding mode function of (37) as follows:
According to Equations (30) and (31), we can obtain the following:
Substitute Equation (39) into Equation (34), since
, ignore
, and get the following:
By combining Equations (39) and (40), and bringing Equation (40) into the first derivative (38) of the sliding mode function, we can obtain the following:
Based on the principle of sliding mode controller, the sliding mode control input has the following form:
where
is the equivalent control items of the system, and
is the nonlinear switching control items of the system.
Make Equation (41) approach zero, to get the following:
Make
, since
is a smooth resistance and has a boundary; then its first and second derivatives are bounded, and the basic running resistance of the train
is bounded, so there is
, and
is the boundary. According to the design principle of sliding mode control, the nonlinear switching term is designed as
, and then the terminal sliding mode control input is as follows:
where
is the control gain.
In addition, since the resistance items of the train are unknown in the actual operation of the train, that is
is unknown, but
is bounded, the actual controller should be changed to the following:
where
;
;
; and
.
Select the Lyapunov function of the following form:
Let
, derivate Equation (46), and bring the controller formula, Equation (45), into it:
The system is stable, but in order to reduce the control burden of switching items and reduce the possibility of chattering, a parameter adaptation mechanism is introduced.
First modify the controller to the following form:
where
are the estimated values of parameters
, respectively.
Construct the following Lyapunov function:
where
Bring Equation (48) into Equation (49), and we can get the following:
To meet the stability rules of the Lyapunov formula, organize Equation (50), and the available parameter adaptation law is as follows:
As showed in Equations (50) and (51), discontinuous switching function , will result in unconnected control signals to keep the train running on the sliding surface. However, due to resistance deviation, measurement deviation, and other factors, the train cannot stay on the sliding surface, which will eventually cause chattering. To avoid chattering, replace signum function , in Equation (47) with saturation function, and is the saturation width.
Generally speaking, the distinguishing feature of typical sliding mode control is the use of a signum function. After the introduction of a saturation function, the controller will no longer be a classic sliding mode controller. According to References [
41,
42], the sat function implies a pseudo sliding mode in fact.
5.3. Disturbance Observer
The introduction of the saturation function will reduce the robustness of the system, but the line interference may be greater during train operation, especially on ramps. Disturbance observer (DOB) has been one of the most widely used robust control tools which can effectively improve robustness and anti-disturbance [
43,
44]. Therefore, a disturbance observer is introduced in this part.
The Equation (48) is introduced into the disturbance observation value and transformed into the following form:
As shown in
Figure 13, the actual model and the nominal model are as follows:
According to the principle of disturbance observer, by considering the input term as the transfer function of
and
, we can get the following:
Furthermore, the expression of the disturbance observation can be obtained as follows:
According to Equation (56), if
, the required disturbance value can be accurately observed, but in the actual system, there may be a deviation between
and
. From Equation (56), the deviation of the estimated value mainly comes from the error of the model:
where
,
is the low pass filter, and
.
According to Equation (57), if there is an error in the model, the obtained observations will also have a certain deviation. The magnitude of the model error determines the magnitude of the deviation of the observations. After introducing the disturbance observer, the gain of the switching function will be weakened to a certain extent. In general, the error of the train model is small, so the disturbance observer can improve the anti-interference ability of the system and reduce the switching gain, to avoid chattering.