Next Article in Journal
Bilinear Quadratic Feedback Control of Modular Multilevel Converters
Previous Article in Journal
Optimizing Lifespan of Circular Products: A Generic Dynamic Programming Approach for Energy-Using Products
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Energy-Efficient Train Driving Based on Optimal Control Theory

by
Wolfram Heineken
1,*,
Marc Richter
1 and
Torsten Birth-Reichert
1,2
1
Fraunhofer Institute for Factory Operation and Automation IFF, Sandtorstraße 22, 39106 Magdeburg, Germany
2
Hochschule für Angewandte Wissenschaften, Berliner Tor 5, 20099 Hamburg, Germany
*
Author to whom correspondence should be addressed.
Energies 2023, 16(18), 6712; https://doi.org/10.3390/en16186712
Submission received: 10 August 2023 / Revised: 11 September 2023 / Accepted: 12 September 2023 / Published: 19 September 2023
(This article belongs to the Section E: Electric Vehicles)

Abstract

:
Efficient train driving plays a vital role in reducing the overall energy consumption in the railway sector. An energy minimising control strategy can be computed using the framework given by optimal control theory; in particular, the Pontryagin maximum principle can be used. Our optimisation approach is based on an algorithm presented by Khmelnitsky that considers electric trains equipped with regenerative braking. A derivation of Khmelnitsky’s theory from a more general formulation of the maximum principle is given in this article, and a complete list of switching cases between different driving regimes is included that is essential for practical application. A number of numerical examples are added to visualise the various switching cases. Energy consumption data from real-life operation of passenger trains are compared to the calculated energy minimum. In the presented study, the optimised strategy was able to save 37 percent of the average energy demand of the train in operation. The sensitivity of the energy consumption to deviations of the train speed from the optimum speed profile is studied in an example. Another example illustrates that the efficiency of regenerative braking has an effect on the optimum speed profile.

1. Introduction

Energy-efficient driving strategies have been shown to substantially lower the total energy consumption of railway trains. With an increased demand for energy saving technologies, numerous research studies have been devoted to this subject in recent years, showing that, in many cases, an average energy saving from 15 to 30% could be achieved by assisted driving [1,2,3,4,5,6,7,8]. A comprehensive review on energy-efficient train control, commonly abbreviated as EETC, is given by Scheepmaker, Goverde, and Kroon [9]. The majority of the mathematical models that have been developed to calculate energy-minimised driving are based on Pontryagin’s maximum principle (PMP), a fundamental theorem of optimal control theory. The maximum principle was formulated in 1961 by Pontryagin, Boltyanskii, Gamkrelidze, and Mishchenko [10,11]. The first application of Pontryagin’s maximum principle to the operation of trains dates back to Japanese research in 1968 (Ichikawa [12]). The model of Ichikawa was restricted to a flat track and linear train resistance depending on speed. Another early contribution is due to Strobel, Horn, and Kosemund [1]. They considered a non-constant altitude and found that the optimal solution is always in one of the drive regimes, namely full traction, partial traction with constant speed, coasting, partial braking with constant speed, and full braking. Based on their model, they set up a first driver advisory system for the Berlin S-Bahn (suburban trains).
Since those pioneering research on this topic, models have gained complexity. Important developments have been: (1) the application to variable altitude and variable speed limit, (2) the inclusion of regenerative braking in electric locomotives or motor coaches, (3) improved modelling for train resistance, tractive, and brake force depending on speed, (4) the energy minimisation with non-zero speed boundary conditions, (5) advanced modelling of efficiency as dependent on speed, tractive, and brake force, (6) the modelling of discrete throttle settings for diesel engines instead of continuous traction, and (7) multiple-train problems.
Regarding the second aspect, some (especially the earlier) authors are dealing only with the mechanical braking of trains. When regenerative braking is considered, most authors base their energy-efficient driving strategy on the exclusive use of the regenerative brake. A few recent contributions cover the combined application of the mechanical and regenerative brake in their optimisation strategy.
We proceed with some important developments for purely mechanical braking. Starting from the Ph.D. thesis of Milroy [13] in 1980, a research team at the University of South Australia introduced the system Metromiser (Howlett and Pudney [2]). The system covered both timetable planning and driver advisory. Metromiser was restricted to constant track gradients during both coasting and braking. The system was first applied on trains in Adelaide (Australia) in 1984, followed by Toronto (Canada), Melbourne (Australia), and Brisbane (Australia). Howlett, Milroy, and Pudney [14] were the first to assume a variable speed limit, but only for the special case of discrete traction control. In 2003, Liu and Golovitcher [15] presented a method covering variable altitude, variable speed limit, continuous traction control, and quadratic train resistance dependence on speed. Except being bound to purely mechanical braking, this algorithm is therefore quite general with respect to the other assumptions. The article contains a complete description of switching between the possible drive regimes. In 2013, Albrecht, Howlett, Pudney, and Vu [16] introduced the driver advisory system Energymiser that was applied by the French railway SNCF on TGV high-speed trains. They proved the uniqueness of the energy optimum being derived. A multi-objective algorithm combining energy minimisation and the punctuality of trains was developed by Aradi, Becsi, and Gaspar [6].
Regenerative braking was first considered in 1985 by Asnis, Dmitruk, and Osmolovskii [17] in the context of train energy minimisation. They assumed a constant altitude and speed limit in their model. Franke, Terwiesch, and Meyer [3] introduced an energy optimisation algorithm for variable altitude, variable speed limit, and purely regenerative braking that was restricted to piecewise constant traction and brake force. They suggested a discrete dynamic programming (DDP) method as a solution. In 2000, Khmelnitsky [18] developed a model that covered variable altitude, variable speed limit, continuous traction control, arbitrary maximum tractive and brake force depending on speed, quadratic train resistance with respect to speed, and exclusively regenerative braking. Regarding these features, they provide broadly applicable approaches with few restrictions. Due to these advantages, the algorithm of Khmelnitsky was chosen as the basis for our research. In a recent article by Ying et al. [19], the EETC problem solved by Khmelnitsky was revisited, with a special focus on solutions reaching the speed limit. The paper of Ying et al. [19] can be especially recommended for its comprehensive illustration of the large variety of cases that can arise during the construction of the optimum solution. The EETC problem with prescribed non-zero speed boundary conditions was solved by Ying et al. [20].
There are relatively few research contributions dealing with the EETC problem using combined mechanical and regenerative braking. Baranov, Meleshin, and Chin’ [21] were the first considering this topic, but an algorithm to solve the problem was left as an open question in their article. Lu et al. [22] studied combined mechanical and regenerative braking, but excluded both variable speed limit and variable track gradients. Zhou et al. [23] studied the synchronisation of accelerating and braking trains, considering both kinds of braking, but again did not include the variable speed limit. Fernández-Rodríguez et al. [24] combined both brake systems in a multi-objective optimisation method, but did not derive a rigorous energy minimum. In the paper of Scheepmaker and Goverde [25], the optimisation problem with combined mechanical and regenerative braking was solved for the first time under general assumptions, i.e., with a speed dependent tractive and brake effort, variable speed limit, and variable altitude. Scheepmaker and Goverde used the Gauss–Radau pseudospectral method implemented in GPOPS [26] to solve the optimal control problem.
Several authors considered advanced efficiency modelling, in particular efficiency depending on speed, tractive, and brake force. This assumption is more realistic but substantially complicates the approach. Su, Tang, and Li [27], Zhu et al. [28], and Su et al. [29] modelled traction efficiency depending on traction force and speed; they apply the energy-distribution-based method (EDBM) to solve the problem. Research on advanced efficiency modelling was also carried out by Franke, Terwiesch, and Meyer [3], Ghaviha et al. [30], Kouzoupis et al. [31], and Feng, Huang, and Lu [32]. Discrete traction control, which is relevant for diesel engines with discrete throttle settings, was considered by Howlett, Milroy, and Pudney [14]; Cao et al. [33]; and Zhang, Cao, and Su [34]. Multiple train problems were often been studied with a focus on timetable design in order to synchronise accelerating and braking trains for the best distribution of regenerative energy. See, for example, Zhou et al. [23] for research on this topic. More references can be found in Scheepmaker and Goverde [25]. Szkopiński and Kochan [35] studied energy-efficient train driving when approaching a train on the track ahead. Su et al. [36] introduced a stabilisation strategy that avoids frequent switching between traction and braking on sections where the energy optimum would actually require a coasting regime. An algorithm for energy-efficient train driving with multiple time targets that can be used for timetable design was developed by Fernández-Rodríguez et al. [37].
Apart from the many research contributions based on an application of Pontryagin’s maximum principle, a variety of other methods exist to solve the EETC problem. In particular, the so-called direct methods have recently gained attention. They are characterised by the fact that the problem is first discretised, usually by dividing the track length into smaller intervals. The resulting nonlinear programming (NLP) problem is then solved by nonlinear programming methods as, for example, the pseudospectral one. Wang et al. [38] were the first to apply the pseudospectral method to an EETC problem. Scheepmaker and Goverde [25] also utilised a pseudospectral approach. Direct methods have the advantage of being very flexible with respect to the problem under consideration. In addition, optimal control solvers such as DIDO (Direct and Indirect Dynamic Optimization) [39,40,41], PSOPT [42,43], or GPOPS-II (General Pseudospectral Optimal Control Software) [26,44] are readily available. The implementation is much easier than a solution of the problem along the lines of Pontryagin’s maximum principle. On the other hand, the pseudospectral method often shows an inaccurate oscillatory behaviour of solutions, and it tends to be time consuming [31]. Recently, Kouzoupis et al. [31] were able to reduce the computing time of a direct method using multiple shooting.
Based on the above-mentioned approach of Khmelnitsky [18], a programme called Leda (low energy driver advisory) was developed at Fraunhofer Institute for Factory Operation and Automation (IFF), Magdeburg, Germany. The code is written in MATLAB and currently features energy minimisation with exclusively regenerative braking, except for the fastest train motion in the case of tight timetables where additional mechanical braking is considered. The choice of Khmelnitsky’s method was mainly motivated by superior accuracy at an acceptable computing time. An extension of the code to fully include mechanical braking into the optimisation is subject to further development, as well as non-zero speed boundary conditions and smooth switching from fast to energy-efficient driving in the case of a train delay. The code was tested offline in numerous cases based on real railway tracks and timetables, and the computing time was substantially reduced by code optimisation. A couple of tests with driver advisory in a real-life operation demonstrated energy savings of about 20% compared to the average of unassisted runs.
The main motivations of this article can be summarised in the following way:
  • It is intended to show a direct derivation of Khmelnitsky’s theory from a more general formulation of Pontryagin’s maximum principle given by Hartl, Sethi, and Vickson [45]. By seeing Khmelnitsky’s theory in this general framework, extensions to other conditions could be elaborated.
  • Our aim is to provide a comprehensive illustration of the behaviour of the method of Khmelnitsky using a couple of numerical examples. We felt that this very useful method would strongly benefit from some illustrative examples that show the behaviour of the trajectory field including kink points as well as the large variety of switching cases from one driving regime to another. We, therefore, put a strong focus on the examples in Section 2.4 and Section 2.5, in which most of the possible switching cases can be studied in detail.
  • Another aim of the article is to estimate the energetic potential of an optimal driving strategy by comparing the energy optimum to energy consumption in a real-life train operation. As a result of a study of 100 train runs on a railway line in Germany, the optimal strategy consumed only 63% of the average energy demand of the unassisted runs in a real operation.
  • The sensitivity to deviations from the optimum speed profile was studied. This will always be of practical relevance, since in a real application the train driver will manually try to follow the suggestions of an advisory system and therefore will never exactly meet the optimum speed profile.
