Two-step Optimization Method of Freight Train Speed Curve Based on Rolling Optimization Algorithm and MPC

1. Introduction

In recent years, the freight transportation has quickly increased [1]. As the freight trains are heavy and long, it is difficult to operate the freight trains smoothly. If the longitudinal impulse is too big, it will shorten the life of train couplers, and even cause fatal accidents such as coupler breakdown and derailment. The operation strategies also influence the train energy consumption. Due to the differences of driving skills and experience of the drivers, the train energy consumption is also different. Even for the same driver, the train operation is also affected by his state, physical and mental, operating conditions, line conditions, and the train condition. Therefore, it is necessary to improve the train operation by providing a smooth and energy-efficient train speed curve, which can be integrated into the Driver Advisory System (DAS) or Automatic Train Operation (ATO) system.

When optimizing the train speed curve, establishing a train dynamic model is an essential step. Single mass model and multi-mass model are both widely studied by scholars. The single mass model is characterized by its comparatively fast computation speed. Howlett used the single mass model to investigate the energy-efficient train control method [2]. Based on the single mass model, Howlett used a new local energy minimization principle to calculate the optimal switching points with steep gradients [3]. Pudovikov simulated the train operation dynamic processes to improve the control effect of the automatic control system for freight train speeds [4]. The single mass model cannot depict the in-train force, while the multi-mass model can consider the in-train force. Consequently, the optimization outputs based on multi-mass model are more applicable to the real-word situation. Gruber introduced the concept of multi-mass model by studying the problem of optimizing the heavy trains towed by multiple locomotives [5]. Chou conducted an optimal control study for the energy saving and speed tracking problems using the LQR (Linear Quadratic Regulator) methodology based on the validated multi-mass model [6]. Zhuan proposed an output regulation method with measurement feedback control of heavy-duty trains based on the multi-mass model, controlling the train speed within a specified speed range [7].

The train speed curves are initially optimized by Pontryagin maximum principle rules and other algorithms. Howlett provided the key equations for optimal control strategies with minimum energy consumption and punctual operation under continuously changing gradients [8]. Khmelnitsky investigated the problem of train energy-saving operation under arbitrary changes of speed and gradient, and realized the corresponding energy-saving control strategy in this context [9]. Martinis proposed a method for generating energy-saving speed curves for freight trains by analyzing data from a real-time vehicle monitoring system [10]. Ko used a dynamic programming method to generate energy-efficient train operating speed curves, and reduced the computation time and speed errors by improving the dynamic programming method [11]. Albrecht proposed a two-layer optimization algorithm, which is able to compute and generate a sequence of energy-saving optimization operations for trains with a minimum number of condition transitions [12]. Guo used the particle swarm algorithm and the NSGA-II (Non-dominated Sorting Genetic Algorithm-II) algorithm to propose a multi-objective optimization algorithm based on the energy consumption and the running time [13]. Yi used the greedy algorithm and improved multi-objective bald eagle search algorithm to solve a multi-objective optimization problem for freight trains [14]. Yang designed a fuzzy adaptive genetic algorithm (FAGA) for multi-objective curve optimization strategy, and used the improved high-order model-free adaptive iterative learning control algorithm to track the optimized speed curve [15].

There are many methods that can be used in freight train control to track the reference speed curve, in which the Model Predictive Control Algorithm considering the train multi-mass model is widely used in solving the train control problem. The MPC algorithm can also be seen as a short-distance optimization method for the train speed curve. Isna Silvia applied the MPC method to predict and adjust the dwell time and the running time of the train [16]. Song designed a constrained MPC closed-loop control algorithm, customized the target curve acting on the MPC according to the operating conditions of the loaded train [17]. Liu solved the optimal control problem based on the MPC framework under the constraints of operating conditions, and analyzed the impact of different weight factors. The simulation results of MPC are compared with fuzzy control and PID (Proportion Integration Differentiation) control, and verified the feasibility of the model predictive controller [18]. A nonlinear MPC algorithm is proposed to optimize train adhesion control by identifying the trajectory parameters online using a dynamic forgetting factor recursive least squares method, and then these parameters are applied to the MPC to achieve optimal adhesion control [19]. In order to solve the optimal operation of heavy-duty trains, the MPC method was used to optimize train performance and operation through the analysis of various penalty factors, weight coefficients, etc., and the MPC reference speed are self-setting step-type curves [20,21,22,23]. The above studies on MPC algorithms usually use a step curve to test the control performance, and it is possible to optimize the train operation by taking advantage of the MPC algorithm. MPC can only be used in short term optimization of the freight train speed curve.

