Optimization for Train Energy Performance

(1)

uptecf14052

Examensarbete 30 hp December 2014

Optimization for Train Energy Performance

Johan Brändström

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Optimization for Train Energy Performace

Johan Brändström

In many studies efforts are made to decrease the energy consumption of trains by optimizing their drive style, e.g. accelerate and brake optimally and regenerate electricity when braking. In other studies the goal is to distribute the run time between stations in an optimal way to decrease the energy consumption, given a relatively simple drive style. In

this report the goal is to combine these two energy saving methods to obtain as low energy consumption as possible. By coupling one software containing a drive style optimizer with another software which by different optimization methods calculates the optimal run time distribution on a given track this is accomplished. The study also contains a comparison between drive styles, with the goal to find a relatively simple but energy efficient drive style. Finally the dependence between run time distribution and energy consumption is further analysed.

The results show that by redistributing the run times the energy consumption can be decreased compared to previously existing time tables. They also show that a relatively simple drive style gives

comparable energy consumption compared to the one obtained using a drive

style optimizer. Finally the results show that the dependence between run time and energy consumption can be approximated with a simple second order equation.

(3)

Populärvetenskaplig sammanfattning

I många studier görs försök att sänka tågs energiförbrukning genom att optimera deras körstil, dvs gasa och bromsa vid rätt tillfällen samt återmata elektricitet vid inbromsning. I andra studier är målet att distribuera gångtiden mellan stationer på ett optimalt sätt för att sänka energiförbrukningen, givet en relativt enkel körstil. I den här rapporten är målet att kombinera dessa två energibesparande metoder för att få en så låg energiförbrukning som möjligt. Detta görs genom att koppla ihop en mjukvara innehållande en körstilsoptimerare med en annan mjukvara vars uppgift är att genom olika optimeringsmetoder räkna ut optimal tidsdistribuering på en given bana. I studien ingår även en jämförelse mellan olika körstilar i syftet att hitta en relativt enkel men energieffektiv körstil. Slutligen analyseras även relationen mellan gångtidsdistribuering och energiförbrukning djupare.

Resultaten visar att genom att distribuera om gångtiderna kan energiförbrukningen sänkas jämfört med tidigare tidtabeller. De visar också att en relativt enkel körstil ger om något högre ändå jämförbar energiför- brukning med den från körstilsoptimeraren. Till slut visar resultaten att sambandet mellan gångtid och energiförbrukning kan approximeras med en enkel andragradsekvation.

(4)

Acknowledgements

I am using this opportunity to express my gratitude to everyone who supported me throughout the course of this master thesis. I am sincerely grateful for the constructive help and path finding thoughts that have helped guiding me to the final results of this report. I express my warm thanks to Astrid Herbst, Johan Lundin, Rasmus Myklebust and Erik Wik for their expertise, support and guidance at Bombardier. Fur- thermore I thank Michel Chapuis, Cecilia Söderberg and Karl-Johan Åhs giving me personal and practical support throughout the whole project. Finally I want to thank my girlfriend Monica Ricão, my parents Bir- gitta Brändström and Tomas Brändström and everyone else that has given me support to finish this project.

Thank you!

Johan Brändström

(5)

Abbreviations and Units

Table 3.1: Abbreviations

Abbreviation Significance

SSH Secure Shell

SET Server name

GUI Graphical User Interface

TEP Train Energy Performance

MF ModeFrontier

DAS Drive assistance system

AOD All-Out Drive style

AODC All-Out Drive style combined with Coasting

GA Genetic Algorithm

DOE Design of Experiments

ULH Uniform Lattice Hypercube

MOGA Multi Objective Genetic Algorithm

RMS Root Mean Square

SVD Singular Value Decomposition

NN Neural Networks

RBF Radial Basis Function

MAE Mean Absolute Error

MRE Mean Relative Error

AIC Akaike Information Criterion

GMTD Generic Metro Track Data

json File name extension of used in- and output files in TEP mat File name extension: Microsoft Access Table

EC Energy Cost = Cost in Euro per kWh ES Energy Saved = net energy saved per km

(8)

Table 3.2: Units

Symbol Unit Significance

rti s run time between station i and i+1

rt_i_min s minimum run time between station i and i+1 rti_max s maximal run time between station i and i+1

rt_i_j s run time between station i and i+1 using design j rttot s total run time between first and last station

vmax m

s maximum allowed velocity on a given track

vbrake m

s velocity of train when starting braking before a station sc m distance from station to point where train starts to coast sb m distance from station to point where train starts to brake

sacc m acceleration distance

sconst m distance running on vmax: sc−sacc

scoast m coasting distance vmax: s_b−sc

stot m total distance between station i and i+1

s_brake m braking distance: stot−sc

ti s arrival/departure time at station i

tacc s acceleration time

tconst s time running on vmax: tc−tacc

tc s time from station to point where train starts to coast

t_b s time when train starts braking

tcoast s coasting time: rti−tc

t_brake s braking time: rti−t_b

rttest s run time calculated in script testing a value of vbreak

rtact s run time calculated in TEP using AODC and coasting point from script

rtdi f f s rti- rtact

∆t s rti_j−rti_j₊₁

Etot kWh Total energy consumption

Ereg kWh Total regenerated energy

Enet kWh Etot−Ereg

Ebas kWh Enetcalculated on GMTD using Baseline d/y km/year The distance a train is projected to run over a year

fi - Prediction value for design i when calculating MAE

y_i - Real value for design i when calculating MAE

(9)

Introduction

4.1 Background information

During the last couple of decades trains have become more and more important in the transportation sector, partly because of their comfort and relatively small pollution contribution compared to other transports.

To be economically comparable to other transports it is important to always increase the efficiency of the trains with regard to energy consumption. This can be done in different ways, e.g. by optimizing drive style using energy saving methods or by distributing the run time between stations optimally. Schedule optimization is of great importance to train operators when minimizing energy consumption. There are many studies on energy optimization with regard to either drive style or time table optimization, but they are rarely connected in one study [4].

An important part of Bombardier is to deliver Metros to customers. In order to be able to make a competitive bid it is important to keep down the energy costs, since they are a big part of the total cost.

