*Article* **Model Predictive Controller Based on Online Obtaining of Softness Factor and Fusion Velocity for Automatic Train Operation**

#### **Longda Wang <sup>1</sup> , Xingcheng Wang 1,\*, Zhao Sheng <sup>2</sup> and Senkui Lu <sup>1</sup>**


Received: 8 February 2020; Accepted: 16 March 2020; Published: 19 March 2020

**Abstract:** This paper develops an improved model predictive controller based on the online obtaining of softness factor and fusion velocity for automatic train operation to enhance the tracking control performance. Specifically, the softness factor of the improved model predictive control algorithm is not a constant, conversely, an improved online adaptive adjusting method for softness factor based on fuzzy satisfaction of system output value and velocity distance trajectory characteristic is adopted, and an improved whale optimization algorithm has been proposed to solve the adjustable parameters; meanwhile, the system output value for automatic train operation is not sampled by a normal speed sensor, on the contrary, an improved online velocity sampled method for the system output value based on a fusion velocity model and an intelligent digital torque sensor is applied. In addition, the two improved strategies proposed take the real-time storage and calculation capacities of the core chip of the controller into account. Therefore, the proposed improved strategies (I) have good performance in tracking precision, (II) are simple and easily conducted, and (III) can ensure the accomplishing of computational tasks in real-time. Finally, to verify the effectiveness of the improved model predictive controller, the Matlab/simulink simulation and hardware-in-the-loop simulation (HILS) are adopted for automatic train operation tracking control, and the tracking control simulation results indicate that the improved model predictive controller has better tracking control effectiveness compared with the existing traditional improved model predictive controller.

**Keywords:** model predictive controller; automatic train operation; softness factor; fusion velocity; online obtaining; hardware-in-the-loop simulation

#### **1. Introduction**

The urban rail transit system with automatic train operation system has the advantages of safety, stability, economy, and comfort, and it has become one of the most popular and efficient means of the urban public transportation [1]. The tracking control functional module makes the velocity trajectory track at the optimal target speed obtained by the upper-layer optimal loop, and according to the appropriate and efficient tracking control algorithm, it is an indispensable crucial system and necessary to ensure optimal safety, comfort, energy-efficiency, punctuality, and parking accuracy for train operation process, which requires the corresponding algorithm to possess good control performance [2]. Therefore, aiming at improving the multi-objective performance index of the train operation process, an automatic train operation system has been developed rapidly and is widely applied in urban rail trains operation [3–5]. Meanwhile, various improved intelligent optimization control algorithms have been proposed and applied for the automatic train operation system [6–8].

In recent years, many improved algorithms have been applied in the automatic train operation tracking control field, such as robust adaptive automatic control, model predictive control, online learning, iterative learning control, matter-element model, etc. [9–13]. An online approximation-based robust adaptive automatic train control method is proposed for the automatic train operation (ATO) system [9]. A fuzzy model predictive control approach is proposed to provide locomotive operation instructions for mainline railways continuously, and extensive simulations show that the proposed approach can provide sufficient solution optimality in reasonable computational time and energy consumption in train operations is reduced [10]. A novel online learning control strategy is proposed to solve the train automatic stop control (TASC) problem [11]. An iterative learning control based on automatic train operation is proposed to deal with the trajectory tracking control problem under certain velocity constrains [12]. Matter-element theory is applied to the established models to optimize speed trajectory for achieving multi-objective optimization, and the relative performance indices weighting is determined in different stages so that the more satisfied decision speeds could be calculated with the goodness evaluation method [13]. The above research can improve the tracking control performance of the traditional control algorithms.

Among numerous algorithms, Model Predictive Control (MPC) is one of the most effective control algorithms, which is characterized by good robustness, fast tracking speed, accurate tracking target speed, etc. [14]. A linear time-varying MPC is used to obtain the power split between the combustion engine and electrical machines and the system operating points at each sample time [15]. A coordinated energy dispatch based on Distributed model predictive control (DMPC) is proposed, and the corresponding simulation results show the effectiveness of the proposed method [16]. A co-design of the self-triggered mechanism and distributed model predictive control (DMPC) is proposed to achieve the cooperative objectives while efficiently exploiting communication network [17]. A model predictive control-based droop current regulator to interface PV in smart dc distribution systems is proposed [18]. From the various model predictive control algorithms, the dynamic matrix control model predictive control (DMC MPC) is an effective algorithm among them due to its characteristics of strong robustness, fast tracking speed, high precision for tracking control, avoiding the parameter identification for the transfer function model, and solving the problem of delay process effectively. A new method that linearizes the RC equivalent circuit model and predicts available battery power according to original Dynamic Matrix Control algorithm is proposed [19]. An application of dynamic matrix control (DMC) to a drum-type boiler–turbine system is proposed [20]. Of particular note is the use of the improved DMC MPC for automatic train operation tracking control scenario [21].

It is necessary to conduct further study on the basis of previous research findings, and key parameters adjusting and improving the sampling accuracy should be taken seriously. A method for calculating the traction characteristics of a traction motor is proposed [22]. A new method to identify the train key design variables against the running performance indicators based on the sensitivity analysis is proposed, which in turn bases itself on simulation-oriented surrogate models [23]. A novel adaptive sampling algorithm for power management in the automated monitoring of the quality of water in an environment is devised and applied [24].

Traditional simulations based on a pure software environment cannot truly reflect the actual automatic train operation process, and the situation representing the actual automatic train operation experiment is difficult to implement because it is expensive, has restricted experimental conditions, high construction difficulties, and high security protection requirements. Hardware-in-the-loop simulation (HILS) is a new simulation technology for solving this difficult issue [25,26]. At present, numerous relative research findings have achieved improvements in the traction control system [27,28].

An improved model predictive controller based on online obtaining of softness factor and fusion velocity for automatic train operation is developed. The following summarizes the main contributions of this paper.


The paper is organized as follows. Section 2 introduces the model predictive controller for automatic train operation tracking control. Section 3 illustrates the improved DMC model predictive controller based on online obtaining of softness factor and fusion velocity developed in this paper. Section 4 provides the Matlab/simulation results and hardware-in-the-loop simulation (HILS) results to illustrate the proposed method. Section 5 concludes this article.

#### **2. Model Predictive Controller for Automatic Train Operation Tracking Control**

#### *2.1. Evaluation Index for Automatic Train Operation Tracking Control*

The integral of time multiplied by the absolute value of error (ITAE) is the frequently used evaluation index for tracking control performance [29]. The specific formula for the evaluation index ITAE is as follows,

$$ITAE = \int t \left| \mathfrak{e} \left( t \right) \right| dt \tag{1}$$

where *t* represents the sample time of control process, and |*e*(*t*)| represents the absolute value of error between target speed and actual tracking control speed.

As automatic train operation has its own unique characteristics and requirements, the multi-objective performance index is more appropriative, and it used universally. The computation model of multi-objective performance index *Pk* is as follows,

$$\begin{cases} \begin{aligned} P\_k &= \sum\_{i=1}^4 \omega\_i \times \frac{f\_i - \min(f\_i)}{\max(f\_i) - \min(f\_i)} \times f\_i \\ & (f\_1, f\_2, f\_3, f\_4) = (\Delta s\_\lambda \Delta t, K\_{\text{Jerk}} E) \\ \Delta v \frac{dv}{ds} &= F\left(u, v\right) - R\left(v, s\right) - B\left(u, v\right) \\ \frac{dt}{ds} &= \frac{1}{v} \\ v\left(s\right) &\leq v\_{\text{lim}}\left(s\right) \\ R\left(v, s\right) &= r(v) + R\_I\left(s\right) \\ \Delta s = \left|s\_z - D\right| &< \Delta s\_{\text{max}} \\ \Delta t = \left|T - T\_r\right| &< \Delta t\_{\text{max}} \\ K\_{I\text{r}k} &= \int \left|\Delta a\right| d\_5 \left|\right| \\ E &= \int \left(Ma - R\right) d\_5 \end{aligned} \end{cases} \tag{2}$$

where *ω<sup>i</sup>* represents the index importance weight factor ( *k* ∑ *i*=1 *ω <sup>i</sup>* = 1), which reflects the relative importance of the *i* th optimization index; *t* represents the actual running time of the train; *s* represents the actual position of the train; *a* represents the actual acceleration of the train; |Δ*a*| represents the actual impingement rate of the train; *M* is the mass of the train; *Ft* (*u*, *v*) and *Br* (*u*, *v*) are the traction force and braking force of the current velocity, respectively; *R* (*v*,*s*) is the resistance of the train determined by the current speed and line position; *sz* is the terminal position; *Tr* is the actual running time; *D* is the actual running distance; *v* (*s*) represents the instantaneous velocity in the position *s*; *T*¯ represents the prospective running time; *vlim* (*s*) represents the upper limit velocity in the position *s*; Δ*s*max represents the allowable maximum parking error; Δ*t*max represents the allowable maximum time error; Δ*s* and Δ*t* represent the actual parking error and time error, respectively; *u* represents the train control quantity; *KJerk* represents the comfort performance index; and *E* represents the energy consumption during the train operation process [2,30].

In addition, security index should be taken into account as well. Traveling over the velocity limit is the main risk and non-negligible factors can cause an unsafe environment. The computation formula of security index *Ksaf e* is as follows,

