Robust energy-efficient train speed profile optimization in a scenario-based position–time–speed network

Abstract

Train speed profile optimization is an efficient approach to reducing energy consumption in urban rail transit systems. Different from most existing studies that assume deterministic parameters as model inputs, this paper proposes a robust energy-efficient train speed profile optimization approach by considering the uncertainty of train modeling parameters. Specifically, we first construct a scenario-based position–time–speed (PTS) network by considering resistance parameters as discrete scenario-based random variables. Then, a percentile reliability model is proposed to generate a robust train speed profile, by which the scenario-based energy consumption is less than the model objective value at \alpha confidence level. To solve the model efficiently, we present several algorithms to eliminate the infeasible nodes and arcs in the PTS network and propose a model reformulation strategy to transform the original model into an equivalent linear programming model. Lastly, on the basis of our field test data collected in Beijing metro Yizhuang line, a series of experiments are conducted to verify the effectiveness of the model and analyze the influences of parameter uncertainties on the generated train speed profile.

Graphical abstract

Keywords

robust train speed profile / percentile reliability model / scenario-based position–time–speed network / mixed-integer programming

1 Introduction

Urban rail transit systems play an irreplaceable role in economic growth and social development due to their characteristics, such as high carrying capacity, fast speed, low carbon emissions, safety, and punctuality. Until January 2020, 43 cities in China have constructed urban rail transit systems. Shanghai owns 22 metro lines with a total length of 705 km, ranking first in the world. Beijing has 18 metro lines, and the total length is approximately 700 km, which is the second in the world. Generally, the metro transit systems in Shanghai and Beijing carry more than 12 million passengers each day, which greatly improve the mobility of urban citizens and alleviate the pressure on the urban road network.

Meanwhile, the energy consumption of urban rail transit systems rises sharply, which has placed tremendous pressures on the rail management governors. For example, the total energy consumption of Beijing metro is more than 20 million kWh per day, accounting for 1% to 2% of the power supply of the entire city. More than half of the energy consumption of Beijing metro is used by the traction operations to move the trains (Gao and Yang, 2019). Therefore, reducing this type of energy consumption in urban rail transit systems has become a key problem that operators are eager to solve.

In practice, the trains of urban rail transit systems are operated by human drivers or using an automatic train control (ATC) system. The ATC system can automatically track a pregiven train speed profile when running between adjacent stations. Theoretically, more than one train speed profile satisfies a given running time. These profiles correspond to different train traction energy consumptions (Howlett and Pudney, 1995). Thus, train speed profile optimization is regarded as an effective method to reduce the energy consumption of trains, and numerous studies have addressed this problem with different mathematical formulations in the past years (Ke et al., 2011; Wang and Goverde, 2019; Wang et al., 2020). A systematic review can be found in Scheepmaker et al. (2017) on energy-efficient train control methodologies.

Nevertheless, in accordance with our field tests in Beijing metro, the resistance parameters of a train operation model are usually affected by many external disturbances, such as the wear and tear of the traction/braking system, weather (e.g., rain, wind, and snow), and air resistances, leading to the deviation between theoretical (based on deterministic formulations) and practical models (Yin et al., 2020). As a result, most existing approaches cannot fully capture the practical requirement to optimize energy utilization. Although several recent studies have considered the uncertainty of the resistance parameters, the parameters are usually modeled as a specific random distribution within a certain boundary, such as Gaussian, uniform, or Pareto distribution (Li et al., 2013a; Bešinović et al., 2013). In this paper, we consider using scenario-based random variables to model the uncertainty of resistance parameters, given that we can easily construct stochastic scenarios with different historical data. In addition, as an effective representation to describe the spatiotemporal movement of commodities, space–time networks aim to integrate physical transportation networks with vehicle time-dependent movements/trajectories, which have been widely used in the studies of traffic problems, such as shortest path (Yang and Zhou, 2014; 2017), vehicle routing (Mahmoudi and Zhou, 2016), train scheduling (Yin et al., 2017b), and high-speed train timetabling (Zhou et al., 2017) problems. Combining the concept of space–time network with parameter uncertainties, we first propose a percentile reliability model (PRM) by constructing a scenario-based space–time network (here, we call it scenario-based position–time–speed (PTS) network) to formulate the train speed profile opti-mization problem with uncertainty resistance parameters.

1.1 Literature review

The energy-efficient train control problem has been a popular research topic in the past two decades, and many studies have proposed various methodologies for the optimization of train speed profiles (Howlett and Pudney, 1995; Howlett, 2000; Liu and Golovitcher, 2003; Wang et al., 2013; Wang and Goverde, 2019). The relevant studies can generally be classified into two categories: Analytical and numerical methods.

Analytical methods can obtain the optimal solution accurately. Nonetheless, the computation process is relatively complicated, and the model structure is usually required to satisfy some good properties. For example, on the basis of Pontryagin’s maximum principle, Howlett and Pudney (1995) proved that an optimal train speed profile exists only when its Hamiltonian function is maximized with a set of boundary conditions. Then, they proved that an optimal energy-efficient speed profile should consist of maximum acceleration, cruising, coasting, and maximum braking phases. Considering various speed limits and grade values, Liu and Golovitcher (2003) addressed the problem of energy-efficient train speed profile designing on a long-distance rail line, in which an analytical model was proposed to solve the switching time among different phases (e.g., accelerating, cruising, coasting, and braking). Su et al. (2013) proposed a bilevel model to optimize the overall energy consumption in consideration of the speed profiles of a single train in every section. In the upper level, they determined the optimal allocation of reserve time for each section. In the second level, with the given trip time for each section, energy-efficient speed profiles were obtained using an efficient algorithm based on Pontry-agin’s maximum principle. Albrecht et al. (2016a; 2016b) addressed the problem of finding an energy-efficient train speed profile on an undulating track with steep grades. The tractive and braking control forces were bounded by nonincreasing speed-dependent magnitude constraints, and the rate of energy caused by frictional resistance was given using a nonnegative strictly convex function with respect to train speed.

Compared with analytical methods, numerical methods, e.g., mathematical programming or metaheuristics, have relatively less requirements for the objective function of the model and can make a better trade-off between the optimization performance and computational efficiency (Yin et al., 2017a). For instance, Ko et al. (2004) proposed a nonlinear mathematical model for optimizing the energy consumption of trains under fixed running time, limited electric motive force, and electric brake force. Then, a dynamic programming (DP) algorithm was designed to solve the model. Considering energy consumption and riding comfort in the objective function, Wang et al. (2013) proposed two mathematical formulations. The first was a nonlinear programming formulation, which was solved using the pseudospectral method; and the second was a mixed-integer linear programming (MILP) formulation, which was solved using a MILP solver. Calderaro et al. (2014) proposed a two-step approach to optimize the train speed profile in consideration of the slopes and curves of the track. First, based on a DP algorithm, a set of semioptimal energy-efficient train speed profiles was generated. Then, a simulation tool was constructed to evaluate the speed profiles on a whole metro network. Owing to the complexity of solving mathematical formulations, existing literature has also proposed a series of metaheuristics to search a good-quality solution efficiently. For example, Kim and Chien (2011) developed a simulated annealing algorithm to search the optimal energy-efficient train speed profile in consideration of track alignment, speed limit, and schedule adherence. Lu et al. (2013) discussed and compared the performance of three solution methodologies, namely, ant colony optimization, a genetic algorithm (GA), and DP, for optimizing energy-efficient train speed profiles. The results demonstrated that DP could obtain the best-quality solution among the three algorithms but with the severest calculation burden. As heuristic algorithms, the ant colony algorithm performed better than GA, especially for cases with longer journey time.

To build accurate and adaptable models, some recent studies have begun to use historical data to study the energy-efficient train speed profile optimization problem. For example, de Martinis and Corman (2018) indicated the hidden potential in large sets of data for improving energy efficiency through novel, data-driven approaches, which demonstrated the necessity of data-driven methods in the study of the optimal energy-efficient train speed profile generating problem. Amrani et al. (2018) combined GA with a data-driven random forest model to optimize the train speed profile. The random forest model, constructed using a large amount of historical data, was used for estimating energy consumption as the fitness function of GA. Huang et al. (2019) proposed a support vector machine regression model and a random forest regression algorithm to identify the importance degree of train speed contributing to the total energy consumption at different positions.

Different from most existing studies on deterministic environments, several researchers have recently addressed the importance of incorporating uncertain factors into the optimization of energy-efficient train speed profiles. For instance, Li et al. (2013b) proposed an expected value model to solve the energy-efficient train speed problem based on the assumption that resistance coefficients are random variables, and the model was then solved using an iterative algorithm. Considering the uncertainty of manual driving parameters, i.e., driver response time and driving behaviors, Sicre et al. (2014) developed a fuzzy rule-based train speed profile optimization model and designed a customized GA. Fernández-Rodríguez et al. (2018) characterized the uncertainty in manual driving by introducing fuzzy numbers (driving behaviors and anticipation/delay time from crisp switching points). A multiobjective optimization, involving running time, energy consumption, and the risk of arrival delay, was designed. Wang et al. (2017) proposed a stochastic constrained shortest path approach to generate the optimal train speed profile, in which the energy consumption of the train on different arcs in the speed–distance network was regarded as a random variable. Given the large number of integer variables, a Lagrangian relaxation heuristic was designed to solve the model in an acceptable time. To improve the on-time arrival of trains, Wang et al. (2020) considered the uncertainty of train traction force and resistance in the trajectory optimization problem. They formulated the problem as a Markov decision process and solved it by combining off-line approximate DP and online search process.