4.2 Objectives

By coupling an optimization software with an energy performance tool developed by Bombardier different parametrizations are investigated in this paper. The two main objectives are schedule improvement and drive style optimization. The schedule improvement aims at lowering the fuel consumption for trains by optimizing the time distribution in the schedule. Desired is to find a drive style that is both easily implemented but also energy efficient. The drive styles are implemented in an energy performance tool, that also calculates the energy consumption between every station. The properties and accuracy of one of the drive styles have not been thoroughly tested and the last part of this project aims at comparing this optimizer with a simpler drive style. Furthermore the dependence of the input parameter run time is investigated with the help of different response surfaces. The accuracy of the algorithms used to create the response surfaces are investigated in different ways.

4.3 Scope

4.3.1 Limitations of the study

The project was limited to make the scope of the project appropriate:

• When comparing drive styles three stations were used although if more time were given seven stations would have been preferable.

• The time of the schedule optimization. A schedule optimization can theoretically run for weeks, which in this project was obviously not possible.

• The input parameter range was restricted. Increasing it drastically increases the calculation time and therefore the range was chosen as wisely as possible to not loose accuracy in the results.

(10)

Theory

5.1 Time table optimization

There are different parameters that can affect a time table, and therefore different ways that a time table optimization can be made. In some studies the goal is to minimize the time table delays [1], while in other different constraints such as station capacities [2] or passenger waiting time [3] are used. Usually the time tables are set using these constraints and the energy optimization lies in the drive style [4]. In this report the focus lies on lowering the energy consumption and optimizing the schedule and the drive style on a single track while outer circumstances such as delays or passenger flow are left out. The energy consumption calculations are made in one software, while the optimization is made in another. The basics of these software and drive style properties are explained in the following sections.

5.2 TEP - Train Energy Performance

The last couple of years Bombardier has developed a software called Train Energy Performance (TEP), which allows a detailed model of a rail vehicle to be simulated under a variety of conditions. The software utilizes detailed loss models for many of the electrical components, and has been validated against real- world energy consumption tests with satisfactory results. The program can be executed using complex input data regarding the train, the track, surrounding environment, and drive style. Many train energy optimization software use simplified models of the train (e.g purely mechanical models and constant train- set efficiency). The models in TEP are much more detailed with regard to the propulsion system losses, thus offering higher accuracy than many other software train models. With this software the models can be made very precise, opening up to improved results in optimizations.

5.3 Drive style optimization

There are different ways the drive style of a train can be optimized to decrease energy consumption. This chapter explains the most basic drive style often called all-out drive style (AOD). Some energy saving techniques are explained that lay the ground for energy efficient driving which is described at the end of the section.

5.3.1 All out Drive Style

The AOD is perhaps not famous for being the most energy efficient drive style. When a vehicle applies AOD it drives at maximum allowed speed all the time, and it also utilizes maximum traction and braking force when accelerating and braking [6]. Consequently, the AOD yields the shortest travel time for any

(11)

sconst= s_b−sacc

and

tconst= tb−tacc

where sbis the point where the train starts to brake at time tb.

• Braking: The last part of the AOD consists of braking until the station is reached. The distance covered under this part is referred to as s_brakewhile the corresponding time is t_brake:

s_brake= stot−s_b and

tbrake= rt−tb

where stotis the total distance between two stations and rt is the run time between two stations.

Figure 5.1: Visualization of a speed profile using an AOD

5.3.2 Energy saving techniques

There are several ways to save energy on drive style. In this section the most commonly used and efficient are explained.

Coasting

Coasting is an energy-saving technique that is commonly combined with the all-out drive style in order to save energy [6]. This drive style is referred to as AODC (All Out Drive style combined with Coasting) in this report. When the train has run with maximum allowed velocity for a while the propulsion equipment is turned off. Coasting makes it possible for the train to lose some speed before it is time to brake, which naturally leads to a lower average speed. With the propulsion equipment turned off, the electrical losses in the propulsion equipment ceases. No energy is thus drawn from the line and ideally no energy is lost. This

(12)

technique is especially efficient when running in "start-stop" traffic, common on tracks with short distances between stations [7]. In figure 5.2 the principle is shown. The coasting time, tcoast, and coasting distance, scoast, can be relatively complicated to calculate since the velocity depends on many different surrounding conditions such as track gradient, curve radius, air resistance etc.

Figure 5.2: A visualization of AODC

Many stations are located higher than the rest of the track, giving the train an advantage when accelerating. If the propulsion system is turned off when it is approaching a station the kinetic energy will be lost due to running resistance instead of by braking.

Feedback

A technique that has been utilized the last couple of years is called feedback and can be used by trains equipped with a feedback system in order to regenerate electricity to the line when braking [7]. The technique is commonly used both in subways and bigger trains. In order to use the technique it is preferable if other trains connected to the same line can use the energy simultaneously, even if it can be saved in conductors to some extent.

Energy optimizer of Drive Style (EODS)

Many studies made on drive style optimization are made on an arbitrary vehicle using simple energy models [8]. A drive style optimizer often minimizes the power, or work done over time [9]. In some studies train energy performance optimizations using the above mentioned energy saving techniques feedback and coasting are made [10].

In TEP an energy optimizer developed by Bombardier using feedback is implemented. The physics behind the models used in TEP are more or less the same as in the studies mentioned above.

(13)

where n is the number of values of the current I_igiven in different time steps.

If the train has to use a lot of traction force, for example in a steep uphill, these currents can increase.

For some optimizations it might be of interest to keep the RMS-currents low in order to decrease the temperature in the engine.

5.4 Optimization algorithms

When solving an optimization problem the goal is to find either the maximum or minimum value of one or more functions. The task can consist of one or more input parameters and the goal is to find the values of these parameters that generate the highest or lowest goal functions. The more parameters and goal functions there are the more complicated the problem becomes. Optimizations with one goal function are called single-objective optimizations, while optimizations with more than one goal function are called multi-objective optimizations [11].

Mathematically the functions that calculate the goal functions can be divided into linear and non-linear.