This article is structured in the following way: In Section 2.1, the model of the train motion is introduced. Section 2.2 is devoted to the Pontryagin maximum principle that provides a couple of necessary conditions the energy minimising solution has to fulfil. In Section 2.3, the maximum principle is applied to the specific minimum energy problem for train motion. This leads to the observation that only four driving regimes, full traction, constant speed, coasting, and full regenerative braking, are feasible for an energy-minimising motion. The regime constant speed can be driven only at two specific velocities or at the speed limit. Section 2.4 and Section 2.5 describe the algorithm of Khmelnitsky [18] for the construction of the minimum energy solution. A complete list of switching cases is given and explained in Section 2.4, and four numerical examples are introduced to illustrate those cases. Some directions are given that lead to a substantial reduction of computing time. Section 3 is devoted to the results. It contains the comparison of the optimum strategy to a real-life train operation in Section 3.1, the sensitivity to deviations from the speed optimum in Section 3.2, and finally a study that shows a slight effect of the regenerative brake efficiency on the optimum speed.

2. Materials and Methods

2.1. Model of Train Motion

The equation of motion of a train with an electric engine and both regenerative and mechanical brakes can be given by
F tr F br F mbr F air F roll F sl = c rot m a ,
where  F tr  is tractive force,  F br  is the force applied by the regenerative brake,  F mbr  is the force of the mechanical brake,  F air  is air drag,  F roll  is rolling friction,  F sl  is the downhill slope force,  c rot  is a coefficient representing rotating masses, m is the mass, and a is the acceleration of the train. This is Newton’s second law of motion with the extension of rotating masses, for example, the wheels of the train, are accounted for by a factor  c rot . A model of this type is commonly used for train motion, see, for example, the review article of Scheepmaker, Goverde, and Kroon [9] or the textbook of Ihme [46]. According to Ihme [46],  c rot = 1.06 1.11  for passenger trains, depending on their length, where longer trains will generally have smaller values of  c rot .
The tractive force  F tr  is limited by both engine power and rail friction (adhesion). According to Fassbinder [47],
F tr F tr , max = min ( P tr , mech / v , μ ad g m loc ) ,
where  P tr , mech  is the mechanical engine power used for traction, v is the velocity of the train,  μ ad  is the adhesion coefficient, g is gravity, and  m loc  is the mass of the locomotive. The adhesion coefficient  μ ad  was obtained experimentally in 1943 by Curtius and Kniffler [48], see, e.g., Schlecht [49], leading to the empirical relation
μ ad = 7.5 v km / h + 44 + 0.161 .
The adhesion coefficient  μ ad  attains a maximum value of 0.331 when  v = 0 .
The regenerative brake force  F br  is, as the tractive force, restricted by engine power and rail adhesion. However, it must be further limited to avoid a derailing of coaches behind the braking locomotive. Thus, for the regenerative brake force there holds
F br F br , max = min ( P br , mech / v , μ ad g m loc , F br , lim ) ,
where  P br , mech  is the mechanical engine power for regenerative braking, and  F br , lim  is the additional limit to avoid derailing. In general,  P tr , mech = P br , mech  will hold. According to [50], the limit force  F br , lim  was recently enlarged in Germany from 150 kN to 240 kN, a value that is also considered in Scheepmaker and Goverde [25].
The first term  P tr , mech / v  in the minimum in (2) is the engine power limit of traction force, and the second term  μ ad g m loc  is the rail friction limit. Likewise, the first term  P br , mech / v  in the minimum in (4) is the engine power limit of the brake force. If  P tr , mech = P br , mech  holds, then both engine power limits are equal. The maximum forces  F tr , max  and  F br , max  are shown in Figure 1.
Remark 1.
The correlation (3) for the adhesion coefficient was obtained experimentally under dry conditions. On wet rails, adhesion can be considered to be about 30% less [49]. The roughness of the rail and wheel, as well as the wheel diameter, will have an influence on adhesion, but this was not studied by Curtius and Kniffler [48]. Strictly speaking, the rail friction limit would be  μ ad g m loc cos ϕ , where ϕ is the slope angle of the track. However, even for extremely steep railway tracks, the value of  cos ϕ  would be very close to one, and the factor  cos ϕ  can therefore be neglected.
Remark 2.
In practical operation, the regenerative braking of trains is not possible when the speed is too small, meaning that for a final halt, the mechanical brake always has to be applied. This was pointed out by Scheepmaker and Goverde [25], who have estimated a minimum speed of 8 km/h for the application of regenerative braking based on data from Netherlands Railways. However, since kinetic energy below 8 km/h is relatively low, we have neglected this consideration in our model.
When a train brakes, using the regenerative brake is clearly advantageous with respect to saving energy. However, the force of the regenerative brake is limited and, especially for long freight trains, much weaker than the force of the mechanical brake. This is due to the fact that regenerative braking applies only to the wheels of the locomotive, while the mechanical brake acts on every wheel of the train. If the time required by the timetable is too short for a certain distance, purely regenerative braking might not be sufficient. However, as Scheepmaker and Goverde [25] have proven by means of optimal control theory, the mechanical brake is always ‘second choice’ in an energy-minimising solution. This means that mechanical braking is only applied when regenerative braking is operating at maximum force, i.e., when  F br = F br , max  holds. In the optimisation method presented here, the mechanical brake is not incorporated into the theory and will only be considered in the calculation of the fastest possible motion. This will be explained in more detail in the following sections.
The mechanical work applied for traction is given by
W tr = s 1 s 2 F tr d s = t 1 t 2 P tr , mech d t ,
where s is track length and t is time. The electric energy required for traction is
E tr = W tr / η tr .
The efficiency  η tr  of electric locomotives is usually given in a range from 83 to 87 percent [25,47,51]. Within this paper, we will use an intermediate value of  η tr = 0.85 . Likewise, the mechanical work applied for regenerative braking is
W br = s 1 s 2 F br d s = t 1 t 2 P br , mech d t ,
and the electric energy returned will be
E br = η br W br .
We use a braking efficiency  η br = η tr = 0.85  according to Fassbinder [47]. The difference  E net = E tr E br  is called the net energy, which should be minimised.
In a real operation, additional energy  E add  will be required that is not directly related to traction or brake, including energy for air compression, air conditioning, lighting and so on. However, since this additional energy is mainly a function of time, it cannot be reduced by a driving strategy when the total time is constant as given by the train’s timetable. Therefore, the additional energy does not enter the net energy minimisation problem and can be excluded here.
For air drag, Ihme [46] recommends the so-called Hanover Formula that goes back to Voß, Gackenholz, and Wiebels [52]. It takes the form
F air = 1 2 ρ air A ref c w v 2 ,
with air density  ρ air = 1.25   kg / m 3 . In this formula,  A ref  is not the cross-sectional area of the train, but a reference area of  10   m 2 . The air drag coefficient  c w  is calculated according to
c w = c w , loc + c w , first + ( n 2 ) c w , middle + c w , last ,
n being the number of coaches. Ihme [46] provides the values  c w , loc = 0.26 c w , first = 0.13 c w , middle = 0.10 c w , last = 0.23  as appropriate for Intercity coaches.
The rolling friction can be modelled by  F roll = c roll m g  with  c roll = 0.0015  according to Ihme [46]. The downhill slope force is given by
F sl = m g d z d s ,
where z is the altitude of the track.
Remark 3.
The models for the air drag and rolling friction that are used here are simplified relations taken from the literature. In particular, the air drag model considered here can be seen as a rather rough estimate since it does not consider the current wind conditions, which will be difficult to model but might have a significant effect. Some authors (as, for example, Scheepmaker and Goverde [25]) use an overall train resistance equation that was obtained experimentally based on a particular train. If reliable data are available for the specific train under consideration, train resistance modelling based on such data could be considered as a more accurate alternative.

2.2. The Maximum Principle

The problem of minimising the net energy of a scheduled train can be formulated as an optimal control problem, as it was proposed by Khmelnitsky [18], who also presented an algorithm to obtain the unique net energy minimum. The algorithm of Khmelnitsky is essentially based on the maximum principle of optimal control theory that was developed by Pontryagin, Boltyanskii, Gamkrelidze, and Mishchenko [10,11]. A variety of formulations of the maximum principle can be found in Hartl, Sethi, and Vickson [45]. While some of the results of the maximum principle for the particular problem of a minimum energy train ride are already given in Khmelnitsky’s paper [18], a direct derivation from the more general formulation of the maximum principle as presented in Hartl, Sethi, and Vickson [45] is not included in Khmelnitsky. We are, therefore, going to show this derivation here. In Section 2.2, the maximum principle is presented, based on the ‘Informal Theorem 4.1’ in [45]. The subsequent Section 2.3 contains the application to the train problem. Equipped with this general framework, extensions of the theory (for example, the inclusion of the mechanical brake) are possible.
We consider the following problem: Let s be the track length coordinate between two stations at  s = 0  and  s = S . The motion of the train is described by the time  t ( s )  it takes for the train to move from the first station to position s. We assume  t ( 0 ) = 0  and a fixed duration  t ( S ) = T  given by the timetable. The velocity v of the train is restricted by a piecewise constant function  v max ( s )  that accounts for speed limits on the track. Tractive and regenerative brake force,  F tr  and  F br , are limited according to (2) and (4), respectively. The optimisation problem is to find a motion, i.e., a function  t ( s ) , that minimises the net energy  E net .
We introduce some notation, closely following Khmelnitsky [18]. Let  E kin = c rot m v 2 / 2  be the kinetic energy of the train. Dividing by  c rot m  leads to a specific kinetic energy
K = E kin c rot m = v 2 2
limited by speed restrictions:  K ( s ) K max ( s ) = ( v max ( s ) ) 2 / 2 . Likewise, a specific potential energy is defined by
P = E pot c rot m = g z c rot .
A total specific mechanical energy is then given by  E = K + P . In the same way, specific forces are defined as
u tr = F tr c rot m ,
u br = F br c rot m ,
g tr = F tr , max c rot m ,
g br = F br , max c rot m ,
w = F air + F roll c rot m .
A total recuperation efficiency is given by  α = η tr η br . The aim of minimising the net energy  E net  is equal to maximising
J = η tr c rot m E net = 0 S u tr ( s ) α u br ( s ) d s .
With this notation, the optimisation problem can be written in the following canonical form:
  • State vector  x = ( E , t ) ;
  • Control vector  u = ( u tr , u br ) ;
  • State differential equations;
    d E d s = u tr u br w
    d t d s = 1 2 K
  • State boundary conditions;
    E ( 0 ) = P ( 0 )
    t ( 0 ) = 0
    E ( S ) = P ( S )
    t ( S ) = T
  • Maximum;
    J = 0 S u tr ( s ) α u br ( s ) d s max
  • State constraint  E K max + P ;
  • Control constraints  0 u tr g tr  and  0 u br g br .