1.2 Focus of this study

As mentioned in the literature review, numerous studies have been conducted on the energy-efficient train speed profile optimization problem, but only a few of them have considered uncertainties in mathematical formulation. All the proposed models are based on the use of an expected value function as an evaluation index. The limitation is that the expected value-based formulation only generates a single expected train profile and cannot fully meet the diversity requirements of rail decision makers. To this end, our study aims to formulate PRM to design a robust optimal energy-efficient train speed profile, in which the objective function is defined using \alpha-critical value of random total energy consumption. To overcome the computational complexity of solving real-world scale instances, a scenario-based PTS network is introduced with several node/arc reduction algorithms. To show the main differences between existing research and our work clearly, a detailed comparison with relevant publications is provided in Table 1.

Specifically, this paper aims to offer the following contributions to the study of the train speed profile optimization problem:

(1) One of the contributions of this paper is the construction of a scenario-based PTS network. The main difference from existing space–time networks (Zhou et al., 2017) is the consideration of scenario-based basic resistance coefficients. We extend space–time networks into a scenario-based PTS network, for the first time, by regarding the basic resistance parameters of the Davis equation as random variables in the network. To handle the uncertainty of train model parameters, which affects not only the energy consumption but also the feasibility of the network arc, a traditional space–time network must construct a number of networks with different parameters, which will lead to a dramatic increase in model solution time. Accordingly, we integrate the networks in different scenarios into a scenario-based PTS network by introducing a series of network node/arc reduction algorithms.

(2) For generating a robust train speed profile, this paper proposes a novel PRM with the constructed scenario-based PTS network. The total energy consumption of the speed profile is identified as a random variable due to the uncertainty of train model parameters. The objective function of the proposed model is formulated through defining \alpha-critical value of random total energy consumption to find a robust train speed profile with the scenario-based energy consumption that is less than the objective value at \alpha confidence level. In addition, for improving the solution efficiency, a linear equivalent formulation of PRM is deduced.

The remainder of this paper is organized as follows. In Section 2, we provide a detailed description of the problem about generating a robust optimal train speed profile with a scenario-based PTS network, including the introduction of the train speed profile optimization problem, the construction and reduction of the PTS network, and the process of building the scenario-based PTS network. In Section 3, PRM is formulated. Then, we present a method to obtain the linear equivalent of the proposed model. In Section 4, on the basis of our field test data in Beijing metro Yizhuang line, we perform a series of experiments to verify the effectiveness of the model and analyze the influence of relevant parameters on the generated optimal speed profile. Lastly, conclusions and further studies are presented in Section 5.

2 Problem description

For the convenience of description, Table 2 presents detailed definitions for all the relevant parameters and decision variables used in the following content of this paper.

Consider a train that runs from one station A to the next station B, in which the distance S and planned travel time T between these two stations are pregiven in accordance with the data from rail managers. The movement of a train (i.e., the train speed–position profile) is characterized by the acceleration/deceleration value (a) along the rail track, which is related to the force F provided by the train (traction force F>\mathrm{0}, braking force ) and the resistance force R of the train. The movement of the train is strictly subjected to a series of constraints, such as the train speed limit constraints, maximum traction/braking force constraints, and constraints related to riding comfort (avoidance of jerk). The dynamic train control model is also affected by the rail track, e.g., the gradient of the track (related to gradient resistance RG).

Theoretically, with the given input data, involving the planned timetable (i.e., the planned running time between the two stations), railway conditions (track slope and speed limitation), vehicle performance (maximum traction force and braking force), and resistances, more than one speed profile satisfies these operation constraints but with different traction energy consumption. An example is illustrated in Fig. 1, in which we plot three profiles with the same travel time yet different energy consumption between two stations with track gradient and speed limits. In particular, the red profile (profile 1) consumes the least energy consumption and involves four processes: Process s_{\mathrm{0}}\to s_{\mathrm{1}} represents the traction operation, in which we have F>\mathrm{0} and a>\mathrm{0}. The cruising operation process is s_{\mathrm{1}}\to s_{\mathrm{2}} with a=\mathrm{0}. s_{\mathrm{2}}\to s_{\mathrm{3}} refers to the coasting process, where F=\mathrm{0} and . Lastly, s_{\mathrm{3}}\to s_{\mathrm{4}} indicates the braking operation, and we have and . Therefore, the objective of the train speed profile optimization problem is to find one unique speed profile that minimizes the train energy consumption while satisfying the operational constraints, given deterministic input data.

However, applying deterministic formulations (Su et al., 2013; Wang et al., 2013; Calderaro et al., 2014) in practical implementations faces a critical issue: Some of the parameters (e.g., resistances and traction/braking efforts) in the train control model cannot be easily determined, and the values of various parameters may vary due to different weather, rail line, and locomotive conditions in accordance with our field tests in Beijing metro. The optimized train speed profile under deterministic formulations may become inefficient or even infeasible due to operational constraints. In this sense, considering the uncertainty of these parameters is essential to generate a robust energy-efficient train speed profile, which can be highly practical in real-world operations. Here, we consider the uncertainty of resistance parameters by introducing the scenario-based PTS network. Next, the detailed construction process of the network is provided.

3 Scenario-based PTS network

3.1 Construction of a PTS network

A PTS network is constructed on the basis of a 3D space with discrete position, time, and speed states (Fig. 2). We first divide the position interval [\mathrm{0},\mathrm{}S] along the railway into H_{\mathrm{1}} discrete segments with the same length \mathrm{\Delta }s, and H_{\mathrm{1}}=\lfloor S/\mathrm{\Delta }s\rfloor. As a result, we can use \{\mathrm{0},\mathrm{}\mathrm{\Delta }s,\mathrm{ }\mathrm{}\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{...},H_{\mathrm{1}}\mathrm{\Delta }s\} to represent the position dimension [\mathrm{0},\mathrm{}S] in the PTS network. Then, the time period [\mathrm{0},\mathrm{}T] of the train running between two adjacent stations is discretized into H_{\mathrm{2}} segments with the same length \mathrm{\Delta }t, where T is the planned running time given by the railway schedule. Set \{\mathrm{0},\mathrm{ }\mathrm{}\mathrm{\Delta }t,\mathrm{ }\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{...},\mathrm{ }\mathrm{}{H}_{\mathrm{2}}\mathrm{\Delta }t\} can be regarded as an approximation of [\mathrm{0},\mathrm{}T], where H_{\mathrm{2}}=\lfloor T/\mathrm{\Delta }t\rfloor. Similarly, candidate speed interval [\mathrm{0},\mathrm{}V] is discretized into \{\mathrm{0},\mathrm{}\mathrm{\Delta }v,\mathrm{2}\mathrm{\Delta }v,\mathrm{}\mathrm{...},\mathrm{ }\mathrm{}{H}_{\mathrm{3}}\mathrm{\Delta }v\}, where H_{\mathrm{3}}=\lfloor V/\mathrm{\Delta }v\rfloor and V is larger than the maximal allowed speed (speed limit) in all sections. As shown in Fig. 2, after the discretization of the three coordinate axes (position, time, and speed), the 3D space is approximately represented by many countable discrete nodes, which are termed set \mathcal{N} in the PTS network. A node (s,\mathrm{}t,\mathrm{}v) in the PTS network represents a kind of state of the train, including position (s), time (t), and speed (v). In Fig. 2, for example, N_{\mathrm{0}}, N_{\mathrm{1}}, N_{\mathrm{2}}, N_{\mathrm{3}}, and N_{\mathrm{4}} are nodes in \mathcal{N}. The connection of two arbitrary nodes (s,\mathrm{}t,\mathrm{}v)\in \mathcal{N} and (s^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{N} (t^{\prime }=t+\mathrm{1}) is regarded as an arc that represents the state transition of a train over time period \mathrm{\Delta }t. All the connections similar to this form the set of arcs \mathcal{A}, i.e.,

\mathcal{A}=\left\{\right(s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime }\left)\right|(s,\mathrm{}t,\mathrm{}v),\mathrm{}(s^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{N},\mathrm{}

(1)t^{\prime }=t+\mathrm{1}\}.

For example, in Fig. 2, the connection between nodes N_{\mathrm{0}} and N_{\mathrm{1}} is an arc (termed arc (\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v)) in the network, which indicates the state transition of the train from position 0 at time 0 with speed 0 to position \mathrm{\Delta }s at time \mathrm{\Delta }t with speed \mathrm{2}\mathrm{\Delta }v. After defining these nodes and arcs, we then use \mathcal{G}(\mathcal{N},\mathrm{}\mathcal{A}) to denote the constructed PTS network.

Remark 3.1: Theoretically, the discretization of position, time, and speed should be carefully designed to model the train movement process properly, while the network size could be relatively large. To balance accuracy and computational intensity to reduce the error caused by the construction of the network, here, we set \mathrm{\Delta }v=\mathrm{2}\mathrm{\Delta }s/\mathrm{\Delta }t on the assumption that each time unit is small enough so that the acceleration/deceleration of the train does not change in time period (t,\mathrm{}t+\mathrm{\Delta }t). Accordingly, if the train speed changes from v_{\mathrm{1}}=h_{\mathrm{1}}\mathrm{\Delta }v=\mathrm{2}{h}_{\mathrm{1}}\mathrm{\Delta }s/\mathrm{\Delta }t to v_{\mathrm{2}}=h_{\mathrm{2}}\mathrm{\Delta }v=\mathrm{2}{h}_{\mathrm{2}}\mathrm{\Delta }s/\mathrm{\Delta }t during time \mathrm{\Delta }t, the running distance of the train is (h_{\mathrm{1}}+h_{\mathrm{2}})\mathrm{\Delta }s, which is a feasible value in the PTS network. A detailed analysis with respect to the discretization of position, time, and speed will be presented in Section 5.

The introduction of the PTS network allows explicitly modeling the train speed profile as a direct route from the origin node to the destination node. As shown in Fig. 2, route N_{\mathrm{0}}\to N_{\mathrm{1}}\to N_{\mathrm{2}}\to N_{\mathrm{3}}\to N_{\mathrm{4}} (the black curve, where N_{\mathrm{0}} and N_{\mathrm{4}} are the origin and destination nodes, respectively) is a 3D speed profile (PTS profile) of the train between stations A and B. Moreover, the PTS network is useful to depict other information of the profile. For example, the projection of the 3D profile in Speed–Time plane (i.e., the blue curve) represents the speed–time profile of the train. The projection in Speed–Position plane (i.e., the red curve) of the 3D speed profile is the speed–position profile. Similarly, the time–position profile of the train is expressed by the projection in Time–Position plane (i.e., the yellow curve).

Given that the goal of this study is to generate an optimal train speed profile with the lowest traction energy consumption, the corresponding traction energy consumption between two nodes in the network can be regarded as the weight of the arc. Next, we will introduce the calculation of arc weight (i.e., the traction energy consumption) between two nodes (s,\mathrm{}t,\mathrm{}v)\in \mathcal{N} and (s^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{N}, where t^{\prime }=t+\mathrm{1}.

First, we consider the traction energy consumption e_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} of the train on arc (s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A}, which can be calculated by multiplying the traction force (when F_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}>\mathrm{0}, F_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is the force provided by the train moving on arc (s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })) and the running distance (s^{\prime }-s), i.e.,