In a linear function the extreme values lie at the end of the input parameter ranges and therefore the extreme values are not very complicated to find. The goal for an optimization problem is to find the global extreme value, but in non-linear functions local extreme values can cause troubles. A good optimization algorithm has to be able to overcome this problem and to be fast at the same time. The following sections go through the most common optimization algorithms solving both single and multi-objective problems.

The structure of the theory is based on the structure in the optimization software ModeFrontier (MF).

5.4.1 Design of Experiments - DOE

The first step before running an optimization is to choose initial designs. This step is very important since it affects the speed of the convergence [12]. The choice of starting points is commonly referred to as Design of Experiments (DOE). There are many different ways to create a DOE and some are more sophisticated than others. A good sampling method manages to avoid clustering in order to make the algorithm test the whole sampling space. If the starting points are clustered the risk of stopping at a local extreme point during the optimization increases. When running the DOE it is also important to have in mind how many starting points that are needed. The desired number of points depends on which algorithms that are used, some need more than others. In the following paragraphs some of the most common sampling methods are explained.

Random Sampling

This method simply generates randomly seeded points over the whole sample space [13] as can be seen in figure 5.3a. The input values X = X1, ... ,XNdo not depend on earlier generated values.

Stratified Sampling

Stratified sampling divides the sample space into equally large sampling spaces, S_i[13]. A random sample, X_i, is then obtained from each section as can be seen in figure 5.3b. The sectioning of the sampling space helps this method avoiding clustering of the input values.

Latin Hypercube Sampling

This method has some similarities with the stratified sampling, but instead of smaller sampling spaces the Latin Hypercube uses columns and rows [13]. Firstly the sampling space is divided into N columns and rows, where N is the number of input values. When the first value is chosen its column and row is stored so that the next value will not share either of them. This way the input value distribution becomes relatively scattered as can be seen in figure 5.3c.

(14)

(a) Example of random sampling

(b) Example of stratified sampling

(c) Example of Latin Hypercube Sampling Figure 5.3: Three different input parameter sampling types

Factorial designs

A factorial design is an example of a statistical design, which means that it is suitable for statistical analysis [14]. It can either be a full factorial which means that every combination of the input parameters are investigated in order to find out how the output responds to changes of the input parameters. A reduced factorial is commonly used when the number of designs becomes too large on a full factorial and the user wants to save time. The number of designs in a full factorial is equal to 2^kwhere k is the number of input parameters.

5.4.2 Hill Climbing Algorithms

A common way of solving non-linear optimization problems is to use a hill climbing method. This method aims at finding the global maximum of a function and is often used to solve single-objective problems [15].

These kinds of algorithms have a tendency to "get stuck" at local extreme points. This can be avoided by only using them on unimodal functions, which means that the function only has one extreme value in the given region. SIMPLEX is a very efficient hill climbing method which has been combined with newer algorithms in order to be used on multi-objective problems.

The SIMPLEX algorithm

A simplex is a geometrical figure that can consist of different number of points [15]. The SIMPLEX algorithm is started by creating a polytope consisting of N+1 vertices, where N is the number of dimensions of the input parameters [16]. This means that the number of input values from the DOE is N+1. The method iterates the function values and moves the vertices of the polyhedron towards the optimal result by comparing the values in the corners. The following example is a 2-dimensional problem, but the method can be used in any number of dimensions. In two dimensions the polytope becomes a triangle.

Firstly the objective functions for the different points are compared in order to find the point with the worst value. For a minimization problem the highest value is the worst while the lowest is worst in a maximization problem. The next step is to reflect this point in the opposite face. In a triangle with the corners A,B and C, assuming that corner A has the worst value, the reflection will become as in figure 5.4.

(15)

Figure 5.4: Reflection of a simplex using 2 dimensions.

Now two things can happen, either the new corner gives a better value than the others, or it is still the worst one. In the first case the next step is simply to reflect the worst point. In the second case this can not be done since it would result in an endless loop, therefore the second worst point is reflected instead. This will continue until a global maximum has been found, resulting as in figure 5.5.

Figure 5.5: Result after 13 reflections using the SIMPLEX method

In the Nelder-Mead method some modifications have been made to make the method more efficient [17]. To find a more precise value something called contraction can be used, as visualized in figure 5.6.

By decreasing the distance between the corners the method approaches the optimal value with a higher accuracy.

(16)

Figure 5.6: Contraction of a simplex in 2 dimensions.

In order to reach the optimal value faster something called expansion can be used, as shown in figure 5.7. If the reflected value is the best value, it is quite probable that the expansion will take the triangle closer to the optimum than the original method.

Figure 5.7: Expansion of a simplex in 2 dimensions.

5.4.3 Multi Objective Algorithms

In most engineering problems there are more than one goal function that has to be optimized. When optimizing with regard to more than one goal function the hill climbing algorithms can not be used, instead so called Genetic Algorithms (GA) are commonly used in this purpose [18]. The GA uses the evolutions theory "survival of the fittest", which means that the weakest species or individuals will extinct by natural selection while the strong ones will pass on their genes to the next generations. Small changes in the gene pool that give advantages will survive while those who give disadvantages will become eliminated by natural selection as can be seen in figure 5.8.

(17)

Figure 5.8: Inheritance and mutation according to the survival of the fittest-theory and as used in genetic algorithms

One general approach that this algorithm uses is to find a set of solutions called Pareto solutions. These solutions are non-dominated with respect to each other and other solutions. This means that none of these solutions is better than all the others with respect to all objectives, there always has to be made sacrifices on one objective to make another better. This information is used in order to find better designs in the next generations.

When the Pareto solutions are found, one way to analyze the result could be to use a weight factor that decides how important the different objectives are. If all the objectives y are multiplied with their corresponding weighting factor w as in equation (5.2), where i is the number of objectives, the problem can be turned into a single-objective optimization.

y=w1∗y1+w2∗y2+...+wi∗yi (5.2) The following commonly used multi-objective algorithms are based on genetic algorithms.

MOGA