$$\begin{aligned} \min & K\_{safc} \\ K\_{safc} &= \sum\_{is=1}^{\text{su}} \text{YS} \left( \text{is} \right) \Big/ \text{sn} \\ \text{YS} \left( \text{is} \right) &= \begin{cases} 0 & v \left( \text{is} \right) > v\_{\text{lim}} \left( \text{is} \right) \\ 1 & v \left( \text{is} \right) \le v\_{\text{lim}} \left( \text{is} \right) \end{cases} \end{aligned} \tag{3}$$

where *is* represents the index of sampling point, *YS* (*is*) represents the security evaluation value of the *is*-th sampling point, and *sn* represents the number of sampling points [1].

#### *2.2. Conventional Dynamic Matrix Control Model Predictive Control*

DMC MPC uses three methods, including the DMC predictive model, rolling optimization, and feedback correction, to control the controlled object [31].

#### 2.2.1. DMC Predictive Model

The DMC predictive model is one of significant models for DMC MPC. The unit step response model reflecting the dynamic performance is adopted as the DMC predictive model for controlled object, and the predictive value of system output is obtained by the step response characteristic for controlled object.

If the model length is *N*, then the *N* sampled values of the controlled object unit step response can be used to describe the dynamic response characteristics of the system. The specific calculation formula for the predictive value of system output is as follows,

$$Y\_p\left(k\right) = \mathcal{Y}\_0\left(k\right) + A\Delta l I\left(k\right) \tag{4}$$

where *Yp* (*k*) = *yp* (*k* + 1|*k*), *yp* (*k* + 2|*k*), ..., *yp* (*k* + *N*|*k*) *<sup>T</sup>* represents the predictive value of system output, *<sup>Y</sup>*<sup>0</sup> (*k*) <sup>=</sup> [*y*<sup>0</sup> (*<sup>k</sup>* <sup>+</sup> <sup>1</sup>|*k*), *<sup>y</sup>*<sup>0</sup> (*<sup>k</sup>* <sup>+</sup> <sup>2</sup>|*k*), ..., *<sup>y</sup>*<sup>0</sup> (*<sup>k</sup>* <sup>+</sup> *<sup>N</sup>*|*k*)]*<sup>T</sup>* represents the predictive value of predictive model, <sup>Δ</sup>*<sup>U</sup>* (*k*) <sup>=</sup> [Δ*<sup>u</sup>* (*k*), <sup>Δ</sup>*<sup>u</sup>* (*<sup>k</sup>* <sup>+</sup> <sup>1</sup>), ..., <sup>Δ</sup>*<sup>u</sup>* (*<sup>k</sup>* <sup>+</sup> *<sup>N</sup>* <sup>−</sup> <sup>1</sup>)]*<sup>T</sup>* represents the incremental sequence for control, and *A* represents the dynamic matrix. The specific dynamic matrix *A* and the specific calculation formula for the element of *Yp* are as follows,

$$A = \begin{bmatrix} a\_1 & 0 & 0 & \cdots & 0 \\ a\_2 & a\_1 & 0 & \cdots & 0 \\ & \vdots & \vdots & \ddots & \vdots \\ a\_N & a\_{N-1} & \cdots & a\_1 \end{bmatrix} \tag{5}$$

$$y\_p(k+i|k) = y\_0(k+i|k) + \sum\_{j=1}^{i} a\_{i-j+1} \Delta u \,(k+j-1) \tag{6}$$

where *i* represents the element index of *Yp*, *i* ∈ {1, 2, ..., *N*}; *k* represents the initial point of DMC predictive model [31,32].

#### 2.2.2. Rolling Optimization

With the aim of avoiding the violent fluctuations in the control process effectively, it is necessary to make the final output value *yf* reach to the reference target value *yr* along the predetermined smooth path by the DMC MPC system, so as to enhance the system robustness. Thus, a popular reference path used in DMC MPC is as follows,

$$y\_f\left(k+i\right) = a^i y\left(k\right) + \left(1 - a^i\right) y\_r\tag{7}$$

where *yf* (*k* + *i*) represents the final output value expected, *α<sup>i</sup>* represents the *i*th softness factor (0 < *α<sup>i</sup>* < 1), *y* (*k*) represents the actual output value of the system, and *yr* represents the reference target value of the system.

The quadratic rolling optimization object of the system is necessary for rolling optimization. If the predictive length is *M* and control length is *L*, in general, *L* ≤ *M* ≤ *N*. The specific quadratic rolling optimization object of system is as follows,

$$\begin{aligned} J &= \left\| Y\_p \left( k \right) - Y\_f \left( k \right) \right\|\_{Q}^{2} + \left\| \Delta \mathcal{U}\_L \left( k \right) \right\|\_{R}^{2} \\ &= \sum\_{i=1}^{M} q\_i \left[ y\_p \left( k + i \middle| k \right) - y\_f \left( k + i \right) \right] + \sum\_{i=1}^{L} r\_i \Delta u \left( k + i - 1 \right) \end{aligned} \tag{8}$$

where *Yf* (*k*) = *yf* (*k* + 1), *yf* (*k* + 2), ..., *yf* (*k* + *M*) *T* represents the control sequence of the system, *R* = *diag*[*r*1,*r*2, ...,*rL*] *<sup>T</sup>* represents the weight coefficient matrix of constraint for error revise, *Q* = *diag*[*q*1, *q*2, ..., *qM*] *<sup>T</sup>* represents the weight coefficient matrix of constraint for error increment revise, and *diag* represents the diagonal matrix.

The necessary condition for obtaining the minimum value of objective function *J* is *<sup>∂</sup><sup>J</sup> <sup>∂</sup>*Δ*UL*(*k*) <sup>=</sup> <sup>0</sup> through extreme value theory under unconstrained conditions. Therefore, the control sequence optimal solution can be obtained by rolling optimization. The specific calculation formula of the control sequence optimal solution is as follows.

$$
\Delta l I\_L \left(k - 1\right) = \left(A^T Q A + R\right)^{-1} A^T Q \left[\Upsilon\_f\left(k\right) - \Upsilon\_0\left(k\right)\right] \tag{9}
$$

Then, the actual control quantity *u*(*k*) can be obtained. The specific calculation formula of the actual control quantity *u*(*k*) is as follows,

$$
\mu\left(k\right) = \mu\left(k-1\right) + \Delta\mu\left(k-1\right) \tag{10}
$$

In the next control period, i.e., the *k* + 1 th control period, the corresponding Δ*u*(*k*) and *u*(*k* + 1) can be obtained by the above way. Thus, it can realize the rolling optimization of the actual control quantity in the iterative control process [21,33].

#### 2.2.3. Feedback Correction

Feedback correction is an important component of DMC MPC; it is used to reduce the influence of system disturbance for the control system, so as to achieve the ideal control effectiveness. The specific calculation formula of the error between actual system output value and the predicted output value in the present control period (the *k* th control period) is as follows.

$$
\epsilon \left( k \right) = y \left( k \right) - y\_{\mathcal{P}} \left( k \right) \tag{11}
$$

After feedback correction calculation, the predicted output value can be corrected to certain extent [21,33,34]. The specific corrected calculation formula is as follows,

$$Y\_{p2}\left(k\right) = Y\_0\left(k\right) + A\Delta l I\left(k-1\right) + HC\tag{12}$$

where *C* represents the error corrected matrix, and its length is *N*; *H* represents the corresponding transformed matrix.

#### *2.3. Fuzzy DMC Model Predictive Controller for Automatic Train Operation*

Aiming at improving the precision of the automatic train operation tracking control for DMC MPC, using the fuzzy model prediction based on the train operation mechanism is a good choice. The slope and velocity error are the most important train operation information. For example, when the train runs in steep uphill and current velocity is far less than target velocity, the conversion degree for train operation is "PB", that is, the maximum extent to drew train is adopted to assist the climb by accelerating or keeping the velocity. In this time, if the maximum traction is used according to the intrinsic DMC MPC, the addition traction incremental quantity is not necessary; otherwise, the appropriate addition traction incremental quantity should be used to correct this error. The fuzzy sets are divided into ['N4',.......,'Z',.......,'P4'] [10,35]. The specific calculation model for fuzzy model prediction is as follows,

$$\begin{array}{l} \mu\_{f\\_p}(k) = \mathbb{C}\_{fuxzy1}\left(\omega\left(k\right), \epsilon\left(k\right)\right) \\ \Delta\mu\_{f\\_p}(k) = \mathbb{C}\_{fuxzy2}\left(\mu\_{f\\_p}(k), \mu\left(k\right)\right) \\ \mu\_c\left(k\right) = \mu\left(k\right) + \Delta\mu\_{f\\_p}\left(k\right) \end{array} \tag{13}$$

where *Cf uzzy*<sup>1</sup> and *Cf uzzy*<sup>2</sup> represent the fuzzy inference functions by using two kinds of fuzzy rules, respectively; *uf* \_*<sup>p</sup>* (*k*) represents the calculated control quantity by using fuzzy rule about slope and velocity error; Δ*uf* \_*<sup>p</sup>* (*k*) represents the calculated control quantity by using fuzzy rule about control quantity calculated by intrinsic DMC MPC and control quantity calculated by fuzzy logic and train operation information; and *uc* (*k*) represents the final calculated control quantity for the automatic train operation tracking control.

The fuzzy rules for fuzzy model prediction and partial membership function for control quantity are shown in Figure 1.