(2)e_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}=\max \left(F_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}\right(s^{\prime }-s),\mathrm{}\mathrm{0}).

In accordance with the force analysis of the train and Newton’s second law, we have

(3)F_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}-R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}-R{G}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}=m{a}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}.

Therefore,

(4)F_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}=m{a}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}+R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}+R{G}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }},

where a_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is the acceleration/deceleration of the train on arc (s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A}, and m is the weight of the train, including onboard passengers. In particular, the value of a_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is denoted by a_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}=(v^{\prime }-v)/\mathrm{\Delta }t assum-ing that the acceleration/deceleration of the train is constant on each arc in this paper. In addition, the basic resistance R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is denoted by mg({\beta }_{\mathrm{1}}{({\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }})}^{\mathrm{2}}+{\beta }_{\mathrm{2}}{\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}+{\beta }_{\mathrm{3}}), where g is the acceleration parameter of gravity, {\beta }_{\mathrm{1}}, {\beta }_{\mathrm{2}}, and {\beta }_{\mathrm{3}} are the parameters of Davis’ formula, and {\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is the average speed of the train on arc (s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A}. For simplicity, {\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is approximated as \overline{v}=(v+v^{\prime })/\mathrm{2}. Hence, the basic resistance R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is reformulated as

(5)R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}=mg({\beta }_{\mathrm{1}}{\left({\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}\right)}^{\mathrm{2}}+{\beta }_{\mathrm{2}}{\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}+{\beta }_{\mathrm{3}})=mg({\beta }_{\mathrm{1}}{(v+v^{\prime })}^{\mathrm{2}}/\mathrm{4}+{\beta }_{\mathrm{2}}(v+v^{\prime })/\mathrm{2}+{\beta }_{\mathrm{3}}).

As mentioned above, the resistance caused by the gradient of the track R{G}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }} is involved, and it is calculated as

(6)R{G}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}=mg\sin {\theta }_{\overline{s}},

where {\theta }_{\overline{s}} is the gradient at position \overline{s} (\overline{s}=\lfloor (s+s^{\prime })/\mathrm{2}\rfloor ).

The basic resistance R{B}_{s,t,v} and the resistance caused by the gradient of the track R{G}_{s,t,v} at node (s,\mathrm{}t,\mathrm{}v) can be calculated, respectively, as

(7)R{B}_{s,t,v}=mg({\beta }_{\mathrm{1}}{\left(v_{s,t,v}\right)}^{\mathrm{2}}+{\beta }_{\mathrm{2}}{v}_{s,t,v}+{\beta }_{\mathrm{3}}),

(8)R{G}_{s,t,v}=mg\sin {\theta }_s,

which will be used in the following subsections.

3.2 Reduction of the PTS network

The train speed profile problem can be satisfactorily transformed into a resource-constraint shortest path problem through discretization of position, time, and speed. However, a critical issue arising from network discretization is the large number of discrete variables. For example, consider a rail track section with a length of 1000 m between two adjacent stations, in which the maximum speed that the train can reach is 20 m/s and the running time of the train is 80 s. If the discrete length of position (\mathrm{\Delta }s), time (\mathrm{\Delta }t), and speed (\mathrm{\Delta }v) are 5, 1, and 1, respectively, the total number of nodes is \mathrm{200}\times \mathrm{80}\times \mathrm{20}= 320000, and the total number of arcs is at most 125 million in the PTS network. Meanwhile, we observe that only a small part of nodes and arcs are feasible in the PTS network due to the unique properties of train operational constraints. For example, the origin and destination nodes of the train in actual operation must be (\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{0}) and (H_{\mathrm{1}},\mathrm{}{H}_{\mathrm{2}},\mathrm{}\mathrm{0}), respectively. The speed of the train at the intermediate nodes between the origin and destination nodes must be greater than 0, given that no parking is allowed as the train runs between two stations. In addition, the speed of all the nodes in the network cannot exceed the speed limit. Therefore, we denote the following two steps of the PTS network reduction strategies: A node-reduction strategy and an arc-reduction strategy.

3.2.1 Node reduction of the PTS network

In our study, the reduction of nodes for the PTS network is conducted via the following two steps: 1) projecting the 3D profile into different coordinate planes and 2) imposing node-cut strategies with practical experiences. We consider the projection region (or feasible region) on Speed–Time, Speed–Position, and Time–Position planes (Fig. 3) and then propose a series of constraints to eliminate the infeasible nodes in the PTS network.

We begin with the feasible region in Speed–Time plane, as shown in Fig. 3(a). Here, V^{\max } is the maximum speed the train can reach, a^{\max } and d^{\max } are the maximum acceleration and deceleration rates of the train, respectively, and T is the total running time of the train. Then, we can obtain the boundaries of the feasible region in Speed–Time plane, given as P_{\mathrm{1}}{P}_{\mathrm{2}}, P_{\mathrm{2}}{P}_{\mathrm{3}}, P_{\mathrm{3}}{P}_{\mathrm{4}}, and P_{\mathrm{4}}{P}_{\mathrm{1}} (P_{\mathrm{1}} to P_{\mathrm{4}} are the vertices of the feasible region). In particular, the slope of curve P_{\mathrm{1}}{P}_{\mathrm{2}} is a^{\max }; thus, curve P_{\mathrm{1}}{P}_{\mathrm{2}} can be presented as v\left(t\right)=a^{\max }t, t\in [\mathrm{0},\mathrm{}{t}_{\mathrm{1}}], t_{\mathrm{1}}=\lfloor \left(V^{\max }/a^{\max }\right)/\mathrm{\Delta }t\rfloor. Curve P_{\mathrm{2}}{P}_{\mathrm{3}} is v\left(t\right)=V^{\max }, t\in [t_{\mathrm{1}},\mathrm{}{t}_{\mathrm{2}}], t_{\mathrm{2}}=H_{\mathrm{2}}-\lfloor \left(V^{\max }/d^{\max }\right)/\mathrm{\Delta }t\rfloor. The slope of curve P_{\mathrm{3}}{P}_{\mathrm{4}} is –d^{\max }, such that the formula of P_{\mathrm{3}}{P}_{\mathrm{4}} is v\left(t\right)=d^{\max }(T-t), t\in [t_{\mathrm{2}},\mathrm{}T], and P_{\mathrm{4}}{P}_{\mathrm{1}} is v\left(t\right)=\mathrm{0}, t\in [\mathrm{0},\mathrm{}T]. As a result, the feasible region in Speed–Time plane is