Remark 4.
The state boundary conditions (22) and (24) are equivalent to zero speed boundary conditions  v ( 0 ) = v ( S ) = 0 .
Remark 5.
Further constraints could be added here, especially  0 t T  and  K 0 . However,  K > 0  will always hold for the energetic optimum solution inside the interval  ( 0 , S ) . The condition  0 t T  follows from (21)–(23). Therefore, conditions  0 t T  and  K 0  need not be set explicitly as a constraint. The simultaneous traction and brake, i.e.,  u tr ( s ) > 0  and  u br ( s ) > 0  for the same s, is also not explicitly prevented by a constraint since such a solution will surely not be energetically optimal.
Remark 6.
In Equation (20), we have excluded mechanical braking. This means that, in accordance with Khmelnitsky [18], we formulate the optimisation problem only for purely regenerative braking. An extension of the presented theory to the case with a simultaneous application of both regenerative and mechanical brake is possible in a straightforward manner, and is, in a slightly different formulation, given by Scheepmaker and Goverde [25].
Introducing the notation
f ( x , u , s ) = u tr u br w 1 / 2 K ,
F ( u ) = u tr + α u br ,
g ( x , u , s ) = u tr u br g tr u tr g br u br ,
h ( x , s ) = K max K ,
we have an optimisation problem of the form given in Hartl, Sethi, and Vickson [45]:
d x d s = f ( x , u , s ) ,
J = 0 S F ( u ( s ) ) d s max ,
g ( x , u , s ) 0 ,
h ( x , s ) 0 .
In Hartl, Sethi, and Vickson [45], the problem is formulated with a time t being the independent variable instead of s. This is due to the fact that many practical optimisation problems are formulated in a time-dependent way. In our case, however, the position-dependent formulation has some advantages, especially since the speed restriction also depends on position, not on time.
The function
H = F + λ T f
with Lagrange multipliers  λ ( s )  is called the Hamiltonian of problem (31)–(34). (Here and in the following,  ( · ) T  means the transposed of a column vector, i.e.,  λ T f  is the scalar product of vectors  λ  and  f .) The Lagrange multipliers  λ ( s )  are also called the costates of the problem. In order to agree with the notation of Khmelnitsky [18], we denote the costates according to
λ ( s ) = ψ ( s ) ψ T ( s ) .
Furthermore, the Lagrangian is introduced by
L = H + μ T g + ν h ,
where  μ ( s )  and  ν ( s )  are also called Lagrange multipliers. The Lagrange multipliers  λ ( s ) μ ( s )  and  ν ( s )  are continuous functions in s, except for positions  s i  where  h ( x , s i ) = 0 .
The ‘Informal Theorem 4.1’ of Hartl, Sethi, and Vickson [45] states the following necessary conditions for  ( x , u )  to be a solution of the optimisation problem (31)–(34):
  • The control vector  u ( s )  maximises the Hamiltonian H pointwise for any  s [ 0 , S ] , i.e.,
    u opt = arg max u H ( x opt , u , λ opt , s )
    is fulfilled. (Here, the subscript ‘opt’ corresponds to the optimal solution.)
  • There holds
    L u = H u + μ T g u = 0 .
  • The costate equation
    d λ d s = L x
    is satisfied.
  • Let  μ i  be the components of vector  μ , and  g i  be the components of vector  g . There holds
    μ i ( s ) 0 and μ i ( s ) g i ( s ) = 0
    for all i and all  s [ 0 , S ] . This is called complementary slackness in Khmelnitsky [18].
  • Another complementary slackness condition,
    ν ( s ) 0 and ν ( s ) h ( s ) = 0 ,
    holds for all  s [ 0 , S ] .
  • The following jump condition is satisfied: At every point  s i  where  λ  is discontinuous, there exists a number  η i  with
    λ ( s i + ) λ ( s i ) = η i h x ,
    H ( s i + ) H ( s i ) = η i h s ,
    η i 0 and η i h ( s i ) = 0 .
    Here, the argument  s i  corresponds to the left-hand limit, and  s i +  to the right-hand limit at  s i . The vector Equation (43) is meant to be component-wise.
The Hamiltonian maximisation condition (38) is called the Pontryagin Maximum Principle, and Equations (39)–(42) are referred to as the Karush–Kuhn–Tucker (KKT) conditions [53,54]. In Hartl, Sethi, and Vickson [45], the maximum principle is presented also for the case of multiple state constraints, i.e., when h is extended to a vector. The maximum principle does not, as we shall see, tell the solution to the energy minimisation problem directly, but it provides essential information such that a construction of the solution becomes possible in an iterative trial-and-error process.

2.3. Application of the Maximum Principle to the Energy Minimisation Problem

From the costate Equation (40), it follows for the second component of  λ
d ψ T d s = L t = 0 ,
since L does not depend explicitly on t. Equation (43) provides, again for the second vector component,
ψ T ( s i + ) ψ T ( s i ) = η i h t = 0 .
Thus, the Lagrange multiplier  ψ T  is a constant. Equation (39) results in
L u tr = 1 + ψ + μ 1 μ 3 = 0 ,
L u br = α ψ + μ 2 μ 4 = 0 .
From the costate Equation (40), first component, we have
d ψ d s = L E = ψ d w d K + ψ T · ( 2 K ) 3 / 2 μ 3 d g tr d K μ 4 d g br d K + ν .
When travelling below the speed limit, i.e.,  h > 0 , then  ν = 0  holds due to (42). (In Khmelnitsky’s notation,  μ 3  and  μ 4  are called  a tr  and  a b , respectively.) We now apply a case distinction to the multiplier  ψ , and evaluate  u tr u br , and the  μ i  using the conditions (38), (41), (48), (49), and the continuity of the multipliers. The results are given in Table 1.
From the values given in Table 1 it follows that, in the interior of the state domain,  μ 3  and  μ 4  can be directly expressed in terms of  ψ :
μ 3 ( s ) = ψ ( s ) 1 if ψ ( s ) 0 0 if ψ ( s ) < 0
μ 4 ( s ) = α ψ ( s ) if ψ ( s ) α 0 if ψ ( s ) > α
The cases indicated in Table 1 are the only driving modi that are possible for the optimum solution of the problem. This means that the interval  [ 0 , S ]  can be completely segmented into subintervals  S i = ( s i 1 , s i )  with  0 = s 0 < s 1 < < s n = S , where every subinterval  S i  corresponds to one of the drive modi, namely full traction, partial traction, coasting, partial regenerative brake, and full regenerative brake. The cases full traction, coasting, and full regenerative brake are regular in the sense that the control variables  u tr  and  u br  are both defined, and therefore the equation of motion is completely provided by the underlying physics. On the contrary, the modi partial traction and partial brake are singular, meaning that here only one control variable is defined. Therefore, they need some additional consideration. We will distinguish the following four cases:
Case 1: Partial traction below speed limit. Let  S i  be an interval of partial traction in the interior of the state domain, i.e., an interval where  ψ ( s ) 1  and  v < v max  hold. It follows that  μ 3 ( s ) μ 4 ( s ) 0  on this interval; see Table 1. The costate Equation (50) reduces to
ψ T = ( 2 K ) 3 / 2 d w d K = 2 3 / 2 R ( K )
with  R ( K ) = K 3 / 2 d w d K  being a known and monotonically increasing function in K. Moreover, the function value  R ( K )  is strictly positive for  K > 0 . Since  ψ T  is constant, we have
K ( s ) R 1 2 3 / 2 ψ T = : K ptr ,
where  R 1  stands for the inverse function of R. This means that, for the optimal solution, partial traction in the interior of the state domain is only possible at a constant speed with  K ( s ) K ptr . From Equation (54) it follows that  ψ T < 0  must hold.
Case 2: Partial regenerative brake below speed limit. Likewise, if  S i  is an interval of partial regenerative brake in the interior of the state domain, then  ψ ( s ) α  and  μ 3 ( s ) μ 4 ( s ) 0  will hold, and the costate equation now reads
ψ T = α ( 2 K ) 3 / 2 d w d K = 2 3 / 2 α R ( K ) .
Since  ψ T  is constant,
K ( s ) R 1 2 3 / 2 ψ T / α = : K pbr
holds, meaning that partial regenerative brake in the interior of the state domain is only possible at a constant speed with  K ( s ) K pbr . (In Khmelnitsky [18], the constants  K ptr  and  K pbr  are called  K s  and  K b s , respectively.)
Case 3: Partial traction on the speed limit. Let  S i  be an open interval of partial traction on the speed limit. Then,  ψ ( s ) 1  for  s S i . Both  u tr  and  μ 3  are unknown, but the equation of motion is entirely determined by the speed limit with  v ( s ) = v max ( s )  and  K ( s ) = K max ( s ) . Since  v max  is piecewise constant, and v can not be discontinuous, v v max , K, and  K max  must be constant in the interval  S i .
Case 4: Partial regenerative brake on the speed limit. Likewise, if  S i  is an open interval of partial regenerative brake on the speed limit, then  ψ ( s ) α  for  s S i . Both  u br  and  μ 4  are unknown in this case, but the equation of motion is again entirely determined by the speed limit with  v ( s ) v max  and  K ( s ) K max . Both  v max  and  K max  must be constant on  S i , using the same argument as in Case 3.
The location of the sub-intervals  S i  corresponding to the drive modi full traction, partial traction, coasting, and partial regenerative brake, and full regenerative brake is not defined a priori. The intervals  S i  will be derived in the algorithm of Khmelnitsky, a process that is shown in the following Section 2.4 and Section 2.5. The process starts with the definition of intervals where constant speed could possibly be driven and proceeds with trials to link these intervals by regular motion trajectories. Finally, if a complete link from the departure point to the arrival point exists that consists of periods with regular motion and admissible constant speed sections, and if this link meets the prescribed total time, then the unique optimal solution is found.
Remark 7.
In recent literature (not in Khmelnitsky), the constant speed driving regimes are often called cruising.