To leverage the advantages of the single mass model and multi-mass model, a Two-Step Optimization Method is proposed in this paper. Firstly, using ROA to generate a speed curve based on a single mass model. Secondly, inspired by the MPC control method, the paper will use the MPC method to re-optimize the train speed curve based on the multi-mass model.

2. The Two-Step Optimization Method

In this paper, a Two-Step Optimization Method is adopted to optimize the freight train speed curve, and the framework of the proposed algorithm is shown in Figure 1. ROA is designed based on the previous work of our group [24], which is on the basis of the single mass model. In this algorithm, the length and weight of the train are taken into account to equalize the values of the line gradient according to the gradient equivalence method. The optimization index includes the train running time and energy consumption.

2.1. Gradient Equivalent Method

When the single mass model is used in the freight train simulation, the gradient resistance force will change suddenly at the slope changing positions. Therefore, a gradient equivalence method is proposed, taking the length and weight of the train into account. Assume that the weight of the train is uniformly distributed along the length of the train, and the train slope resistance is distributed along the length of the train.

Define the equivalent gradient $i_c\left(x\right)$ when the position of the train head is x, which can be calculated by

$$i_c\left(x\right)={\displaystyle \frac{{\sum }_{k=1}^n{m}_k{l}_k{i}_k\left(x\right)}{LM}},$$

(1)

where n is the number of all locomotive/wagon in the whole train, $m_k$ is the weight of the kth locomotive/wagon, $l_k$ is the length of the kth locomotive/wagon, $i_k\left(x\right)$ is the original gradient value where the kth locomotive/wagon is located, L is the total length of the train, and M is the total weight of the train. A gradient comparison is given by Figure 2, where the sub-figure above gives the original gradient and the sub-figure below gives the equivalent gradient of a demo line.

When the cruising speed is given, the line gradients can be classified into four types as shown in Table 1, where $i_s$ is defined as the gradient on which the train can keep the cruise speed with maximum traction force. Similarly, $i_x$ is defined as the gradient on which the train can keep the cruise speed with maximum braking force. Since there is no steep downhill section in the selected simulation line, it will not be discussed in the subsequent research.

According to Equation (1), calculate the equivalent gradient at each position and combine the successive gradients with the same type as a equivalent slope. The equivalent slopes of the line can be obtained as Figure 3. In order to find the optimal operation solution, set three slopes as a group, six slope combinations can be obtained as [a b a], [a c a], [b a c], [b a b], [c a b], [c a c]. As can be seen from the Figure 3, the slope combinations include [a b a], [b a c], [a c a].

2.2. Rolling Optimization Algorithm (ROA) for Freight Train Speed Curve

A Rolling Optimization Algorithm is proposed to find the optimal freight train speed curve. The flowchart of the rolling optimization algorithm is shown in Figure 4. The cruise speed is calculated according to the train information, and the switching points of the optimal operation strategies are optimized by the bisection method. Based on each cruise speed, the speed curve is calculated in a rolling strategy. Firstly, equivalent gradients of the line are calculated by Equation (1). Secondly, take three gradients as a group in each rolling step and calculate the optimal operation strategies according to the optimization rules. Roll forward to the next three gradients to continue the calculation, and stop the process until reaching the end of the line.

2.2.1. Train Speed Curve Optimization for Each Gradient Group

The ROA takes punctuality and energy saving as the optimization objectives, and establishes the optimization objective function as shown in

$$J=\sum _{i=1}^N\alpha ·|T_s-T_p|+\beta ·E_{S_s,S_e}\to min,$$

(2)

where $T_s$ is the simulated train running time in the section, $T_p$ is the planned running time in the section, $E_{S_s,S_e}$ is the energy consumption of the train running in the section $[S_s,S_e]$, $S_s$ is the starting point of the current train operating section, $S_e$ is the end point of the current train operating section, $\alpha$ and $\beta$ are the weights of the objective function, N is the slope combination number.

Considering the time constraints and throttle constraints of the train during operation, the constraint analysis is given. Regarding the time constraints, the optimized train running time will have a certain error with the planned running time. Define the maximum deviation value as $t_{max}$, and then the section running time constraints shown in (3). In this way, the whole process running time constraints can also be defined as

$$\begin{array}{c}\hfill |T_s^i-T_p^i|\le t_{max}\\ \hfill |\sum _{i=1}^N(T_s^i-T_p^i)|\le T_{max},\end{array}$$

(3)

where $T_{max}$ is the maximum error value of the whole running time.

Regarding the throttle constraints, it can be expressed by:

$$G_{min}\le G\le G_{max},$$

(4)

where G is the current throttle of the train, $G_{max}$ is the maximum throttle that the train is allowed, $G_{min}$ is the minimum throttle that the train is allowed.