MOGA is a typical example of a GA and works in the following way: When the initial points have been chosen in the DOE the function values are calculated. Each individual in the population is then assigned a rank depending on how good the function value is. An individual which has no other individual that is better with regard to all the objectives, thus being a Pareto solution, is given rank 1. Furthermore the individuals with worse value than all the individuals with rank 1 (with respect to all objectives), but better than the rest of the individuals, is assigned rank 2. Then the rank decreases with the number of solutions that beat a certain individual, which can be seen in figure 5.9. Individuals with the same rank have at least one objective in which their value is better than the other with the same rank.

The next step is for the most fit individuals to pass on their genes two by two to the next generation. The method used for this is called inheritence. When the fittest individuals are chosen the rest of the population

(18)

is filled up with mutations. By stochastically modifying some parameters on some individuals a so called mutation method will make it possible for them to go on to the next generation, even if they are not the fittest ones. This is where GA is different from other algorithms, making it powerful to avoid local extreme points.

Figure 5.9: Rankning method used in multi objective genetic algorithms

MOGA-II

With the help of two different methods an improved version of MOGA called MOGA-II has been developed [19]. One of the methods is called elitism and is used to make the convergence more robust. Every iteration step the goal function value from the most fit individual is saved, and the method continues until the elitism is increasing. The other method is called crossover and is used to find a faster convergence. In directional crossover the fitness value of different individuals are compared in order to evaluate a direction of improvement. The parameters of the individuals will be changed randomly within the direction of improvement to create a new generation. The number of generations decides how long the goal function will converge. There is no rule of thumb showing how many generations that should be used, but the more generations that are used the better the result. The number of generations is simply restricted by the calculation time.

(19)

Figure 5.10: Directional crossover between the individuals Ind_i, Ind_jand Ind_k

5.5 Response surface method - RSM

A response surface is an approximation of the goal function in the whole solution space. They are commonly used to give the user an idea of how the goal function responds to changes in the input parameters.

To create a response surface first a couple of designs are chosen to calculate the goal function in those values.

An algorithm then approximates the values of the rest of the solution space.

5.5.1 Response surface calculation models

In ModeFrontier various models can be chosen to create a response surface. These can be divided into two groups; approximating and interpolating RSM:s, where the interpolating RSM:s go through the training points, while the approximating RSM:s do not. Some are more advanced than other but they all have the same purpose; to create a response surface given a certain number of points in the design space. In the following paragraphs some of the most commonly used RSM:s are explained.

Polynomial Singular Value Decomposition

One of the most basic models is the Polynomial Singular Value Decomposition (SVD). This model creates a fitting polynomial by minimizing the squares of the error predictions [20]. In ModeFrontier the degree of the polynomials varies between 1 and 10. This is not the most accurate response surface model, but can generate a reliable guess of the behaviour of the response which can be used to detect specific regions of interest.

Neural Networks

Neural Networks (NN) is an advanced method for creating response surfaces using the brain structure as inspiration. This method can approximate many different types of functions with desired accuracy, making it very efficient if the function is not known. The algorithms used in the model are based on layers, where the input parameters create an input layer as a starting layer [21]. By using a non-linear function, a so-called hidden layer is created and the final output-layer found by interpolation of the hidden layer.

Radial basis function

The radial basis function (RBF) is an example of an interpolating RSM [20]. The RBF works well when working with scattered data, meaning that the data is not sampled on a regular grid. The RBF interpolates

(20)

the function using the input data and a certain radial function. There are five different radial functions available in ModeFrontier. These are Gaussians (G), Duchon’s Polyharmonic Splines (PS), Hardy’s Multi- Quadrics (MQ), Inverse MultiQuardics (IMQ) and Wendland’s Compactly Supported C²(W2). All of these are commonly used in the field and could well work for the kind of problems that this report solves except the W2 that can not be used with more than 5 dimensions.

5.5.2 RSM validation methods

There are several ways to validate if a RSM is a good approximation of the solution to a problem. In the following subsections some constants that can give information about the error are described as well as the more complex constant Akaike Information Criterion (AIC).

Mean absolute error

One of the most simple but efficient ways is to calculate the mean absolute error (MAE). The error, or the difference between the real values and the calculated ones, is also commonly referred to as the residual.

The absolute error is calculated as in equation 5.3.

e_i = |f_i−y_i| (5.3)

Here f_iis the function value of design i and y_ithe true value of design i. The absolute mean error is obtained by summing up the mean errors and dividing the result by the number of designs as in equation 5.4.

MAE= ¹ n

∑

n i=1

|fi−yi| = ¹ n

∑

n i=1

|ei| (5.4)

Mean relative error

Another validation tool that can give additional information about the accuracy of a model is the mean relative error, which is calculated as in equation 5.5.

MRE= ¹⁰⁰ n

∑

n i=1

|f_i−y_i| yi

= ¹⁰⁰ n

∑

n i=1

|e_i|

yi (5.5)

Basically the absolute error is divided by the function value in every point and the total sum is multiplied by 100 to obtain the answer in percent.

Akaike Information Criterion

One way to measure the quality of RSM is to use the Akaike information criterion (AIC). The AIC does not give any information about how well the model fits the data, instead it estimates how much information that is lost when generating the model [22]. AIC is calculated as in equation 5.6.

AIC=2k−2 ln L (5.6)

(21)

Over-fitted functions

When making a polynomial-fitted curve a higher polynomial usually gives a smaller value of the MAE for the training points [23]. However when the polynomial increases the risk for over-fitting occurs, which means that the fitted curve starts to prioritize the accuracy in the training points rather then learning and generalize from them, which results in less accuracy in the new predicted data. If a function is over-fitted this will be noticed in the AIC.

(22)

Methods

The first part of the project consisted of setting up the environment which is the foundation to both the time table optimization and the drive style comparison. The train energy performance was optimized by coupling TEP to the optimization software ModeFrontier (MF). TEP was installed on a Linux server while MF was run from a workstation using Windows. By using a secure shell (SSH) MF was connected to the Linux server. After coupling the software, optimizations using different set-ups and parametrizations were made and different drive styles compared. Lastly different response surfaces were made in order to find an input parameter dependence, whereupon the response surface methods were examined.