(9)\{\begin{array}{l}v\left(t\right)\le a^{\max }t\hfill \\ v\left(t\right)\le V^{\max }\hfill \\ v\left(t\right)\le d^{\max }(T-t)\hfill \\ v\left(t\right)\ge \mathrm{0}\hfill \end{array}.

As shown in Fig. 3(b), curves P_{\mathrm{1}}{P}_{\mathrm{2}}, P_{\mathrm{2}}{P}_{\mathrm{3}}, P_{\mathrm{3}}{P}_{\mathrm{4}}, and P_{\mathrm{4}}{P}_{\mathrm{1}} are the boundaries of the feasible region in Speed–Position plane, and S is the length between two stations. The vertices of the region are P_{\mathrm{1}}(\mathrm{0},\mathrm{}\mathrm{0}), P_{\mathrm{2}}(s_{\mathrm{1}},\mathrm{}{V}^{\max }), P_{\mathrm{3}}(s_{\mathrm{2}},\mathrm{}{V}^{\max }), and P_{\mathrm{4}}(S,\mathrm{}\mathrm{0}), where s_{\mathrm{1}}={\left(V^{\max }\right)}^{\mathrm{2}}/\left(\mathrm{2}{a}^{\max }\right) and s_{\mathrm{2}}=S-{\left(V^{\max }\right)}^{\mathrm{2}}/\left(\mathrm{2}{a}^{\max }\right). Curves P_{\mathrm{1}}{P}_{\mathrm{2}}, P_{\mathrm{2}}{P}_{\mathrm{3}}, and P_{\mathrm{3}}{P}_{\mathrm{4}} are respectively given as v\left(s\right)=\sqrt{\mathrm{2}{a}^{\max }s}, s\in[\mathrm{0},\mathrm{}{s}_{\mathrm{1}}], v\left(s\right)=V^{\max }, s\in [s_{\mathrm{1}},\mathrm{}{s}_{\mathrm{2}}], and v\left(s\right)=\sqrt{\mathrm{2}{d}^{\max }(S-s)}, s\in [s_{\mathrm{2}},\mathrm{}S]. Curve P_{\mathrm{4}}{P}_{\mathrm{1}} can be represented by v\left(s\right)=\mathrm{0}, s\in [\mathrm{0},\mathrm{}S]. Therefore, the feasible region in Speed–Position plane can be represented by

(10)\{\begin{array}{l}v\left(s\right)\le \sqrt{\mathrm{2}{a}^{\max }s}\hfill \\ v\left(s\right)\le V^{\max }\hfill \\ v\left(s\right)\le \sqrt{\mathrm{2}{d}^{\max }(S-s)}\hfill \\ v\left(s\right)\ge \mathrm{0}\hfill \end{array}.

Similarly, the vertices of the region in Fig. 3(c) are P_{\mathrm{1}}(\mathrm{0},\mathrm{}\mathrm{0}), P_{\mathrm{2}}(s_{\mathrm{1}},\mathrm{}{t}_{\mathrm{1}}), P_{\mathrm{3}}(S,\mathrm{}{t}_{\mathrm{3}}), P_{\mathrm{4}}(S,\mathrm{}T), and P_{\mathrm{5}}(\mathrm{0},\mathrm{}T), where t_{\mathrm{3}}=(S-s_{\mathrm{1}})/V^{\max }. The boundaries of the feasible region in Time–Position plane consist of curves P_{\mathrm{1}}{P}_{\mathrm{2}}, P_{\mathrm{2}}{P}_{\mathrm{3}}, P_{\mathrm{3}}{P}_{\mathrm{4}}, P_{\mathrm{4}}{P}_{\mathrm{5}}, and P_{\mathrm{5}}{P}_{\mathrm{1}}, which are denoted, respec-tively, by t\left(s\right)=\sqrt{\mathrm{2}s/a^{\max }}, s\in [\mathrm{0},\mathrm{}{s}_{\mathrm{1}}], t(s)=\mathrm{1}/\left(V^{\max }s\right)+V^{\max }/\left(\mathrm{2}{a}^{\max }\right), s\in [s_{\mathrm{1}},\mathrm{}S] (the slope is \mathrm{1}/V^{\max }), s\left(t\right)=S, t\in [\mathrm{0},\mathrm{}T], t\left(s\right)=T, s\in [\mathrm{0},\mathrm{}S], and s\left(t\right)=0, t\in [\mathrm{0},\mathrm{}T], where s_{\mathrm{1}}={\left(V^{\max }\right)}^{\mathrm{2}}/\left(\mathrm{2}{a}^{\max }\right). Then, the feasible region in Time–Position plane is represented by

(11)\{\begin{array}{l}t\left(s\right)\ge \sqrt{\mathrm{2}s/a^{\max }}\hfill \\ t(s)\ge \mathrm{1}/\left(V^{\max }s\right)+V^{\max }/\left(\mathrm{2}{a}^{\max }\right)\hfill \\ s\left(t\right)\le S\hfill \\ t\left(s\right)\le T\hfill \\ s\left(t\right)\ge \mathrm{0}\hfill \end{array}.

In accordance with the formulated constraints, a node-reduction algorithm is designed to identify feasible nodes in the PTS network. The details about the algorithm are presented in Algorithm 1 (Table 3).

3.2.2 Arc reduction of the PTS network

After the reduction of nodes in the PTS network, we consider reducing the size of arcs in the network to eliminate infeasible state transitions. Here, an algorithm (Algorithm 2, Table 4) based on DP is designed to generate the feasible arc set {\mathcal{A}}^{*} in the PTS network. The idea is that we first divide the nodes in the PTS network into three categories: The origin node where only acceleration is allowed, the destination node where only deceleration is allowed, and other nodes with no limits. In addition, considering that the transition of train state is strictly subject to the maximum train acceleration a^{\max } and train deceleration d^{\max }, we go through the entire time horizon from t=\mathrm{0} to t=T with the time step length of \mathrm{\Delta }t for traversal searching of the follower nodes (s,\mathrm{}t,\mathrm{}v) that the current node (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}}) can reach. If it is reachable between nodes (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}}) and (s,\mathrm{}t,\mathrm{}v), arc (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}},\mathrm{}s,t,\mathrm{}v) is an feasible arc in the PTS network, which is then added into feasible arc set {\mathcal{A}}^{*}.

To illustrate the process of DP clearly, we consider Fig. 4 as an example.

Example 3.1: In this example, our algorithm starts with the origin node N_{\mathrm{0}}=(\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{0}). When regarding origin node N_{\mathrm{0}} as the current node, only the acceleration process is involved because the train must be accelerated from the station. Then, the minimum speed {\stackrel{⌣}{v}}_{\mathrm{0},\mathrm{0},\mathrm{0}} is \mathrm{\Delta }v (i.e., the minimum speed unit). The maximum speed {\widehat{v}}_{\mathrm{0},\mathrm{0},\mathrm{0}} of the follower nodes is calculated as

(12){\widehat{v}}_{\mathrm{0},\mathrm{0},\mathrm{0}}=\min ({\widehat{a}}_{\mathrm{0},\mathrm{0},\mathrm{0}}\mathrm{\Delta }t,\mathrm{}{\widehat{V}}_s),

where {\widehat{V}}_s is the speed limit at position s, and

(13){\widehat{a}}_{\mathrm{0},\mathrm{0},\mathrm{0}}=(m{a}^{\max }-R{B}_{\mathrm{0},\mathrm{0},\mathrm{0}}-R{G}_{\mathrm{0},\mathrm{0},\mathrm{0}})/m.

The yellow nodes in the figure are the ones with the speed value in [{\stackrel{⌣}{v}}_{\mathrm{0},\mathrm{0},\mathrm{0}},\mathrm{}{\widehat{v}}_{\mathrm{0},\mathrm{0},\mathrm{0}}]=[\mathrm{\Delta }v,\mathrm{}\mathrm{2}\mathrm{\Delta }v]. Then, the reachable node (s,\mathrm{}t,\mathrm{}v) from the current node (s_{\mathrm{0}},\mathrm{\hspace{0.17em}}{t}_{\mathrm{0}},\mathrm{\hspace{0.17em}}{v}_{\mathrm{0}}) is constrained by the possible travel distance of the train from node (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}}) to node (s,\mathrm{}t,\mathrm{}v), i.e.,

(14)s=\lfloor \left(v\mathrm{\Delta }v/\mathrm{2}\right)\mathrm{\Delta }t/\mathrm{\Delta }s\rfloor .

In Fig. 4, N_{\mathrm{1}}=(\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v) is the reachable node of the current node N_{\mathrm{0}}, as a result of which, (\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{2}\mathrm{\Delta }s,\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v) is the feasible arc of the PTS network (i.e., (\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v)\in {\mathcal{A}}^{*}).