**Figure 1.** The fuzzy rules for fuzzy model prediction and partial membership function for control quantity. (**a**) Fuzzy rule for train operation information. (**b**) Fuzzy rule for control quantity prediction. (**c**) Partial membership functions for control quantity.

Fuzzy dynamic matrix control model predictive control (Fuzzy DMC MPC) is a control method that considers the step response characteristics and fuzzy logic for the train operation mechanism of the control object. The fuzzy DMC model predictive controller is widely used in automatic train operation due to its characteristics of simple design scheme and high tracking precision. The fuzzy DMC model prediction controller is mainly composed of four function chips (Fuzzy model prediction function chip, DMC model prediction function chip, rolling optimization function chip, and feedback correction function chip), and it is used to realize four function modules of Fuzzy DMC MPC (Fuzzy model prediction, DMC model prediction, rolling optimization, and feedback correction). The schematic diagram of the Fuzzy DMC model predictive controller for automatic train operation is shown in Figure 2.

**Figure 2.** Schematic diagram of the Fuzzy DMC model predictive controller for automatic train operation.

As can be seen from Figure 3, it is impossible to obtain the actual output value (real-time velocity) for automatic train operation tracking control system. In addition, as can be seen from Formula (7), the real-time softness factor is also an important factor that impacted the multi-objective performance index for automatic train operation tracking control. Therefore, it is necessary to improve the real-time velocity sampling accuracy and softness factor accuracy for automatic train operation tracking control as much as possible.

#### **3. Model Predictive Controller Based on Online Obtaining of Softness Factor and Fusion Velocity**

#### *3.1. Fusion Velocity Computation Model and Corrected Model Based on Online Obtaining*

#### 3.1.1. Fusion Velocity Computation Model Based on Online Obtaining

According to the multi-objective performance index for automatic train operation tracking control, the fusion velocity model based on online obtaining is necessary to take into account the energy consumption, running time, comfort, and parking accuracy. Taking into account the sampling effect, hardware technology (storage and computing ability), funds, space, and other factors, three kinds of velocity sampling sources are sufficient (motor speed, motor torque and train instantaneous displacement) and are selected and synthesized. The fusion velocity computation model based on online obtaining is established as follows,

$$\begin{cases} \boldsymbol{\upsilon}\_{\rm is,\upsilon} = \boldsymbol{n}\_{\rm is} \times \boldsymbol{t} \boldsymbol{r}\_{\rm ntr} \times \eta\_{\rm g} \times \eta\_{\rm is,T} \times \eta\_{\rm is,ict} \\\ \boldsymbol{\upsilon}\_{\rm is,F} = \boldsymbol{\upsilon}\_{\rm is-1,a} + \boldsymbol{F}\_{\rm is} - \boldsymbol{w}\_{\rm is} / \boldsymbol{M} \cdot \Delta t \\\ \boldsymbol{\upsilon}\_{\rm is,s} = \Delta \boldsymbol{s} / \Delta t = \boldsymbol{s}\_{\rm is} - \boldsymbol{s}\_{\rm is-1} / \Delta t \end{cases} \tag{14}$$

$$
\upsilon\_{\rm is,a} = \lambda\_{\rm is,v} \times \upsilon\_{\rm is,v} + \lambda\_{\rm is,F} \times \upsilon\_{\rm is,F} + \lambda\_{\rm is,s} \times \upsilon\_{\rm is,s} \tag{15}
$$