Therefore, the optimization objective function and constraints of the ROA are as follows:

$$\begin{array}{c}\hfill J=\sum _{i=1}^N\alpha ·|T_s^i-T_p^i|+\beta ·E_{S_s,S_e}\to min\\ \hfill s.t.\left\{\begin{array}{c}\hfill |T_s^i-T_p^i|\le t_{max}\hfill \\ \hfill |\sum _{i=1}^N(T_s^i-T_p^i)|\le T_{max}\hfill \\ \hfill G_{min}\le G\le G_{max}\hfill \end{array}\right..\end{array}$$

(5)

The main optimization rules of the ROA are given as follows. When a train enters an uphill section, it will search for an optimal switching point of the traction condition in advance. So that when the train pulls out of the uphill section, the running speed will not be much lower than the cruise speed. When a train enters a downhill section, it will search for an optimal switching point of the coast condition in advance. So that when the train pulls out of the downhill section, the running speed will not be much higher than the cruise speed. The comparison of different traction switching points for the uphill section are shown in the Figure 5. The comparison of different coast switching points for the downhill section are shown in the Figure 6. The $v_{pro}$ is the limit speed.

Combined with the optimization objective function, there exists an optimal position of the switching point, which makes the objective function minimum. Therefore, the trisection method is used to find the optimal position of switching points in uphill and downhill sections, and the specific steps are as follows.

Step 1: Set the boundaries of the switching point $[S{P}_l,S{P}_r]$ as the start and end points of the gentle gradient section. Calculate the corresponding objective functions $[O{F}_l,O{F}_r]$ of the switching point boundaries.

Step 2: Divide $[S{P}_l,S{P}_r]$ into three equal parts and compute two trisections $S{P}_{mid1}$ and $S{P}_{mid2}$, where $S{P}_{mid1}$ is close to $S{P}_l$ and $S{P}_{mid2}$ is close to $S{P}_r$. Calculate the corresponding objective function values $O{F}_{mid1}$ and $O{F}_{mid2}$.

Step 3: If $O{F}_{mid1}>O{F}_{mid2}$, update the left array boundary to $S{P}_l$ according to the (6). If $O{F}_{mid1}<O{F}_{mid2}$, update the right array boundary to $S{P}_r$ according to the (7).

$$S{P}_l=S{P}_{mid1}=S{P}_l+{\displaystyle \frac{S{P}_r-S{P}_l}{3}}.$$

(6)

$$S{P}_r=S{P}_{mid2}=S{P}_l+{\displaystyle \frac{2\times (S{P}_r-S{P}_l)}{3}}.$$

(7)

Step 4: Repeat the above steps, and each step will narrow the interval down to two-thirds of its original length, until find the optimal switching point.

2.2.2. Speed Curves Connecting Method

According to the gradient equivalence method, combining three continuous slopes as a group. The ROA is carried out on the basis of this combination to obtain the speed curve of each group. And then the optimized speed curves are connected to obtain the train running speed curve of the whole line. Figure 7 is the speed curves rolling connecting method of some groups.

2.2.3. Cruise Speed Updating Method

Firstly, the initial cruise speed $V_{initial}$ is calculated based on the start and end points of the train section and the planned running time, the calculation formula is shown in

$$V_{initial}={\displaystyle \frac{S_{ea}-S_{sa}}{T_{pa}}},$$

(8)

where $T_{pa}$ is the planned running time in the whole train operating line, $S_{sa}$ is the starting point of the whole train operating line, $S_{ea}$ is the end point of the whole train operating line.

When the optimized train running time does not satisfy the running time constraints, the cruise speed needs to be adjusted with the following rule:

Step 1: Given the difference threshold $T_{op}$ between the optimized train running time and the planned running time, compare the $T_{sa}$ and $T_{pa}$ in the section, where $T_{sa}$ is the simulation running time of the whole train operating line, and let $T=|T_{sa}-T_{pa}|$. If $T<T_{op}$, the running time constraint is satisfied, and then the optimization is finished; Otherwise, the cruise speed needs to be adjusted and re-optimized.

Step 2: Obtain the adjustment value $\Delta V$ of cruise speed according to

$$\Delta V=\left\{\begin{array}{c}\hfill 1,T>2T_{op}\hfill \\ \hfill {\displaystyle \frac{T_{op}}{T}},T_{op}<T\le 2T_{op}\hfill \end{array}\right..$$

(9)

Step 3: When adjusting the speed of the cruise, if $T_{sa}>T_{pa}$, which indicates that the current cruise speed is comparatively low and the train is running slow, so let $V_{new}=V+\Delta V$. If $T_{sa}<T_{pa}$, which indicates that the the current cruise speed is comparatively high and the train is running faster, then let $V_{new}=V-\Delta V$. Where $V_{new}$ is the adjusted cruise speed.

Step 4: Use the $V_{new}$ to re-optimize the speed curve until the running time constraint is satisfied.

2.3. MPC Algorithm for Freight Train Speed Curve Re-Optimization

The principle of the MPC algorithm can be specifically described as follows. At each sampling moment, an open-loop optimization problem is solved in real time for a finite period of time based on the current system state information. Then, the first few groups of data in the generated control sequence are applied to the controlled model as the control output. At the next sample moment, the optimization problem is solved again using the updated system state information as the initial condition to predict the future state of the system repeatedly [25]. Thus, the main advantage of the MPC algorithm is that it balances the effects of ideal optimization and practical uncertainty on the future finite time horizon, because the rolling optimization strategies that always establishes new optimization objectives on a practical basis. This is more practical and effective than the traditional optimal control established under ideal conditions. The principle of the MPC algorithm is shown in Figure 8.

2.3.1. The Multi-Mass Train Model

The multi-mass model takes into account not only the length of the train but also the coupler force, as shown in Figure 9. The coupler force is a direct index for evaluating the longitudinal impulse. For long-distance (heavy-duty) freight trains, the operation regulation recommends that the maximum coupler pulling force should be less than 2000 kN, and the maximum coupler pressing force should be less than 2250 kN [26]. The subsequent constraints on the coupler force will consider this recommendation.

The coupler can be regarded as a spring with damping characteristics, and then the coupler force is calculated by

$$f_i=k_i(x_i-x_{i+1}-l_c)+d_i(v_i-v_{i+1}),$$

(10)

where $k_i$ and $d_i$ are the coefficients related to the elastic and damping force of the coupler, $x_i$ and $v_i$ are the displacement and speed of the ith locomotive/wagon respectively, $l_c$ is the distance of two neighboring couplers when they are in nature state. In the subsequent calculations, we only consider the elastic force to simplify the simulation. Based on the multi-mass model, the force analysis for a single locomotive/wagon is shown in Figure 10.

For the ith locomotive/wagon, $u_i$ is the traction or braking force, $f_{i-1}$ is the front coupler force, $f_i$ is the latter coupler force, and $m_i$ is the mass of the ith locomotive/wagon.

For the first locomotive, it has no coupler force in front of it, and for the last wagon, it has no coupler force behind it. Based on the force analysis of the locomotive and wagons, the whole train motion equation is defined by

$$\left\{\begin{array}{cc}\hfill m_1{\ddot{x}}_1=\hfill & \hfill u_1-f_1-w_1\hfill \\ & \hfill ⋮\hfill \\ \hfill m_i{\ddot{x}}_i=\hfill & \hfill u_i-f_i+f_{i-1}-w_i\hfill \\ & \hfill ⋮\hfill \\ \hfill m_n{\ddot{x}}_n=\hfill & \hfill u_n+f_{n-1}-w_n\hfill \end{array}\right.,$$

(11)

where $w_i$ is the total resistance force of the ith locomotive/wagon, which is

$$w_i=w_0+w_g,$$

(12)

where $w_g$ is the additional resistance on slopes. The basic resistance force $w_0$ is defend by

$$w_0=a+bv+c{v}^2.$$

(13)

In order to facilitate the subsequent design of the control algorithm, the train multi-mass model is simplified. Considering the train basic resistance force $w_0$ and unit slope additional resistance force $w_g$, define $\delta u_i=(a_i+c_i{v}_i^2)m_i+w_{g_i}(2\le i\le n-1)$ and ${\overline{u}}_i=u_i-\delta u_i$. Define $x={[v_1,v_2,\dots v_n,x_1,x_2,\dots x_n]}^T$ as the state vectors and $\overline{u}={[{\overline{u}}_1,{\overline{u}}_2,\dots {\overline{u}}_n]}^T$ as the decision vectors, and then the state equation of the linearized system is

$$\dot{x}=Ax+B\overline{u},$$

(14)

where

$$A=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ I_{n\times n}& 0_{n\times n}\end{array}\right],$$

$$A_{11}=\left[\begin{array}{cccc}-b_1& 0& \cdots & 0\\ 0& -b_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & -b_n\end{array}\right],$$

$$A_{12}=\left[\begin{array}{ccccccc}-{\displaystyle \frac{k}{m_1}}& {\displaystyle \frac{k}{m_1}}& 0& \cdots & 0& 0& 0\\ {\displaystyle \frac{k}{m_2}}& -{\displaystyle \frac{2k}{m_2}}& {\displaystyle \frac{k}{m_2}}& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & {\displaystyle \frac{k}{m_{n-1}}}& -{\displaystyle \frac{2k}{m_{n-1}}}& {\displaystyle \frac{k}{m_{n-1}}}\\ 0& 0& 0& \cdots & 0& {\displaystyle \frac{k}{m_n}}& -{\displaystyle \frac{k}{m_n}}\end{array}\right],$$

$$B=\left[\begin{array}{c}{B}_1\\ 0_{n\times n}\end{array}\right],$$

$$B_1=\left[\begin{array}{cccc}{\displaystyle \frac{1}{m_1}}& 0& \cdots & 0\\ 0& {\displaystyle \frac{1}{m_2}}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\displaystyle \frac{1}{m_n}}\end{array}\right].$$

The model after linearization is discretized in the subsequent simulation application, and the discretized model can be expressed as

$$x(k_c+1|k_c)=Ax\left(k_c\right)+B\overline{u}\left(k_c\right).$$

(15)

2.3.2. MPC Framework

According to the state Equation (15) of the multi-mass model of the freight train, assuming that the train operation state variable $x\left(k_c\right)$ is measurable at the moment $k_c$, the derivation of the state equation leads to the predicted state variable $\widehat{x}(k_c+i|k_c)$ is calculated as follows:

$$\begin{array}{cc}\hfill \widehat{x}(k_c+1|k_c)& \hfill =Ax\left(k_c\right)+B\overline{u}\left(k_c\right)\hfill \\ \hfill \widehat{x}(k_c+2|k_c)& \hfill =A\widehat{x}(k_c+1|k_c)+B\overline{u}(k_c+1|k_c)\hfill \\ \hfill & \hfill =A^{2}x\left(k_c\right)+AB\overline{u}\left(k_c\right)+B\overline{u}(k_c+1)\hfill \\ \hfill & \hfill \cdots \hfill \\ \hfill \widehat{x}(k_c+N_p|k_c)& \hfill =A^{N_p}x\left(k_c\right)+A^{N_{p-1}}B\overline{u}\left(k_c\right)+\hfill \\ & \hfill A^{N_{p-2}}B\overline{u}(k_c+1)+\cdots +\hfill \\ & \hfill A^{N_p-N_c}B\overline{u}(k_c+N_c-1)\hfill \end{array}.$$

(16)

Generalizing the above equation yields

$$X=Gx\left(k_c\right)+\Gamma U,$$

(17)

where

$$X={\left[\begin{array}{c}\widehat{x}(k_c+1|k_c),\widehat{x}(k_c+2|k_c),\cdots ,\widehat{x}(k_c+N_p|k_c)\end{array}\right]}^T,$$

$$U={\left[\begin{array}{c}\overline{u}\left(k_c\right),\overline{u}(k_c+1),\cdots ,\overline{u}(k_c+N_p)\end{array}\right]}^T,$$

$$G={\left[\begin{array}{c}A,A^2,\cdots ,A^{N_p}\end{array}\right]}^T,$$

$$\Gamma =\left[\begin{array}{cccc}B& 0& \cdots & 0\\ AB& B& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ A^{N_p-1}B& A^{N_p-2}B& \cdots & A^{N_p-N_c}B\end{array}\right],$$

where X is the train state at each predicted moment, U is the corresponding control input, and G, $\Gamma$ are the coefficient matrices. The predicted state X can be predicted by (17), and the control sequence U within the prediction range can be obtained.

2.3.3. Optimization Objective Function

In freight train operation, energy saving, punctuality and coupler force are the key factors to ensure the safe operation of the train, the operation optimization objective function of the freight train is constructed by

$$J={\int }_{t_0}^{t_1}\left[\sum _{i=1}^{n-1}{k}_a{f}_i^2+\sum _{i=1}^n{k}_b{u}_i^2+\sum _{i=1}^n{k}_d{\left(v_i-v_r\right)}^2\right]dt,$$

(18)

where $v_r$ is the reference speed, $k_a$, $k_b$, $k_d$ are the weighting coefficients for the coupler force, energy consumption and speed tracking error respectively, and $t_0$, $t_1$ represent the start and end time of prediction period respectively.

Based on (14), the above objective function is transformed into

$$J={\int }_{t_0}^{t_1}\left(x^{T}Lx+{\overline{u}}^{T}P\overline{u}+Q_1^{T}x+Q_2^T\overline{u}\right)dt\to min,$$

(19)

where

$$L=\left[\begin{array}{cc}{L}_{11}& 0\\ 0& L_{22}\end{array}\right],$$

$$L_{11}=diag\left\{k_d,\cdots ,k_d\right\},$$

$$L_{22}=\left[\begin{array}{ccccc}{k}_a{k}_1^2& -k_a& \cdots & 0& 0\\ -k_a& 2k_a{k}_2^2& -k_a& \cdots & 0\\ 0& -k_a& 2k_a{k}_3^2& -k_a& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \\ 0& 0& \cdots & -k_a& k_a{k}_{n-1}^2\end{array}\right],$$

$$P=k_b{I}_{n\times n},$$

$$Q_1^T=-2k_d{v}_r\left[1,\cdots ,1,0,\cdots ,0\right],$$

$$Q_2^T=2k_b\left[\delta u_1,\cdots ,\delta u_n\right].$$

$Q_1$ and $Q_2$ are related to the reference velocity $v_r$, which is determined based on the train position, while $\delta u_i$ depends on the train position and speed, so a computational method for handling these variables based on the train position and speed is given. Firstly, the speed in the predicted state $\widehat{x}(k_c+i|k_c)$ is denoted as $v(k_c+i|k_c)$, and the position of the train at the sample moment $k_c$ is denoted as $s\left(k_c\right)$. Then the position of the train at the sample moment $k_c+i$ can be expressed as

$$s(k_c+i|k_c)=s\left(k_c\right)+\sum _{j=1}^{i}T\left[{\displaystyle \frac{1}{2}}(v(k_c+j-1|k_c)+v(k_c+j|k_c))\right],$$

(20)

where T is the sample time.

From the above equation, the predictions of $Q_1$ and $Q_2$ can be expressed respectively as a function of s by

$$\begin{array}{c}\hfill \overline{Q_1}\left(s\right)={\left[Q_1{\left(k\right)}^T,\cdots ,Q_1{(k+N_p)}^T\right]}^T\\ \hfill \overline{Q_2}\left(s\right)={\left[Q_2{\left(k\right)}^T,\cdots ,Q_2{(k+N_p)}^T\right]}^T\end{array}.$$

(21)

Therefore, the train operation optimization objective function under the MPC framework can be expressed as

$$J=U^{T}HU+2U^{T}f\to min,$$

(22)

where

$$H={\Gamma }^T\overline{L}\Gamma +\overline{P},$$

$$f={\Gamma }^T\overline{L}Gx\left(k_c\right)+{\displaystyle \frac{1}{2}}({\Gamma }^T\overline{Q_1}+\overline{Q_2}),$$

$$\overline{L}=diag\left\{L,\cdots ,L\right\},$$

$$\overline{P}=diag\left\{P,\cdots ,P\right\}.$$

2.3.4. Train Operational Constraints

The freight train operation should consider the constraints of the coupler force, the traction and braking force of the locomotive. Therefore, the constraints of the freight train are given as follows:

$$\begin{array}{cc}\hfill f^{min}\le f_{i{n}_k}\le f^{max},& \hfill k=1,2,\cdots ,n-1\hfill \\ \hfill u_i^{min}\le \overline{u_i}\le u_i^{max},& \hfill i=1,2,\cdots ,n\hfill \end{array},$$

(23)

where $f_{i{n}_k}$ is the kth locomotive/wagon coupler force, $f^{min}$ and $f^{max}$ are lower and upper bounds of coupler force, $\overline{u_i}$ is the ith traction and braking force of locomotive/wagon, $u_i^{min}$ and $u_i^{max}$ are lower and upper bounds of ith traction and braking force of locomotive/wagon.

The above constraints can be expressed in terms of state quantities:

$$X^{min}\le X_{in}\le X^{max},$$

(24)

$$\left[\begin{array}{c}{I}_{(n{N}_c\times n{N}_c)}\\ -I_{(n{N}_c\times n{N}_c)}\end{array}\right]U\le \left[\begin{array}{c}{U}_{(n{N}_c\times 1)}^{max}\\ -U_{(n{N}_c\times 1)}^{min}\end{array}\right].$$

(25)

$X^{max}$, $X^{min}$ and $U^{max}$, $U^{min}$ are upper and lower bounds of $X_{in}$ and U respectively.

Firstly, according to the calculation method of train coupler force, a matrix Z is constructed to extract the relative displacement between neighboring compartments in the state variable, where $Z={\left[\begin{array}{cc}{0}_{(n-1)\times n}& I_{(n-1)\times (n-1)}\end{array}\right]}_{(n-1)\times (2n-1)}$, Then $x_{in}\left(k\right)=Zx\left(k\right)$. Secondly, according to the derivation process in (16), the same reasoning can be used to express the calculation of the coupler force as

$$X_{in}=G_{x}x\left(k_c\right)+{\Gamma }_{x}U,$$

(26)

where

$$G_x={\left[\begin{array}{c}ZA,Z{A}^2,Z{A}^3,\cdots ,Z{A}^{N_c}\end{array}\right]}^T,$$

$${\Gamma }_x=\left[\begin{array}{cccc}ZB& 0& \cdots & 0\\ ZAB& ZB& \cdots & 0\\ Z{A}^{2}B& ZAB& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ Z{A}^{N_c-1}B& Z{A}^{N_c-2}B& \cdots & ZB\end{array}\right].$$

Thus, the constraints on the coupler force can be expressed as follows:

$$\left[\begin{array}{c}{\Gamma }_x\\ -{\Gamma }_x\end{array}\right]U\le \left[\begin{array}{c}{X}^{max}-G_{x}x\left(k_c\right)\\ -X^{min}+G_{x}x\left(k_c\right)\end{array}\right].$$

(27)

By combining the train control force constraint (25) and the train coupler force constraint (27), the constraints for this optimization function can be obtained

$$RU\le \sigma ,$$

(28)

where

$$R=\left[\begin{array}{c}{I}_{(n{N}_c\times n{N}_c)}\\ -I_{(n{N}_c\times n{N}_c)}\\ {\Gamma }_x\\ -{\Gamma }_x\end{array}\right],$$

$$\sigma =\left[\begin{array}{c}{U}_{(n{N}_c\times 1)}^{max}\\ -U_{(n{N}_c\times 1)}^{min}\\ X^{max}-{\Gamma }_{x}x\left(k_c\right)\\ -X^{min}+{\Gamma }_{x}x\left(k_c\right)\end{array}\right].$$

Finally, the overall optimization problem is formulated by the optimization objective function (22) and the constraints (28). This is a typical Quadratic Programming problem, which will be solved subsequently using the active set method.

3. Simulation Results

The simulation data in this paper, such as gradient and speed limit, are all the actual data of Dalailong railway in China. The simulated train is composed of a ${\mathrm{DF}}_{8\mathrm{B}}$ diesel locomotive and 42 C80 wagons, with a total length of 525.1 m and a total weight of 4338 t. Part of the railway from Longkougang Station to Changyibei Station is considered in the simulation, with a length of 133.4 km. This simulation was performed in Matlab.

Table 2 gives the train operation results of the MPC algorithm with different objective function weighting coefficients, while the prediction step and optimization step of the MPC algorithm are fixed $(N_p=6,N_c=1)$ when changing the weighting coefficients of the objective function. In the Table 2, $f_{max}$ and $f_{min}$ denote the maximum pulling coupler force and the maximum pressing coupler force for all locomotive/wagon during simulation, respectively, and $\Delta T$ denotes the operation time difference between the MPC algorithm and the ROA. Figure 11 and Figure 12 give the maximum pulling coupler force and the maximum pressing coupler force for all locomotive/wagon during simulation with different weighting coefficients. When $k_a$ = 0.1, $k_b$ = 10, $k_d$ = 10000, $\Delta T$ is 20 s, the maximum pulling and pressing coupler forces are within the recommended range, and the train speed is within the limit. Concurrently, the maximum pulling and pressing coupler forces, under the other objective function weighting coefficients, fall within the recommended limits.

Figure 13 and Figure 14 give the Two-step Optimization Method simulation speed curves for different values of $k_a$ and $k_d$ when $k_b$ = 10, $N_p$ = 6 and $N_c$ = 1. As can be seen from Figure 13, the simulation speed curves with different objective function weighting coefficients have little change, and the speed curves are smoother compared to the ROA simulation curve. From Figure 14, it can be seen that when the train is operated through the speed limit section, none of them exceeds the speed limit value. The coupler force becomes larger as $k_a$ decreases or $k_d$ increases.

Table 3 gives the train operation results of the MPC algorithm with different prediction step($N_p$) and optimization step($N_c$), when the weight coefficients of the objective function are fixed to $k_a$ = 0.1, $k_b$ = 10, $k_d$ = 10000. Figure 15 gives the maximum pulling coupler force and the maximum pressing coupler force for all locomotive/wagon during simulation with different $N_p$ and $N_c$. The optimization result is $N_p$ = 6, $N_c$ = 1.

Figure 16 and Figure 17 give the Two-step Optimization Method simulation speed curves for different values of $N_p$ and $N_c$ when $k_a$ = 0.1, $k_b$ = 10, $k_d$ = 10000. From Figure 16, it can be seen that when the values of $N_p$ and $N_c$ are changed, there are some differences in the generated speed curves, but they are smoother compared to the ROA simulation curve. And it can be seen from Figure 17 that when the train is operated through the speed limit section, the speed limit value is not exceeded whatever the values of $N_p$ and $N_c$. When changing the value of $N_p$ and $N_c$, the train speed curves change obviously, while the change in coupler force is not significant. When $N_p$ increases from 3 to 6 with fixed $N_c$, the simulation speed curves closer to ROA simulation curve and more in line with speed limits. Meanwhile, with the decrease of $N_c$ while the $N_p$ is fixed, the simulation speed gradually changes towards the speed limit range, and when $N_p$ = 6, $N_c$ = 1 the simulation results achieve the reasonable performance. Therefore, the train will have better operation with longer prediction step and shorter optimization step.

Summarizing the above simulation results, when $k_a$ = 0.1, $k_b$ = 10, $k_d$ = 10000, $N_p$ = 6 and $N_c$ = 1, the maximum pulling coupler force is 34.4% of the pulling coupler force limit value, and the maximum pressing coupler force is 35.7% of the pressing coupler force limit value. It can be shown that the optimized train operation is smoother.

Figure 16. Simulation speed curves with different $N_p$ and $N_c$

Figure 16. Simulation speed curves with different $N_p$ and $N_c$

Figure 17. Partially enlarged simulation speed curves with different $N_p$ and $N_c$

Figure 17. Partially enlarged simulation speed curves with different $N_p$ and $N_c$

Figure 18 gives the simulation results when $k_a$ = 0.1, $k_b$ = 10, $k_d$ = 10000, $N_p$ = 6 and $N_c$ = 1. In the first subfigure, the red line is the speed limit curve, the purple line is the reference speed curve, the green line is the actual speed curve, and the black line is the simulation curve of MPC algorithm. The actual speed curve is plotted from the data in the LKJ(Train Operation Monitoring Equipment). In the second subfigure, the blue line represents the maximum pulling coupler force and the red line represents the maximum pressing coupler force in all locomotive/wagons at each sampling time. The third subfigure represents the line gradient curve. The horizontal coordinate indicates the distance. From the first subfigure, it can be seen that the speed curve obtained from the Two-step Optimization Method has a smaller range of speed fluctuations and is smoother. In the second subfigure, the coupler force variations of the MPC method are small and within safe limits.

Figure 19 gives the simulation results of train control strategy. The first subfigure represents the train control strategy of the MPC algorithm. The second subfigure represents that of the Rolling Optimization Algorithm. The third subfigure represents the gradient curve. The vertical coordinate values in the first two subfigures represent the throttles, where $Gear=0$ means that the train has neither traction nor braking force, $Gear>0$ means that the train is in traction, and $Gear<0$ means that the train is braking. The horizontal coordinate indicates the distance. It can be seen that the control strategies generated by the MPC algorithm and the Rolling Optimization Algorithm are similar, and when the train travels between 60 km and 72 km, the throttles generated by the MPC algorithm change less. Hence, the control strategy generated after the Two-step Optimization Method is easier to follow for the drivers during actual operation.

Figure 18. Train speed curves and coupler force comparisons

Figure 18. Train speed curves and coupler force comparisons

Figure 19. Control strategy comparison chart

Figure 19. Control strategy comparison chart

The train fuel consumption based on ROA simulation is 8.09 kg/10⁴ t·km, the train fuel consumption based on MPC algorithm simulation is 6.71 kg/10⁴ t·km, and the actual train fuel consumption calculated by LKJ data is 8.49 kg/10⁴ t·km. Therefore, the Two-step Optimization Method can reduce energy consumption by 20.9% and 17.1% compared to actual operation and ROA simulation respectively.

4. Discussion

In this paper, a Two-step Optimization Method of freight train speed curve is proposed. The Rolling Optimization Algorithm is used to optimize the speed curve, and the optimized speed curve is used as the reference curve of MPC algorithm for the second optimization step. The ROA takes into account the gradient equivalence method and the single-mass train model, while the MPC algorithm employs a multi-mass train model to consider in-train forces and establish an objective function. The optimal train operation results formed by different parameters are calculated by adjusting the weight coefficients of the objective function and the prediction step and optimization step of the MPC algorithm.

Our simulation results, based on data from China’s Dalailong railway, demonstrate that the speed curves generated by the Two-step Optimization Method exhibit superior smoothness and are more aligned with actual operational constraints. Compared to actual operations, this method achieves a remarkable 20.9% reduction in energy consumption, and a 17.1% reduction when compared to the ROA simulation. The optimized throttle sequences require less frequent adjustments, providing a more practical guide for drivers during actual operations.

Looking ahead, we plan to further test and refine the train suggested speed curve optimization program based on the outcomes of the Two-step Optimization Method. The ultimate goal is to integrate this optimization algorithm into the Driver Advisory System (DAS) to offer drivers a more effective reference, thereby reducing train energy consumption and enhancing the overall efficiency of freight train operations.