In the following sections details about the set-ups, parametrizations and response surface methods are explained.

6.1 TEP set-up

TEP can be run in two different modes, either using a Graphical User Interface (GUI) or in batch mode, a mode where no GUI is needed. When coupling TEP to ModeFrontier, TEP must be run in batch mode.

The optimizations were made in MF while the energy consumption calculations were made in TEP. The optimization process can be seen in figure 6.1 where the in- and output values were saved in a file with the file name extension "json".

(23)

In TEP the drive style properties of the train can be chosen. The ones provided are AOD, AODC and one drive style optimizer called Drive Assistance System (DAS). The details of these drive styles are explained in section 6.2.

6.2 Parametrization

The two drive styles that were used in this project were DAS and AODC. The drive style optimizer DAS chooses where the traction- and braking force is applied in order for the train to drive in an energy efficient manner. This is done by testing all the available trajectories in every point of the track and choosing the ones with least consumption for every point. When using DAS it is possible to manually set the time table before running a simulation.

When applying DAS in TEP the user can choose a certain time table that the train shall follow. To help creating a speed profile that results in a certain time, a script that calculates braking points, s_b, is implemented in DAS. The script uses the value of maximum allowed braking to create a curve that shows when the train has to start braking in order to stop at the next station, as shown in figure 6.2.

Figure 6.2: Graphical explanation of braking point calculation in DAS.

When using AODC, TEP is missing the explicit manual setting of time table. In section 6.2.2 it is described how the time table indirectly was set in order to make the comparison between drive styles. Re- gardless of drive style and set-up the following properties were used in every optimization:

• For all optimizations the genetic algorithm MOGA-II was used.

• The run time:

rti=ti+1−ti

which is the travelling time between station i and i+1, leaving from station i at ti and arriving at station i+_{1 at t}_i+1, was used as an input parameter.

• The net energy consumption:

Enet=Etot−Ereg

(24)

where Etotis the total consumed energy and Eregis the regenerated energy.

• Track data was taken from a generic metro network. This track data is referred to as GMTD - Generic Metro Track Data. In figure 6.3 the altitude of the whole track can be seen. For all the optimizations parts of this track were used.

Figure 6.3: Altitude as a function of distance for the GMTD.

6.2.1 Optimization set-ups

Test set-up (*)

Before making a final set-up that could be used both for the drive style comparison and time table optimization a test set-up was made, partly to test the connection between the two software TEP and MF, but also to see that the optimization went in the right direction. By letting the train run between two stations, the fact that a longer run time should result in a lower energy consumption was tested. The following list summarizes the set-up while it is visualized in figure 6.4.

• Goal function: Minimize Enet.

• Input parameter: rt .

(25)

Figure 6.4: Visualization of DAS test set-up, running between two stations.

Final set-up (**)

This parametrization was partly implemented to optimize the distribution of the run times and compare the result from the run time distribution optimization with results from an earlier optimized time table, which will be referred to as "baseline". The parametrization was also used for the DAS and AODC comparison.

In this set-up stations were added and hence stop times. The value of the stop times were set to constant values. The total time was set to a constant value and the run time distribution was optimized in order to minimize the Enet. The following list summarizes the set-up while it is visualized in figure 6.5.

• Goal function: Minimize Enet.

• Input parameters: rt_i, where 1≤i≤24 for the baseline comparison and 1≤i≤2 for the drive style comparison.

• Constraints:

– Constant total time, rttot.

– Minimum value for each run time, rt_i_min. – Maximal value for each run time, rt_i_max.

Figure 6.5: Visualization of the final DAS set-up, running between three stations.

6.2.2 AODC set-up

The set-up (**) in section 6.2.1 was used for the AODC optimizations. Running TEP, rti can not be set explicitly using AODC. This was solved by setting all track- and train data to constant values with two

(26)

exceptions. While letting scbe the parameter that regulated rt_ifor the given train- and track data, vmaxwas used as an input parameter. The reason for letting vmax be an input variable was simply because it can be of interest to know how fast the train should go to consume as little energy as possible. The value of sc, that would generate a certain rtifor given input parameters, was found by using a Matlab script (coast script).

The following parts of rtiare used in the coast script:

• Acceleration time, t_acc: To calculate this time TEP was run separately with given track and train properties. The time to reach the given maximum velocity, as in figure 6.6, was measured and manually inserted in the coast script.

• Time with constant velocity, tconst: tconst = ^s_v^const

max

where sconst =sc−sacc, where saccis the acceleration distance. The maximum allowed velocity on the track was modified to one value between every station to make the calculations easier.

• Coasting time:

tcoast=t_b−tc

where t_bis the time when the train starts to brake and tcis the time the train starts to coast.

• Braking time:

tbrake=rti−tb

(27)

In order to calculate tcoastand t_brakefor given vmax and train- and track parameters, the coast script utilizes the DAS-script containing the braking velocity calculations explained earlier in this section. In the DAS-script specific track- and train data are inserted into a mat-file, a file that can contain both numbers, tables and arrays. With the help of this the DAS-script calculates tcoastand t_brakefor a given value of v_brake, which is used as input to the DAS-script. The coast script starts setting v_braketo zero, calculating the ’tested’

value, rttest, as follows:

rttest =tacc+tconst+tcoast+tbrake

If rttest<rtithe coast script ends with an error message. If rttest >rtithen vbrakeis increased by a small value, which in this case was set to e=0.1, and the coast script calculates the total run time with the new value of vbrake. The steps of the DAS-script are summarized in figure 6.7.

Figure 6.7: Flow diagram summarizing the steps in the script calculating coasting point.

After calculating sc, the speed profile generated by the DAS-script was compared to one received from running TEP using the resulting sc, the used value of vmax and train- and track parameters. Furthermore rt_i was compared to rtact, the actual received run time from TEP. For many combinations of maximum velocities and run times no coasting point could be found that resulted in the desired run time. This is due to the short distance between stations on the used track. If the distance is short there is less room for coasting. To open up for more possible combinations the train was scheduled only to stop every second station.

(28)

6.3 ModeFrontier set-up

