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

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.


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][4][5]. Meanwhile, various improved intelligent optimization control algorithms have been proposed and applied for the automatic train operation system [6][7][8].
(I) An improved whale optimization algorithm (IWOA) based on the Tchebycheff decomposition method, convergence factor nonlinear decline, and genetic evolution measurement is proposed to solve the optimization of the softness factor adaptive adjusting parameters appropriately. (II) Aiming at improving the tracking control performance for automatic train operation, an improved model predictive controller based on online obtaining of the softness factor and fusion velocity is developed for automatic train operation tracking control effectively. (III) To further verify the effectiveness of the developed model predictive control controller, the scenario about rail transit line No.12 in Dalian, China is chosen for simulation test. The results of the Matlab simulation and hardware-in-the-loop simulation (HILS) show that the tracking controller proposed in this paper has good tracking control performance.
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.

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, (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 P k is as follows, where ω i represents the index importance weight factor ( k ∑ i=1 ω i = 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; F t (u, v) and B r (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; s z is the terminal position; T r 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; v lim (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; K Jerk 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 K sa f e is as follows, 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].

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].

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, where Y p (k) = y p (k + 1|k) , y p (k + 2|k) , ..., y p (k + N|k) T represents the predictive value of system output, Y 0 (k) = [y 0 (k + 1|k) , y 0 (k + 2|k) , ..., y 0 (k + N|k)] T represents the predictive value of predictive model, ∆U (k) = [∆u (k) , ∆u (k + 1) , ..., ∆u (k + N − 1)] T 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 Y p are as follows, where i represents the element index of Y p , i ∈ {1, 2, ..., N}; k represents the initial point of DMC predictive model [31,32].

Rolling Optimization
With the aim of avoiding the violent fluctuations in the control process effectively, it is necessary to make the final output value y f reach to the reference target value y r 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, where y f (k + i) represents the final output value expected, α i represents the ith softness factor (0 < α i < 1), y (k) represents the actual output value of the system, and y r 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, represents the control sequence of the system, R = diag[r 1 , r 2 , ..., r L ] T represents the weight coefficient matrix of constraint for error revise, Q = diag[q 1 , q 2 , ..., q M ] T 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 ∂J ∂∆U L (k) = 0 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.
Then, the actual control quantity u(k) can be obtained. The specific calculation formula of the actual control quantity u(k) is as follows, 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].

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.
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, where C represents the error corrected matrix, and its length is N; H represents the corresponding transformed matrix.

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, (13) where C f uzzy1 and C f uzzy2 represent the fuzzy inference functions by using two kinds of fuzzy rules, respectively; u f _p (k) represents the calculated control quantity by using fuzzy rule about slope and velocity error; ∆u f _p (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 u c (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. 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. 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.

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, where v is,a represents the final calculated velocity by speed analyzer ultimately of the i-th sampling point; v is,v represents the velocity calculated based on actual motor speed sampled of the i-th sampling point; v is,F represents the velocity calculated based on actual motor torque sampled of the i-th sampling point; v is,s represents the velocity calculated based on actual train instantaneous displacement sampled of the i-th sampling point; n is represents the actual motor speed sampled by speed sensors of the i-th sampling point; tr ntv represents the transmission ratio of the motor speed to train velocity; η g represents the degree of tooth engagement between the gears; η is,T 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; F is represents the force calculated based on actual motor torque sampled of the i-th sampling point; F is = η is,F × T is / R mr ; T is represents the actual motor torque sampled by torquemeter of the i-th sampling point; η is,F represents the force comprehensive transmission efficiency of the i-th sampling point; R mr represents the radius of motor rotor; w is (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; s is 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, where j = (1, 2, · · ·, m),i = (1, 2, · · ·, n); m represents index number; n represents the number of samples; x ij represents the j-th processed index value of the i-th sample after normalization; x ij 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 x ij 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, where x ij 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, where p ij represents the j-th index weight value of i-th sample in j-th index; e j 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, where d j represents the j-th entropy redundancy value of j-th index, it indicates the difference degree of this index; e j represents the j-th index entropy value; λ j represents the j-th weight value.

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, where v is,x represents a velocity sampling source x from v is,v , v is,F , v is,s using to synthetic final calculated velocity v is,a ; v is,re f represents the reference velocity (auxiliary corrective velocity sampling source) based on actual sampling data sampled by specific auxiliary sensor; ∆v is,c represents the actual error value between v is,x and v is,re f ; ∆v is,p represents the maximum permit satisfied error value between v is,x and v is,re f ; v is,c represents the final corrected value calculated by Formula (24) when correctness condition (Formula (23)) is not satisfied. Reference velocity v is,re f 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.

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, 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; α µy (y (k) , y r ) represents the real-time softness factor calculated based on tracking control condition of the present system output y (k); λ αTs and λ α µy represent the fusion weights of α Ts and α µy (y (k) , y r ), respectively; and λ αTs + λ α µy = 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, where si represents the subinterval index si ∈ {1, 2, ..., si max }; si_max represents the number of subintervals; S si represents the starting position of the si-th subinterval; S si_max+1 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 +1 = α Ts,si_ max ; S 1 and S 2 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 α µy 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, where s 1 and s 2 represent the blur width, which can indicate the requirement of designer, if s 1 = s 2 = 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. 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 α µy 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 α µy 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 α µy for the system dynamic response, the specific exponential calculation formula for softness factor α µy by fuzzy satisfaction degree µ y(k) of system output y(k) is as follows, where α max represents the maximum value of softness factor µ y(k) ; b represents the gain coefficient; it determines the shape of the softening factor α µy (y (k) , y r ) function curve. The corresponding diagram for softness factor α µy of the fuzzy satisfaction degree calculation µ y(k) , and softness factor α µy of system output y(k) are shown in Figures 4 and 5.

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, where x represents the solution vector; F(x) represents the target vector; Ω represents feasible solution space of x; g ig (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, where z * represents the reference point, (z * i = min{ f i (x)|x ∈ Ω}, i = 1, ..., m), which is the optimal solution of each objective function at present; λ i is the weight of the i th objective, 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, where X * (t) represents the optimal position vector obtained by the current optimization; D p = |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]; P s represents the probability of surrounding prey of humpback whales, P s ∈ [0, 1]; the probability of spiral hunting is 1 − p s ; 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; T max is the maximum number of iterations; a represents convergence factor; A and C represent the correlation coefficients respectively; r 1 and r 2 are random numbers, The specific calculation formulas for convergence factor a, correlation coefficients A, and C is as follows, 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, where X rand 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, where p a (t) represents the behavioral selection probability of the convergence factor a, p a (t) ∈ [0, 1]; P a represents the probability of cosine decline of the convergence factor a, P a ∈ [0, 1]; and the probability of chaotic random is 1 − P a . 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 * 1 , z * 2 , ..., z * m ), and z * j = min( f j (x)), j = 1, 2, ..., m, m represents the number of objectives, a uniformly distributed weight vector set λ 0 is generated, and If the current iteration number is greater than 1, the weight λ t,i for solution x t,i are need to be recalculated. According to the literature [40], in the t-th iteration, the specific calculation formula for weight λ t,i,k of the k-th optimization index of the i-th individual (solution) x t,i in the population is as follows, where i ∈ {1, 2, ..., Nw}, k ∈ {1, 2, ..., m}.
For any solution target z c = (z c 1 , z c 2 , ..., z c m ) of Pareto front, its weight vector is Because the Pareto front is not easily available, it is replaced by the nearest solution target Z re f 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 N A.
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.
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.
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, α A (z) = 1; otherwise, α A (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 (P k , ITAE max(ITAE) , K sa f e max(K sa f 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, 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. If the hypervolume indicator I * H (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 = 1 1+e −x is as follows, where P m is the mutation probability value; p m_ max is the maximum mutation probability; ap is the shape factor for the sigmoid function of mutation probability; N s is the demarcation point of the whale population; f it(x) is the normalized value of fitness function value f it(x) of whale individual x in the whale population. 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, where x c is the solution after multimodal crossover operation; x p1 and x p2 are two parent solutions for blended crossover; pc is the behavioral selection probability about multimodal crossover operation; pc b and pc u are the crossover probabilities for blended crossover and unimodal normal distribution crossover, respectively; x p is the midpoint for µ parent solutions for unimodal normal distribution crossover; d is the differential vector; e ic is the ic th orthogonal basis; ε and ν ic are the random numbers obey normal distribution N 0, σ 1 2 and N 0, σ 2 2 .

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 P s is 0.6; population size Nw is 50; allowed limit archive size N A is 30; probability of cosine decline of the convergence factor P a 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 N s is 0.8; probability of blended crossover pc b and unimodal normal distribution crossover pc u 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.  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 I * H (A(tp)) , A (op) and A (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.

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. 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 v is,a and softness factor α for rolling optimization and feedback correction based on three kinds of velocity sampling sources (v is,v , v is,F , and v is,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. 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.

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. 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 1 = s 2 = 0.05 km/h; maximum and minimum value of design expectation y max = y r + 0.05 km/h and y min = y r − 0.05 km/h; connected length S 1 = S 2 = 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 λ α µy , λ αTs = 0.582, 0.592, 0.573 and λ α µy = 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 14. The optimization results for the softness factor adaptive adjusting parameters optimization of each optimization algorithm. 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 is the volume of selected objective space Ω S by experience, and it is 0.16 ( 3 5 × 2 3 × 2 5 ) in this paper.

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. 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.   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.

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.  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).    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.

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.
(I) Aiming at improving the efficiency of the whale optimization algorithm based on the Tchebycheff decomposition method, the strategy of cosine decline combined with chaotic random method for convergence factor nonlinear decline is proposed, so as to obtain more satisfactory softness factor adaptive adjusting parameters for tracking control. (II) Not only is an improved online adaptive adjusting method for softness factor based on fuzzy satisfaction of system output value and velocity distance trajectory characteristic adopted, but also a fusion velocity model and a corrected model of real-time sampling for automatic train operation tracking control are adopted. Thus, compared with traditional improved model predictive controller, the improved model predictive controller developed in this paper based on online obtaining of softness factor and fusion velocity could enable the train in the optimal working state as much as possible, so as to obtain a more ideal tracking control result with more satisfactory performance indexes, including energy saving, punctuality, parking precision and comfort, ITAE, and security index. (III) For any tracking control system, the accomplishing capacity of computational tasks real-time is significant important. The only purpose of the improved strategies is the online obtaining of optimal softness factor and fusion velocity, so as to obtain the more reasonable real-time control quantity u(k), and enable the automatic train operation tracking control system robustness and rapidity as much as possible. Thus, the quantity of additional computational tasks is not very large. In addition, some complex computational tasks for adjustable parameters optimization for softness factor adaptive adjusting and setting synthetic weight of the velocity sampled according to entropy weight method are achieved offline. Then, the advanced urban rail vehicle velocity monitoring device and the additional function chip of online obtaining of softness factor and fusion velocity are applied to improved the accomplishing capacity of computational tasks real-time is significant important. Finally, some complex functions such as logarithm or exponential function are avoided in online computation. Thus, it has advantage of simple and easily conducted.
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.

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

Abbreviations
The following abbreviations are used in this manuscript.

ATO
automatic train operation HILS hardware-in-the-loop simulation DMC MPC dynamic matrix control model predictive control DMC MPC I fuzzy DMC MPC based on online obtaining of softness factor which softness factor adaptive adjusting parameters optimization using MOEA/D DMC MPC II fuzzy DMC MPC based on online obtaining of softness factor which softness factor adaptive adjusting parameters optimization using IWOA DMC MPC III fuzzy DMC MPC based on online obtaining of fusion velocity DMC MPC IV fuzzy DMC MPC based on online obtaining of softness factor and fusion velocity which softness factor adaptive adjusting parameters optimization using dMOPSO DMC MPC V fuzzy DMC MPC based on online obtaining of softness factor and fusion velocity which softness factor adaptive adjusting parameters optimization using IWOA