Next, regarding the destination nodes of the feasible arcs found in the previous step as the current nodes (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}}), we can calculate the maximum train acceleration {\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}} (Eq. (13)) and deceleration {\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}} at node (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}})

(15){\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}=(m{d}^{\max }+R{B}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}+R{G}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}})/m.

Then, the possible speed of the follower nodes is in [{\stackrel{⌣}{v}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}},\mathrm{}{\widehat{v}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}], where

(16){\widehat{v}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}=\min (v_{\mathrm{0}}+{\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}\mathrm{\Delta }t,\mathrm{}{\widehat{V}}_s),

(17){\stackrel{⌣}{v}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}=\max (v_{\mathrm{0}}-{\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}\mathrm{\Delta }t,\mathrm{}\mathrm{\Delta }v).

In addition, the constraint of the reachable node (s,\mathrm{}t,\mathrm{}v) from the current node (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}}) based on the possible travel distance over \mathrm{\Delta }t is

(18)s=s_{\mathrm{0}}+\lfloor ((v_{\mathrm{0}}+v)\mathrm{\Delta }v/\mathrm{2})\mathrm{\Delta }t/\mathrm{\Delta }s\rfloor .

For example, let N_{\mathrm{1}} be the current node, and the possible speed of the follower nodes is in [\mathrm{\Delta }v,\mathrm{}\mathrm{4}\mathrm{\Delta }v]. Nodes N_{\mathrm{2}}=(\mathrm{7}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{3}\mathrm{\Delta }v), N_{\mathrm{3}}=(\mathrm{6}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v), and N_{\mathrm{4}}=(\mathrm{5}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{\Delta }v) are the reachable nodes. Moreover,

(\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v,\mathrm{}\mathrm{7}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{3}\mathrm{\Delta }v),\mathrm{}(\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v,\mathrm{}\mathrm{6}\mathrm{\Delta }s,\mathrm{}

\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v),\mathrm{}(\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v,\mathrm{}\mathrm{5}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{\Delta }v)\in {\mathcal{A}}^{*}.

The calculation process is terminated until the destination node of the train is reached. In Fig. 4, the destination node is N_{\mathrm{5}}=(\mathrm{8}\mathrm{\Delta }s,\mathrm{}\mathrm{3}\mathrm{\Delta }t,\mathrm{}\mathrm{0}), and arc (\mathrm{6}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v,\mathrm{}\mathrm{8}\mathrm{\Delta }s,\mathrm{3}\mathrm{\Delta }t,\mathrm{}\mathrm{0}) is the feasible arc, given that node (\mathrm{6}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{2}\mathrm{\Delta }v) satisfies constraints (19) and (20):

(19)s_{\mathrm{0}}+\lfloor \left(v_{\mathrm{0}}\mathrm{\Delta }v/\mathrm{2}\right)\mathrm{\Delta }t/\mathrm{\Delta }s\rfloor =H_{\mathrm{1}},

(20)a_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}},H_{\mathrm{1}},H_{\mathrm{2}},\mathrm{0}}\ge –{\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}.

Lastly, we remove the disconnected arcs (the destination node of this type of arc is not the origin node of any feasible arc in {\mathcal{A}}^{*}) from the feasible arc {\mathcal{A}}^{*}, such as the red dotted arrows in the figure. Consequently, in this example, we can obtain a subset of arcs in \mathcal{A}, given by

{\mathcal{A}}^{*}=\left\{\right(\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v),\mathrm{}(\mathrm{2}\mathrm{\Delta }s,\mathrm{}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v,\mathrm{}\mathrm{6}\mathrm{\Delta }s,\mathrm{}

(21)\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v),\mathrm{}(\mathrm{6}\mathrm{\Delta }s,\mathrm{}\mathrm{2}\mathrm{\Delta }t,\mathrm{}\mathrm{2}\mathrm{\Delta }v,\mathrm{}\mathrm{8}\mathrm{\Delta }s,\mathrm{}\mathrm{3}\mathrm{\Delta }t,\mathrm{}\mathrm{0}\left)\right\}.

From Example 3.1, the number of arcs can be reduced markedly by eliminating infeasible arcs in accordance with our DP algorithm. For simplicity, we next use \mathcal{A} and \mathcal{N} to represent the arcs and nodes after reduction in the following content of this paper.

3.3 Construction of a scenario-based PTS network

Resistance parameters are practically significant to model the movement of trains, but the values of these parameters are usually affected by the external and internal factors of the train. For example, the aerodynamic resistance parameter ({\beta }_{\mathrm{3}} in this paper) of a new British diesel multiple-unit train is 0.026, whereas it is 0.0574 for an old one (Li et al., 2013a). Therefore, several studies have considered the uncertainty of resistance parameters. In this paper, we discretize uncertain resistance parameters into random variables on the basis of scenarios. Here, we select the values of uncertain resistance parameters in several key scenarios as the realization values of these random variables, i.e., {\beta }_{\mathrm{1}}^{\omega }, {\beta }_{\mathrm{2}}^{\omega }, and {\beta }_{\mathrm{3}}^{\omega }, where \omega is the index of scenario (\omega \in \Omega). The probability of each scenario \omega is P_{\omega }.

We consider the following two aspects of modifications caused by parameter uncertainty: 1) the traction energy consumption on each arc (s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A} and 2) the generation of feasible arcs in the PTS network (including the generation of feasible nodes). Thus, we propose the following strategies to construct a scenario-based PTS to handle parameter uncertainties.

The energy consumption on arc (s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime }) under scenario \omega can be derived from Eqs. (2)–(6), i.e.,

(22)e_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}^{\omega }=\max \left(\right(m{a}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}+R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}^{\omega }+R{G}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}\left)\right(s^{\prime }-s)\mathrm{\Delta }s,\mathrm{}\mathrm{0}),

where R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}^{\omega } is the basic resistance received by the train under scenario \omega and can be calculated as

(23)R{B}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}^{\omega }=mg({\beta }_{\mathrm{1}}^{\omega }{\left({\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}\right)}^{\mathrm{2}}+{\beta }_{\mathrm{2}}^{\omega }{\overline{v}}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}+{\beta }_{\mathrm{3}}^{\omega }).

Next, we consider the influence of parameter uncertainties on the generation of feasible arcs. Given that the resistance parameters {\beta }_{\mathrm{1}}^{\omega }, {\beta }_{\mathrm{2}}^{\omega }, and {\beta }_{\mathrm{3}}^{\omega } directly affect the calculation of maximum train acceleration {\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}} and deceleration {\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}} at the current node (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}}) under each scenario, constructing a total of \left|\Omega \right| PTS networks by using {\beta }_{\mathrm{1}}^{\omega }, {\beta }_{\mathrm{2}}^{\omega }, and {\beta }_{\mathrm{3}}^{\omega } as input to Algorithm 2, i.e., one scenario corresponds to one PTS network, is theoretically required. However, this requirement will dramatically increase the space and time cost for solving the problem. Consequently, to improve the solution efficiency, we pro-pose to aggregate these uncertain scenarios and construct only one scenario-based PTS network by revising Steps 5, 10, and 16 in Algorithm 2. Specifically, the calculation of {\widehat{a}}_{s,t,v}, {\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}, and e_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}},s^{\prime },t^{\prime },v^{\prime }}^{\omega }, (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}},\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A}, in Algorithm 2 is updated by considering the uncertain values of {\beta }_{\mathrm{1}}^{\omega }, {\beta }_{\mathrm{2}}^{\omega }, and {\beta }_{\mathrm{3}}^{\omega }, which are given by

(24){\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}={\max }_{\omega \in \Omega }\left({\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }\right),

(25){\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}={\max }_{\omega \in \Omega }({\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }),

(26)e_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}},s^{\prime },t^{\prime },v^{\prime }}^{\omega }=\{\begin{array}{ll}{e}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}},s^{\prime },t^{\prime },v^{\prime }}^{\omega },\hfill & –{\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }\le a_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}},s^{\prime },t^{\prime },v^{\prime }}\le {\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }\hfill \\ \\ M,\hfill & \mathrm{otherwise}\hfill \end{array},

where {\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega } and {\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega } represent the maximum acceler-ation and deceleration at node (s_{\mathrm{0}},\mathrm{}{t}_{\mathrm{0}},\mathrm{}{v}_{\mathrm{0}})under scenario \omega, respectively, and they are calculated as

(27){\widehat{a}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }=(m{a}^{\max }-R{B}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }-R{G}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}})/m,

(28){\widehat{d}}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }=(m{d}^{\max }+R{B}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}}^{\omega }+R{G}_{s_{\mathrm{0}},t_{\mathrm{0}},v_{\mathrm{0}}})/m.

4 Mathematical model

In this paper, we translate the energy-efficient train speed profile optimization problem into a constrained shortest path problem by introducing a scenario-based PTS network. The path is expressed by a series of binary decision variables x_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}, (s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A},