The optimization described in the previous chapter takes place in MF. All the commands in MF are con- trolled by nodes that each have different properties. One node was used for the coupling by using SSH- information. This way the input- and output files were sent back and forth between the Linux server and the workstation. By running the terminal on the server through the SSH node TEP was started in batch mode. This way the energy consumption was calculated for every set of input parameters. In the following subsection the node configuration in MF is explained.

6.3.1 Node configuration using DAS

The basic structure of the node configuration using DAS is shown in Figure 6.8. The nodes were used as follows:

Figure 6.8: ModeFrontier Structure using DAS.

• DOE: Sets Design of Experiment, throughout the whole project ULH (Uniform Lattice Hypercube) was used.

• MOGA-II: Node that chooses optimization algorithm. For all optimizations MOGA-II was used.

Number of individuals in every generation was chosen by using the rule of thumb: 2*(Number of input parameters)*(number of output parameters) [User manual MF].

(29)

• Input file: A node that parses the input file that is send to the Linux server.

• SSH node: Contains SET-server name and login information.

• File transfer: Node that copies the output file from the SSH-server and sends it to a folder where it can be found by the energy template node.

• Energy template: Searches the output file for DAS errors and Enet. This is done by setting a number tag that has a relative position to a certain string.

• If-node: Node that discards designs that do not fulfil the constraints.

• Warnings: Node that informs the if-node if the design should be discarded or not.

• Calc energy: Calculates Enetin kWh.

• Net energy output nodes: Transfer the values of the energy to the next nodes.

• Min energy: Node that decides whether the goal function should be maximized or minimized.

6.3.2 Node configuration using all-out drive style

In figure 6.9 the ModeFrontier structure using AODC is displayed. The structure is basically the same as when using DAS with a few exceptions, listed as follows:

Figure 6.9: ModeFrontier Structure using AODC.

• v max: In this set-up arrival times and departure times are not needed since the time table is set by vmax and sc. Therefore a vmax node is used instead of stop time node.

• Matlab node: Instead of the calculator node a Matlab-node is added because the DAS-script is a Matlab-script.

• Coasting point and v maxms: The TEP input parameters are coasting point and vmaxinstead of arrival time and departure time.

(30)

• If node: Discards the designs found in Matlab-script with combinations of vmax and total time that result in an infeasible design.

• Rt min and Rt sec nodes: Nodes that were implemented in order to compare the desired run times with the resulting ones using the given coasting point and maximum velocity.

When both the drive styles had been implemented the modified track data shown in figure 6.10 was used for the comparison between drive styles. The drive styles were compared with regard to Enet, speed profile and run time distribution.

Figure 6.10: Original and modified altitude as a function of distance from first station on GMTD.

6.4 Time table optimization

When optimizing the GMTD-time table, set-up (**) in section 6.2.1 and all the track data in section 6.2 were used. Part of the reason for using DAS was that the results from this optimization were compared with the baseline, for which Enetwas received by using DAS as well. In the baseline, the time table was optimized by manually distributing the different rtibetween stations.

6.5 RMS-impact investigation

The impact of what would happen if using RMS-currents as a constraint was investigated. No extra con-

(31)

in section 6.2.1. Normally a full factorial could be used in this purpose, but for the scope of this project the optimization designs were used.

The MOGA-II optimization results were used to create the RSM-surface. Different RSM:s were made in order to find one that correctly describes how Enetdepends on the distribution of the run time. To examine if any polynomial dependence existed 1st, 2nd and 3rd order SVD polynomials were calculated. These were then compared to the more advanced methods Neural Networks and Radial basis functions.

6.6.1 RSM-validation

In ModeFrontier there are various ways to validate a RSM-model. In this project the absolute mean error between the original design data and the RSM:s were taken in regard. In addition to this the AIC was calculated since it is a good way to validate statistical models. Firstly this was done using only the original data points from ULH and MOGA-II optimization. An additional test of the models was made by adding randomly selected data points and investigate if the absolute mean error changed.

(32)

Results

7.1 Optimization results using DAS

In this section the results from set-up (*) and (**), explained in section 6.2, that were used to optimize GMTD time table are shown. The used drive style throughout the whole section is DAS.

7.1.1 Test set-up using DAS

The goal with the test set-up was to find out if a higher value of the run time resulted in a lower Enet as according to theory. In the first run with this set-up the resolution of the run time,∆t, defined as follows:

∆t = rti_j−rt_i_j₊₁

was one second, where rt_i_j is the run time using design j. Figure 7.1 shows the result from the first test set-up run where the train runs between two stations. The figure shows that for most values of the run time, Enetdecreases with higher run time, with a few exceptions. These exceptions were explained by some properties of DAS, a small change of input parameter can sometimes result in unexpected changes in the output.

Figure 7.1: Test set-up showing Enetas a function of run time using DAS and two stations with variable

(33)

Figure 7.2: Test set-up showing Enetas a function of run time using DAS and two stations with variable run time and six seconds resolution.

7.1.2 Final set-up using DAS

In figure 7.3 the result from the final set-up using three stations is shown. Unlike the set-up (*), the distribution of the run time is of interest. The optimal time distribution according to these results is when rt1= 105 seconds and rt2= 65 seconds.

Figure 7.3: Final DAS set-up showing Enetas a function of run times for three stations using GMTD.

7.2 Optimization preparation using AODC

When creating a method that properly makes the train run between two stations on a desired time using the AODC the speed profile obtained from the DAS-script calculating coasting point was compared to the

(34)

speed profile generated by TEP using the same input parameters.

At first the script generated an error message for most of the combinations of run time and maximum allowed velocity. A few combinations of run time and maximum velocity resulted in a speed profile, but when comparing it with the one obtained from TEP they did not match as can be seen in figure 7.4. The error appeared when the altitude gradient was negative at the point when the train started braking. This was solved by modifying the track data, as seen in section 6.3.2, so that the gradient was never negative on the parts of the track where the train coasts and brakes.

Figure 7.4: Speed profile using AODC and original track data on two stations of the GMTD.