2.4. The Algorithm of Khmelnitsky for Fixed  ψ T

In this section, we describe the algorithm of finding a solution to the optimisation problem under the assumption that the Lagrange multiplier  ψ T < 0  is a given number. Then,  K ptr  and  K pbr  are defined according to (54) and (56), respectively. Following Khmelnitsky [18], we define intervals with possible constant speed that correspond to the four cases studied in Section 2.3. In the case of constant speed,  d K / d s = 0  holds, and the state differential Equation (20) takes the form
u tr u br = d P d s + w .
Definition 1.
A sub-interval I of  [ 0 , S ]  with
0 d P ( s ) d s + w ( K ptr ) g tr ( K ptr ) and K ptr < K max ( s )
for all  s I  that cannot be extended, i.e., any enlargement of the interval would violate (58), is called a PT-interval. PT stands for partial traction.
Definition 2.
A sub-interval I of  [ 0 , S ]  with
g br ( K pbr ) d P ( s ) d s + w ( K pbr ) 0 and K pbr < K max ( s )
for all  s I  that cannot be extended is called a PB-interval. PB stands for partial regenerative brake.
Definition 3.
A sub-interval I of  [ 0 , S ]  with
0 d P ( s ) d s + w ( K max ( s ) ) g tr ( K max ( s ) ) and K max ( s ) = const .
for all  s I  that cannot be extended is called a PT-SL-interval. PT-SL stands for partial traction on speed limit.
Definition 4.
A sub-interval I of  [ 0 , S ]  with
g br ( K max ( s ) ) d P ( s ) d s + w ( K max ( s ) ) 0 and K max ( s ) = const .
for all  s I  that cannot be extended is called a PB-SL-interval. PB-SL stands for partial regenerative brake on speed limit.
Intervals of type PT, PB, PT-SL, and PB-SL are intervals where a constant speed motion of the optimum solution would be allowed by the control constraints. Those intervals are summarised under the name pcs-intervals, meaning ‘possible constant speed’. (In Khmelnitsky [18], PT is called minor grade, PB is called steep fall, PT-SL is called minor grade imitation, and PB-SL is called steep fall imitation.) Due to their definition, pcs-intervals will never intersect. In this paper, the start point at  s = 0 , the pcs-intervals, and the stop point at  s = S  are summarised under the name ports. Additional speed limit ports might be introduced, as will be explained later. All ports are numbered in the order of increasing s.
The optimal solution is found by connecting ports by trajectories of regular motion, namely full traction, coasting, or full regenerative brake. Below the speed limit, regular motion is governed by the differential equations
d K d s = d P d s + u w ,
d ψ d s = ψ d w d K + ψ T · ( 2 K ) 3 / 2 μ 3 d g tr d K μ 4 d g br d K .
with
u = g tr if ψ > 1 0 if α ψ 1 g br if ψ < α
and  μ 3  and  μ 4  as defined in (51) and (52), respectively. Equation (62) follows from (20) and Table 1, while (63) results from (50). Since u is discontinuous at  ψ = α  and  ψ = 1 , the trajectory of K will have a kink point whenever  ψ  crosses these values.
It follows from the maximum principle that if the optimal solution touches the speed limit at a position  s i , then  h ( s i ) = 0  and  η i  might be positive by Equation (45). This means for the first component of Equation (43),
ψ ( s i + ) ψ ( s i ) = η i h E = η i h K = η i 0 ,
i.e.,  ψ  might have a positive jump whenever the speed limit is attained. The height of this jump, however, is not defined in the maximum principle. It must be found by a one-dimensional search.
If one tries to connect the ports A and B by a trajectory of regular motion, and this trajectory violates the speed limit, then a connection from A to B is not possible by this trajectory. In this case, one looks for a connection to the speed limit itself. If a trajectory is found that touches the speed limit in a single point, a new port is inserted there. It is called a speed limit touchpoint (SL-TP). Ports are now renumbered to be in the order of increasing s again. If the new speed limit touchpoint is located between the begin and the end position of a pcs-interval, then numbering is conducted such that the speed limit touchpoint comes first. If the speed limit touchpoint has number i, the  ψ -value of the incoming trajectory is stored as  ψ i . At the speed limit touchpoint,  ψ  is allowed to have a positive jump. A connection of ports in order to construct the optimal solution is only allowed in the direction of increasing port numbers.
When a trajectory of regular motion leaves a port, it is called a take off, and when it arrives at a port, it is called a landing. When trying to connect a port A with a port B by a trajectory of regular motion, there is always one of the variables s, K, and  ψ  not fixed at both take off and landing. These variables can be adjusted to make the connection possible. The following cases of take off and landing can exist:
  • (T1): Take off from the start point at  s = 0  with  K = 0 . The value of  ψ  is not fixed, but must be greater than 1 since full traction is applied.
  • (T2): Take off from interval of type PT with  K = K ptr  and  ψ = 1 ± ϵ , with some small  ϵ > 0 . Since  d K / d s  is discontinuous at  ψ = 1 , the trajectories of K will leave in different directions depending on the choice of  ψ  slightly above or below 1, so both must be checked. The take-off position s is not fixed.
  • (T3): Take off from the interval of type PB with  K = K pbr  and  ψ = α ± ϵ , again with some small  ϵ > 0 . The take-off position s is not fixed.
  • (T4): Take off from interior of a PT-SL interval with  K = K max  and  ψ = 1 . The take-off position s is not fixed.
  • (T5): Take off from the end of a PT-SL interval with  K = K max  and  ψ 1 , since a jump in  ψ  is allowed here.
  • (T6): Take off from interior of a PB-SL interval with  K = K max  and  ψ = α . The take-off position s is not fixed.
  • (T7): Take off from the end of a PB-SL interval with  K = K max  and  ψ α , since a jump in  ψ  is allowed here.
  • (T8): Take off from an SL-TP with number i: start with  K = K max  and  ψ ψ i , since a jump in  ψ  is allowed here.
  • (L1): Landing on PT interval with  K = K ptr  and  ψ = 1 . The landing position s is not fixed.
  • (L2): Landing on PB interval with  K = K pbr  and  ψ = α . The landing position s is not fixed.
  • (L3): Landing on the start of a PT-SL interval with  K = K max  and  ψ < 1 . Then, a new speed limit touchpoint is inserted at the landing position, is connected with the PT-SL interval, and port renumbering is conducted such that the new speed limit touchpoint comes before the PT-SL interval. If the new speed limit touchpoint has number i, the  ψ -value of the incoming trajectory is stored as  ψ i .
  • (L4): Landing on the start of a PT-SL interval with  K = K max  and  ψ = 1 . Since any jump of  ψ  here would lead to full traction, it would violate the speed limit. Therefore, no new speed limit touchpoint is inserted.
  • (L5): Landing on the interior of a PT-SL interval with  K = K max  and  ψ = 1 . The landing position s is not fixed.
  • (L6): Landing on the start of a PB-SL interval with  K = K max  and  ψ < α . Then, a new speed limit touchpoint is inserted at the landing position, is connected with the PB-SL interval, and port renumbering is done such that the new speed limit touchpoint comes before the PB-SL interval. If the new speed limit touchpoint has number i, the  ψ -value of the incoming trajectory is stored as  ψ i .
  • (L7): Landing on the start of a PB-SL interval with  K = K max  and  ψ = α . Since any jump of  ψ  here would lead to coasting, it would violate the speed limit. Therefore, no new speed limit touchpoint is inserted.
  • (L8): Landing on the interior of a PB-SL interval with  K = K max  and  ψ = α . The landing position s is not fixed.
  • (L9): Landing on an SL-TP with  K = K max . The  ψ -value of the incoming trajectory is stored, and  ψ  is allowed to jump here.
  • (L10): Landing on the end point at  s = S  with  K = 0 . The value of  ψ  is not fixed.
The construction of an optimal solution is best illustrated using a numerical example.
Example 1.
A train is driven from Station A at  s = 0  to Station B at  s = S = 20   km . The altitude is given by  z = 40   m · sin ( s / km ) + 100   m ; see Figure 2. Parameters are set according to Table 2, no speed limit is assumed, and we set  ψ T = 1 .
Remark 8.
The sine-shaped altitude is only assumed here as an academic example. In fact, railway lines are typically built with a rather piecewise-linear altitude profile, since this minimises the slope. Two examples based on an existing railway line will be shown later in Section 3.2 (Example 5) and Section 3.3 (Example 6).
Remark 9.
As mentioned above,  ψ T  is only required to be negative, and any negative value of  ψ T  would result in an energy minimisation problem.  ψ T  is only set to  1  as an example; this choice is not mandatory here.
Equations (54) and (56) lead to  K ptr = 536   m 2 / s 2  and  K pbr = 665   m 2 / s 2 . Figure 3 shows the pcs-intervals calculated with Equations (58) to (61), together with the kinetic energy of the fastest motion. In our example, this leads to the exclusion of the first PT-interval since it cannot be reached.
The first step in connecting ports would be the take off from the start point at  s = 0 . This is the take-off case (T1). The trajectories of K and  ψ  are calculated according to Equations (62) and (63). Different values of  ψ  at  s = 0  lead to different trajectories that are shown in Figure 4 and Figure 5. When the  ψ -trajectory crosses the value 1, full traction changes to coasting, and the corresponding K-trajectory has a kink point. When the  ψ -trajectory crosses the value  α = η tr η br = 0.7225 , coasting changes to full regenerative brake, and the corresponding K-trajectory has a kink point again. The blue line marks the only trajectory that would land on the PB-interval with number 2 (landing case (L2)). It is found by an iterative bisection of the  ψ -values at  s = 0 .
Figure 6 and Figure 7 show the take off from the PT-interval number 3. This is the take-off case (T2). On the interval, the trajectories show an unstable behaviour with respect to  ψ . If  ψ  is slightly above the value of 1, both K and  ψ  will move upwards. If  ψ  is slightly below the value of 1, both K and  ψ  will move downwards. Note that  ψ  takes off tangentially on the entire interval while K does so only from the ends of the interval. Close to the right end of the interval, the upwards moving K- and  ψ -trajectories will soon turn downwards, and by that change the driving modus from full traction to coasting. This is a typical picture for the take off from PT- to PB-intervals.
Figure 8 and Figure 9 show the K- and  ψ -trajectories when trying to connect intervals 4 and 8. The K-trajectories have a kink point when the driving modus changes from coasting to full traction or full regenerative brake, corresponding to  ψ  crossing the values 1 or  α . Both the K- and  ψ -trajectories are able to cross pcs-intervals, but the trajectory field often splits at pcs-intervals, as here at interval 7, where a shadowed region lies behind that cannot be reached by the trajectories.
Remark 10.
There is great potential for speeding up the algorithm by the detection of those shadowed regions. For example, one can conclude immediately from Figure 4 that interval 3 cannot be reached from starting point 1, and from Figure 8 that interval 9 cannot be reached from interval 4. When solving the optimisation problem, it will always pay off to invest in good statistics, showing which connections should be checked and which ones can safely be excluded.
Example 2.
We consider Example 1, but now with a speed limit according to Table 3.
The multiplier  ψ T  is again set to  1 . Figure 10, Figure 11 and Figure 12 show that all types of pcs-intervals occur. The speed limit touchpoints 4, 10, and 14 have been inserted during the run of the algorithm.
Table 4 shows the possible port connections in Example 2.
It was shown by Khmelnitsky [18] that there always exists a unique connection from the start point (1) to the stop point (here, 17). In Example 2, this is the connection  1 2 4 5 9 10 11 13 14 15 17 . The total time  T ˜  required on this connection is not known a priori, but can be calculated from the K-curve. Since  K = v 2 / 2  is the specific kinetic energy,
T ˜ = 0 S d s 2 K
holds.
Remark 11.
Here, we distinguish between the scheduled time T and the time  T ˜  that is evaluated from the algorithm. The final goal of the algorithm is to match  T ˜  to the prescribed T by a variation of  ψ T . This will be explained in Section 2.5.
Example 3.
We consider Example 2, but now with  ψ T = 0.3 . Trajectories of K and ψ are illustrated in Figure 13 and Figure 14, respectively.
The connections shown in Figure 13 and Figure 14 are listed in Table 5.

