Multi-Objective Optimization and Weight Selection Method for Heavy Haul Trains Trajectory

Energy-saving, punctuality and smoothness are difficult to achieve in the operation of heavy haul train. This paper proposes a multi-objective method for speed curve optimization. Firstly, a quadratic programming optimization model is established, using train kinetic energy and train force as independent variables and several goals of energy conservation, smoothness, and target speed tracking. Under the supplied weight vector, the model can realize the ideal speed curve under unbounded time. Secondly, in order to address the demand for punctuality, this paper develops a weight selection algorithm to optimal weight vector for any given trip time. Finally, a time-bounded multi-objective speed curve optimization model is proposed, which achieves energy conservation, punctuality and smoothness train operation. The simulation indicates that the proposed algorithm outperforms human driving in multi-objectives. And the designed time bounded weight selection algorithm is more efficient than the conventional linear approaching weight search method.


I. INTRODUCTION
Large and heavy cargo has received much attention as it is cost-effective transportation with reliable security and a high level of carrying capacity [1][2][3]. The most common goal in railway speed curve optimization challenges is to maintain the operation safe. Furthermore, energy conservation and punctuality are two crucial but opposing optimization goals [4]. In addition, a heavy haul train entails many wagons carrying huge loads [5]. Coupler breakage and train derailment may even occur in extreme instances [6]. Thus, the additional goal of smoothness is addressed in the heavy haul train speed optimization. Overall, the heavy haul train speed optimization problem is converted to a multi-objective optimization problem, including safety, energy-saving, punctuality and smoothness.
Heavy haul train has a mass and hence has a large inertia. It also determines that its control force is incapable of mutation. The general environment of heavy haul train is the continuous undulating ramp, resulting in a complex train-slope strongly coupled relationship and a challenging energy-saving control method. Heavy haul railway is responsible for transporting products, and the journey distances are typically longer than standard railway, making it difficult to effectuate punctual operation. The complexity of optimization is much higher than that of other types of trains. Currently, heavy haul train is mainly operated based on operation diagrams completed by the experience of excellent drivers. [7,8] This method cannot find the optimal trajectory that meets multiple objectives and doesn't schedule for the unplanned event. For this reason, realtime planning is also one of the requirements for the highcomplexity planning algorithm of heavy-haul train [9]. Consequently, global optimization of heavy-haul trains is difficult to achieve with an efficient algorithm.
Recent research studies show that the optimization method of heavy haul trains can be divided into the analytical method, heuristic algorithm method, and numerical optimization method [10]. The analytical method is based on Pontryagin's Maximum Principle (PMP), in which the differential equation of the co-state variable is calculated and the relation between optimal control modes and the co-state variable is explicit, which includes full power, partial power, coasting, partial braking, full braking [11][12][13]. The optimal control for a specific journey on a non-steep track is unique [14,15]. But the solving process of the analytical method is a difficulty that hinders his further development [16,17]. The pseudo-spectral method is based on analytical model via a numerical method [18][19][20]. In comparison to the analytical method, it enhances computation efficiency while causing an accept an acceptable loss, and it is stated that this provides an excellent offline solution method. But this method has computational efficiency limitations due to the rationale for the highcomplexity optimization problem of heavy-haul trains. Alternately, the heuristic algorithm is based on historical data to train the control output results. No longer adopt complex heavy haul train models, but the algorithm uses data-driven programming to fit its nonlinear coupled relationship [21][22][23]. Nevertheless, demand for a large amount of support data as a learning sample is tough to accomplish. Besides this method may not be able to find the optimum solution.
The numerical optimization method refers to discretizing the dynamic system into a problem with finite variables, and then directly solving the problem through the mathematical programming method [24]. Many scholars hold the view that the speed trajectory planning of the heavy haul train is associated with quadratic programming problems [25][26][27][28], which equips with mature solver with high computational efficiency [29]. This method obtains an approximate optimal solution at a shorter time scale. Therefore, it has the potential for real-time optimization in engineering.
In the multi-objective optimization problem, it is generally difficult for multiple optimization objectives to be their absolute minimums at the same time. The Pareto optimum reaches the relative optimum under mutual constraints. In essence, multi-objective optimization problems intend to solve Pareto optimum boundary. The general approach of multiobjective optimization algorithms is to convert multiobjectives into single-objective solutions in a certain way. Representative methods are based on aggregation selection, target-vector approaches, and lexicographic ordering, etc. [30][31][32][33].Fuzzy multi-objective method is proposed in hypersonic vehicle [34] engaging with the deep neural network but demands for powerful computing capacity. Ant colony algorithm and particle swarm algorithm are also used in multiobjective optimization problems, in which each solution is treated as an individual to iterate toward the optimal solution [35][36][37]. Due to the characteristic of multi-objective optimization of heavy haul trains, few studies have instigated systematic methods for weight selection, but based on empirical values [38,39]. Multi-objective problems are widely investigated in the operation system of the metro train [40][41][42][43]. The objective optimization function is cogitated as a weighted function of energy consumption and stability and formed mixed-integer linear programming through a piecewise affine function [44]. Model [45] uses a three-objective dynamic algorithm with fuzzy parameters, which detect the delays to balance energy consumption and the risk of delays in arrival.
A longitudinal dynamical model[ [46] and the quadratic programming algorithm [3,47,48]are used to solve the multiobjective optimization problem, which has good real-time performance.
Train operation at particular periods are required to finish within a restricting time due to the constraints of the train control system. As a result, the optimization model should be able to adjust the weight in order to accomplish time-bounded operation. The multi-objective constraint is essentially a broad constraint that does not impose a precise time limitation. There is no obvious weight selection technique in existing research, and it is largely reliant on engineering debugging. Finally, we should examine the distribution trend of the weight space and implement the weight matching algorithm.
In this research, we build a multi-objective trajectory optimization model based on QP and propose a weight selection algorithm for given trip time. Finally, the train's energy-saving, punctuality, and smoothness performance can be realized.
The remainder of this paper is organized as follows. In Section 2-Model Statement, we define the basic heavy haul train motion model, transform it into a quadratic programming problem. In Section 3-Weight Selection Algorithm we analyze the impact of weight vector and propose weight selection methods. Section 4-Simulation Verification presents simulation results, numerical examples with actual data. We conclude this paper in Section 5-Conclusion.

II. MODEL STATEMENT
The train motion model is subjected to traction force, braking force, basic running resistance, and additional line resistance during running. The braking force consists of pneumatic braking and electric braking. These independent braking forces are considered separately for the sake of making distinctions in braking mode. The dynamic model of the train can be described in the following Equation (1).
Where v , s , t are the speed, position and running time of the train. M is the total train mass, and γ is the rotary mass coefficient. Ft(s) is the train traction force, Fd(s) is the train electric braking force and Fm(s) is the train pneumatic braking force, they are all functions determined by s. R(s, v) is the environmental influences determined by train speed v and train position s. In this model, traction and braking force are the decision-making variables, which are scheduled to be optimized.
It can be seen from the dynamic model of the train is composed of continuous, non-linearized strongly coupled differential equations. This form of the equation is not conducive to the solving process and application. This paper discretizes each variable based on distance. The discrete step length is set to Δs, dividing the whole process into N segments. To achieve the linearization of the original problem, kinetic energy Ek is a substitute for velocity vk T 12 2 [ , ,......, ] Simultaneously, all control variables can be stacked into an (N×1) matrix according to the discrete step: The operation safety must be prioritized, which is ensured through speed limit overspeed kinetic energy protection Epro,k for the train.
The kinetic energy flow of the train can be expressed by a linear equation.
The mentioned train model is widely used in various types of trains according to the different parameters. Nevertheless, the composition of heavy haul trains resistance is unique. It should be noticed that the environmental impact R(s, v) consists of the basic running resistance R0, the curve additional resistance Rr, and the ramp additional resistance Rg.
Where r refers to the unit basic resistance. mj is the mass of the j th vehicle. Cj and ij respectively represent the radius of the curve and the slope at the location of the j th vehicle, and C is a constant value of radius. Train operation is constrained by the vehicle characteristics, shown in Figure 1. The maximum values Ft max,k Fd max,k of train control variables are derived from continuous profiles determined by the speed v, which meets the following conditions: The vehicle characteristics constraints are transformed into multiple inequality constraints with kinetic energy Ek. As a result, the train is constrained by inequations of the maximum traction and brake force-fitting curves Ft1, Ft2, Fd1, Fd2, Fd3 determined by kinetic energy and the linear approximation of the original problem is realized.
Furthermore, the concern of heavy-haul train operation includes smoothness of control force, energy consumption and target speed tracking.
(1)The energy consumption of whole train operation is: Where ηt and ηd are the efficiency under traction and electric braking conditions, the constant α is the regenerative energy utilization coefficient.
(2) The target speed tracking item is proposed to maintain at the destined cruising speed. By adjusting the destined kinetic energy Edes,k, the cruising speed of train operation is achieved.
(3)The smoothness of the train during operation is expressed by the force changing rate between adjacent discrete steps, which is expressed as:

A. Application of Quadratic Programming
Based on the mentioned train operation optimization objectives with allocated weights to respective target, the multi-objective optimization function of train operation optimization can be described as below: Where w1, w2 and w3 are the non-negative weight values corresponding to the energy efficiency, target speed tracking and smoothness respectively. This paper uses QP for numerical calculations. Quadratic programming refers to optimization problems in the form of Equation 15.
During the operation of the train, the variables that can be controlled by the driver are traction force Ft, electric braking force Fd, and pneumatic braking force Fm. Additionally, the train speed is the key status indicator, which is represented by kinetic energy E. Therefore, the optimization variable vector follows the form of the following equation.
The objective function Equation 14 is transformed into its equivalent form Equation 17 to meet the requirements of the quadratic programming problem. The highest-order term of the target variable of the equation is quadratic, conforming the matrix in QP as Equation 18.
Then the degree 1 term of the objective function can be expressed by 4N-dimensional vector f. 22 To conform to the equality constraint form of QP Ax=b, the basic resistance (Equation 8) and kinetic energy flow (Equation 6) are expressed in the following way. Basic resistance is the only resistance determined by speed.
The matrix form of the above formulas are represented as follows: After these definitions and derivations, the multi-objective optimization model of the heavy haul train can be transformed into the standard QP form. With the protection of the speed limit, the train can reach the destination consistently, so this QP problem should be solved with a numerical solution. But the discrete step cannot be too large, otherwise, in specific conditions, it may lead to discontinuities, which signifies no solution. Therefore, the discrete step length in this article is set to 100m, relative to the hundreds of kilometers of freight line, which can certainly guarantee the astringency of complicated conditions in heavy haul trains. The QP problems have many mature solvers such as quadratic programming in MATLAB and OSQP [49] which uses a specialized ADMM-based firstorder method with custom sparse linear algebra routines. These provide a sufficient foundation for solving the above problem.

B. Trip Time Calculation
The target speed tracking term, which is a compromise to the QP problem, indirectly reflects the punctuality of heavy haul train operation. Despite the fact that computation efficiency improves, this algorithm is unable to meet the train's punctuality requirement. The most intuitive way to evaluate the speed profile is to use trip time. As a result, the trip time of the QP solution must be investigated.
Equation 26 shows the algorithm for trip time. After solving J through quadratic programming, the corresponding solution variable x can be obtained, which is the data source for calculating trip time t through the derived form of the kinetic energy equation. The significance of the weight value is to balance the three goals, and adjusting the magnitude of the weight value cannot change the optimize the solved solution. There are just two degrees of freedom, despite the fact that there are three objectives. w3 is set to 1 by default for ease of presentation, and variables w1 and w2 can be considered independent variables.
The formula above is given by a function T, where the operation time is a function of the weight vector with QP, to make mathematical representation easier. 12 ( , ) t

A. Linear Approaching Weight Search Method
To realize the punctual trip time of trajectory. The adjustment of the weight is indispensable, and the flowchart below Figure  2 shows a matching method based on linear approaching, in which we adjust the target speed tracking-related weight w2 to implement the on-time requirements. The linear approaches method's major characteristics are built on feedback adjustments depending on each successive simulation outcome. Although the linear approach method is simple to implement in practice, it fails to account for all of the multi-objective influencing factors because on-time adjustment is achieved solely by increasing a single target weight rather than systematically distributing weights among target speed tracking, energy efficiency, and smoothness. Furthermore, because the approach uses the QP solution process and the implicit time calculation function T in each iteration, each iteration requires a lengthy weight matching calculation time.

B. Time Bounded Weight Selection Algorithm
For the sake of time-bounded, the optimization model above needs to consider the impact of different combinations on solution results. J refers to the minimum function of optimization result with the weight vector as the independent variable. After the quadratic programming model is established, the running time, energy consumption and train smoothness can be regarded as functions of w1 and w2.
In effect, the weight vector has complex nonlinear factors affecting multi-objective optimization function, which can be observed from Figure 3. This also leads to the inability to select effective weights using the linear approaching method in the previous chapter. However, each target amount can't be deduced from the overall target distribution. Therefore, it is necessary to analyze the weight distribution of each specific target to realize the ontime weight selection of the speed trajectory. 5×103 kWh and 6.6×10 3 s. As the respective weighting factors increase, the indicator they represent also has positive correlation changes and that they do not represent turns on the contrary. Overall, the change of the variable relative to the energy consumption factor is stable and relatively linear, but the influence on the trip time factor changes relatively very sharply, especially in the case of variable in the range of [0.2×10 -4 , 0.9×10 -4 ], which is caused by a deviation between the trip time and the target speed tracking in the objective function. Figure 4(c) represents the smoothness of the train. The smoothness weight exists as the denominator in substitution variables, and, logically, its optimal region is near the origin. Meanwhile, the distribution of smoothness is largely affected by the route condition. However, the algorithm is unsuitable for solving smoothness-inclined, whose operation sequence violates common sense. Therefore, the weight distribution at the origin should be appropriately discarded. In the range shown in the figure, the variable w2 has a relatively more important impact on train smoothness.
Weight distribution makes it difficult to portray human subjectivity. As a result, determining artificially optimal weights and ideal trajectory is quite challenging. There is no question that a weight selection method should be proposed. Where εt is the maximum allowable time deviation. The weight selection algorithm for given trip time is proposed in Figure 5. First, offline data of indicators with weight are traversed as the independent variable. The initial interval of the weight value set S1 can be delineated by setting the allowable time deviation range, so as to narrow down the feasible weight range. Afterward, we search weight combinations with lower energy by the smallest amount of energy consumption Emin in the feasible area to obtain the set S2. Finally, a feasible solution S3 is obtained through the constraint of smoothness control Js of the train. On account of the feasible solution is based on the minimum value of each solution set in feasible weight range, therefore, this method should have a solution if the target time is reasonable. As the temporal deviation requirement grows, the possible solution set shrinks and approaches a point set.
In the operation of heavy haul trains, standard operation time tsta is the most intuitive and basic constraint. in the flowchart, by setting scheduled trip time, the majority of weight combinations can be filtered in time bounded filtration, framing a smaller traversal range, improving traversal algorithm efficiency. By creating a set of feasible weight vectors (w1, w2), we will search for the optimal weight combination in the reaming weight selection module.

V. SIMULATION VERIFICATION
The practical train HXD1 is used as a simulation model in this paper to verify the feasibility and weight selection methods. The detailed train parameter is shown in Table 1.

A. QP and Pseudo-spectral optimization comparison
The pseudo-spectral method is based on an analytical model. It is claimed that this method gives an outstanding offline solution method for subway and high-speed trains trajectories. This method, as well as the quadratic programming method, can be used for comparative analysis.  As a result, it has been determined that pseudo-spectral is not suitable for heavy haul trains. Under the same train running time limitation, the pseudo-spectral approach can only superimpose a single goal of energy. The pseudo-spectral method has defect in adapting complex ramps, hence the approach stubbornly accepts acceleration, cruising, coasting and braking, similar to the subway train. The pseudo-spectral method outperforms QP in energy optimization. Nevertheless, frequent control force adjustment is nearly hard to execute. At the transition of some ramps, there is also noticeable jitter in the control force.
Some QP optimization features are similarly comparable to the pseudo-spectral, which proves that quadratic programming optimization is beneficial in terms of energy-saving and punctuality. QP's computing efficiency is substantially higher than that of the pseudo-spectral, so it has better application prospects for heavy haul train.

B. Actual Operating Data and Simulation Comparison
To verify the feasibility and engineering practicability of the quadratic programming algorithm, the comparison between actual operating data and simulation results is indispensable. Thanks to the time bounded weight matching algorithm, the final weight value is set as [2.841×10 3 , 2.7×10 -4 , 1], whose trip time is slightly modified from the actual operation time. The interval of operation comparison is set between 621.3 to 760.9 kilometers, as the longest start-stop cycle in actual operation, to make comparison in the same experimental environment. In this case, as described in Figure 8 the train can make a good trade-off between energy conservation, punctuality and smoothness. Since the data recorded by the train only has discrete and irregularly distributed speed information, the linear interpolation method is used in solving the additional information of train operation. The existence of errors is possible, especially in the control force, where sudden changes in force affect experimental results. The summary of QP simulation is shown in Table Ⅲ.   Figure 8, it can be found that the optimal trajectory calculated by the quadratic programming has better adaptability to the route. For example, in Scenario a, the starting phase: the actual operation is under braking condition at 1km, whose trip time is around 1563.3s. In contrast, the QP optimization curve completely adopts the traction condition with trip time of 1168.83s, which meets the requirements of starting phase of trains operation. QP optimization method is also able to select operation mode similar to the excellent manual operation, as illustrated in Scenario b: the algorithm has foreseen the descending ramp ahead, so it starts coasting before the ramp to ensure that the speed is maintained within the corresponding range and save energy as much as possible. Meanwhile, the actual average operating speed reached 83.98km/h, higher than that of QP simulation of 80.72km/h, signifying more energy consumption in the related interval. At 742 km, as presented in Scenario 3: the train is in downhill but about to face a very large uphill, so QP optimization curve is in coasting mode in advance. Simultaneously, it is ensured that the train does not exceed the speed limit. Manually driven train is relatively lack of foresight into the future route condition, which just maintains a specific speed, even allows longer time for braking mode, signifying extra energy consumption.
Overall, manual operation and QP algorithm can both complete transportation trip with guaranteed safety. But the former partially ignores the factor of energy, punctuality and smoothness, while the quadratic programming algorithm can make a good optimization of the trajectory speed nowadays.

C. Weight Selection Comparison
In order to further verify the effectiveness and superiority of the weight selection algorithm, we compare it with the linear approaching algorithm, which adjusts individual weight through the time difference between target and current weight as presented in Figure 10. Since the linear approaching method can only improve single weight value, we first set other substitution weight as the target ideal value for comparison and verification in the same conditions.  The results of Figure 10 show that the number of iterations calculated by the linear approach method is significantly larger. Each iteration will perform at least one QP solving process, which increases the calculation time. Although its final trip time of linear approaching may be more accurate through multiple iterations, the w1 of weight vector is provided in advance (for the sake of the same experimental environment). In practice, w1 needs to be adjusted with linear approaching, and the calculation process is far more complicated than the experiment.
In view of the fact that the time-priority weight selection algorithm driven by data, it can solve the optimal weight of the time bounded through one step by narrowing down the accessible weight area, and the final time errors are 1.64s, 1.23s and -0.43s as shown in Table Ⅳ. On the contrary, the linear approximation method is not efficient as the superior. For the convenience of display, the time deviation of more than 500 seconds is ignored in 140km and 70km intervals, deviation more than 150 seconds in 30km interval either.
From this comparison, it can be confirmed that time bounded weight selection algorithm has high computational efficiency, and it is able to select a multi-objective adaptation weight on the basis of guaranteeing punctual trip time.

VI. CONCLUSION
In this paper, the optimization of trajectory speed curve for heavy haul trains is studied. The multi-objective optimization model is presented in terms of consideration of energy consumption, punctuality and smooth changing of control force. The solving process of this problems is based on the QP algorithm. Due to the multi-objective nature of heavy haul train optimization, the distribution of various indicators with weight is presented and we proposed a time bounded weight selection algorithm of high computational efficiency, which is able to consider the overall multi-factors. Compared with the actual operation, it can be concluded that the QP method can be more adapted to real route condition and time bouned weight section method has the ability to calculate rapidly punctual adaptive weights.