where *vis*,*<sup>a</sup>* represents the final calculated velocity by speed analyzer ultimately of the *i*-th sampling point; *vis*,*<sup>v</sup>* represents the velocity calculated based on actual motor speed sampled of the *i*-th sampling point; *vis*,*<sup>F</sup>* represents the velocity calculated based on actual motor torque sampled of the *i*-th sampling point; *vis*,*<sup>s</sup>* represents the velocity calculated based on actual train instantaneous displacement sampled of the *i*-th sampling point; *nis* represents the actual motor speed sampled by speed sensors of the *i*-th sampling point; *trntv* represents the transmission ratio of the motor speed to train velocity; *η<sup>g</sup>* represents the degree of tooth engagement between the gears; *ηis*,*<sup>T</sup>* represents the speed transmission efficiency of the *i*-th sampling point; *ηis*,*itc* represents the efficiency for the train to overcome idling, taxiing, and creep sliding of the *i*-th sampling point; *Fis* represents the force calculated based on actual motor torque sampled of the *<sup>i</sup>*-th sampling point; *Fis* <sup>=</sup> *<sup>η</sup>is*,*<sup>F</sup>* <sup>×</sup> *Tis*/*Rmr*; *Tis* represents the actual motor torque sampled by torquemeter of the *i*-th sampling point; *ηis*,*<sup>F</sup>* represents the force comprehensive transmission efficiency of the *i*-th sampling point; *Rmr* represents the radius of motor rotor; *wis*(*v*,*s*) represents the actual resistance of the *i*-th sampling point (*v*,*s*); Δ*t* represents the sampling time-interval, Δ*t* is 500 μs in this paper; Δ*s* represents the sampling displacement-interval; *sis* represents the actual train instantaneous displacement sampled by displacement pickup of the *i*-th sampling point; and *λis* = {*λis*,*v*, *λis*,*F*, *λis*,*s*} represents the synthetic weight of the velocity sampled by different ways.

Synthetic weight is vital for real-time sampling precision in the tracking control process, the importance of each velocity sampling sources needs to be considered, so as to give the appropriate synthetic weight. The synthetic weight indicates the importance of the real-time velocity obtained by different speed sampling sources. Yet, the selection of the synthetic weight by traditional methods lacks the specific theoretical basis, so there is certain subjective limitation in actual applied. As automatic train operation tracking control is an extremely complex issue, there is a trajectory characteristic for automatic train operation tracking control curve dominated by velocity target curve, train parameters, line conditions and running requirements, and the traditional methods for setting synthetic weight based on experience empower is subjective and blind, so it is necessary to be improved. In this paper, an synthetic weight empower using entropy method is applied for automatic train operation tracking control. First, the whole tracking control curve is divided by a position according to trajectory characteristic and line conditions; second, a large number of real-time data, including velocity, force, and position information for the whole tracking control process, should be sampled to prepare for calculation; finally, the entropy method is used to calculate the synthetic weight of each divided subinterval of tracking control curve.

Entropy is a measure of uncertainty for information calculation. The entropy weight method is utilizes the entropy characteristics and assigns a weight to each index in an event by calculating the entropy value. The entropy weight method is an objective weight empower method, because it simply depends on the discreteness of data itself. The specific steps of computational process for entropy weight method are as follows.

A certain number of samples (as many as possible) must be collected to prepare for the calculation, and their index values also needed to be recorded.

To eliminate the negative influences caused by the difference between units and magnitude orders, the index values must be normalized. The calculation formulas for the normalization can be expressed as follows,

$$\mathbf{x}\_{ij}^{\prime} = \frac{\mathbf{x}\_{ij} - \min\left\{\mathbf{x}\_{ij\prime}, \mathbf{x}\_{2j\prime}, \dots, \mathbf{x}\_{nj}\right\}}{\max\left\{\mathbf{x}\_{ij\prime}, \mathbf{x}\_{2j\prime}, \dots, \mathbf{x}\_{nj\prime}\right\} - \min\left\{\mathbf{x}\_{ij\prime}, \mathbf{x}\_{2j\prime}, \dots, \mathbf{x}\_{nj\prime}\right\}} \tag{16}$$

$$\mathbf{x}\_{ij}^{'} = \frac{\max\left\{\mathbf{x}\_{i\bar{j}\prime}, \mathbf{x}\_{2\bar{j}\prime}, \dots, \mathbf{x}\_{n\bar{j}}\right\} - \mathbf{x}\_{i\bar{j}}}{\max\left\{\mathbf{x}\_{i\bar{j}\prime}, \mathbf{x}\_{2\bar{j}\prime}, \dots, \mathbf{x}\_{n\bar{j}}\right\} - \min\left\{\mathbf{x}\_{i\bar{j}\prime}, \mathbf{x}\_{2\bar{j}\prime}, \dots, \mathbf{x}\_{n\bar{j}}\right\}} \tag{17}$$

where *j* = (1, 2, ···, *m*),*i* = (1, 2, ···, *n*); *m* represents index number; *n* represents the number of samples; *xij* represents the *j*-th processed index value of the *i*-th sample after normalization; *xij* represents the *j*-th original index value of the *i*-th sample before normalization; *max* and *min*, respectively, represent the maximum and minimum values of the array. If index value *xij* is a positive number, Formula (16) is used to normalize; otherwise, Formula (17) is used to normalize.

The normalized index values are necessary to filtered out the zero value further, so as to avoid illegal logarithmic function (ln(0)) in next subsequent calculation processes. The specific formula for filtered out zero value is as follows,

$$
\Delta \mathbf{x}\_{\text{ij}}^{\prime\prime} = \Delta z + (1 - 2 \times \Delta z) \times \mathbf{x}\_{\text{ij}}^{\prime} \tag{18}
$$

where *xij* represents the *j*-th processed index value of *i*-th sample after filtering out the zero value; Δ*z* represents the tiny value reasonable; Δ*z* is 0.01 in this paper.

Then, the entropy values of each index values are necessary to be calculated. The specific formulas for calculating the entropy values are as follows,

$$p\_{ij} = \frac{\mathbf{x}\_{ij}^{\prime\prime}}{\sum\_{i=1}^{n} \mathbf{x}\_{ij}^{\prime\prime}} \tag{19}$$

$$\mathcal{e}\_{\vec{j}} = -k \times \sum\_{i=1}^{n} \left( p\_{\vec{i}\vec{j}} \times \ln \left( p\_{\vec{i}\vec{j}} \right) \right) \tag{20}$$

where *pij* represents the *j*-th index weight value of *i*-th sample in *j*-th index; *ej* represents the *j*-th index entropy value; *k* represents the entropy coefficient, the value is the reciprocal of ln (*n*), *k* = 1 ) ln (*n*).

Finally, the index weight value could be calculated. The specific formula for calculating the index weight values is as follows,

$$d\_j = 1 - \varepsilon\_j \tag{21}$$

$$
\lambda\_{\bar{j}} = \frac{d\_{\bar{j}}}{\sum\_{j=1}^{m} d\_{\bar{j}}} \tag{22}
$$

where *dj* represents the *j*-th entropy redundancy value of *j*-th index, it indicates the difference degree of this index; *ej* represents the *j*-th index entropy value; *λ<sup>j</sup>* represents the *j*-th weight value.

#### 3.1.2. Fusion Velocity Corrected Model Based on Online Obtaining

If the factor of the velocity sampling sources sampled inaccurately is not considered, the improvement effect for automatic train operation tracking control will inevitably be restricted. To improve the real-time sampling velocity precision, a corrected method of real-time sampling velocity for automatic train operation tracking control using auxiliary corrective velocity sampling source is popularly applied in various types of urban rail vehicles. The specific evaluation and corrected formulas are as follows,

$$
\Delta v\_{\rm is,c} = \left| v\_{\rm is,x} - v\_{\rm is,ref} \right| \le \Delta v\_{\rm is,p} \tag{23}
$$

$$
\Delta v\_{\rm is,c} = \frac{\Delta v\_{\rm is,p}}{\Delta v\_{\rm is,c}} \times v\_{\rm is,x} + \frac{\Delta v\_{\rm is,c} - \Delta v\_{\rm is,p}}{\Delta v\_{\rm is,c}} \times v\_{\rm is,ref} \tag{24}
$$

where *vis*,*<sup>x</sup>* represents a velocity sampling source *x* from *vis*,*v*, *vis*,*F*, *vis*,*<sup>s</sup>* using to synthetic final calculated velocity *vis*,*a*; *vis*,*ref* represents the reference velocity (auxiliary corrective velocity sampling source) based on actual sampling data sampled by specific auxiliary sensor; Δ*vis*,*<sup>c</sup>* represents the actual error value between *vis*,*<sup>x</sup>* and *vis*,*ref* ; Δ*vis*,*<sup>p</sup>* represents the maximum permit satisfied error value between *vis*,*<sup>x</sup>* and *vis*,*ref* ; *vis*,*<sup>c</sup>* represents the final corrected value calculated by Formula (24) when correctness condition (Formula (23)) is not satisfied.

Reference velocity *vis*,*ref* will exert a measure of influence over the velocity-corrected effect. Thus, the choice of auxiliary corrective velocity sampling source is significant. The specific gear speed on the vehicle wheel side of the gear box is a good choice.

#### *3.2. Softness Factor Adaptive Adjusting Model Based on Online Obtaining*

The softness factor is a key parameter for DMC MPC; it plays plays an important role in balancing the degree of robustness and rapidity for the DMC MPC tracking control system. If softness factor *α* is chosen as a larger value, the system will have slower response speed and stronger robustness; by contrast, if softness factor *α* is chosen as a smaller value, the system will have faster response speed and worse robustness [36]. Thus, both response speed and robustness must be taken into account for softness factor *α* setting.

Considering the trajectory characteristic and tracking control condition for automatic train operation, the softness factor adaptive adjusting model based on online obtaining is established as follows,

$$
\mathfrak{a} = \lambda\_{\mathfrak{a}Ts} \times \mathfrak{a}\_{Ts} \begin{pmatrix} \mathfrak{s} \end{pmatrix} + \lambda\_{\mathfrak{a}\_{\mathfrak{y}\mathfrak{y}}} \times \mathfrak{a}\_{\mathfrak{y}\mathfrak{y}} \begin{pmatrix} \mathfrak{y} \begin{pmatrix} k \end{pmatrix} \end{pmatrix} , \mathfrak{y}\_r \tag{25}
$$

where *α* represents the final calculated real-time softness factor; *αTs* represents the real-time softness factor calculated based on the trajectory characteristic of the present position; *αμ<sup>y</sup>* (*y* (*k*), *yr*) represents the real-time softness factor calculated based on tracking control condition of the present system output *y* (*k*); *λαTs* and *λαμ<sup>y</sup>* represent the fusion weights of *αTs* and *αμ<sup>y</sup>* (*y* (*k*), *yr*), respectively; and *λαTs* + *λαμ<sup>y</sup>* = 1.

The whole tracking control curve must be divided into several different types of subintervals by position according to the trajectory characteristic and line conditions. The specific types of subintervals are described as follows.

Type 1: The vibrating area nearby inflection point of tracking control curve.

In this area, there is the strong velocity fluctuation in the velocity trajectory. Thus, aiming at improving the system robustness as much as possible, softness factor *αTs* should be an appropriate larger value at the cost of reduce acceptable system rapidity.

Type 2: The smooth area of tracking control curve.

In this area, there is no obvious velocity fluctuation in the velocity trajectory. Thus, aiming at improving the system rapidity as much as possible, softness factor *αTs* should be chosen as an appropriate smaller value at the cost of reduce acceptable system robustness.

Type 3: The connected area in the middle of smooth area and the vibrating area of the tracking control curve.

In this area, the system rapidity and rapidity are taken into account for softness factor *αTs* setting. Thus, softness factor *αTs* should be choose a appropriate intermediate value.

In addition, although in the same type of subintervals, the softening factor *αTs* almost varies because of the different intensity degrees of velocity fluctuation. The specific calculation formula for softness factor *αTs* based on trajectory characteristic of the present position is described as follows,

$$\mathbf{a}\_{Ts}\left(\mathbf{s}\right) = \begin{cases} \mathbf{a}\_{\mathrm{Tr},si} + \left(\mathbf{a}\_{\mathrm{Tr},si} - \mathbf{a}\_{\mathrm{Tr},si-1}\right) \times \frac{\left(\left(\mathbf{S}\_{si} + \mathbf{S}\_{1}\right) - s\right)}{\mathbf{S}\_{1} + \mathbf{S}\_{2}} & s < \mathbf{S}\_{si} + \mathbf{S}\_{1} \\\mathbf{a}\_{\mathrm{Tr},si} & \mathbf{S}\_{si} + \mathbf{S}\_{1} \le s \le \mathbf{S}\_{si+1} - \mathbf{S}\_{2} \\\mathbf{a}\_{\mathrm{Tr},si} + \left(\mathbf{a}\_{\mathrm{Tr},si+1} - \mathbf{a}\_{\mathrm{Tr},si}\right) \times \frac{\left(s - \left(\mathbf{S}\_{si} - \mathbf{S}\_{2}\right)\right)}{\mathbf{S}\_{1} + \mathbf{S}\_{2}} & s > \mathbf{S}\_{si+1} - \mathbf{S}\_{2} \end{cases} \tag{26}$$

where *si* represents the subinterval index *si* ∈ {1, 2, ...,*simax*}; *si*\_*max* represents the number of subintervals; *Ssi* represents the starting position of the *si*-th subinterval; *Ssi*\_*max*+<sup>1</sup> represents the terminal position of tracking control curve, it is a target parking position; *αTr*,*si* represents the reference value of soften factor in the *si*-th subinterval; *αTs*,0 = *αTs*,1, *αTs*,*si*\_ max <sup>+</sup><sup>1</sup> = *αTs*,*si*\_ max; *S*<sup>1</sup> and *S*<sup>2</sup> represent the connected length, in the connected area; the softness factor *αTs* is reduced or increased linearly and smoothly, so as to avoid the instability of tracking control system.

Aiming at solving this control problem with fuzzy characteristic, an fuzzy adaptive adjusting method for online obtaining softness factor *αμ<sup>y</sup>* is applied. First of all, the satisfaction degree of control is defined, so as to the automatic train operation tracking control problem can be transformed into an optimization decision-making problem by fuzzy reasoning; then, the corresponding real-time parameters of the controller are adjusted online to meet the requirements of the system control quality, so as to achieve the purpose of system optimization control. The specific calculation formula for fuzzy satisfaction degree *μy*(*k*) of system output *y*(*k*) is as follows,

$$\mu\_{y(k)} = \begin{cases} 0 & y \left( k \right) < y\_{\text{min}} - s\_1 \\ 1 + \frac{y(k) - y\_{\text{min}}}{s\_1} & y\_{\text{min}} - s\_1 \le y \left( k \right) < y\_{\text{min}} \\ 1 & y\_{\text{min}} \le y \left( k \right) < y\_{\text{max}} \\ 1 + \frac{y(k) - y\_{\text{max}}}{s\_2} & y\_{\text{max}} \le y \left( k \right) < y\_{\text{max}} + s\_2 \\ 0 & y\_{\text{min}} - s\_1 \ge y \left( k \right) \end{cases} \tag{27}$$

where *s*<sup>1</sup> and *s*<sup>2</sup> represent the blur width, which can indicate the requirement of designer, if *s*<sup>1</sup> = *s*<sup>2</sup> = 0, the requirements for the control system are strict, and the automatic train operation tracking control is not so, and this represents a combination of the practical situation; *y*max and *y*min represent the maximum and minimum value of design expectation, respectively, if *y*max = *y*min, it will be shown as trigonometric membership function; otherwise, it will be shown as trapezoid membership function. The corresponding diagram for fuzzy satisfaction degree calculation *μy*(*k*) of system output *y*(*k*) is shown in Figure 3.

**Figure 3.** Diagram for fuzzy satisfaction degree calculation *μy*(*k*) of system output *y*(*k*).

The error between the output value and the reference target value of system (i.e., fuzzy satisfaction degree *μy*(*k*) of system output *y*(*k*)) should also be considered. If the fuzzy satisfaction degree *μy*(*k*) is larger, it can indicate that the error between the output value and the reference target value of system is smaller, at this time, there is a small overshoot of the system and softness factor *αμ<sup>y</sup>* so that an appropriate lager value needs to chosen to increase the system rapidity; by contrast, there is an obvious overshoot of the system and softness factor *αμ<sup>y</sup>* so that an appropriate small value needs to be chosen to reduce the system rapidity to ensure system robustness [36]. According to the influence of softness factor *αμ<sup>y</sup>* for the system dynamic response, the specific exponential calculation formula for softness factor *αμ<sup>y</sup>* by fuzzy satisfaction degree *μy*(*k*) of system output *y*(*k*) is as follows,

$$\mathfrak{a}\_{\mu y}\left(y\left(k\right),y\_{r}\right) = \mathfrak{a}\_{\text{max}} + \left(\mathfrak{a}\_{\text{max}} \times e^{-\left(b \times \mu\_{y(k)}\right)} - \mathfrak{a}\_{\text{max}}\right) \times \mu\_{y(k)}\tag{28}$$

where *α*max represents the maximum value of softness factor *μy*(*k*); *b* represents the gain coefficient; it determines the shape of the softening factor *αμ<sup>y</sup>* (*y* (*k*), *yr*) function curve.

The corresponding diagram for softness factor *αμ<sup>y</sup>* of the fuzzy satisfaction degree calculation *μy*(*k*), and softness factor *αμ<sup>y</sup>* of system output *y*(*k*) are shown in Figures 4 and 5.

**Figure 4.** Diagram for fuzzy satisfaction degree calculation *μy*(*k*) of fuzzy satisfaction degree calculation *μy*(*k*).

**Figure 5.** Diagram for fuzzy satisfaction degree calculation *μy*(*k*) of system output *y*(*k*).

*3.3. Improved Whale Optimization Algorithm for Softness Factor Adaptive Adjusting Parameters Optimization*

Optimization algorithms are used to obtain a set of adjustable parameters for the satisfactory tracking control effect in actual automatic train operation scenarios. The specific softness factor adaptive adjusting parameters optimization model is as follows,

$$\begin{array}{ll}\min & F(\mathbf{x}) = (P\_{\mathbf{k}\prime} \frac{ITAE}{\max(ITAE)}, \frac{K\_{safe}}{\max\{K\_{safe}\}})\\\text{s.t.} & \mathbf{x} = (\lambda\_{dTs}, \lambda\_{a\mu y}, \alpha\_{T\prime}, \alpha\_{\max}, b) \\ & g\_{i\mathbf{j}\prime}(\mathbf{x}) \le 0, \ i\mathbf{g} = 1, 2, \cdots, n\mathbf{g} \\ & \mathbf{x} \in \Omega^{\prime\prime} \end{array} \tag{29}$$

where *x* represents the solution vector; *F*(*x*) represents the target vector; Ω represents feasible solution space of *x*; *gig*(*x*) represents the *ig*-th equality or inequality constraint for automatic train operation tracking control problem, *ng* represents the number of equalities and inequality constraints; the five adjustable parameters (*λαTs*, *λαμy*, *αTr*, *α*max, *b*) are decision variables.

Objective decomposition is an effective method to solve the multi-objective optimization problems. The Tchebycheff decomposition method is selected in this paper among many objective decomposition methods [37]. The specific calculation formula for the aggregate function value of the Tchebycheff decomposition method is as follows,

$$\begin{aligned} \text{minim} & \mathcal{G}^{tc}(\mathbf{x}|\lambda, z^{\*}) = \max\_{1 \le i \le m} \{\lambda\_i |f\_i(\mathbf{x}) - z\_i^{\*}|\} \\ \text{s.t.} & \mathbf{x} \in \Omega^{\prime \prime} \end{aligned} \tag{30}$$

where *z*∗ represents the reference point, (*z*∗ *<sup>i</sup>* = min{ *fi*(*x*)|*x* ∈ Ω}, *i* = 1, ..., *m*), which is the optimal solution of each objective function at present; *<sup>λ</sup><sup>i</sup>* is the weight of the *<sup>i</sup>* th objective, *<sup>m</sup>* ∑ *i*=1 *λ<sup>i</sup>* = 1.

Whale optimization algorithm with strong global optimization ability is chosen in this paper. Whale optimization algorithm (WOA) is a new metaheuristic optimization algorithm learned from whale predatory behavior. There are two operators (position update and prey searching) in the computation process of the whale optimization algorithm [38]. The specific calculation formula for the position update of the basic whale optimization algorithm is as follows,

$$X(t+1) = \begin{cases} X^\*(t) - A \cdot D & p < P\_s \\\ X^\*(t) + D\_p \cdot e^{Bl} \cdot \cos(2\pi l) & p \ge P\_s \end{cases} \tag{31}$$

where *X*∗(*t*) represents the optimal position vector obtained by the current optimization; *Dp* = |*X*∗(*t*) − *X*(*t*)| represents the distance between humpback whales and their prey; *p* represents the behavioral selection probability of humpback whales, *p* ∈ [0, 1]; *Ps* represents the probability of surrounding prey of humpback whales, *Ps* ∈ [0, 1]; the probability of spiral hunting is 1 − *ps*; *B* represents a constant, which is used to define the shape of spiral; *l* represents the random number in (−1, 1); *t* is the current iteration number; *Tmax* is the maximum number of iterations; *a* represents convergence factor; *A* and *C* represent the correlation coefficients respectively; *r*<sup>1</sup> and *r*<sup>2</sup> are random numbers, *r*<sup>1</sup> ∈ [0, 1], *r*<sup>2</sup> ∈ [0, 1].

The specific calculation formulas for convergence factor *a*, correlation coefficients *A*, and *C* is as follows,

$$a = 2 - 2 \times \text{t} / T\_{\text{max}} \tag{32}$$

$$A = 2a \times r\_1 - a \tag{33}$$

$$\mathbf{C} = \mathbf{2} \times \mathbf{r}\_2 \tag{34}$$

After the position updated, prey searching is implemented by means of random individual positions. The specific calculation formula for prey searching of the basic whale optimization algorithm is as follows,

$$D = \left| CX\_{rand} - X(t) \right| \tag{35}$$

$$X(t+1) = X\_{rand} - A \cdot D \tag{36}$$

where *Xrand* is the position vector of randomly selected whales. If *A* ≥ 1, a search leader individual is randomly selected, and the position of other whales is updated based on the whale position of the leader individual, so as to guide the whales to leave the prey and find a more suitable prey to enhance the global search ability of the algorithm.

The relatively fixed method of linear decline of convergence factor *a* will reduce the population diversity maintenance ability, so that the algorithm can easy to fall into local convergence in the late iteration. Aiming at solving this problem, the strategy of cosine decline combined with chaotic random method for convergence factor nonlinear decline is proposed in this paper. The specific calculation formula is as follows,

$$\begin{cases} a = 2 \cdot \cos\left(\frac{\pi}{2} \cdot \frac{t}{T\_{\text{max}}}\right) & p\_a(t) < P\_a \\ a = 2 \cdot rand^2 \cdot \sin\left(\pi \cdot rand\right) & p\_a(t) > P\_a \end{cases} \tag{37}$$

where *pa* (*t*) represents the behavioral selection probability of the convergence factor *a*, *pa* (*t*) ∈ [0, 1]; *Pa* represents the probability of cosine decline of the convergence factor *a*, *Pa* ∈ [0, 1]; and the probability of chaotic random is 1 − *Pa*.

Compared with the linear decline strategy, the decline rate of the convergence factor is significantly different in the whole iteration cycle caused by the nonlinear decline strategy with a certain degree of chaos uncertainty for convergence factor *a*, and it is helpful for maintaining the population diversity, thus the algorithm global convergence performance will be improved [39].

According to the Tchebycheff decomposition method, the aggregate function value is the fitness index for the multi-objective optimization algorithm. After the computation process of each iteration, the newly generated non-dominated solutions of the current population are put into the elite archive. The archive must kept within a certain size by some elite individuals with small differences from other elite individuals, so as to avoid computational burden of the algorithm. According to the updating rules of the whale optimization algorithm, the reference point *z*∗ plays an important role in guiding the direction of global convergence, and a certain degree of local convergence due to this fixed foraging behavior. Meanwhile, the selector, crossover, and mutation of the genetic algorithm can generate a large number of new solutions with great differences for the whale optimization algorithm based on evolutionary processes, so as to further improve the global convergence performance due to more powerful population diversity maintenance ability.

The specific steps of improved whale optimization algorithm proposed in this paper are as follows. Step 1: Initialization.

Initialize the whale population (the size is *Nw*), and the Tchebycheff aggregation function values of each whale individual are calculated.

Step 2: Iterative computations.

the *Archive* is obtained;

*Archive* = ∅, the reference point *z*∗ = (*z*∗ <sup>1</sup>, *z*<sup>∗</sup> <sup>2</sup>, ..., *z*<sup>∗</sup> *<sup>m</sup>*), and *z*<sup>∗</sup> *<sup>j</sup>* = min(*fj*(*x*)), *j* = 1, 2, ..., *m*, *m* represents the number of objectives, a uniformly distributed weight vector set *λ*<sup>0</sup> is generated, and *λ*1=*λ*0.

If the current iteration number is greater than 1, the weight *λt*,*<sup>i</sup>* for solution *xt*,*<sup>i</sup>* are need to be recalculated. According to the literature [40], in the *t*-th iteration, the specific calculation formula for weight *λt*,*i*,*<sup>k</sup>* of the *k*-th optimization index of the *i*-th individual (solution) *xt*,*<sup>i</sup>* in the population is as follows,

$$\lambda^{t\dot{t},\dot{k}} = \frac{1}{f(\mathbf{x}^{t\dot{t}})^k - Z^{ref}k} \left( \sum\_{i\mathbf{k}=1}^m \frac{1}{f(\mathbf{x}^{t\dot{t}})^{i\dot{k}} - Z^{ref,\dot{\mathbf{k}}}} \right)^{-1} \tag{38}$$

where *i* ∈ {1, 2, ..., *Nw*}, *k* ∈ {1, 2, ..., *m*}.

For any solution target *z<sup>c</sup>* = (*z<sup>c</sup>* <sup>1</sup>, *<sup>z</sup><sup>c</sup>* <sup>2</sup>, ..., *<sup>z</sup><sup>c</sup> <sup>m</sup>*) of Pareto front, its weight vector is 1 *f*(*x*)−*Z<sup>c</sup> m* ∑ *ik*=1 1 *<sup>f</sup>*(*x*)−*Zc*,*ik* <sup>−</sup><sup>1</sup> . Because the Pareto front is not easily available, it is replaced by the nearest solution target *Zref* in *Archive*;

The strategy of cosine decline combined with chaotic random method is used to calculate convergence factor *a*;

The updating rules of the whale optimization algorithm are used to update each individual whale. Step 3: The archive and genetic evolution mechanism.

The Pareto front of the current whale population is obtained, and it is used to expand the *Archive*. Some elite individuals with small differences from other elite individuals are deleted, until the size of *Archive* is not exceed allowed limit archive size *NA*.

Three operators (selection, crossover, and mutation) of the genetic algorithm are applied in the whole whale population, so as to further improve the population diversity maintenance ability.

Step 4: Termination judging.

The hypervolume indicator is chosen as termination judging indicator. Using the hypervolume indicator of the dominated portion of the objective space as a measure for the quality of Pareto set approximations is appropriate and effective. The specific formula for hypervolume indicator for objective vector set *A* with reference point (0,0,...,0) is as follows.

$$I^\*\_{\
u}(A) = \int\_{(0,0,\dots,0)}^{(1,1,\dots,1)} a\_A(z)dz\tag{39}$$

where *A* is any objective vector set in objective space Ω; if there is an objective vector *a*, Pareto superior *z* and *a* belongs to *A*, *α<sup>A</sup>* (*z*) = 1; otherwise, *α<sup>A</sup>* (*z*) = 0 [41].

Thus, the hypervolume indicator indicates the dominated portion in objective space Ω. Generally, Hypervolume as Klee's Measure Problem (HKMP) is an effective calculation method for hypervolume indicator [42]. As only 3 objects (*Pk*, *ITAE* max(*ITAE*), *Ksaf e* max(*Ksaf e*) ) must be taken into account, the calculation method by using equivalent volume model for the volume of irregular objects can be used. The specific formulas for the hypervolume indicator for 3 objects by using equivalent volume model is as follows,

$$\operatorname{I}^\*\!\_{H}(A) = \frac{\sum\_{ia=1}^{na} \sum\_{ib=1}^{nb} \sum\_{ic=1}^{nc} a\_A \left( cp(ia, ib, ic) \right)}{na \times nb \times nc} \tag{40}$$

where *na*, *nb*, and *nc* are the split numbers for each normalization objective domain (0,1); *cp*(*ia*, *ib*, *ic*) represents the central point of the (*ia*, *ib*, *ic*)th cube of normalization objective space.

The schematic diagram of the hypervolume indicator for 3 objects by using equivalent volume model is shown in Figure 6.

**Figure 6.** Schematic diagram of hypervolume indicator for 3 objects by using equivalent volume model.

If the hypervolume indicator *I*∗ *<sup>H</sup>* (*A*) is reached beforehand, an unchanged number of hypervolume indicators *nH* will be resetted; otherwise, make *nH* = *nH* + 1.

If the maximum unchanged number of the hypervolume indicator *nHmax* is reached, the calculation will be terminated; otherwise, return Step 2.

The flowchart of improved whale optimization algorithm proposed in this paper is shown in Figure 7.

Aiming at improving the global searching ability, the evolution law is must be considered in the computation process. The excellent individuals should have a small mutation probability, so that they can accumulate optimization results effectively, and the poor individuals should choose a large mutation probability, which can be fully eliminated, so as to enhance the capacity of exploration [43]. The mutation probability calculation formula based on sigmoid function *y* = <sup>1</sup> <sup>1</sup>+*e*−*<sup>x</sup>* is as follows,

$$p\_{\rm mf} = p\_{\rm m\\_max} \frac{1}{1 + e^{-ap\left(N\_s - fit(x)^\prime\right)}} \tag{41}$$

where *Pm* is the mutation probability value; *pm*\_ max is the maximum mutation probability; *ap* is the shape factor for the sigmoid function of mutation probability; *Ns* is the demarcation point of the whale population; *fit*(*x*) is the normalized value of fitness function value *fit*(*x*) of whale individual *x* in the whale population.

**Figure 7.** The flowchart of the improved whale optimization algorithm proposed in this paper.

This mutation probability calculation method has certain fairness, and the whale individuals have appropriate mutation probability according to fitness function value, so as to prevent the population controlled by advantage individuals and persist evolution opportunity for disadvantaged individuals.

A multimodal crossover method is conducive to finding a more satisfactory optimal solution for the complex optimization issue [44,45]. The multimodal crossover combining popular blended crossover and unimodal normal distribution crossover is applied in this paper. The specific calculation formula for multimodal crossover is as follows,

$$\begin{aligned} \mathbf{x}\_{\varepsilon} &= r\_{p1} \times \mathbf{x}\_{p1} + \left(1 - r\_{p1}\right) \times \mathbf{x}\_{p2} & p\varepsilon &< p\mathbf{c}\_{b} \\ \mathbf{x}\_{\varepsilon} &= \mathbf{x}\_{p} + \varepsilon d + \sum\_{\substack{i = 1, i \neq p}}^{\mu} \nu\_{ii}\mathbf{c}\_{i\varepsilon} & p\mathbf{c}\_{b} &< p\mathbf{c} < \left(p\mathbf{c}\_{b} + p\mathbf{c}\_{u}\right) \end{aligned} \tag{42}$$

where *xc* is the solution after multimodal crossover operation; *xp*<sup>1</sup> and *xp*<sup>2</sup> are two parent solutions for blended crossover; *pc* is the behavioral selection probability about multimodal crossover operation; *pcb* and *pcu* are the crossover probabilities for blended crossover and unimodal normal distribution crossover, respectively; *xp* is the midpoint for *μ* parent solutions for unimodal normal distribution crossover; *d* is the differential vector; *eic* is the *ic* th orthogonal basis; *ε* and *νic* are the random numbers obey normal distribution *N* - 0, *σ*<sup>1</sup> 2 and *N* - 0, *σ*<sup>2</sup> 2 .

#### *3.4. Performance Analysis of Optimization Algorithms Based on Standard Test Functions*

Aiming at verifying the effectiveness of IWOA proposed in this paper, standard test functions (ZDT1, ZDT3, and DTLZ2) are selected as optimization objects, and multi-objective particle swarm optimization based on decomposition (dMOPSO) [46] and multi-objective evolutionary algorithm based on decomposition (MOEA/D) [47] are selected as contrasted optimization algorithms. The improved whale optimization algorithm parameters are set as follows; maximum number of iterations is 100; probability of surrounding prey of humpback whales *Ps* is 0.6; population size *Nw* is 50; allowed limit archive size *NA* is 30; probability of cosine decline of the convergence factor *Pa* is 0.9; shape constant of spiral *b* is 1; shape factor for sigmoid function of mutation probability *ap* is 4.9; demarcation point of whale population *Ns* is 0.8; probability of blended crossover *pcb* and unimodal normal distribution crossover *pcu* are 0.3 and 0.3, respectively; selection probability is 0.5; split number for each normalization objective *na*, *nb*, and *nc* are all 50; the maximum unchanged number of hypervolume indicator *nHmax* is 25. The Matlab/simulink platform is used for verifying, and the Matlab/simulink revision and the computer processor type are 2016b, MathWorks and CPU Core i9-7920X @ 2.9GHZ. The specific optimization results (approximate Pareto solution set) for test functions of each optimization algorithms are shown in Tables 1 and 2 and Figure 8.

**Figure 8.** The optimization results for test functions of each optimization algorithms. (**a**) Optimization results for ZDT1. (**b**) Optimization results for ZDT3. (**c**) Optimization results for DTLZ2.


**Table 1.** The hypervolume indicator ratio *<sup>I</sup>*<sup>∗</sup> *<sup>H</sup>* (*A*(*op*)) *<sup>I</sup>*<sup>∗</sup> *<sup>H</sup>* (*A*(*tp*)) for test functions of each optimization algorithm.

**Table 2.** The computation time for test functions of each optimization algorithm.


As can be seen from Tables 1 and 2, compared with other optimization algorithms (dMOPSO and MOEA/D), IWOA has been improved to a considerable extent, not only in the computation speed, but also in the optimization effect for approximate Pareto solution set reflected by the hypervolume indicator ratio *<sup>I</sup>*<sup>∗</sup> *<sup>H</sup>*(*A*(*op*)) *<sup>I</sup>*∗*H*(*A*(*tp*)) , *<sup>A</sup>* (*op*) and *<sup>A</sup>* (*tp*) are the approximate Pareto solution sets obtained by optimization algorithms and real Pareto front set, respectively. According to Figure 9, compared with other optimization algorithms (dMOPSO and MOEA/D), only the better approximate Pareto solution sets closer to the Pareto fronts of ZDT1, ZDT3, and DTLZ2 were found using IWOA, but also the distribution for obtained approximate Pareto solution set was more evenly. This indicates that the IWOA proposed in this paper has better optimization effectiveness.

#### *3.5. Improved DMC Model Predictive Controller and Hardware-In-The-Loop Simulation Platform*

Based on the traditional Fuzzy DMC model predictive controller, a new functional module is necessary to be realized the function of online obtaining of softness factor and fusion velocity. The schematic diagram of improved DMC model predictive controller for automatic train operation designed in this paper is shown in Figure 9.

**Figure 9.** The schematic diagram of improved DMC model predictive controller for automatic train operation designed in this paper.

In Figure 9, the improved DMC model predictive controller could provide control commands for the corresponding equipments in real-time using fuzzy DMC MPC based on online obtaining of softness factor and fusion velocity, enabling the the urban rail vehicle to track the target velocity trajectory; the Speed and softness factor analyzer could provide the precision instantaneous fusion velocity *vis*,*<sup>a</sup>* and softness factor *α* for rolling optimization and feedback correction based on three kinds of velocity sampling sources (*vis*,*v*, *vis*,*F*, and *vis*,*s*) and a set of softness factor *α* adjustable parameters. The intelligent digital torque sensor, gear speed sensor, and displacement pickup are data acquisition devices. The physical diagram of the intelligent digital torque sensor is shown in Figure 10.

**Figure 10.** The physical diagram of the intelligent digital torque sensor.

In Figure 10, the acquisition equipment is an intelligent digital torque sensor of model No. JN338; the fixed bracket is made of iron and has two functions of fixing and preventing jitter; the power source supplies electricity to torque sensor of model No. JN338; the communication module transmits the real-time value of motor speed and torque using the 485 communication protocol; the sampling data connectors are connectors for sampling data (motor speed and torque); and can be connected to communication module, controller, or other equipment; rotation shaft of intelligent digital torque sensor must be connected to rotation shaft of traction motor and load (breaker or rheostat box).

To more effectively test the performance of the tracking control algorithm in actual automatic train operation tracking control scenarios, the dSPACE HILS technology is adopted. In this way, the optimization algorithm or control algorithm needed to be verified is written into the chip of the optimizer or controller. The structure diagram of HILS platform used in this paper for automatic train operation tracking control scenario, and the physical diagram of controller cabinet and simulation cabinet for HILS platform are shown in Figures 11 and 12.

In Figure 11, the Displacement generator is used to generate the train instantaneous displacement, so that the static (no actual displacement) HILS platform can truly reflect the actual automatic train operation tracking control scenario; various sensors are used to feed electrical waves of sampling sources back to the Controller in real-time; the Conditioning circuit can regulate electrical signals properly for the Tracking controller appropriately; the Motor controller could provide electrical control commands for Traction moto and other corresponding equipments of Electrical loop in real-time using a proper electrical control algorithm. DC power source, Converter system, Traction motor, Digital rheostat box, and Gear box are simulation electric hardware equipments.

In Figure 12, the 'train controller cabinet' and 'train emulator cabinet' are the vital equipments for automatic train operation HILS, except for the controller and emulator, the conditioning circuit, signal processing unit, and corresponding switch groups are included. The 'emulator' provides some correlative simulation environments for the automatic train operation HILS, the related models included such as accurate braking model, traction transformer model, running line model, velocity fluctuation model, etc. The 'conditioning circuit' can regulate electrical signals properly for 'Controller'. The 'signal processing unit' can regulate net signals properly for 'Optimizer' appropriately.

**Figure 11.** The structure diagram of the hardware-in-the-loop simulation (HILS) platform used in this paper for automatic train operation tracking control scenario.

**Figure 12.** The physical diagram of controller cabinet and simulation cabinet for HILS platform.

#### **4. Simulation for Automatic Train Operation Tracking Control Scenario**

#### *4.1. Data and Parameters for Automatic Train Operation Tracking Control Scenario*

The automatic train operation tracking control scenario for rail transit line No.12 in Dalian, China is chosen as the experimental simulation object. Rail transit line No.12 is a significant urban rail transit line with 40.38 kilometers from Hekou station to Lvshun New Port. The running simulation line of scenario about rail transit line No.12 is from Lvshun New Port to Tieshan Town, there are three long steep ramps and three velocity limit subintervals in running interval. The main parameters of the automatic train operation tracking control scenario are shown in Table 3, and the target velocity trajectory, slopes, and limited velocity curves for automatic train operation are shown in Figure 13.


**Table 3.** The main parameters of the automatic train operation tracking control scenario.

The basic DMC MPC parameters are set as follows; sampling time is 500 μ*s*; model length *N* is 60; control length is 15; predictive length is 15. The addition adjustable parameters for online obtaining of softness factor are set by practical experience as follows; blur width *s*<sup>1</sup> = *s*<sup>2</sup> = 0.05 km/h; maximum and minimum value of design expectation *y*max = *yr* + 0.05 km/h and *y*min = *yr* − 0.05 km/h; connected length *S*<sup>1</sup> = *S*<sup>2</sup> = 0.4 m. Considering the online real-time calculation efficiency and tracking control effect of the improved DMC MPC proposed in this paper, the following parameters are given based on the relevant scientific literature, field experience, and simulation results of multiple experiments. The Matlab/simulink simulation platform is used for softness factor adaptive

adjusting parameters optimization and the experience parameters are as follows; maximum allowed multi-objective performance index, *ITAE* index, and security index are 0.8, 750, and 6%, respectively; ideal multi-objective performance index, *ITAE* index, and security index are 0.2, 250, and 3%, respectively. The chosen addition adjustable parameters for online obtaining of softness factor obtained by each optimization algorithms (improved whale optimization algorithm, MOEA/D, dMOPSO) are as follows; fusion weights of *λαTs* and *λαμ<sup>y</sup>* , *λαTs* = 0.582, 0.592, 0.573 and *λαμ<sup>y</sup>* = 0.418, 0.408, 0.427; maximum value of softness factor *αmax* is 0.939, 0.941, and 0.930; gain coefficient *b* is 0.909, 0.911, and 0.913. The subinterval range obtained by practical experience, reference value of soften factor obtained by each optimization algorithms, synthetic weight of fusion velocity obtained by entropy weight method are shown in Table 4. The specific optimization results (approximate Pareto solution set) for softness factor adaptive adjusting parameters optimization of each optimization algorithms are shown in Table 5 and Figure 14.

**Figure 13.** The target velocity trajectory, slopes, and limited velocity curves for automatic train operation tracking control scenario. (**a**) Target velocity trajectory. (**b**) Slopes and limited velocity curves.

**Table 4.** The optimization results for subinterval range, reference value of soften factor *αTr*, and synthetic weight of fusion velocity *λis*.


**Figure 14.** The optimization results for the softness factor adaptive adjusting parameters optimization of each optimization algorithm.


**Table 5.** The hypervolume indicator ratio *<sup>I</sup>*<sup>∗</sup> *<sup>H</sup>* (*A*(*op*)) *<sup>V</sup>*(Ω*S*) and computation time for softness factor adaptive adjusting parameters optimization of each optimization algorithms.

In Figure 14, the approximate Pareto solution sets have been obtained by optimization algorithms (IWOA, dMOPSO, and MOEA/D), and one of the approximate Pareto solutions has been chosen for automatic train operation tracking control scenario. As can be seen from Figure 14, compared with other optimization algorithms (dMOPSO and MOEA/D), the wider dominated portion for approximate Pareto solution sets obtained by using IWOA. This indicates that the IWOA proposed in this paper has better optimization effectiveness. As can be seen from Table 5, compared with other optimization algorithms (dMOPSO and MOEA/D), IWOA has been improved to a considerable extent not only in the computation speed but also in the optimization effect for approximate Pareto solution set reflected by the hypervolume indicator ratio *<sup>I</sup>*<sup>∗</sup> *<sup>H</sup>*(*A*(*op*)) *<sup>V</sup>*(Ω*S*) . *<sup>V</sup>* (Ω*S*) is the volume of selected objective space Ω*<sup>S</sup>* by experience, and it is 0.16 ( <sup>3</sup> <sup>5</sup> <sup>×</sup> <sup>2</sup> <sup>3</sup> <sup>×</sup> <sup>2</sup> <sup>5</sup> ) in this paper.

#### *4.2. Matlab/simulink Simulation Results for Automatic Train Operation Tracking Control Scenario*

The sampling process is not necessary in the Matlab/simulink simulation platform. According to the automatic train operation tracking control scenario of rail transit line No.12 in Dalian, China, the Matlab/simulink simulation results are obtained by using the fuzzy DMC MPC based on online obtaining of softness factor *α*, and the softness factor adaptive adjusting parameters optimization using IWOA proposed in this paper (denoted as DMC MPC II); the fuzzy DMC MPC is based on online obtaining of softness factor *α*, which uses softness factor adaptive adjusting parameters optimization using MOEA/D (denoted as DMC MPC I) and traditional fuzzy DMC MPC. The specific configuration of the Matlab/simulink platform used in this paper is described as follows; the Matlab/simulink revision is 2016b, MathWorks; the type of computer processor is CPU Core i9-7920X @ 2.9GHZ. The specific Matlab/simulink results are shown in Figures 15–18 and Tables 6 and 7.

**Figure 15.** The Matlab/simulink velocity trajectory curves of different DMC MPC algorithms for automatic train operation tracking control scenario.

As can be seen from Tables 6 and 7, the tracking control results obtained by the DMC MPC II are superior to that of DMC MPC I and traditional fuzzy DMC MPC, and four indexes of multi-objective performance index (energy saving, punctuality, parking precision, and comfort), ITAE index, and security index for automatic train operation have been improved considerably. As can be seen from the Figures 15 and 16, the DMC MPC II can make the tracking control curves closer to target curves, so as to obtain the ideal tracking control results as smooth as possible. As can be seen from the six enlarged areas of the velocity trajectory curves in Figures 15 and 16, the velocity fluctuation degree is weaker and the velocity trajectory is closer to the target by using DMC MPC II. As can be seen from the one enlarged areas of time distance curves of Figure 16, compared with DMC MPC I and traditional

fuzzy DMC MPC, the time traceability of DMC MPC II is more powerful, so as to obtained the time distance curve closer to target. As can be seen from Figure 17, the smaller velocity error effect can be obtained by using DMC MPC II. As can be seen from Figure 18, the more ideal parking results can be obtained by using DMC MPC II; both the distance and time errors of parking are reduced to certain extent, so as to improve the punctuality and fixed position effect.

**Figure 16.** The Matlab/simulink time traceability curves of different DMC MPC algorithms for automatic train operation tracking control scenario. (**a**) Time–velocity curves. (**b**) Time–distance curves.

**Figure 17.** The Matlab/simulink time velocity error curves of different DMC MPC algorithms for automatic train operation tracking control scenario.

**Table 6.** The Matlab/simulink tracking control results of energy saving, punctuality, parking precision, and comfort for automatic train operation.


**Table 7.** The Matlab/simulink tracking control results of multi-objective performance index, ITAE index, and security index for automatic train operation.


Compared with traditional fuzzy DMC MPC and DMC MPC I, DMC MPC II has several obvious superiorities in the matlab/simulation environment. However, as there is no hardware equipment in the actual automatic train operation tracking control scenario in matlab/simulation environment, the effectiveness of DMC MPC proposed in this paper must be further tested and verified.

**Figure 18.** The Matlab/simulink parking error curves of different DMC MPC algorithms for automatic train operation tracking control scenario. (**a**) Distance–velocity curves. (**b**) Time–velocity curves. (**c**) Time–distance curves.

#### *4.3. HILS Results for Automatic Train Operation Tracking Control Scenario*

In this way, sampling accuracy must be taken into account. To further verify the effectiveness of the algorithm, according to the automatic train operation tracking control scenario of rail transit line No.12 in Dalian, China, the HILS results are obtained by using the fuzzy DMC MPC based on

online obtaining of softness factor *α* and fusion velocity, which uses the softness factor adaptive adjusting parameters optimization using IWOA proposed in this paper (denoted as DMC MPC V), the fuzzy DMC MPC based on online obtaining of softness factor *α* and fusion velocity, which softness factor adaptive adjusting parameters optimization using dMOPSO (denoted as DMC MPC IV) and the fuzzy DMC MPC based on online obtaining of fusion velocity (denoted as DMC MPC III). The specific configuration of the automatic train operation HILS platform used in this paper is described as follows; the Matlab/simulink revision is "2016b, MathWorks"; the type of computer processor is "CPU Core i9-7920X @ 2.9GHZ"; the core chip of "Tracking controller" and "Motor optimizer" is "TMS320F28335"; the simulation software of "dSPACE emulator" is dSPACE control desk (revision is control desk 6.1); the communication protocol of the HILS platform is MVB (multifunction vehicle bus); the fuzzy PID (proportion integration differentiation) algorithm is adopted as motor control algorithm; vehicle velocity proportion is (0.83 × 400 rad/min)/(80 km/h). The specific HILS results are shown in Figures 19–22 and Tables 8 and 9.

**Table 8.** The HILS tracking control results of energy saving, punctuality, parking precision, and comfort for automatic train operation.


**Table 9.** The HILS tracking control results of multi-objective performance index, ITAE index, and security index for automatic train operation.


According to the HILS results of different algorithms from Tables 8 and 9, compared with DMC MPC III and DMC MPC IV, DMC MPC V has an obvious performance improvement effectiveness, the multi-objective performance index (energy saving, punctuality, parking precision, and comfort) of the tracking control trajectory has been improved considerably; meanwhile, the ITAE index and security index have also been reduced considerably. In Figure 19, during the automatic train operation tracking control experiment simulation, all the pilot lights and buttons are in normal. As can be seen from Figures 19 and 20, the DMC MPC V can bring the tracking control curves closer to the target curves, so as to obtain the ideal tracking control results as smooth as possible. As can be seen from the six enlarged areas of velocity trajectory curves of Figures 19 and 20, the velocity trajectory curves obtained by DMC MPC V were smoother; compared with the traditional improved tracking control algorithm (DMC MPC), DMC MPC V enables the train to be in the optimal working state as much as possible, so as to reduce the velocity fluctuation degree and obtain more ideal tracking control results. As can be seen from the one enlarged areas of the time–distance curves in Figure 20, compared with DMC MPC III and DMC MPC IV, the time traceability of DMC MPC V is more powerful, so as to obtain a time–distance curve closer to target. As can be seen from Figure 21, the velocity error obtained by using DMC MPC V is smaller. As can be seen from Figure 22, the more ideal parking point (parking time and position) can be obtained by using DMC MPC V, its parking point is closer to prospective parking point (180 s and 2940 m).

**Figure 19.** The HILS velocity trajectory curves of different DMC MPC algorithms for automatic train operation tracking control scenario.

**Figure 20.** The HILS time traceability curves of different DMC MPC algorithms for automatic train operation tracking control scenario. (**a**) Time–velocity curves. (**b**) Time–distance curves.

**Figure 21.** The HILS time–velocity error curves of different DMC MPC algorithms for automatic train operation tracking control scenario.

**Figure 22.** The HILS parking error curves of different DMC MPC algorithms for automatic train operation tracking control scenario. (**a**) Distance–velocity curves. (**b**) Time–velocity curves. (**c**) Time–distance curves.

The above HILS results show that DMC MPC V is a tracking control algorithm with good practical tracking control effect for automatic train operation tracking control scenario.

#### **5. Conclusions**

Tracking control optimization for automatic train operation is a sophisticated optimization problem, and the model predictive controller is widely used to solve this problem due to its advantages of strong robustness and good performance in tracking speed and tracking precision. Aiming at obtaining a more ideal tracking control performance for automatic train operation, an improved model predictive control algorithm and corresponding controller based on online obtaining of softness factor and fusion velocity for automatic train operation are proposed and developed, and an improved whale optimization algorithm based on Tchebycheff decomposition method was proposed for softness factor adaptive adjusting parameters optimization. This clearly shows that model predictive controller has been improved to a considerable extent in tracking control optimization for automatic train operation not only in the pure software scenarios but also in the hardware-in-the-loop simulation scenarios (the

*ITAE* index obtained by DMC MPC II is 16.3% and 31.2% lower than that of DMC MPC I and Fuzzy DMC MPC, and that obtained by the DMC MPC V is 7.3% and 16.1% lower than that of DMC MPC IV and DMC MPC III). The specific advantages are described below.


According to the Matlab/simulink results and ATO HILS results (compare with the other DMC MPC algorithms for comparison), the improved DMC MPC based on online obtaining of softness factor and fusion velocity proposed in this paper has better tracking control performance, so it can obtain more ideal tracking control results. In actuality, it can also be designed by other ways for online obtaining of softness factor and fusion velocity. It is important to note that, as the computational tasks such as various matrix computation of basic DMC MPC are onerous, the designed scheme should be simple and appropriate caused by the limited computation margin in real-time.

**Author Contributions:** The work presented here was performed in collaboration among all authors. L.W. designed, analyzed, and wrote the paper, and completed the simulation experiment. X.W. guided the full text and provided simulation conditions. Z.S. conceived the idea and involved simulation experiment. S.L. involved simulation experiment and analyzed the data. All authors have contributed to and approved the manuscript.

**Funding:** This research was funded by the Nature Science Foundation of China (grant number 60574018).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript.



#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