2.5. The Algorithm of Khmelnitsky: Variation of the Multiplier  ψ T

Khmelnitsky [18] has shown that the time  T ˜  between two stops of the train is a strictly monotonously increasing function of the Lagrange multiplier  ψ T . Moreover,  ψ T < 0  and
lim ψ T 0 T ˜ =
hold. Therefore, whenever the scheduled time T is possible to be driven on a track section by a particular train, it can be approached by the optimisation algorithm by adjusting  ψ T  with iterative bisection. The complete algorithm for the minimum energy operation would read as follows:
Consider a track section with stops at  s = 0  and  s = S  to be driven in a scheduled duration T.
  • Step 1: Calculate the fastest possible motion on the track section with purely regenerative braking. If the time  T ˜  needed for that is larger than the scheduled time T, add mechanical braking and leave the algorithm. If  T ˜ < T , choose an arbitrary negative value for  ψ T  and proceed with Step 2 in the algorithm.
  • Step 2: Calculate  K ptr  by (54),  K pbr  by (56), and evaluate the pcs-intervals PT, PB, PT-SL, PB-SL by (58)–(61). Skip pcs-intervals that cannot be reached by the fastest motion. Numerate the ports, i.e., the remaining pcs-intervals, the start at  s = 0 , and the stop at  s = S , in the order of ascending s.
  • Step 3: Try to connect ports by regular motion with Equations (62) and (63). Use one of the take-off cases (T1)–(T8) and one of the landing cases (L1)–(L10). Add speed limit touchpoints if necessary, according to the instructions given above. Step 3 is complete when a connection from the start point at  s = 0  to the stop point  s = S  has been found.
  • Step 4: Calculate  T ˜  according to (66). If  T ˜  is sufficiently close to T, the algorithm is successfully completed. If not, adjust  ψ T  and proceed with Step 2.
If n is the number of ports then a maximum of  n 2 = ( n 2 n ) / 2  possible connections has to be checked, unless the start–stop connection is found earlier. This means that the number of possible connections grows quadratically with n. The algorithm can be considered as a search tree, with port connections being the branches of the tree that need to be checked. Therefore, it is crucial to follow a clever search strategy to keep computing time at an acceptable level. We mention four important measures that dramatically shortened computing time when the code was developed:
Parallel path exclusion. Khmelnitsky [18] has proven that any two ports can only be connected by at most one path. Therefore, it is wise to exclude all parallel paths in the search tree. For example, if a connection  1 2 3  is established, the parallel direct link  1 3  is not possible and need not be checked.
Start from the treetop. When choosing the next connection to check for in Step 3, start at the highest port number that is connected to port 1 and not a known dead end. This strategy usually leads to an early discovery of the start–stop connection.
Look out for shadowed ports. If ports lie in shadow, exclude impossible connections. See the remark in the discussion of Example 1 in Section 2.3.
Estimate  ψ T  by interpolation. The adjustment of  ψ T  in Step 4 can be sped up using interpolation techniques.
Example 4.
Example 4 is equal to Examples 2 and 3, except that  ψ T  is not prescribed, and T is set to 16 min.
Applying the algorithm,  ψ T  converges to  0.5237 . The final K- and  ψ -trajectories are shown in Figure 15 and Figure 16.
Connections and types of take off and landing are given for the final value  ψ T = 0.5237  in Table 6.
In Figure 17 and Figure 18, speed and electric energy are shown for Example 4 at  ψ T = 0.5237 .

3. Results and Discussion

3.1. Estimating the Potential of Efficient Train Driving

In order to compare the optimal driving performance to the energy demand of trains in a real-life operation, a study was carried out covering a total of 100 regional train runs on four track sections in Germany. The trains were drawn by electric locomotives equipped with regenerative braking. Time, velocity, energy supply from the catenary, and energy return by regenerative braking have been recorded. Trains with five, six, or seven coaches were included in the study. Train runs affected by speed restriction due to signalling in the interior of the track section have been excluded.
The total energy demand is defined by  E tot = E tr E br + E add , where  E tr  is the electric energy used for traction,  E br  is the electric energy returned by regenerative braking, and  E add  is an additional energy demand that is used, e.g., for air-conditioning, pressurised air, ventilation, etc. Recorded data are the energy supply  E tr + E add  and the energy return  E br . Both traction energy  E tr  and brake energy  E br  can also be computed from the measured velocity using the model given in Section 2.1, and it is interesting to compare the measured energy quantities with the ones calculated from the model. However, a direct comparison of energy return with the regenerative brake energy was difficult since the energy return turned out to be generally far less than the potential of the brake energy calculated from the model. This indicates that in many cases the regenerative brake was little used by the train drivers. A comparison of energy supply, on the other hand, had to deal with the problem that the additional energy  E add  was not recorded separately from traction energy  E tr . Therefore, we have fitted the unknown additional energy such that the measured energy supply  E tr + E add  comes close to the energy supply derived from the model. This leads to the estimate of  E add = 254   kW . The comparison of energy supply is shown in Figure 19.
Result. The energy supply was fitted to lie close to the equality line. The graph still shows that the model and measurement do not perfectly agree. This could be due to the following reasons:
  • The additional energy demand  E add  might differ from one train to another.
  • Wind conditions are not accounted for.
  • The air drag modelled by the ‘Hanover formula’ (9), (10) is only a rough approximation.
  • The engine efficiency is simplified in the model to be a constant; speed and force dependent efficiency is not considered.
In the following study, we will compare the measured energy demand to the energy optimum derived by the algorithm. In this study, it would not be appropriate to compare the energy demand of a particular train run to the energy optimum based on the duration given by the timetable. Trains that are delayed need to drive faster; thus, the shorter time they have available needs to be considered in the optimum calculation. On the other hand, trains that arrive early are judged to unnecessarily waste energy by driving too fast. Therefore, for each train run, a reference time is considered that allows a fair comparison to the optimum. Let
  • t start TT  be the departure time according to the timetable;
  • t stop TT  be the arrival time according to the timetable;
  • t start rec  be the recorded departure time, and
  • t stop rec  be the recorded arrival time of the train run.
Then, the reference time duration is defined by  T = max ( t stop TT , t stop rec ) t start rec , and the energy optimum is calculated with respect to this reference time. Figure 20, Figure 21, Figure 22 and Figure 23 display the total energy  E tot  over the reference time T. The stars in the figures are recorded measurements, while the lines indicate the calculated energy optimum.
Result. Figure 20, Figure 21, Figure 22 and Figure 23 show that the recorded total energy spreads over a rather wide range. A factor of approximately two lies between the minimum and maximum energy consumption on each track section, even for cases with a quite similar reference time. Contrary to the expectation, the recorded values do not show an increase in energy demand with train length, but this might be due to the small size of the sample and the wide spread of values. For many train runs, the measured energy is far above the corresponding optimum line, in some extremal cases by a factor of three.
Remark 12.
In one case shown in Figure 20, a seven-coach train was in fact better than the optimum. This is not a contradiction to the optimality property. As was pointed out, the optimal control theory ensures that the algorithm returns the energy minimum, but this only holds for the given train motion model. This model is only a simplification, and some assumptions might not have been met during the train drive experiment. For example, wind speed is neither measured nor considered in the model, and the additional energy requirement  E add  is only averaged for a lack of individual measurements.
If, for the i-th train run, the recorded total energy is denoted by  E tot , i rec , and the corresponding optimum energy is  E tot , i opt , then the energy ratio  η i = E tot , i opt / E tot , i rec  can be evaluated. Within the concept of the Physical Optimum (PhO) [55,56],  η i  is called a PhO factor that quantitatively evaluates the efficiency with respect to a feasible and most efficient reference process. This energy ratio is shown in Figure 24 for each of the 100 train runs of the study. A total energy ratio  η tot = i E tot , i opt / i E tot , i rec = 0.627  is obtained by summing over all train runs studied. This can be seen as an estimate of the potential that energy optimised driving would have. Based on the results of the present study, an amount of 37.3 percent of energy could be saved by following an energy minimising driving strategy. One should, however, bear in mind two limitations: First, due to model simplifications, the computed driving strategy might differ from the real optimum under the current conditions. Secondly, even if an optimum strategy is displayed to the driver by an assistance system, it can only be approached, and interactions with other traffic will sometimes not allow one to exactly follow the suggestions. However, since the theoretical optimum turned out to be substantially lower than the measured energy demand in operation, we think that there is still a large potential for energy saving by the optimal control-based assistance for train drivers.

3.2. Sensitivity to Deviations from the Optimum Driving Profile

