Abstract
Automatic heavy-haul train (HHT) operation technology has recently received considerable attention in the field of rail transportation. In this paper, a discrete-time-based mathematical formulation is proposed to address the speed profile optimization problem in order to ensure the safe, efficient, and economical operation of heavy-haul trains (HHTs). Due to the presence of long and steep downgrades (LSDs) on some heavy-haul lines, the brake forces of the HHT are typically jointly determined by air braking and electric braking. The time characteristics of the air braking, such as the command delay and the change process caused by the air pressure, are taken into account, and then formulas are presented to calculate the air brake force. In addition, the influence of the neutral section on the control of the electric braking is considered via space-based constraints. The resulting problem is a nonlinear optimal control problem. To achieve linearization, auxiliary 0-1 binary variables and the big-M approach are introduced to transform the nonlinear constraints regarding slope, curve, neutral section, air brake force, and air-filled time into linear constraints. Moreover, piecewise affine (PWA) functions are used to approximate the basic resistance of the HHT. Finally, a mixed integer linear programming (MILP) model is developed, which can be solved by CPLEX. The experiments are carried out using data from a heavy-haul railway line in China, and the results show that the proposed approach is effective and flexible.
1. Introduction
Heavy-haul transportation is of great importance for improving transport efficiency, increasing economic profits, and lowering transport costs. It has developed rapidly in many countries. Specifically, the axle load of heavy-haul trains (HHTs) has increased significantly to improve the transport capacity, which in turn presents higher requirements for the drivers of HHTs. Moreover, the increase in weight and size of HHTs is also limited by the infrastructure and the hauling capacity of the locomotives. In particular, the application of the automatic train operation (ATO) system for HHTs [1–3] has attracted increasing attention from researchers in recent years, compared with the traditional manual operation.
With the increasing demand for iron ore year by year in Australia, the AutoHaul system [2] developed by Rio Tinto has been applied to achieve autonomous heavy-haul train (HHT) operation in 2018. There are three main benefits to this endeavour in preference to manual operation. The first is to increase productivity by shortening cycle times and eliminating stops for driver changes. Furthermore, operating costs would be reduced with fewer drivers and lower fuel consumption. The third is to reduce the possibility of driver error, thereby improving the safety of HHT operations. However, the implementation of the ATO system in Chinese HHTs is a challenging open issue due to the differences in HHT speed control methods between Australia and China. Specifically, HHTs in Australian railways have achieved automatic operation through the electronically controlled pneumatic (ECP) braking mode, which allows all wagons to receive control commands consistently [4]. HHTs in China are still equipped with conventional air brake systems, where control commands are delivered with a transmission delay depending on the speed of the air wave and the length of the HHT.
Although there has been much research on the application of ATO in Chinese urban rail transit [5–7]. Researchers have mainly focused on the following two critical issues: the speed profile optimization that can be achieved off-line and the train speed controller design with real-time performance. In order to achieve improved performance of train operation, there are still many challenges in the operation of HHTs. Passenger trains, such as urban rail trains and high-speed trains, are mainly controlled by traction and electric braking, which has higher control accuracy and faster response time than air braking; in particular, air braking is only applied when the electric brake force is too small to stop the train at low speed without considering the effect on the next regime. For HHTs, due to their large tractive mass and many sections with long and steep downgrades (LSDs), air braking and electric braking should be applied collaboratively to keep the speed within the given limits, i.e., speed limit and minimum release speed. In particular, we note that the air braking should be applied intermittently for the operation of HHTs on LSDs, as shown in Figure 1. This strategy is also known as the cyclic air braking, where the driver should consider when to apply or release the air brakes and the amount of pressure reduction, taking into account the effect on the next brake application or release, until the trains have left the LSDs [8]. In addition, electric braking prevents the sharp increase in the speed of HHTs and makes the air-filled times of the train pipes more sufficient. It is noted that the electric brake force should be equal to zero when HHTs are running in a neutral section.
Train speed profile optimization [9, 10] is commonly formulated as an optimal control problem aimed at enhancing train operation performance by optimizing the outputs like regime sequences and switching points. This study centers on the optimization of strategies aimed at ensuring the safe, efficient, and economical operation of HHTs on LSDs. While various approaches like pontryagin maximum principle (PMP) [11], approximate dynamic programming (ADP) [12, 13], quadratic programming (QP) [14, 15], mixed integer linear programming (MILP) [16], and heuristic algorithms (HAs) [17–19] have been used to develop HHT operation strategies, there exists a notable gap in the discussions pertaining to critical aspects. Specifically, the time characteristics of air braking, and a constraint within engineering applications where power supply is interrupted while travelling through the neutral section, have not been adequately addressed. Furthermore, the cycle air braking strategy was not taken into consideration in [14, 15], warranting further investigation. In addition, the complexity of line conditions makes it challenging for the PMP to yield optimal solutions, while the results obtained from both ADP and HAs are approximate solutions. On the other hand, mathematical programming methods like QP and MILP benefit from mature commercial solvers that ensure the identification of globally optimal solutions [20, 21]. This study was conducted building upon the aforementioned discussion. The main contributions of this paper are summarized as follows:(1)A formula is developed to calculate air brake force, considering time characteristics for the first time. Subsequently, a discrete-time-based MILP model is presented to optimize the HHT speed profile. In this model, the total running time is divided into multiple intervals.(2)Real-world line conditions and operational constraints, including slopes, neutral sections, and air-filled times, are considered. To facilitate approximate linearization, the big-M approach and binary variables are introduced.(3)A hybrid scheme that integrates coarse-grained and fine-grained models is devised to strike a balance between computation efficiency and accuracy. Then this complex challenge is tackled by using the well-known solver CPLEX. Moreover, a series of experiments are conducted to demonstrate the solutions achieved.
The remainder of this paper is organized as follows. Section 2 delves into the existing literature concerning the optimization of HHT operation strategies. Section 3 outlines the model of the HHT running on LSDs integrating a neutral section. Moving forward, the discrete-time-based MILP approach is given in Section 4. The effectiveness and flexibility of the proposed approach are validated through simulation experiments in Section 5. Finally, conclusions are drawn in Section 6.
2. Literature Review
The optimization of train operation has been aptly described as an optimal control problem, as noted in [5]. Various optimization methods have been employed to improve the operation performance of long HHTs.
In early years, the researchers from the Scheduling and Control Group (SCG) of the University of South Australia developed the FreightMiser system to reduce fuel consumption for long-haul trains, and optimal speed profiles for relatively simple scenarios, such as isolated steep uphill or downhill sections, were obtained using PMP [22]. Subsequently, another PMP-based approach was proposed in [11] for the purpose of energy saving, which considered the operation requirements of a Chinese HHT. However, the pursuit of rigorous mathematical solutions led to model simplification, ignoring factors such as changing slopes and air braking characteristics. With the improvement of computing power, recent years have witnessed a multistage decision process to determine the control sequence for HHTs running on LSDs, where lookup-table-based [12] and value-iteration-based [13] ADP were used to yield safety-centric, cost-effective, and operationally efficient driving strategies. The advent of deep Q-network (DQN) [23] mitigated the challenge of high-dimensional spaces by approximating action-value functions with neural networks. Furthermore, a QP-based approach was developed to address the multiobjective trajectory optimization problem in the context of HHTs. This approach encompassed considerations of key constraints, including line resistance, neutral section, as well as traction and electric braking, as highlighted in [14]. Moreover, a purposeful selection scheme was designed to obtain the optimal weighting factors [15]. Notably, these investigations were conducted on railway scenarios characterized by less intricate line conditions and the absence of LSDs. In addition, HAs find predominant application in resolving the HHT operation strategy on LSDs due to their less stringent model requirements. For instance, to ensure safety and improve efficiency during operations, methods like particle swarm optimization (PSO) [17], genetic algorithm (GA) [18], and artificial bee colony (ABC) [19] have been embraced. These algorithms enable the identification of optimized switching points for air brake application and release. In the realm of optimizing train speed profiles, MILP has established itself as a widely recognized and effective method. Its utility spans across urban rail trains [24–26], high-speed trains [27, 28], and medium-speed maglev trains [29]. Notably, it has recently been employed in the domain of HHT [16], where the linearization description of the air brake force and the air-filled time constraint is not well treated and still full of challenges.
However, in the aforementioned studies, when addressing air braking, the force was calculated using empirical formulas on the assumption of instantaneous characteristics due to the complex working principle of the air brake system, and the delay of air brake application or release, as well as the progression of air brake force alteration were overlooked. Moreover, real constraints such as the neutral section were not given due consideration. To highlight achievements and gaps in existing research, a brief comparison of the most relevant literature is presented in Table 1.
3. Model Formulation
3.1. Assumptions
In this paper, we make the following assumptions:(1)The pressure reduction for brake control typically ranges from 50 to 140 kPa in intervals of 10 kPa. On lines with LSDs, a pressure reduction of 50 kPa can be employed to control the speed of the HHT. It is important to note that routine stops and stops resulting from suboptimal control are not taken into account.(2)The output for electric brake control is expressed as a ratio of the characteristic curve.
3.2. Notations
For a better understanding of the paper, Table 2 lists the notations used in this section.
3.3. Longitudinal Dynamics of the HHT
In this paper, the HHT is regarded as a rigid multimass model when running on the line, as in [30], to reduce unnecessary complexity.
When the HHT with cars runs on LSDs, , is the number of locomotives and is the number of wagons. , is the index of cars and is the mass of car . The continuous-time model of the HHT moving along the track is expressed as follows:where is the position of the HHT, and are the mass and speed of the HHT, is the time of the HHT running, is the electric brake force provided by the locomotives, is the air brake force acting on all cars, and and are the basic and line resistances.
3.3.1. Resistance
The basic resistance [31], related to the speed of the train, can be expressed as follows:where [32], , , and are positive coefficients related to the type of car. denotes the gravity constant. As the number of locomotives is much smaller than that of wagons, it can be assumed that the unit basic resistance of a locomotive is the same as that of a wagon.
The line resistances [31] composed by the slope and curve resistances can be computed as follows:where is the empirical coefficient.
For the position of the car , the gradient is measured by and the curve resistance can be calculated by the curve radius , which can be described by the following piecewise linear functions:where is the position of car , is the physical length between cars and , and ; are the starting and ending positions of slopes, are the gradient of slopes; are the starting and ending positions of curves, and are the radius of curves.
3.3.2. Electric Brake Force
The electric brake force [31] can be formulated as follows:where is the output ratio and is the maximum electric brake force that can be achieved according to the characteristics of a locomotive.
3.3.3. Air Brake Force
The air brake force [31] is generated by compressed air, and it can be calculated as follows:where for car , is the pressure acting on each brake shoe, and is the brake shoe friction coefficient. is the number of brake shoes equipped on a car.where is the diameter of the brake cylinder, is the transmission efficiency of the braking device, is the braking leverage, and is the number of brake cylinders equipped on a car.
In equation (7a), is a constant for a given brake system. Therefore, the description of the brake cylinder pressure is the key to calculating the air brake force. Based on the function proposed in [33] and test data, the mathematical formulations are as follows:where and are the pressure of the brake cylinder of car during air brake application and release, and are the binary variables indicating the application and release of the air brake, denotes the amount of pressure reduction, and are the timing of the air brake application and release, and and are the delay times for car at which the pressure begins to rise and fall, depending on the propagation speed of the wave. and are the times required for the pressure to change to the maximum value and zero, and and are the functions that describe the processes of air brake application (pressure rise) and air brake release (pressure fall) for each car.
3.4. Optimization Problem
3.4.1. Objective
Given a specific running time, the pursuit of enhanced operation efficiency necessitates a higher average speed, thereby maximizing the covered distance. Besides, to mitigate the wear of brake shoes caused by elevated temperatures and reduce maintenance cost, the duration of air brake application must be minimized.
The running distance and the air brake application time can be calculated by
The objective function can be written as follows:where and are weight coefficients, . is the running distance when the HHT runs at maximum speed (see below “speed limit”). is the given running time of the HHT.
3.4.2. Constraints
Train dynamic constraint: The movement of the HHT should follow equations (1a) and (1b).
Speed constraints: For safety, trains cannot run above the speed limit or below the minimum release speed . When the speed is lower than , the in-train force will increase and the coupler may even break.
Electric brake constraint: Generally, , that is, the output of the electric brake force should not be greater than the maximum electric brake force. So in this paper, the following set should be satisfied.
Air brake constraints: “” denotes application and “” denotes release, the relationship between them can be described as follows:
The rated pressure of the brake pipe is set at 500 or 600 kPa, with a chosen output of 50 kPa when the HHT is navigating on lines with LSDs.
Besides, to allow sufficient time for the train pipe to be filled, a release interval between adjacent air brake applications must be at least as long as the air-filled time.where is the air brake release time, and is the air-filled time.
Neutral section constraint: For the neutral section of the electric railway, the output of the electric brake force equates to zero, which can be rewritten as follows:where and are the starting position and ending position of a neutral section, respectively.
4. Solution Approach
The given running time is divided into intervals. , and in each interval , values of the speed limit, the slope resistance, and the curve resistance are constant; we also assume that the traction or braking force is taken as constant. Note that and , and are the beginning and ending times.
4.1. Decision Variables
Table 3 presents the decision variables for optimizing the speed profile of the HHT.
4.2. MILP Approach
Inspired by [24], the optimization problem described in Section 3 was solved using the MILP approach. Specifically, the trajectory optimization problem of passenger trains was studied in [24], and a discrete-space-based model was obtained, with kinetic energy per unit mass and time as state variables, and position as the independent variable. In this paper, a discrete-time-based model has been established, with speed and position as state variables, and time as the independent variable. In addition, the modelling of the air brake force and the air-filled time of the train pipe are the focus, and the neutral section, the gradient and the curve radius are all piecewise functions related to the position. Binary variables must therefore be introduced to achieve linearization.
PWA functions are used to realize the linearization, and the nonlinear function , equation (2), can be written as a piecewise linear function of ,where the values of , , , , and are determined by the fitting process.
Besides, , , and can be written as the form: , , and . So defining , , , and , the model in equation (1a) can be rewritten as follows:
Then, the zero-order holder and the trapezoidal integration rule are used to obtain the discrete-time dynamic model:where , , and ; , , and are the initial speeds of the HHT.
According to the properties introduced in [34], the PWA model in equation (19a) can be transformed as follows:
Here, the auxiliary logical variables , , and are defined to transform the model
Therefore, the following linear constraints should be satisfied:where and are the minimum and maximum electric brake force, and are the minimum and maximum air brake force, and are the minimum and maximum line resistance, and notably, is a positive small number.
For equations (4a) and (4b), binary variables are introduced to describe these piecewise linear functions:where is a positive large number, , .
Based on equations (6)–(8c), the total air brake force of the HHT can be obtained and expressed as a function of time for the specific pressure reduction.where and are the functions of air brake force and time during brake application and release, respectively. Here, when , the value range of is from 0.75 to 0.84, which is not large. So, in equation (7b), the average speed is used for the calculation. Then, we define and ; equation (24) can be rewritten as . Next, general recursive expressions are proposed to calculate the air brake force, taking into account the time characteristics and the effects of historical regimes:where and are number of the selected values used in force calculation for air brake application and release. , , , , , ; , , , , , ; , .
Take air brake application for example, as shown in Figure 2, , : if , ; if , ; if , ; if , .
Here, , , . Also, these logical conditions should be rewritten as the following linear constraints:
For equation (11), that is,where and can be calculated based on the minimum release speed and the speed limit . Two binary variables and are also introduced to indicate the timing of air brake application and release, then ensure that the release time is not less than the air-filled time given in equation (15),
Thus, logical condition, , is expressed as follows:where is a positive large number, which is introduced to ensure when , , and . In addition, the binary variable is introduced to define the neutral section,Then equation (5) can be rewritten as follows:
Also, three binary variables are introduced to transform equation (30),where is a positive large number. For equation (31), defining , and the following linear constraints should be satisfied:where and are the minimum and maximum values determined by the control rules.
Finally, the objective function in equation (10) can be rewritten as follows:
The constraints consist of those related to train dynamics (19b), (20), (22a), (22b), and (22c), line slope and curve (23a) and (23b), air brake force (24), (25a), (25b), and (26), and train operation (27), (28a), (28b), and (29)−(33).
Above all, the speed profile optimization problem of the HHT can be transformed into a MILP model. In this paper, the MILP model is solved using CPLEX.
4.3. Balancing the Efficiency and Accuracy
However, the selection of the time interval significantly affects the computation efficiency and accuracy. Smaller intervals enhance solution accuracy at the cost of increased computation time, while larger intervals yield the reverse effect. Specifically, smaller intervals could enable a more detailed consideration of air braking characteristics. As illustrated in Figure 2, assuming , additional data points from the air brake force curve are employed for modelling. Therefore, this paper proposes a hybrid scheme that combines coarse-grained and fine-grained models to optimize operation strategies. Firstly, the speed profile of HHT running is acquired using a larger time interval, constituting a coarse-grained model. Then, leveraging the calculated solution, a fine-grained model is constructed with a shorter time interval. This sequential refinement leads to improved solutions without significantly prolonging computational time.
Assuming the initial solution with is . Then, we define a search scope, , and are positive integers. So when , it is known that and . Finally, the solution will be obtained.
5. Case Study
In this section, we employ a 10, 988-tonne HHT with 116 wagons as the subject of our numerical experiments. The partial line information of the actual railway is shown in Figure 3. Other relevant details are given in Tables 4 and 5. The experiments are then carried out on a computer with 3.20 GHz Intel i7 CPU and 16 GB RAM, and CPLEX 12.9.0 is used to solve the proposed model. Also, the computation time is denoted by .
Remark 1. For the sake of simplicity, we define two schemes for the following discussion. Scheme I entails the use of a coarse-grained model or a fine-grained model, while Scheme II involves the combination of coarse-grained and fine-grained models.
Remark 2. Due to the gradual variation in slope gradient and the limited impact of curve resistance, we opt to calculate the line resistance for the entire HHT based on the gradient and curve radius of the foremost position, the first car.
Remark 3. Our primary focus is on scenarios, where the HHT traverses a neutral section with an ongoing air brake application regime, necessitating adherence to the relevant constraint, i.e., .
We then delve into an analysis of the impact of linearization approximation, presenting the average errors between the approximated and true values of both basic resistance and the air brake force in Table 6. Here, we define that , where and denote the approximate value and the true value, respectively.
From Table 6, it becomes evident that the linearization methods proposed in this paper yield results that closely approximate the true values. The average approximation errors for basic resistance and air brake force are 1% and 8.9%, respectively. Notably, the 8.9% error is obtained at a time interval of 5 seconds during the air brake application, as an illustration.
5.1. Case 1: Effectiveness Verification
In this case, and , and the initial speeds are defined as 13.89 m/s and 19.44 m/s, respectively. Initially, Scheme I is employed to elucidate the impact of the time interval on the HHT operation performance and computation time. In this context, the time intervals are set to 30, 20, 10, and 5 seconds, respectively. Scheme II is then introduced to demonstrate that there is an improvement in efficiency. The experimental results are shown in Tables 7 and 8.
Noticeably, in Scheme I, as the time interval decreases, the number of intervals increases, resulting in performance improvements, including increased running distances and reduced air brake application times. This implies that the use of the fine-grained model in conjunction with MILP leads to improved performance. However, this comes at the cost of escalated computation times, surging from 11.29 and 6.55 seconds to over 5000 seconds. Transitioning to Scheme II brings about 91.59% ( is 13.89 m/s) and 91.67% ( is 19.44 m/s) reduction in computation time compared to Scheme I when the time interval is 10 seconds, while preserving performance levels. In particular, when the time interval is 5 seconds, Scheme I struggles to provide solutions within the 5000-second time; with Scheme II, the average time taken is 515.89 seconds. To strike a balanced compromise between operation performance and computation time, a time interval of 5 seconds proves to be optimal for further discussions in Scheme II.
Moving on, Figure 4 shows the relationship between the number of intervals, the objective function, and the computation time for the two initial speeds. Applying Scheme I with a finely tuned time interval results in a marginal improvement in the objective function relative to the substantial increase in computation time. Also, the application of Scheme II improves this situation.
5.2. Case 2: Flexibility Verification
In this case, to clarify the flexibility of Scheme II, we apply it at a time interval of 5 seconds and investigate the performance of different weighting coefficients.
Table 9 presents the optimization results of variable weight coefficients when the initial speeds are set to 13.89 m/s and 19.44 m/s, the number of intervals . Solutions are successfully obtained, averaging around 604.38 and 541.54 seconds for the respective initial speeds. These values markedly exceed those recorded in the scenario with a time interval of 10 seconds. This escalation can be attributed to the significant increase in computational complexity as the number of intervals augments and the search space broadens. Furthermore, for identical initial speeds, as increases, the air brake application time when changes significantly compared to ; as increases, the running distance is when changes significantly compared to .
Figure 5 illustrates the speed profiles and corresponding forces as aligned with the results presented in Table 9.
In Figure 5, the HHT operates within confines of the speed limit and the minimum release speed. It can be seen that distinct initial regimes manifest with varying initial speeds. However, the influence of the initial speed on the speed profile is progressively mitigated through the coordinated interplay of air braking and electric braking. This trend arises due to the presence of a fixed neutral section around the 20.7 km mark. In accordance with the requirements of operation safety, the HHT control must respect the air brake application regime when traversing the neutral section.
Furthermore, Table 10 provides insights into air brake release times. Notably, a minimum of 180 seconds is allocated for air brake release processes, in compliance with the air-filled time mandate of the train pipe. Noteworthy is the prolonged period required for the initial release of air braking at = 13.89 m/s, serving the purpose of strategic operation adjustment.
On the basis of the above experiments, it is shown that the proposed MILP approach and Scheme II of combining coarse and fine-grained models can address the problem of speed profile optimization well.
6. Conclusions
In this paper, the speed profile optimization problem has been investigated for the safe, efficient, and economical operation of heavy-haul trains (HHTs) on a railway with long and steep downgrades (LSDs). The optimization problem is reformulated as a discrete-time-based mixed integer linear programming (MILP) model, where a set of practical constraints is considered. Specifically, a number of mathematical formulas have been introduced for the air brake application and release process to obtain a more accurate air brake model. Also, for the neutral section, the regime, i.e., air brake application, adapted to the LSDs has been adopted. Moreover, Scheme II, which combines the coarse-grained model and the fine-grained model, has been presented to balance the efficiency and accuracy of the computation. The computational results suggest that the optimal speed profiles of HHTs can be achieved via the mathematical model and solution approach presented.
The proposed method can be used to obtain the speed profile for automatic train operation, which helps to make the operation strategy more realistic, reduce the work intensity of the drivers, and improve the operation speed of HHTs. While the HHT is passing through the neutral section, where the coasting should be followed, we can modify the constraints to ensure the safe and efficient operation.
We focus mainly on the 10, 000-tonne trains with a single locomotive that make up a large proportion in heavy-haul railway. For the application in 20, 000-tonne trains with “1 + 1” formation, which relies on wireless remote multitraction synchronization control, the model in this paper is no longer applicable. In future work, we will extend this research to more complex train formations, such as 30, 000-tonne trains, and we will consider the characteristics of traction, electric braking and air braking, operation rules, and the coupler force that limit train operation. An efficient algorithm will then be developed to solve the optimization model for long train formations.
Data Availability
The data supporting the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Fundamental Research Funds for the Central Universities 2023JBZY017 and Scientific and Technological Innovation Program of Shuohuang Railway under Grant nos. 2019-148 and 2021-234, Beijing Laboratory of Urban Rail Transit.