(29)x_{s\mathrm{,}t\mathrm{,}v\mathrm{,}s\text{'}\mathrm{,}t\text{'}\mathrm{,}v\text{'}}\mathrm{=}\{\begin{array}{ll}\mathrm{1},\hfill & \mathrm{arc}\mathrm{}(s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\mathrm{is}\mathrm{\hspace{0.17em}}\mathrm{selected}\mathrm{\hspace{0.17em}}\hfill \\ \hfill & \mathrm{in}\mathrm{\hspace{0.17em}}\mathrm{the}\mathrm{\hspace{0.17em}}\mathrm{shortest}\mathrm{\hspace{0.17em}}\mathrm{path}\hfill \\ \mathrm{0},\hfill & \mathrm{otherwise}\hfill \end{array}.

After the definition of the decision variables, the corresponding energy consumption Z\left(X\right) is introduced as the evaluation index to evaluate the generated path X (X=\left\{x_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}\right|(s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A}\}). As mentioned before, the energy consumption of an arc becomes a scenario-based random variable because the basic resistance parameters are regarded as discrete random variables. Therefore, the scenario-based total energy consumption of a path X over scenario \omega is defined as Z(X,\mathrm{}\omega ), \omega \in \Omega.

4.1 Mixed-integer programming formulation

Considering the randomness of Z(X,\mathrm{}\omega ), \omega \in \Omega, we adopt a stochastic programming approach to formulate the energy-efficient train speed profile problem with the aid of the scenario-based PTS network. Instead of using a traditional expected value of energy consumption as the objective function (Li et al., 2013b; Wang et al., 2017), we adopt a more practical approach that aims to find a robust energy-efficient train speed profile. In particular, the robustness is captured using PRM with confidence level \alpha. That is, our objective is to find a train speed profile in the scenario-based PTS network such that its percentile reliability is not lower than \alpha. Next, a detailed formulation process of the model is provided.

First, the energy consumption of path X over scenario \omega \in \Omega, i.e., Z(X,\mathrm{}\omega ), is given as follows:

Z(X,\mathrm{}\omega )={\displaystyle \sum _{(s,t,v,s^{\prime },t^{\prime },v^{\prime })\in \mathcal{A}}{e}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}^{\omega }}{x}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }},\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}

(30)\omega \in \mathrm{\Omega }.

Then, to deal with the randomness of Z(X,\mathrm{}\omega ), \alpha-critical value Z{\left(X\right)}_{\inf }\left(\alpha \right) of Z\left(X\right) is introduced on the basis of the critical value definition given in Liu (2002)

(31)Z{\left(X\right)}_{\inf }\left(\alpha \right)=\inf \left\{\phi \right|\mathrm{Pr}\left\{Z\right(X,\mathrm{}\omega )\le \phi \}\ge \alpha \}.

Using \alpha-critical value Z{\left(X\right)}_{\inf }\left(\alpha \right) of Z\left(X\right), we formulate PRM for the energy-efficient train speed profile optimization problem, which is given as follows:

(32)\min \phi

(33)s.t.\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{Pr}\left\{Z\right(X,\mathrm{}\omega )\le \phi \}\ge \alpha ,

(34){\displaystyle \sum _{(s^{\prime },t^{\prime },v^{\prime })\in \mathcal{N}}{x}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}}-{\displaystyle \sum _{(s^{\prime },t^{\prime },v^{\prime })\in \mathcal{N}}{x}_{s^{\prime },t^{\prime },v^{\prime },s,t,v}}=\{\begin{array}{l}\begin{array}{ll}\mathrm{1},\hfill & \mathrm{}(s,\mathrm{}t,\mathrm{}v)=(\mathrm{0},\mathrm{}\mathrm{0},\mathrm{}\mathrm{0})\hfill \end{array}\hfill \\ \begin{array}{ll}-\mathrm{1},\hfill & (s,\mathrm{}t,\mathrm{}v)=(H_{\mathrm{1}},\mathrm{}{H}_{\mathrm{2}},\mathrm{}\mathrm{0})\hfill \end{array}\hfill \\ \begin{array}{ll}\mathrm{0},\hfill & \mathrm{otherwise}\hfill \end{array}\hfill \end{array},

{\displaystyle \sum _{(s^{\prime },t^{\prime },v^{\prime })\in \mathcal{N}}{x}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}}{a}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}

(35)-{\displaystyle \sum _{(s^{\prime },t^{\prime },v^{\prime })\in \mathcal{N}}{x}_{s^{\prime },t^{\prime },v^{\prime },s,t,v}}{a}_{s^{\prime },t^{\prime },v^{\prime },s,t,v}\mathrm{\Delta }t\le \mathit{\gamma },

(36)x\in \{\mathrm{0},\mathrm{}\mathrm{1}\},

(37)\phi \in R^{+}.

In the above formulation, Formulas (32) and (33) aim to minimize the critical value of the random return Z(X,\mathrm{}\omega ), which implies that the scenario-based returns are smaller than the optimal objective at least with a probability confidence level \alpha. Formula (34) is the flow balance constraint, which ensures that the feasible solution is an acyclic path from the origin node to the destination node in the PTS network. Inequality (35) indicates that the acceleration variation between two connected nodes cannot exceed \mathit{\gamma } to reduce the jerk of the train and meet the comfort requirements of passengers (Yin et al., 2014).

The proposed model (Formulas (32)–(37)) is typically nonlinear and nonconvex, which will probably lead to the increase in complexity in the solution process. Next, we provide a method to transform it into a linear form through introducing two types of decision variables according to Yang and Feng (2007).

Proposition 4.1: PRM is equivalent to the following linear programming model (denoted as LPRM), in which M is a sufficiently large number:

(38)\min L

(39)\mathit{s.t.}{\displaystyle \sum _{(s,t,v,s^{\prime },t^{\prime },v^{\prime })\in \mathcal{A}}{e}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}^{\omega }}{x}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}\le L+y_{\omega }M,\omega \in \mathrm{\Omega },

(40){\displaystyle \sum _{\omega \in \Omega }{P}_{\omega }}{y}_{\omega }\le \mathrm{1}-\alpha ,

(41)L\in R^{+},

(42)x_{s,t,v,s^{\prime },t^{\prime },v^{\prime }},y_{\omega }\in \{\mathrm{0},\mathrm{}\mathrm{1}\}.

Proof: Let (X^{*},\mathrm{}{\phi }^{*}) be the optimal solution for PRM, and set Z(X^{*},\mathrm{}\mathrm{1})\le Z(X^{*},\mathrm{}\mathrm{2})\le \mathrm{...}\le Z(X^{*},\mathrm{}\left|\Omega \right|). Then, we have {\phi }^{*}=Z(X^{*},\mathrm{}{k}^{\prime }), where k^{\prime }=\min \{\widehat{k}{\sum }_{\omega =\mathrm{1}}^{\omega =k}{\displaystyle P_{\omega }\ge \alpha }\}. Next, let y_{\omega }=\mathrm{1} when \omega \ge k^{\prime } and y_{\omega }=\mathrm{0} when , then constraints (39) and (40) hold. Therefore, ({\phi }^{*},\mathrm{}{X}^{*},\mathrm{}Y) is a feasible solution for LPRM (Y=\{y_{\mathrm{1}},\mathrm{}{y}_{\mathrm{2}},\mathrm{}\mathrm{...},\mathrm{}{y}_{\left|\Omega \right|}\}). Denote the optimal solution for LPRM by (L^{**},\mathrm{}{X}^{**},\mathrm{}{Y}^{**}); thus, we have

(43){\phi }^{*}\ge L^{**}.

Here, we set B=\left\{\omega \right|y_{\omega }=\mathrm{0}\} and B^c=\left\{\omega \right|y_{\omega }=\mathrm{1}\}. In accordance with Formula (40),

(44)\mathrm{Pr}\left\{B^c\right\}={\displaystyle \sum _{\omega \in \Omega }{P}_{\omega }}{y}_{\omega }\le \mathrm{1}-\alpha .

Hence,

(45)\mathrm{Pr}\left\{B\right\}=\mathrm{1}-\mathrm{Pr}\left\{B^c\right\}\ge \alpha .

Then, we have

(46)\mathrm{Pr}\{{\displaystyle \sum _{(s,t,v,s^{\prime },t^{\prime },v^{\prime })\in {\mathcal{A}}^{*}}{e}_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}^{\omega }}{x}_{ij}\le L^{**}\}\ge \mathrm{Pr}\left\{B\right\}\ge \alpha ,

which indicates that (L^{**},\mathrm{}{X}^{**},\mathrm{}{Y}^{**}) is a feasible solution of PRM. Accordingly,

(47){\phi }^{*}\le L^{**}.

From Formulas (43) and (47), we have {\phi }^{*} equals to L^{**}, indicating that PRM is equivalent to LPRM.

4.2 Complexity analysis of the proposed formulation