In this section, it is studied how deviations from the optimum driving profile will increase the energy requirement. The effect is shown using the following example.
Example 5.
We consider the railway section between Dresden-Neustadt and Riesa. Both stations are located on the Dresden–Leipzig line in Germany. The altitude of the section is shown in Figure 25. Train parameters are the same as in Examples 1–4, see Table 2. The only difference to Examples 1–4 is that now an additional energy requirement  E add = 250   kW  is assumed. The total time T is set to 22 min. Speed limits are not assumed in this example.
After calculation of the optimum driving profile, seven neighbouring trajectories are constructed, called variant 1 to variant 7. Every variant is driven with the same train parameters on the same track, and the total time needed is equal for every variant. Variants 1–3 are too fast at the beginning. Then, at a certain position, the speed of each variant is corrected such that it will still arrive at the scheduled total time of 22 min. On the other hand, variants 4–7 are too slow at the beginning. Variant 4 manages to arrive at the prescribed time by completely avoiding the coasting regime. Variants 5–7 accelerate at certain positions and then switch to coasting. In order not to overload the figures, we first show the comparison of variant 1–3 to the optimum in Figure 26, Figure 27 and Figure 28 and proceed then with variant 4–7 versus the optimum in Figure 29, Figure 30 and Figure 31.
Result. The total electric energy demand  E tot  is compared
As could be expected, the optimum driving strategy indeed leads to a minimum value of  E tot  at arrival. However, most variants studied here require only slightly more total electric energy. An exception is variant 4 that has no coasting section included and needs substantially more energy. The values of  E tot  and percentages are given in Table 7.
Remark 13.
In fact it would be easy to construct driving strategies that are far worse than the ones considered here, for example driving with frequent switching between acceleration and braking, but this was not the focus of our study. It was our intention to compare the optimum to some variants that are still quite reasonable. Most of the variants perform quite well, since from the end point of their cruising period on, they are constructed such that they actually continue in an optimum way. The ‘good news’ from this study is that, whenever a deviation from optimum is noticed early enough, counteraction by following a recalculated optimum profile will keep the error small.

3.3. Influence of the Total Recuperation Efficiency  α  on the Optimum Driving Profile

As the total recuperation efficiency  α = η tr η br  is directly included in the maximisation problem, Equation (19), it can be expected that a variation of  α  leads to a certain change in the optimum driving strategy. In this section, the effect of  α  on the optimum speed profile is briefly illustrated in a numerical example.
Example 6.
Example 6 is equal to Example 5 (Dresden-Neustadt—Riesa) except for the following changes: (1) The total time T is set to 23 min. (2) A speed limit is considered according to Table 8. (3) The efficiency  η br  of the regenerative brake will be varied.
The optimum train driving profiles are calculated for three values of  η br . The resulting curves for speed v and total electric energy  E tot  are shown in Figure 32 and Figure 33, respectively.
Result. A small effect of different recuperation efficiency  α  on speed is visible during the coasting periods. The position where coasting starts is slightly different. The effect on total energy is much more pronounced due to the different amount of recuperation energy.
Example 6 is only meant to show that a certain (but small) effect on speed could be observed. Braking regimes lead to different energy losses when the recuperation efficiency changes. This can result in different optimum speed profiles. A deeper analysis of this subject is beyond the focus of this study.

4. Conclusions

Within this paper, energy-efficient train driving has been studied for the case of an exclusive usage of the regenerative brake in electric trains. The code Leda written at Fraunhofer Institute Magdeburg is based on the algorithm of Khmelnitsky that constructs the unique minimum energy solution. A derivation of the statements given by Khmelnitsky from a more general formulation of Pontryagin’s maximum principle is presented. In addition to the theory in Khmelnitsky’s article, a complete list of switching cases was provided and illustrated by a number of numerical examples. A comparison to energy consumption data in a real operation showed that the energy minimising strategy was able to save, on average, about 37% of energy. A study of some sub-optimum driving profiles showed, in most cases, only moderate energy loss when deviations are corrected using a recalculated optimum. A slight influence of regenerative brake efficiency on the optimum speed profile was detected in a numerical example.
Extensions of the code to include mechanical braking, non-zero speed boundary conditions, and dynamic response to train delays are subject to further research.

Author Contributions

Writing—original draft preparation, W.H.; software, W.H.; project administration, M.R. and T.B.-R.; funding acquisition, M.R. and T.B.-R. All authors have read and agreed to the published version of the manuscript.

Funding

The authors gratefully acknowledge external funding by the Federal Office for Economic Affairs and Export Control (Bundesamt für Wirtschaft und Ausfuhrkontrolle) of Germany within the pilot programme ‘Einsparzähler’.

Data Availability Statement

External data that have been processed cannot be provided nor referenced for confidentiality reasons.

Acknowledgments

The authors would like to thank the team of traveltainer GmbH & Co. KG, Altenbeken, Germany for their valuable collaboration and support during this research project.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
DDPDiscrete dynamic programming
EETCEnergy-efficient train control
KKTKarush–Kuhn–Tucker (conditions)
NLPNonlinear programming
PBPartial regenerative braking
PB-SLPartial regenerative braking at speed limit
pcs-intervalPossible constant speed interval
PMPPontryagin maximum principle
PTPartial traction
PT-SLPartial traction at speed limit
SL-TPSpeed limit touchpoint
SNCFSociété nationale des chemins de fer français (French national railways)
TGVTrain à grande vitesse (French high-speed train)