Figure 7.5 shows a part of the speed profile when coasting obtained from the comparison between TEP and the DAS-script using the modified track data. As can be seen in the figure the results are improved, but two differences between the speed profiles still exist. Firstly the slope of the curves differ, and secondly they do not ’turn’ simultaneously. The track data showed that the points where the gradient changed was the same as when the gradient on the track changed. TEP responded to the change in gradient a bit later than the script, which depended on the train length. In TEP the actual train length was used while the script approximated the train as a point.

(35)

Figure 7.5: Speed profile using AODC and modified track data on two stations of the GMTD.

By changing the train length to the lowest allowed in TEP and running the calculations again the speed profiles in figure 7.6 were obtained. The speed profiles now react to the change of altitude gradient simultaneously. Even if a closer look can tell that the slope on some places still differed, the agreement was regarded as sufficient for the scope of the project and the set-up was used in the next AODC calculations.

Figure 7.6: Speed profile using AOD, modified track data, shortened train length and two stations of the GMTD.

(36)

Due to the slope difference rtactdiffered from rt_i. This difference:

rt_{di f f} =rt_i−rtact

was measured by running a simulation using the methods and track data described in section 6.2.2. Since the resolution of the input parameters was six seconds it would be desired for the difference between desired and actual run times to be less than this. In figure 7.7 the difference between rt1and rtactis plotted, while figure 7.8 shows the corresponding information for rt2. The plots show that rtdi f f is less than six seconds for all combinations of rt1and vmax, with a few exceptions. They also show that the error increases with increasing rti.

For one value of rti, vmax differed for every design resulting in different speed profiles. This explains why the time difference is not the same for every value of rti.

Figure 7.7: rt_{di f f} as a function of rt1using GMTD.

(37)

7.3 Comparison between drive styles

The comparison between drive styles was made using set-up (**) in section 6.2.1 and the modified track data described in section 6.2.2. Figure 7.9 shows the resulting speed profile for the four best designs using DAS and three stations, while figure 7.10 shows the corresponding plot using all-out drive style combined with coasting. The first observation is that none of the AODC designs gave as good results as the best DAS-results. The best DAS-design had about 1kWh less energy consumption than the best AODC-result, approximately 5 percent less. The second fact is that the speed profiles have many similarities. DAS chooses to drive with high acceleration until it reaches a high velocity, similar to the all-out drive style. Both drive styles then maintain a high velocity until the train is close to the station where it has to brake.

Figure 7.9: Speed profile from the four best design results using the final DAS set-up on three stations of the GMTD.

Figure 7.10: Speed profile from the four best design results using the final AODC set-up on three stations of the GMTD.

(38)

Figure 7.11 was added to compare how the drive styles distributed the time between stations. No certain pattern for the run time distribution using DAS was observed, while there is a tendency for AODC to distribute 160-170 seconds to rt1and the rest to rt2.

Figure 7.11: Comparison of run time distribution between DAS and AODC on three stations of the GMTD.

7.3.1 RMS-impact

In the comparison between drive styles a possible impact of RMS-currents was investigated using set-up (**) in section 6.2.1. The RMS-currents were measured in the same optimization as earlier in this section. In figure 7.12 the RMS-values both using DAS and AODC are displayed. The figure shows that the trend for both drive styles is that a lower value of Enetresults in lower RMS-value. It also shows that the RMS-values using DAS are lower than when using coasting.

(39)

Figure 7.12: Comparison of RMS-currents between DAS and AODC running three stations of the GMTD.

The upper line shows a linear regression of the RMS for AODC-results while the lower for the DAS-results.

7.4 Time table comparison

The time table comparison was made as described in 6.4. In table 7.1 the results from the final optimization of the time table of a GMTD using a ULH as DOE are shown. In the best design the distribution of run times theoretically saves 1,25 percent energy compared to the baseline. The improvement is calculated the following way:

Improvement (%) =(Ebas−Enet(optimization))/(Ebas) where Ebasis the net energy calculated using baseline.

Table 7.1: Table comparing the eight best results of the MOGA-II optimization of the GMTD schedule using DAS with original time table.

ID Enet(kWh) Improvement (Percent)

1788 125,56 1,25

1834 125,63 1,20

1392 125,73 1,12

1752 125,85 1,03

1562 125,86 1,02

1861 125,86 1,02

1647 125,88 1,01

1751 125,89 1,00

Baseline 127,16 0

(40)

In figure 7.13 the run time distribution for the three best designs are compared to the baseline. The figure shows that the run time distribution does not differ much, in the best time table about 50 percent of the run times are exactly the same as in the manually optimized one, while the rest of the run times do not differ more than six seconds, or one time step.

Figure 7.13: Comparison between baseline and optimized time table using MF.

In figure 7.14 all the design values from the MOGA-II optimization are shown as well as the baseline Enet value. The figure shows that Enet decreased more rapidly during the first 1000 designs of the optimization and then slowed down. The best design was found after 1788 designs and after 1900 designs the optimization was aborted.

(41)

Figure 7.14: Resulting Enetas a function of design ID on MOGA-II optimization of run time distribution on the GMTD time table.

Economical saving results

In order to calculate the amount of money that can be saved on one train driving for one year, with the used optimization tool, equation 7.1 was used. EC is the energy cost/km, ES the Enet saving/km for the operator and d/Y the distance a train covers in one year. The baseline Enet/km was 4,75 kWh/km while the optimized was 4,69 kWh/km. Therefore the used value of the Enetsaving per km is 4,75-4,69 = 0,06kWh.

Saved money per year per train = EC * ES * d/y = 0,1*0,06*150000≈900Euro/year/train (7.1)

7.5 RSM-analysis

In this section the resulting input parameter dependence according to the different RSM:s are displayed as well as the resulting MAE and MRE for every response surface and their AIC-values.

7.5.1 Input parameter dependence according to the response surfaces

The RSM comparison was made as described in section 6.6. In table 7.2 the six best designs from this optimization are displayed.

(42)

Table 7.2: Run time distribution of the six designs resulting in the lowest energy consumptions from optimization of the seven first stations on GMTD.

ID Run time 1 Run time 2 Run time 3 Run time 4 Run time 5 Run time 6 Enet(kWh)