The variables in the two proposed models fall into two categories. The first type refers to the binary variables, i.e., X=\left\{x_{s,t,v,s^{\prime },t^{\prime },v^{\prime }}\right|(s,\mathrm{}t,\mathrm{}v,\mathrm{}{s}^{\prime },\mathrm{}{t}^{\prime },\mathrm{}{v}^{\prime })\in \mathcal{A}\}, which determines the train speed profile in the PTS network, and Y=\{y_{\mathrm{1}},\mathrm{}{y}_{\mathrm{2}},\mathrm{}\mathrm{...},\mathrm{}{y}_{\left|\Omega \right|}\}, which linearizes PRM into LPRM. The second type refers to the continuous variables, i.e., L and \phi, which are the critical values of the scenario-based variable Z(X,\mathrm{}\omega ), \omega \in \Omega, involved in LPRM. Evidently, PRM is a nonlinear model because Constraint (33) is a set of nonlinear inequalities. By contrast, Constraints (34), (35), (39), and (40) are linear equalities or inequalities. Therefore, LPRM is a MILP model, which can be solved using existing solvers, e.g., CPLEX.

Next, the complexity of the proposed model will be discussed. Table 5 lists the total number at most of the variables or constraints involved in PRM and LPRM, respectively, where |\cdot | represents the cardinality of a set. For example, in accordance with the definition, the total number at most of the binary variable X is \left|\mathcal{A}\right| in both models, where \left|\mathcal{A}\right| is the number of the arcs in the scenario-based PTS network. The complexity of these models mainly depends on the size of the network and the number of scenarios, as shown in Table 5, which confirms the importance of network reduction and scenario-based network construction.

5 Numerical experiments

In this section, a series of numerical experiments based on real-world data in Beijing metro is conducted to verify the effectiveness and performance of the proposed approaches. The experiments are conducted using a laptop with 3.79 GHz CPU and 8 GB memory, and the models are solved using CPLEX 12.8.

Here, we consider a metro railway segment on the Yizhuang line (from Rongjingdongjie to Wanyuanjie station) of Beijing metro (Fig. 5 on the left). The segment length is 1280 m, and different positions in this segment have distinct speed limits in accordance with real-world operational restrictions. Specifically, the speed limit of sections [0 m, 130 m] and [1096 m, 1280 m] is 57.6 km/h (16 m/s), that of section [131 m, 950 m] is 79.2 km/h (22 m/s), and that of section [951 m, 1095 m] decreases linearly between 79.2 km/h (22 m/s) and 57.6 km/h (16 m/s), as shown in Fig. 5 on the right. The track gradient is also presented in this figure. The total weight of the train is set to 270 t, consisting of six vehicles (each with a weight of 40 t) and a number of passengers (approximated as 30 t). The basic resistance parameter values in different scenarios are presented in Table 6.

In the next subsections, we conduct three groups of experiments. The first group of experiments is performed to analyze the effects of different scenarios and confidence levels \alpha on the optimal speed profile. The second group of experiments is considered to test the effects of discretization of the PTS network on the performance of the optimal train speed profile. The last group of experiments tests the robust energy-efficient train speed profile with different travel times.

5.1 Experiments on different scenarios and confidence levels

With the purpose of studying the influence of different scenarios and confidence levels on the optimal speed profile, the following experiments are performed on the PTS network with the discretization of \mathrm{\Delta }t=\mathrm{10}, \mathrm{\Delta }s=\mathrm{2.5}, and \mathrm{\Delta }v=\mathrm{0.5}. The total running time of the train between two stations is set to 90 s.

First, to show the influence of different scenarios (uncertain parameters) on the optimal speed profile, we obtain the optimal speed profiles in each scenario (Scenarios 1, 3, 6, 8, 10 in Table 6), and the results are shown in Fig. 6.

From the figure, the optimal speed profiles are relatively different in diverse scenarios. This result indicates that the uncertainty of the basic resistance parameters may influence the optimization of the train speed profile, which shows the necessity of generating a robust train speed profile. Then, to verify the robustness of the designed method, we conduct the following set of experiments, in which five scenarios (Scenarios 1, 3, 6, 8, 10 in Table 6) are considered and the confidence level \alpha is set between 0 and 1.

Table 7 lists the energy consumption in each scenario (using the optimized train speed profile as input), the objective value of the model, and the probability of scenarios that the energy consumptions do not exceed the objective value. As shown in Table 7, the energy consumption of all scenarios is no more than the objective value of the model with confidence level \alpha, i.e., the proposed method is able to obtain a robust energy-efficient train speed profile despite the uncertainty of basic resistance parameters. The decision maker can set the confidence level \alpha in accordance with the actual situation to obtain an optimal robust train speed profile based on all scenarios.

For further analysis, we conduct the following two sets of experiments with different values of confidence level \alpha (from 1.0 to 0.1) under 5 scenarios (Scenarios 1, 3, 6, 8, 10 in Table 6) and 10 scenarios (Scenarios 1–10 in Table 6), respectively.

The results of the two sets of experiments are shown in Table 8, where the objective value of the model, the expected energy consumption (EEC) of the optimal solution, and the computation time (CT) in different experiments are reported. All the instances can be solved efficiently within several seconds. The objective value of the model decreases as the confidence level decreases in these two sets of experiments, which is consistent with the actual situation. An interesting phenomenon is that the EEC first decreases and then increases as the confidence level decreases. The minimum EEC occurs when the confidence level is between 0.4 and 0.6. This phenomenon indicates that the decision maker should properly set the confidence level in accordance with the random distribution of parameters to reduce the energy consumption of trains.

To study the solutions further, Fig. 7 presents the optimal speed profiles obtained from the experiments with 5 scenarios and different confidence levels. From the figure, the following conclusions are drawn: 1) these speed curves are all composed of acceleration, cruising, coasting, and braking operation sequences but with different transition points, which is consistent with the conclusions of existing studies that these sequences lead to the minimum energy consumption (Su et al., 2013); and 2) the optimal speed profiles obtained with distinct confidence levels are also different. In addition, the optimal speed profile generated using the traditional expected value model is also shown in Fig. 7 (the black curve), which is similar to those when the confidence levels are 0.4 and 0.6. That is, the optimal solution of the proposed model is similar to that of the expected value model if the confidence level is set between 0.4 and 0.6.

5.2 Experiments on different network discretizations

In our formulation, the discretization of the PTS network is directly related to the quality of the solution and computation time. To explore the effects of the PTS network discretization, a set of experiments is conducted as follows. The experiments consider five scenarios (Scenarios 1, 3, 6, 8, 10 in Table 6) with different values of \mathrm{\Delta }s, \mathrm{\Delta }t, and \mathrm{\Delta }v. The other parameters are the same as those in the former experiments. To compare the performance of different instances, we also denote a new parameter \mathrm{\Delta }a=\mathrm{\Delta }v/\mathrm{\Delta }t to indicate the discretization of the PTS network.

As shown in Table 9, the discretization of the PTS network has a direct impact on EEC and CT. An evident tendency is that smaller \mathrm{\Delta }a reduces the energy consumption on account of more choices for acceleration value. When the values of \mathrm{\Delta }a are the same, the solution with a smaller length of \mathrm{\Delta }v often has less energy consumption. With respect to CT, it is directly related to the scale of the PTS network, given that a larger PTS network has more integer variables that increase the difficulty for solving the MIP model.

From the results, we have the following observations in designing the PTS network: 1) a smaller value of \mathrm{\Delta }a is required to generate a better quality of energy-efficient train speed profile; 2) the CT increases exponentially with a larger PTS network; and 3) the solution quality and CT can be compromised properly by setting appropriate values of \mathrm{\Delta }s, \mathrm{\Delta }t, and \mathrm{\Delta }v. For example, in Table 9, network with \mathrm{\Delta }s=\mathrm{2.5}, \mathrm{\Delta }t=\mathrm{10}, \mathrm{\Delta }v=\mathrm{0.5} (bold in the table) is the best choice in consideration of solution quality and CT.

5.3 Experiments on different total running times

As mentioned above in this paper, the total running time of the train between two adjacent stations is a major factor that influences the speed profile. Therefore, we conduct the following experiments by adjusting the total running time between two stations. The results are shown in Table 10 and Fig. 8. Table 10 lists the EEC of the optimal speed profile, the number of nodes and arcs of the reduced PTS network, and the CT of the experiments (with total running time varying from 86 s to 96 s). As shown in the table, the running time of the train has a considerable effect on the optimal train speed profile. The EEC decreases as the total running time increases, which is consistent with the actual operation. This phenomenon is also a suggestion for rail managers that a small increase in running time (e.g., 2 s) is seamless for passengers but can significantly reduce the energy consumption (by nearly 10%). On the contrary, the scale of the PTS network (number of nodes and arcs) increases when the total running time increases, which leads to a significant increase in CT.

In Fig. 8, the acceleration process of all speed profiles in the first 30 s (position 1–200 m) is the same. Longer total running time increases the time of coasting process. The optimal speed profile with longer total running time starts the breaking process later.

6 Conclusions