References

  1. Strobel, H.; Horn, P.; Kosemund, M. Contribution to optimum computer-aided control of train operation. In Proceedings of the 2nd IFAC/IFIP/IFORS Symposium on Traffic Control and Transportation Systems, Monte Carlo, Monaco, 16–21 September 1974; pp. 377–387. [Google Scholar]
  2. Howlett, P.G.; Pudney, P.J. Energy-Efficient Train Control; Springer: London, UK, 1995. [Google Scholar]
  3. Franke, R.; Terwiesch, P.; Meyer, M. An algorithm for the optimal control of the driving of trains. In Proceedings of the 39th IEEE Conference on Decision and Control, Sydney, Australia, 12–15 December 2000; Volume 3. [Google Scholar]
  4. Albrecht, T. Ein Beitrag zur Nutzbarmachung Genetischer Algorithmen für die Optimale Steuerung und Planung eines Flexiblen Stadtschnellbahnbetriebes (A Contribution to the Utilization of Genetic Algorithms for Optimal Control and Planning of a Flexible Urban Rail Transport). Ph.D. Thesis, Dresden University of Technology, Dresden, Germany, 2005. [Google Scholar]
  5. Howlett, P.G.; Pudney, P.J.; Vu, X. Freightmiser: An energy-efficient application of the train control problem. In Proceedings of the 30th Conference of Australian Institutes of Transport Research (CAITR), Perth, Australia, 10–12 December 2008. [Google Scholar]
  6. Aradi, S.; Becsi, T.; Gaspar, P. A predictive optimization method for energy-optimal speed profile generation for trains. In Proceedings of the 2013 IEEE 14th International Symposium on Computational Intelligence and Informatics (CINTI), Budapest, Hungary, 19–21 November 2013. [Google Scholar]
  7. DAS on-Board Algorithms. Technical Report D6.3, ON-TIME: Optimal Networks for Train Integration Management across Europe. 2014. Available online: https://cordis.europa.eu/project/id/285243 (accessed on 9 September 2023).
  8. Albrecht, A.R.; Howlett, P.G.; Pudney, P.J.; Vu, X.; Zhou, P. Energy-efficient train control: The two-train separation problem on level track. J. Rail Transp. Plan. Manag. 2015, 5, 163–182. [Google Scholar] [CrossRef]
  9. Scheepmaker, G.; Goverde, R.; Kroon, L. Review of energy-efficient train control and timetabling. Eur. J. Oper. Res. 2017, 257, 355–376. [Google Scholar] [CrossRef]
  10. Pontryagin, L.; Boltyanskii, V.; Gamkrelidze, R.; Mishchenko, E. Matematicheskaya Teoriya Optimal’nykh Prozessov; Fizmatgiz: Moscow, Russia, 1961. [Google Scholar]
  11. Pontryagin, L.; Boltyanskii, V.; Gamkrelidze, R.; Mishchenko, E. The Mathematical Theory of Optimal Processes; John Wiley and Sons (Interscience Publishers): New York, NY, USA, 1962. [Google Scholar]
  12. Ichikawa, K. Application of optimization theory for bounded state variable problems to the operation of train. Bull. Jpn. Soc. Mech. Eng. 1968, 11, 857–865. [Google Scholar] [CrossRef]
  13. Milroy, I.P. Aspects of Automatic Train Control. Ph.D. Thesis, Loughborough University, Loughborough, UK, 1980. [Google Scholar]
  14. Howlett, P.G.; Milroy, I.P.; Pudney, P.J. Energy-efficient train control. Control Eng. Pract. 1994, 2, 193–200. [Google Scholar] [CrossRef]
  15. Liu, R.R.; Golovitcher, I.M. Energy-efficient operation of rail vehicles. Transp. Res. Part A Policy Pract. 2003, 37, 917–932. [Google Scholar] [CrossRef]
  16. Albrecht, A.R.; Howlett, P.G.; Pudney, P.J.; Vu, X. Energy-efficient train control: From local convexity to global optimization and uniqueness. Automatica 2003, 49, 3072–3078. [Google Scholar] [CrossRef]
  17. Asnis, I.A.; Dmitruk, A.V.; Osmolovskii, N.P. Solution of the problem of the energetically optimal control of the motion of a train by the maximum principle. USSR Comput. Math. Math. Phys. 1985, 25, 37–44. [Google Scholar] [CrossRef]
  18. Khmelnitsky, E. On an Optimal Control Problem of Train Operation. IEEE Trans. Autom. Control 2000, 45, 1257–1266. [Google Scholar] [CrossRef]
  19. Ying, P.; Zeng, X.; Song, H.; Shen, T.; Yuan, T. Energy-efficient train operation with steep track and speed limits: A novel Pontryagin’s maximum principle-based approach for adjoint variable discontinuity cases. IET Intell. Transp. Syst. 2021, 15, 1183–1202. [Google Scholar] [CrossRef]
  20. Ying, P.; Zeng, X.; Shen, T.; Wang, Y.; Ma, Z.; Wu, Y. Partial Train Speed Trajectory Optimization. In Proceedings of the 2022 IEEE 25th International Conference on Intelligent Transportation Systems (ITSC), Macau, China, 8–12 October 2022. [Google Scholar]
  21. Baranov, L.A.; Meleshin, I.S.; Chin’, L.M. Optimal control of a subway train with regard to the criteria of minimum energy consumption. Russ. Electr. Eng. 2011, 82, 405–410. [Google Scholar] [CrossRef]
  22. Lu, S.; Weston, P.; Hillmansen, S.; Gooi, H.; Roberts, C. Increasing the regenerative braking energy for railway vehicles. IEEE Trans. Intell. Transp. Syst. 2014, 15, 2506–2515. [Google Scholar] [CrossRef]
  23. Zhou, Y.; Bai, Y.; Li, J.; Mao, B.; Li, T. Integrated optimization on train control and timetable to minimize net energy consumption of metro lines. J. Adv. Transp. 2018, 2018, 7905820. [Google Scholar] [CrossRef]
  24. Fernández-Rodríguez, A.; Fernández-Cardador, A.; Cucala, A. Real time eco-driving of high speed trains by simulation-based dynamic multi-objective optimization. Simul. Model. Pract. Theory 2018, 84, 50–68. [Google Scholar] [CrossRef]
  25. Scheepmaker, G.; Goverde, R. Energy-efficient train control using nonlinear bounded regenerative braking. Transp. Res. Part C Emerg. Technol. 2020, 121, 102852. [Google Scholar] [CrossRef]
  26. Rao, A.; Benson, D.; Darby, C.; Patterson, M.; Francolin, C.; Sanders, I.; Huntington, G. Algorithm 902: GPOPS, a MATLAB software for solving multiple-phase optimal control problems using the Gauss pseudospectral method. ACM Trans. Math. Softw. 2010, 37, 22:1–22:39. [Google Scholar] [CrossRef]
  27. Su, S.; Tang, T.; Li, X. Driving strategy optimization for trains in subway systems. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit 2018, 232, 369–383. [Google Scholar] [CrossRef]
  28. Zhu, Q.; Su, S.; Tang, T.; Liu, W.; Zhang, Z.; Tian, Q. An eco-driving algorithm for trains through distributing energy: A Q-Learning approach. ISA Trans. 2022, 122, 24–37. [Google Scholar] [CrossRef]
  29. Su, S.; Zhu, Q.; Liu, J.; Tang, T.; Wei, Q.; Cao, Y. A Data-Driven Iterative Learning Approach for Optimizing the Train Control Strategy. IEEE Trans. Ind. Inform. 2023, 19, 7885–7893. [Google Scholar] [CrossRef]
  30. Ghaviha, N.; Bohlin, M.; Holmberg, C.; Dahlquist, E.; Skoglund, R.; Jonasson, D. A driver advisory system with dynamic losses for passenger electric multiple units. Transp. Res. Part C Emerg. Technol. 2017, 85, 111–130. [Google Scholar] [CrossRef]
  31. Kouzoupis, D.; Pendharkar, I.; Frey, J.; Diehl, M.; Corman, F. Direct multiple shooting for computationally efficient train trajectory optimization. Transp. Res. Part C Emerg. Technol. 2023, 152, 104170. [Google Scholar] [CrossRef]
  32. Feng, M.; Huang, Y.; Lu, S. Eco-driving Strategy Optimization for High-speed Railways Considering Dynamic Traction System Efficiency. IEEE Trans. Transp. Electrif. 2023. early access. [Google Scholar] [CrossRef]
  33. Cao, Y.; Zhang, Z.; Cheng, F.; Su, S. Trajectory Optimization for High-Speed Trains via a Mixed Integer Linear Programming Approach. IEEE Trans. Intell. Transp. Syst. 2022, 23, 17666–17676. [Google Scholar] [CrossRef]
  34. Zhang, Z.; Cao, Y.; Su, S. Energy-Efficient Driving Strategy for High-Speed Trains with Considering the Checkpoints. Chin. J. Electron. 2023, in press. [Google Scholar]
  35. Szkopiński, J.; Kochan, A. Energy Efficiency and Smooth Running of a Train on the Route While Approaching Another Train. Energies 2021, 14, 7593. [Google Scholar] [CrossRef]
  36. Su, S.; Tang, T.; Chen, L.; Liu, B. Energy-efficient train control in urban rail transit systems. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit 2015, 229, 446–454. [Google Scholar] [CrossRef]
  37. Fernández-Rodríguez, A.; Su, S.; Fernández-Cardador, A.; Cucala, A.P.; Cao, Y. A multi-objective algorithm for train driving energy reduction with multiple time targets. Eng. Optim. 2021, 53, 719–734. [Google Scholar] [CrossRef]
  38. Wang, Y.; Schutter, B.; van den Boom, T.; Ning, B. Optimal trajectory planning for trains—A pseudospectral method and a mixed integer linear programming approach. Transp. Res. Part C Emerg. Technol. 2013, 29, 97–114. [Google Scholar] [CrossRef]
  39. DIDO Optimal Control Software. Elissar Global. Available online: https://www.mathworks.com/products/connections/product_detail/dido.html (accessed on 9 September 2023).
  40. Ross, I.; Karpenko, M. A review of pseudospectral optimal control: From theory to flight. Annu. Rev. Control 2012, 36, 182–197. [Google Scholar] [CrossRef]
  41. Ross, I. A Primer on Pontryagin’s Principle in Optimal Control; Collegiate Publishers: San Francisco, CA, USA, 2015. [Google Scholar]
  42. PSOPT Optimal Control Software. Available online: https://www.psopt.net/ (accessed on 9 September 2023).
  43. Becerra, V. Solving complex optimal control problems at no cost with psopt. In Proceedings of the 2010 IEEE International Symposium on Computer-Aided Control System Design, Yokohama, Japan, 8–10 September 2010; pp. 1391–1396. [Google Scholar]
  44. GPOPS-II: Next-Generation Optimal Control Software. Available online: http://gpops2.com/ (accessed on 9 September 2023).
  45. Hartl, R.; Sethi, S.; Vickson, R. A survey of the maximum principles for optimal control problems with state constraints. SIAM Rev. 1995, 37, 181–218. [Google Scholar] [CrossRef]
  46. Ihme, J. Schienenfahrzeugtechnik; Springer Fachmedien: Wiesbaden, Germany, 2016. [Google Scholar]
  47. Fassbinder, S. Wie Energie-Effizient ist der Bahnverkehr Wirklich? Technical Report, Deutsches Kupferinstitut, Düsseldorf. 2010. Available online: https://de.slideshare.net/LeonardoENERGYDeutschland/wie-energieeffizient-ist-der-bahnverkehr-wirklich (accessed on 9 September 2023).
  48. Curtius, E.; Kniffler, A. Neue Erkenntnisse über Haftung zwischen Treibrad und Schiene. Elektrische Bahnen 1944, 20, 25–30. [Google Scholar]
  49. Schlecht, B. Maschinenelemente 2; Pearson Studium: Munich, Germany, 2010. [Google Scholar]
  50. Wie “grün” ist der Schienenverkehr? Technical Report, Deutsches Kupferinstitut, Düsseldorf. 2020. Available online: https://www.kupferinstitut.de/wp-content/uploads/2020/02/BahnEffizienz.pdf (accessed on 9 September 2023).
  51. Frilli, A.; Meli, E.; Nocciolini, D.; Pugi, L.; Rindi, A. Energetic optimization of regenerative braking for high speed railway systems. Energy Convers. Manag. 2016, 129, 200–215. [Google Scholar] [CrossRef]
  52. Voß, G.; Gackenholz, L.; Wiebels, R. Eine neue Formel (Hannoversche Formel) zur Bestimmung des Luftwiderstandes spurgebundener Fahrzeuge. Glas. Annalen–Zeitschrift Für Eisenbahnwesen Und Verkehrstechnik 1972, 96, 166–171. [Google Scholar]
  53. Karush, W. Minima of Functions of Several Variables with Inequalities as Side Constraints. Master’s Thesis, University of Chicago, Chicago, IL, USA, 1939. [Google Scholar]
  54. Kuhn, H.; Tucker, A. Nonlinear Programming. In Proceedings of the 2nd Berkeley Symposium on Mathematical Statististics and Probability, Berkeley, CA, USA, 31 July–12 August 1950; pp. 481–492. [Google Scholar]
  55. Volta, D. Das Physikalische Optimum als Basis von Systematiken zur Steigerung der Energie-und Stoffeffizienz von Produktionsprozessen. Ph.D. Thesis, Clausthal University of Technology, Clausthal-Zellerfeld, Germany, 2014. [Google Scholar]
  56. Eggers, N.; Böttger, J.; Kerpen, L.; Sankol, B.; Birth, T. Refining VDI guideline 4663 to evaluate the efficiency of a power-to-gas process by employing limit-oriented indicators. Energy Effic. 2021, 14, 73. [Google Scholar] [CrossRef]