55 105 68 87 89 83 50 25,6

113 105 62 87 95 83 50 25,6

116 105 62 87 89 89 50 25,6

72 105 68 87 89 89 44 25,7

75 105 68 87 95 77 50 25,8

Figure 7.15 shows how Enet depends on the run times according to the RSM made using a first order SVD. The marked points show Enetusing the run time distribution obtained from the best design in in table 7.2 (ID 55). The figure shows that for some run times Enetincreases with increasing run time, while for other it decreases.

The corresponding results from the second order SVD are shown in figure 7.16. The figure shows that the run times from ID 55 all are close to the minimum of the graph.

Finally the third order SVD results are shown in figure 7.17. In these results it is not that clear that the run time distribution from ID 55 is close to the minimum of the run time dependences.

(43)

Figure 7.16: Input parameter dependance according to the 2nd order SVD using optimized time table on the seven first stations of GMTD.

Figure 7.17: Input parameter dependance according to the 3rd order SVD using optimized time table on the seven first stations of GMTD.

(44)

7.5.2 RSM evaluation

Figure 7.18 shows the MAE and MRE (the residuals) between the response surface obtained from the 1st order SVD and the designs from the optimization using MOGA-II. It also shows the residuals between the response surface and randomly generated designs. The figure shows that the residuals for the randomly generated designs do not differ much from the residuals for the optimization designs.

Figure 7.19 shows the same results for the 2nd order SVD, while figure 7.20 shows that the residuals between randomly generated designs and 3rd order SVD are significantly larger than between the optimization designs and the 3rd order SVD.

The biggest residual difference was observed when comparing designs with NN as can be seen in figure 7.21.

Figure 7.18: Relative (in kWh) and absolute (in %) residuals between RSM points and calculated designs using 1st order SVD. The black line separates the original and random designs.

(45)

Figure 7.19: Relative (in kWh) and absolute (in %) residuals between RSM points and calculated designs using 2nd order SVD. The black line separates the original and random designs.

Figure 7.20: Relative (in kWh) and absolute (in %) residuals between RSM points and calculated designs using 3rd order SVD. The black line separates the original and random designs.

(46)

Figure 7.21: Relative (in kWh) and absolute (in %) residuals between RSM points and calculated designs using NN. The black line separates the original and random designs.

Since the RSM for RBF goes through the training points the MAE is zero for the designs used to calculate the response surface. This and all the previous results regarding MAE are summarized in figure 7.22 where the MAE between optimization designs and response surfaces, before and after adding the randomly generated designs, are shown. The figure shows that the MAE is larger for lower polynomial SVD:s and that NN and RBF have the lowest values. It also shows that the residuals on 1st and 2nd order SVD:s were less affected by adding random designs than 3rd order SVD, NN and RBF. Table 7.3 shows a comparison of the AIC-values for all the models before and after adding extra designs.

(47)

Table 7.3: AIC-values of different RSM:s before and after adding randomly generated designs

RSM AIC before adding random designs AIC after adding random designs

1st order SVD -2,4 -5,9

2nd order SVD -83 -103

3rd order SVD -56 -20

NN -155 -31

RBF -6956 -187

When plotting the AIC-results in figure 7.23 a logarithmic scale was used since the range of the values was very wide, and since the AIC-values are negative log(-AIC) was plotted. The figure shows that the value for the 1st and 2nd order polynomial SVD:s increased when adding the randomly selected points, while value using 3rd order SVD, NN and RBF decreased. Since a low value of the AIC implicates a good model, a high value of log(-AIC) would implicate a good model, which means that the RBF and 2nd order SVD have the best best fitted models according to the AIC, although adding randomly generated designs had a more negative effect on them than the 2nd order SVD.

Figure 7.23: log(-AIC)for the different RSM-models before and after adding random designs

(48)

Discussion

8.1 Drive style comparison

A comparison between the drive style optimizer DAS and the simpler drive style AODC has been made.

This comparison was partly made to cross-check the performance of DAS but also to see if a less complicated drive style can match the energy consumption achieved using a more complex optimizer. The drive style as computing using a DAS is hard to implement exactly for a driver and could only be used either as a reference or be implemented with a computer.

The results showed that the implementation of the drive styles still can be improved using DAS. Small changes of the input-parameters resulted in unexpected changes of the output that in some cases led to unexpected errors affecting the results. This is something already noticed by other users and will be further investigated. Furthermore the coast script did not generate the same time as the time received from TEP (rtact was not equal to rt_i). This was due to the slope of the speed profile calculated with the script that differed from the one calculated in TEP.

When calculating the run time using coasting, the DAS-script, used to calculate coasting time, generated an error message for many of the designs and the calculations were stopped. This happened if the train was coasting on a part of the track where the slope was negative. This was solved by modifying the track data so that the gradients were positive or zero on the parts where the train coasted. The comparison between drive styles was still possible to do and the modification did not affect the results.

The results of the comparison showed that DAS gave approximately 5 percent lower energy consumption than the AODC. It was also observed that the drive styles had many similarities regarding the speed profile. By further investigating different parametrizations there is a good potential that the energy consumption using AODC could be decreased further. Instead of using an advanced drive style such as DAS, AODC could be used which is easy to use even for a human driver.

8.1.1 RMS-currents

While comparing drive styles a comparison between the RMS-currents was made. Since the RMS-currents are correlated to the temperature in the engine it is of interest to keep them low. If one of the drive styles would result in high RMS-currents a constraint could be of interest. The results showed that a lower energy consumption in general resulted in lower RMS-currents, regardless of drive style. They also showed that the RMS-currents when using DAS in general were lower than when using AODC. The last part could be explained by the energy consumption being lower for most of the designs using DAS. If a lower value of the energy consumption would have been found using coasting the RMS-current would have decreased according to the trend. Early in the project adding a constraint on the RMS-currents was suggested, but since no maximum value could be found it was chosen only to monitor the RMS-values for each design.

For the results in this project a constraint would not have changed the results since the lowest energy

Optimization for Train Energy Performance

Examensarbete 30 hp December 2014