This work studied the energy-efficient train speed profile optimization problem under uncertain parameters. Particularly, we addressed the uncertainty of basic resistance parameters, which are critical factors affecting the optimized train speed profile. To handle the effects of these uncertainties, we constructed a scenario-based PTS network and proposed a series of node and arc reduction strategies. With the PTS network, the problem was transformed into a shortest path problem with the lowest energy consumption. Moreover, PRM was proposed, in which the basic resistance parameters were regarded as discrete scenario-based random variables. The proposed nonlinear model was equivalently reformulated into a linear solvable model. Lastly, we conducted a series of experiments based on real-world data in Beijing metro and analyzed the influence of relevant parameters on the optimal train speed profile. The results showed that the proposed approach was able to obtain a robust energy-efficient train speed profile. To balance energy saving and robustness, the decision maker should properly set the confidence level in accordance with the random distribution of the uncertain parameters. The results provided a suggestion for rail managers that a small increase in running time (e.g., 2 s) was seamless for passengers but could significantly reduce the energy consumption (by nearly 10%).

Further research will focus on the following three major aspects. (1) In this paper, we assume that the basic resistance parameters depend on one another on the basis of the actual situation without proving it. In the future, we will focus on analyzing the underlying relationships among these different parameters. (2) Another research direction is to construct stochastic programming models for multiple trains in consideration of the regenerative energy caused by breaking trains, in which the regenerative energy can be fed back into the overhead contact line and immediately used by accelerating trains located in the same electricity supply interval. The uncertainties may also involve the train arrival and departure times. And (3) the train and locomotive schedules also have a direct influence on the energy consumption of urban rail transit systems. Therefore, our future research will consider train schedule and locomotive circulation to improve the energy efficiency of urban rail transit systems further.

References

[1] Albrecht A, Howlett P, Pudney P, Vu X, Zhou P (2016a). The key principles of optimal train control—part 1: Formulation of the model, strategies of optimal type, evolutionary lines, location of optimal switching points. Transportation Research Part B: Methodological, 94: 482–508

[2] Albrecht A, Howlett P, Pudney P, Vu X, Zhou P (2016b). The key principles of optimal train control—part 2: Existence of an optimal strategy, the local energy minimization principle, uniqueness, computational techniques. Transportation Research Part B: Methodological, 94: 509–538

[3] Amrani A, Hamida A, Liu T, Langlois O (2018). Train speed profiles optimization using a genetic algorithm based on a random-forest model to estimate energy consumption. In: Proceedings of 7th Transport Research Arena. Vienna, hal-01767006

[4] Bešinović N, Quaglietta E, Goverde R M P (2013). A simulation-based optimization approach for the calibration of dynamic train speed profiles. Journal of Rail Transport Planning and Management, 3(4): 126–136

[5] Calderaro V, Galdi V, Graber G, Piccolo A, Cogliano D (2014). An algorithm to optimize speed profiles of the metro vehicles for minimizing energy consumption. In: International Symposium on Power Electronics, Electrical Drives, Automation and Motion. Ischia:IEEE, 813–819

[6] de Martinis V, Corman F (2018). Data-driven perspectives for energy efficient operations in railway systems: Current practices and future opportunities. Transportation Research Part C: Emerging Technologies, 95: 679–697

[7] Fernández-Rodríguez A, Fernández-Cardador A, Cucala A P (2018). Balancing energy consumption and risk of delay in high speed trains: A three-objective real-time eco-driving algorithm with fuzzy parameters. Transportation Research Part C: Emerging Technologies, 95: 652–678

[8] Gao Z, Yang L (2019). Energy-saving operation approaches for urban rail transit systems. Frontiers of Engineering Management, 6(2): 139–151

[9] Howlett P G (2000). The optimal control of a train. Annals of Operations Research, 98(1/4): 65–87

[10] Howlett P G, Pudney P J (1995). Energy-efficient Train Control: Advances in Industrial Control. London: Springer

[11] Huang K, Wu J, Yang X, Gao Z, Liu F, Zhu Y (2019). Discrete train speed profile optimization for urban rail transit: A data-driven model and integrated algorithms based on machine learning. Journal of Advanced Transportation, (4): 7258986

[12] Ke B R, Lin C L, Lai C W (2011). Optimization of train-speed trajectory and control for mass rapid transit systems. Control Engineering Practice, 19(7): 675–687

[13] Kim K, Chien S I (2011). Optimal train operation for minimum energy consumption considering track alignment, speed limit, and schedule adherence. Journal of Transportation Engineering, 137(9): 665–674

[14] Ko H, Koseki T, Miyatake M (2004). Application of dynamic programming to the optimization of the running profile of a train. WIT Transactions on the Built Environment, 74: 103–112

[15] Li X, Gao Z, Sun W (2013a). Existence of an optimal strategy for stochastic train energy-efficient operation problem. Soft Computing, 17(4): 651–657

[16] Li X, Li L, Gao Z, Tang T, Su S (2013b). Train energy-efficient operation with stochastic resistance coefficient. International Journal of Innovative Computing, Information and Control, 9(8): 3471–3483

[17] Liu B (2002). Random fuzzy variables. In: Liu B, ed. Theory and Practice of Uncertain Programming. Heidelberg: Physica-Verlag, Springer, 295–308

[18] Liu R, Golovitcher I (2003). Energy-efficient operation of rail vehicles. Transportation Research Part A: Policy and Practice, 37(10): 917–932

[19] Lu S, Hillmansen S, Ho T K, Roberts C (2013). Single-train trajectory optimization. IEEE Transactions on Intelligent Transportation Systems, 14(2): 743–750

[20] Mahmoudi M, Zhou X (2016). Finding optimal solutions for vehicle routing problem with pickup and delivery services with time windows: A dynamic programming approach based on state–space–time network representations. Transportation Research Part B: Methodological, 89: 19–42

[21] Scheepmaker G M, Goverde R M P, Kroon L G (2017). Review of energy-efficient train control and timetabling. European Journal of Operational Research, 257(2): 355–376

[22] Shangguan W, Yan X H, Cai B G, Wang J (2015). Multiobjective optimization for train speed trajectory in CTCS high-speed railway with hybrid evolutionary algorithm. IEEE Transactions on Intelligent Transportation Systems, 16(4): 2215–2225

[23] Sicre C, Cucala A P, Fernández-Cardador A (2014). Real time regulation of efficient driving of high speed trains based on a genetic algorithm and a fuzzy model of manual driving. Engineering Applications of Artificial Intelligence, 29: 79–92

[24] Su S, Li X, Tang T, Gao Z (2013). A subway train timetable optimization approach based on energy-efficient operation strategy. IEEE Transactions on Intelligent Transportation Systems, 14(2): 883–893

[25] Wang L, Yang L, Gao Z, Huang Y (2017). Robust train speed trajectory optimization: A stochastic constrained shortest path approach. Frontiers of Engineering Management, 4(4): 408–417

[26] Wang P, Goverde R M P (2016). Multiple-phase train trajectory optimization with signaling and operational constraints. Transportation Research Part C: Emerging Technologies, 69: 255–275

[27] Wang P, Goverde R M P (2019). Multi-train trajectory optimization for energy-efficient timetabling. European Journal of Operational Research, 272(2): 621–635

[28] Wang P, Trivella A, Goverde R M P, Corman F (2020). Train trajectory optimization for improved on-time arrival under parametric uncertainty. Transportation Research Part C: Emerging Technologies, 119: 102680

[29] Wang Y, de Schutter B, van den Boom T J J, Ning B (2013). Optimal trajectory planning for trains: A pseudospectral method and a mixed integer linear programming approach. Transportation Research Part C: Emerging Technologies, 29: 97–114

[30] Yang L, Feng Y (2007). A bicriteria solid transportation problem with fixed charge under stochastic environment. Applied Mathematical Modelling, 31(12): 2668–2683

[31] Yang L, Zhou X (2014). Constraint reformulation and a Lagrangian relaxation-based solution algorithm for a least expected time path problem. Transportation Research Part B: Methodological, 59(1): 22–44

[32] Yang L, Zhou X (2017). Optimizing on-time arrival probability and percentile travel time for elementary path finding in time-dependent transportation networks: Linear mixed integer programming reformulations. Transportation Research Part B: Methodological, 96: 68–91

[33] Yin J, Chen D, Li L (2014). Intelligent train operation algorithms for subway by expert system and reinforcement learning. IEEE Transactions on Intelligent Transportation Systems, 15(6): 2561–2571

[34] Yin J, Su S, Xun J, Tang T, Liu R (2020). Data-driven approaches for modeling train control models: Comparison and case studies. ISA Transactions, 98: 349–363

[35] Yin J, Tang T, Yang L, Xun J, Huang Y, Gao Z (2017a). Research and development of automatic train operation for railway transportation systems: A survey. Transportation Research Part C: Emerging Technologies, 85: 548–572

[36] Yin J, Yang L, Tang T, Gao Z, Ran B (2017b). Dynamic passenger demand oriented metro train scheduling with energy-efficiency and waiting time minimization: Mixed-integer linear programming approaches. Transportation Research Part B: Methodological, 97: 182–213

[37] Zhou L, Tong L, Chen J, Tang J, Zhou X (2017). Joint optimization of high-speed train timetables and speed profiles: A unified modeling approach using space-time-speed grid networks. Transportation Research Part B: Methodological, 97: 157–181

RIGHTS & PERMISSIONS

Higher Education Press