Figure 1. Maximum force for traction and regenerative brake for an electric locomotive with mass  m loc = 84   t  and  P tr , mech = P br , mech = 5.6   MW .
Figure 1. Maximum force for traction and regenerative brake for an electric locomotive with mass  m loc = 84   t  and  P tr , mech = P br , mech = 5.6   MW .
Energies 16 06712 g001
Figure 2. Track altitude in Example 1.
Figure 2. Track altitude in Example 1.
Energies 16 06712 g002
Figure 3. Example 1: pcs-intervals PT and PB, fastest motion, and numbering.
Figure 3. Example 1: pcs-intervals PT and PB, fastest motion, and numbering.
Energies 16 06712 g003
Figure 4. Example 1: K-trajectories starting from  s = 0  (port 1) for various values of  ψ  at  s = 0 . Trajectories change from full traction to coasting and then to full regenerative braking. The blue trajectory connects port 1 with port 2.
Figure 4. Example 1: K-trajectories starting from  s = 0  (port 1) for various values of  ψ  at  s = 0 . Trajectories change from full traction to coasting and then to full regenerative braking. The blue trajectory connects port 1 with port 2.
Energies 16 06712 g004
Figure 5. Example 1:  ψ -trajectories starting from  s = 0  (port 1) for various values of  ψ  at  s = 0 . Trajectories change from full traction to coasting and then to full regenerative braking. The blue trajectory connects port 1 with port 2.
Figure 5. Example 1:  ψ -trajectories starting from  s = 0  (port 1) for various values of  ψ  at  s = 0 . Trajectories change from full traction to coasting and then to full regenerative braking. The blue trajectory connects port 1 with port 2.
Energies 16 06712 g005
Figure 6. Example 1: K-trajectories starting from PT-interval with number 3.
Figure 6. Example 1: K-trajectories starting from PT-interval with number 3.
Energies 16 06712 g006
Figure 7. Example 1:  ψ -trajectories starting from PT-interval with number 3.
Figure 7. Example 1:  ψ -trajectories starting from PT-interval with number 3.
Energies 16 06712 g007
Figure 8. Example 1: K-trajectories starting from interval 4, heading for interval 8.
Figure 8. Example 1: K-trajectories starting from interval 4, heading for interval 8.
Energies 16 06712 g008
Figure 9. Example 1:  ψ -trajectories starting from interval 4, heading for interval 8.
Figure 9. Example 1:  ψ -trajectories starting from interval 4, heading for interval 8.
Energies 16 06712 g009
Figure 10. Example 2: K-trajectories.
Figure 10. Example 2: K-trajectories.
Energies 16 06712 g010
Figure 11. Example 2:  ψ -trajectories.
Figure 11. Example 2:  ψ -trajectories.
Energies 16 06712 g011
Figure 12. Example 2: Zoom into  ψ -trajectories.
Figure 12. Example 2: Zoom into  ψ -trajectories.
Energies 16 06712 g012
Figure 13. Example 3: K-trajectories.
Figure 13. Example 3: K-trajectories.
Energies 16 06712 g013
Figure 14. Example 3:  ψ -trajectories.
Figure 14. Example 3:  ψ -trajectories.
Energies 16 06712 g014
Figure 15. Example 4: K-trajectories.
Figure 15. Example 4: K-trajectories.
Energies 16 06712 g015
Figure 16. Example 4:  ψ -trajectories.
Figure 16. Example 4:  ψ -trajectories.
Energies 16 06712 g016
Figure 17. Example 4: train speed.
Figure 17. Example 4: train speed.
Energies 16 06712 g017
Figure 18. Example 4: net electric energy.
Figure 18. Example 4: net electric energy.
Energies 16 06712 g018
Figure 19. Comparison of energy supply  E tr + E add ; measurement vs. model.
Figure 19. Comparison of energy supply  E tr + E add ; measurement vs. model.
Energies 16 06712 g019
Figure 20. Total energy demand  E tot  over reference time T for track in Section 1.
Figure 20. Total energy demand  E tot  over reference time T for track in Section 1.
Energies 16 06712 g020
Figure 21. Total energy demand  E tot  over reference time T for track Section 2.
Figure 21. Total energy demand  E tot  over reference time T for track Section 2.
Energies 16 06712 g021
Figure 22. Total energy demand  E tot  over reference time T for track in Section 3.
Figure 22. Total energy demand  E tot  over reference time T for track in Section 3.
Energies 16 06712 g022
Figure 23. Total energy demand  E tot  over reference time T for track in Section 4.
Figure 23. Total energy demand  E tot  over reference time T for track in Section 4.
Energies 16 06712 g023
Figure 24. Energy ratio  η i  for each run, shown in descending order.
Figure 24. Energy ratio  η i  for each run, shown in descending order.
Energies 16 06712 g024
Figure 25. Example 5: track altitude z.
Figure 25. Example 5: track altitude z.
Energies 16 06712 g025
Figure 26. Example 5: train speed v (optimum and variants 1–3).
Figure 26. Example 5: train speed v (optimum and variants 1–3).
Energies 16 06712 g026
Figure 27. Example 5: total electric energy  E tot  (optimum and variants 1–3).
Figure 27. Example 5: total electric energy  E tot  (optimum and variants 1–3).
Energies 16 06712 g027
Figure 28. Example 5: total electric energy  E tot  (optimum and variants 1–3). Zoom into Figure 27.
Figure 28. Example 5: total electric energy  E tot  (optimum and variants 1–3). Zoom into Figure 27.
Energies 16 06712 g028
Figure 29. Example 5: train speed v (optimum and variants 4–7).
Figure 29. Example 5: train speed v (optimum and variants 4–7).
Energies 16 06712 g029
Figure 30. Example 5: total electric energy  E tot  (optimum and variants 4–7).
Figure 30. Example 5: total electric energy  E tot  (optimum and variants 4–7).
Energies 16 06712 g030
Figure 31. Example 5: total electric energy  E tot  (optimum and variants 4–7). Zoom into Figure 30.
Figure 31. Example 5: total electric energy  E tot  (optimum and variants 4–7). Zoom into Figure 30.
Energies 16 06712 g031
Figure 32. Example 6: train speed v for different regenerative brake efficiency  η br .
Figure 32. Example 6: train speed v for different regenerative brake efficiency  η br .
Energies 16 06712 g032
Figure 33. Example 6: total electric energy  E tot  for different regenerative brake efficiency  η br .
Figure 33. Example 6: total electric energy  E tot  for different regenerative brake efficiency  η br .
Energies 16 06712 g033
Table 1. Case distinction for  ψ . ‘cont.’ means that the continuity of the Lagrange multipliers was exploited. The corresponding value is only valid in the interior of the state domain, i.e., when  v < v max .
Table 1. Case distinction for  ψ . ‘cont.’ means that the continuity of the Lagrange multipliers was exploited. The corresponding value is only valid in the interior of the state domain, i.e., when  v < v max .
Case   u tr =   u br =   μ 1 =   μ 2 =   μ 3 =   μ 4 =
  ψ > 1   g tr 00   ψ α   ψ 1 0
(full traction)(38)(38)(41)(49)(48)(41)
  ψ = 1 00   1 α 00
(partial traction) (38)cont.(49)cont.(41)
  α < ψ < 1 00   1 ψ   ψ α 00
(coasting)(38)(38)(48)(49)(41)(41)
  ψ = α 0   1 α 000
(partial reg. brake)(38) (48)cont.(41)cont.
  ψ < α 0   g br   1 ψ 00   α ψ
(full reg. brake)(38)(38)(48)(41)(41)(49)
Table 2. Parameters of Example 1.
Table 2. Parameters of Example 1.
ParameterSymbolValue
Number of coachesn6
Mass of locomotive   m loc   84   t
Total mass of trainm   414   t
Engine efficiency   η tr   0.85
Efficiency of regenerative brake   η br   0.85
Max. mechanical power for traction   P tr , mech   5.6   MW
Max. mechanical power for regenerative brake   P br , mech   5.6   MW
Max. force of regenerative brake   F br , lim   240   kN
Air drag coefficient for locomotive   c w , loc   0.26
Air drag coefficient for first coach   c w , first   0.13
Air drag coefficient for middle coaches   c w , middle   0.10
Air drag coefficient for last coach   c w , last   0.23
Rolling friction coefficient   c roll   0.0015
Coefficient accounting for rotating masses   c rot   1.08
Table 3. Speed limit in Example 2.
Table 3. Speed limit in Example 2.
From Position s [km]To Position s [km]max. Speed  v max  [km/h]
05.5160
5.57.0110
7.09.6150
9.612.0105
12.020.0140
Table 4. Port connections in Example 2.
Table 4. Port connections in Example 2.
Port ConnectionTake Off and Landing TypeRemark
  1 2 (T1), (L2)K has a kink point when full traction changes to coasting.
  2 4 5 (T3), (L3)New SL-TP 4 included, connected to 5,  ψ  jumps at 4.
  5 6 (T5), (L5) ψ  jumps at end of 5.
  5 9 (T5), (L4) ψ  jumps at end of 5. Since  ψ = 1  at the beginning of interval 9, no new speed limit touchpoint is included.
  6 7 (T5), (L1) ψ  jumps at end of 6.
  9 10 11 (T7), (L3)A connection with constant speed, but with a jump of  ψ  at new SL-TP 10, which is connected to 11.
  11 12 (T5), (L1) ψ  jumps at end of 11.
  11 13 (T5), (L2) ψ  jumps at end of 11. Same K-trajectory as  11 12  at the beginning, but then change to coasting.
  13 14 (T3), (L9)Landing at newly included SL-TP 14.
  14 15 (T8), (L1)Take off from SL-TP 14 with  ψ -jump.
  15 16 (T2), (L1)Full traction.
  15 17 (T2), (L10)Coasting to stop.
Table 5. Port connections in Example 3.
Table 5. Port connections in Example 3.
Port ConnectionTake Off and Landing TypeRemark
  1 2 (T1), (L1)
  1 3 (T1), (L2)K-trajectories of  1 2  and  1 3  coincide at the beginning.
  3 4 (T3), (L9)New SL-TP 4 is inside the s-range of interval 5. SL-TP 4 needs to come first in numbering.
  4 5 (T8), (L1) ψ  jumps at SL-TP 4.
  5 6 (T2), (L2)
  6 7 (T3), (L9)SL-TP 7 is inserted, lying exactly between intervals 6 and 8.
  7 8 (T8), (L1) ψ  jumps at SL-TP 7.
  8 9 (T2), (L2)
  8 10 (T2), (L1)Trajectory is close to speed limit but does not touch.
  10 11 (T2), (L10)Coasting to stop.
Table 6. Port connections in Example 4 at final value  ψ T = 0.5237 .
Table 6. Port connections in Example 4 at final value  ψ T = 0.5237 .
Port ConnectionTake Off and Landing TypeRemark
  1 2 (T1), (L2)
  2 3 (T3), (L9)New SL-TP 3 is inside the s-range of interval 4. The SL-TP 3 needs to come first in numbering.
  3 4 (T8), (L1) ψ  jumps at SL-TP 3.
  4 5 (T2), (L2)
  4 6 (T2), (L7)Landing on start of interval 6 with  ψ = α . Any jump of  ψ  here would lead to coasting and violate speed limit. Therefore, no new speed limit touchpoint is inserted.
  6 7 (T7), (L1) ψ  jumps at end of PB-SL 6.
  7 8 (T2), (L2)
  7 9 (T2), (L9)New SL-TP 9 is exactly between intervals 8 and 10.
  9 10 (T8), (L1)Take off from SL-TP 9.
  10 11 (T2), (L10)Coasting and finally full regenerative braking to stop.
Table 7. Electric energy demand at arrival for optimum strategy and variants 1–7.
Table 7. Electric energy demand at arrival for optimum strategy and variants 1–7.
CaseTotal Electric Energy  E tot  at ArrivalPercent Regarding Optimum
Optimum353.0 kWh100%
Variant 1363.1 kWh102.9%
Variant 2372.1 kWh105.4%
Variant 3377.2 kWh106.9%
Variant 4415.0 kWh117.6%
Variant 5367.1 kWh104.0%
Variant 6366.5 kWh103.8%
Variant 7368.9 kWh104.5%
Table 8. Speed limit in Example 6.
Table 8. Speed limit in Example 6.
From Position s [km]To Position s [km]max. Speed  v max  [km/h]
010.9160
10.912.9130
12.936.2160
36.238.7130
38.750.2120
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Heineken, W.; Richter, M.; Birth-Reichert, T. Energy-Efficient Train Driving Based on Optimal Control Theory. Energies 2023, 16, 6712. https://doi.org/10.3390/en16186712

AMA Style

Heineken W, Richter M, Birth-Reichert T. Energy-Efficient Train Driving Based on Optimal Control Theory. Energies. 2023; 16(18):6712. https://doi.org/10.3390/en16186712

Chicago/Turabian Style

Heineken, Wolfram, Marc Richter, and Torsten Birth-Reichert. 2023. "Energy-Efficient Train Driving Based on Optimal Control Theory" Energies 16, no. 18: 6712. https://doi.org/10.3390/en16186712

APA Style

Heineken, W., Richter, M., & Birth-Reichert, T. (2023). Energy-Efficient Train Driving Based on Optimal Control Theory. Energies, 16(18), 6712. https://doi.org/10.3390/en16186712

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop