David Fournier
email: david.l.l.fournier@gmail.com

Thierry Martinez
email: thierry.martinez@inria.fr

François Fages
email: francois.fages@inria.fr

Denis Mulard
email: denis.mulard@ge.com

Metro Energy Optimization through Rescheduling: Mathematical Models and Heuristic Algorithm Compared to MILP and CMA-ES $

Keywords: Energy Optimization, Timetable Optimization, Regenerative Braking, Mass Rapid Transit, Greedy Heuristics

The use of regenerative braking is a key factor to reduce the energy consumption of a metro line. In the case where no device can store the energy produced during braking, only the metros that are accelerating at the same time can benefit from it. Maximizing the power transfers between accelerating and braking metros thus provides a simple strategy to benefit from regenerative energy without any other hardware device. In this paper, we use a mathematical timetable model to classify various metro energy optimization problems studied in the literature and prove their NP-hardness by polynomial reductions of SAT. We then focus on the problem of minimizing the global energy consumption of a metro timetable by modifying the dwell times in stations. We present a greedy heuristic algorithm which aims at locally synchronizing braking trains along the timetable with accelerating trains in their time neighbourhood, using a nonlinear approximation of energy transfers. On a benchmark of the litterature composed of six small size timetables, we show that our greedy heuristics performs better than CPLEX using a MILP formulation of the problem with a linear approximation of the objective function. We also show that it runs ten times faster than a state-of-the-art evolutionary algorithm, called the covariance matrix adaptation evolution strategy (CMA-ES), using the same non-linear objective function on these small size instances. On real data leading to 10000 decision variables on which both MILP and CMA-ES do not provide solutions, our dedicated algorithm computes solutions with a reduction of energy consumption ranging from 5% to 9%.

Introduction

Energy consumption is a major issue for the future and has been the subject of increasing research activities over the last years. The total energy consumed in transportation is estimated to represent 27% of the world energy production [START_REF] Agency | Energy balance for world[END_REF]. In 2006, the London Underground consumed 1173 GW.h [START_REF] O'toole | Environment report[END_REF], representing 2.8% of the Great London total electricity consumption [START_REF] Rose | Sub-national electricity consumption statistics and household energy distribution analysis for 2010[END_REF]. The transportation field is thus of particular importance and industrial companies try to optimize energy by different means, in particular in mass rapid transit such as metros.

Nowadays, almost all metros have energy regenerative braking systems. These devices turn the electric motors into generators during braking, and produce electric energy. It has been shown that the raw energy discount provided by this technology is about 16.5% [START_REF] Greatbanks | Review of the discount for using regenerative braking[END_REF]. Some metros can directly use their own regenerative energy. Super capacitors allow much faster loads and unloads compared to classical batteries, and a metro equipped with super capacitors is able to collect the energy during braking, and give it back to the engine for its own accelerations [START_REF] Ayad | Voltage regulated hybrid DC power source using supercapacitors as energy storage device[END_REF]. The metros that cannot store their own regenerative energy can return it to the DC electrical network, but with important losses on long distances. The electrical substations (ESS) are the devices that convert AC to DC to feed the metro line. Some ESS are revertible and can convert the regenerative power of trains to AC power. In this case, the energy regenerated by metros can be used in other parts of the metro line without important loss, or be sold back to the electricity provider. However, super capacitors, as well as revertible ESSs, are expensive equipment to buy and maintain, and may not be economically justified.

Another way to use the regenerative energy returned to the DC line, without requiring any extra equipment, is to synchronize the braking of a metro with the acceleration of another metro in its close neighbourhood on the DC line. This synchronization can be done, for instance, by modifying the departure and arrival times of the metros in the stations, in order to shift the acceleration phases to the deceleration phases of some other metros in their neighbourhood. The objective function may then combine two criteria: the minimization of power peaks, and the minimization of the global energy consumption during the day. A power peak occurs when too many metros are accelerating at the same time, which may cause two problems. First, a metro network is sized with some maximum power capacity. When the demand exceeds that threshold, a control system prevents the destruction of the electric equipment by a momentary shut down of a part of the network, dropping the quality of service. Second, more than 80% of metro companies pay their energy provider for a certain amount of energy per given time period, typically around 15 minutes [START_REF] Of | Reducing energy consumption in underground systems -an important contribution to protecting the environment[END_REF]. When the energy consumed exceeds a certain limit, the metro company pays a fine to the electricity provider.

Albrecht has shown in [START_REF] Albrecht | Reducing power peaks and energy consumption in rail transit systems by simultaneous metro running time control[END_REF] that it is possible to reduce power peaks by utilizing the reserve time of metros, i.e. the remaining time that a metro has to finish its journey without disturbing the network. It is however tricky to use the reserve time for energy optimization reasons since it is primarily used for traffic regulation. Moreover, this optimization is done by modifying metro interstation times, which may be difficult to implement in a real-time application. Nevertheless, the implementation of this method using a genetic algorithm showed good results in [START_REF] Albrecht | Reducing power peaks and energy consumption in rail transit systems by simultaneous metro running time control[END_REF]. Kim et al. have proposed in [START_REF] Kim | A mathematical approach for reducing the maximum traction energy: The case of Korean MRT trains[END_REF] to optimize the metro departure times in terminals instead of reserve times. They have partially solved a simplified model of this problem using MILP. However, their approximation is not precise enough for a real application since regular timetables are typically second-accurate, whereas their model has a precision of 15 seconds. In [START_REF] Chen | Optimization of an MRT train schedule: Reducing maximum traction power by using genetic algorithms[END_REF], Chen et al. have described a precise electrical network simulator, which leads to an accurate evaluation of the metro power demands. They have managed to reduce the maximum power peak using a genetic algorithm that chooses between only either short or long stopping times for metros in each station.

Concerning the minimization of the total energy consumption on the line, Nasri et al. have shown in [START_REF] Nasri | Timetable optimization for maximum usage of regenerative energy of braking in electrical railway systems[END_REF] that it is possible to decrease the energy consumption by modifying the dwell times, i.e. the stopping times in the stations. Their model uses an exact energy function and an accurate discretization of time at the scale of one second. This work provides an interesting proof of concept but has only be tested on a pilot system, consisting of 4 trains and 4 stations. On a more realistic size problem, Peña et al. [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF] have proposed a MILP model to optimize a night shift timetable for the metro of Madrid by maximizing the overlap of braking and acceleration phases of different metros. Their measurements on field have shown that the optimized timetable could reduce the global energy consumption by 3%.

In this paper, we focus on the problem of minimizing the global energy consumption of a metro timetable by modifying the dwell times in stations. We present a heuristic algorithm which gives better results than general purpose optimization methods. In particular, we compare our algorithm with MILP and with a state-of-the-art evolutionary algorithm called the covariance matrix adaptation evolution strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]. The MILP formulation of Peña et al. [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF] maximizes the overlapping time between accelerating and braking phases and not directly the global energy consumption. Also, to limit the number of binary variables, their model considers only one-to-one pairing between braking and accelerating metros, whereas our objective function takes into account the redistribution of the regenerative energy to several accelerating metros. We show that the MILP approximations of the objective function lead to solutions of lesser quality, eventhough CPLEX [START_REF] Ilog | V12.1: User's manual for CPLEX[END_REF] is able to prove the optimality of the computed solutions. Furthermore, this approach does not scale-up to real data which lead to instances of 10000 decision variables by running out of memory, while handled with our heuristics in 20 minutes. CMA-ES using the same non-linear approximation as us does not scale up on real data either, failing at evaluating the initial population objective functions. On small size instances, our heuristic algorithm also performs better than CMA-ES mainly for the implementation reason that our algorithm computes the objective function incrementally .

The rest of the paper is organised as follows. Section 2 presents a mathematical model of metro timetables taking into account both time constraints and energy optimization objectives. Section 3 uses this model to classify various related problems of the literature, and prove their NP-hardness by polynomial reductions of SAT. Then in Section 4 we model the instant power demand as an electrical network which simulates the energy consumption of the metro line, and describe in Section 5 a power flow approximation based on a distribution matrix. We propose an algorithm to compute this power flow and use it as an objective function for CMA-ES and in our heuristics. This formulation is also used to derive the linear approximation of the objective function for MILP. Section 6 describes our heuristics for minimizing the global energy consumption of a metro line by modifying solely the dwell times in stations. Section 7 compares the results of this heuristics on a benchmark of six small size instances, and shows that it gives better results than MILP and CMA-ES in both computation time and quality of the solutions. The set of these instances is available at http://lifeware.inria.fr/wiki/COR14/Bench. Then in Section 8 on real data, we show that the heuristics is able to reschedule two full timetables containing respectively 9585 and 7679 variables in 20 minutes. On these two examples, the two other methods fail at giving a solution within 30 minutes, CPLEX running out of memory and CMA-ES failing at computing its first iteration, while the heuristics computes solutions which decrease the total amount of energy consumed by respectively 5.15% and 7.54%. We also show that it is possible to save up to 8.91% energy by increasing the tolerance on trip times and headways.

Timetable Model

In this section, we define a generic mathematical model of metro timetabling which will be used to:

1. define and classify different metro energy optimization problems from the literature (Section 3), 2. define the related MILP models, 3. define the instant power demand function (Section 4 and 5), 4. and define our greedy heuristic algorithm (Section 6).

For the sake of simplicity, we assume that the metro line is composed of N stations, S = {S 1 , ..., S N } and the metro timetable is represented as a sequence of M trips, with M being pair, T = (T 1 , ..., T M 2 , T M 2 +1 , ..., T M ). We assume that the trips T 1 , ..., T M 2 cross all stations in the upstream sequence, (S 1 , ..., S N ), and the trips T M 2 +1 , ..., T M cross all stations in the downstream sequence, (S N , ..., S 1 ). We note S t (s) the s th station crossed by the trip T t , i.e. station S s if 1

≤ t ≤ M 2 or station S N -s+1 if M 2 + 1 ≤ t ≤ M .
We assume that physical trains and crew have been allocated to trips beforehand so that the model of timetable described here does not consider metros depot movements, crew rostering, nor turnaround manoeuvres.

Variables

The variables are the dates of the departure and arrival times in stations for each trip, the dates of the starting of the braking phase and of the ending of the acceleration phase at each station for each trip. We consider that the time domain I = {0, 1, ..., I EN D } is discrete, with a precision of 1 second.

• d t,s ∈ I is the departure time of the trip T t , with 1 ≤ t ≤ M , at station S t (s), with 1 ≤ s ≤ N -1,
• a t,s ∈ I is the arrival time of the trip T t , with 1 ≤ t ≤ M , at station S t (s), with 2 ≤ s ≤ N ,

• d acc t,s ∈ I is the ending time of the acceleration phase of the trip T t , with 1 ≤ t ≤ M , leaving station S t (s), with 1 ≤ s ≤ N -1,

• a brk t,s ∈ I is the beginning time of the braking phase of the trip T t , with 1 ≤ t ≤ M , arriving at station S t (s), with 2 ≤ s ≤ N .

Auxiliary Variables

To simplify the equations, it is useful to introduce the following auxiliary variables as functions of the variables. These variables have the dimension of a time in I:

• the interstation time for a given trip,

int t,s = a t,s+1 -d t,s 1 ≤ t ≤ M, 1 ≤ s ≤ N -1,
• the dwell time, or stopping time, in a station for a given trip,

dwe t,s = d t,s -a t,s 1 ≤ t ≤ M, 2 ≤ s ≤ N -1,
• the total trip time,

trt t = int t,1 + N -1 k=2 (dwe t,k + int t,k ) = a t,N -d t,1 1 ≤ t ≤ M,
• the time interval, called headway, between two successive trips running in the same direction in a given station,

hdw t,s = d t,s -d t-1,s t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N.
• The duration of the acceleration phase for a given trip to leave a given station, acc t,s = d acc t,s -d t,s 1 ≤ t ≤ M, 1 ≤ s ≤ N -1, It is worth noticing that every acceleration phase occurs right after a dwell time. Shifting the starting time of an acceleration phase d t,s without modifying its length acc t,s is thus equivalent to modify the length of the adjacent dwell time dwe t,s .

• the duration of the braking phase for a given trip before arriving at a given stations, brk t,s = a t,s -a brk t,s

1 ≤ t ≤ M, 2 ≤ s ≤ N.

Constraints

Each trip t must leave its departure terminal within some bounds set by the metro company according to the quality of service provided to passengers,

d t,1 ≤ d t,1 ≤ d t,1 1 ≤ t ≤ M. (1)
The interstation times are bound according to the different speeds -e.g. economical, nominal or full throttle -a metro can take,

int t,s ≤ int t,s ≤ int t,s 1 ≤ t ≤ M, 1 ≤ s ≤ N -1, (2)
The dwell times are bound according to a minimum quality of service for the passengers,

dwe t,s ≤ dwe t,s ≤ dwe t,s 1 ≤ t ≤ M, 2 ≤ s ≤ N -1, (3)
The acceleration and braking phases are also bound according to the possible interstation times for a given station :

acc t,s ≤ acc t,s ≤ acc t,s 1 ≤ t ≤ M, 1 ≤ s ≤ N -1, ( 4
) brk t,s ≤ brk t,s ≤ brk t,s 1 ≤ t ≤ M, 2 ≤ s ≤ N, (5)
The global trip time is also bound to ensure the feasibility of the timetable and the quality of service to the passengers,

trt t ≤ trt t ≤ trt t 1 ≤ t ≤ M, (6)
Finally, the headways are bound according to security requirements and quality of service,

hdw t,s ≤ hdw t,s ≤ hdw t,s t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N. (7)
It is worth noticing that the headways must be checked at each station. The tolerances on headways differ according to the hour of the day. During peak hours, the headways are small and the tolerance for stretching them is tight. During off peak hours, the tolerances are higher and more important modifications of the timetable are possible, leading to potentially greater energy savings.

Objective Function

All the constraints of the timetable model shown up to now are pretty trivial bound linear constraints. The difficulty of energy optimization problems lies in the definition of the objective function, and more precisely of the instant power demand P i at each time.

Given an appropriate definition of P i , the objective can then be to minimize either the global energy consumption G T T of all trips,

G T T = I EN D i=0 P i , (8)
or the maximum power peak

P P T T = max i∈I P i . (9)
It is worth noticing that the objective function P P T T is widely used in the literature [START_REF] Albrecht | Reducing power peaks and energy consumption in rail transit systems by simultaneous metro running time control[END_REF][START_REF] Kim | A mathematical approach for reducing the maximum traction energy: The case of Korean MRT trains[END_REF][START_REF] Chen | Optimization of an MRT train schedule: Reducing maximum traction power by using genetic algorithms[END_REF][START_REF] Sansó | Trains scheduling desynchronization and power peak optimization in a subway system[END_REF][START_REF] Kim | A model and approaches for synchronized energy saving in timetabling[END_REF], but that a more realistic function would be the number of times that the power exceeds a certain threshold P M AX ,

CP T T = card(i | P i > P M AX ). ( 10
)
This function is more in accordance with the system of fines paid to the electricity provider when the quota is exceeded. The most accurate evaluations of P i are obtained by electrical simulators. In Section 4 we describe an instant power model, from which we derive two approximations: a non-linear approximation based on a power flow used in our algorithm and the linear approximation used in [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF]. Before that, the timetable model can already be used to classify several variants of the problem studied in the literature.

Classification and Complexity of Metro Timetabling Energy Optimization Problems

Problem Classification

As mentioned in the introduction, research is active for optimizing energy in the field of railways and some attempts to classify the studied problems have been made. Xun et al. [START_REF] Xun | Coordinated train control in a fully automatic operation system for reducing energy consumption[END_REF] proposed a classification of the methods used to solve the problems and Li et al. [START_REF] Li | An energy-efficient scheduling and speed control approach for metro rail operations[END_REF] listed without apparent classification some papers in the literature, detailing the decision variables or the algorithms used. To date, there is no formal classification of metro energy optimization problems. We propose a classification based on the previous timetable by a triple (G/P P/CP, dep/dwe/int, sim/nonlin/lin) denoting the choice of the objective function (G, P P or CP ), the decision variables (departure times, dwell times, interstation times or any combination Problem Equations References (P P, dwe, sim)

3, 6, 7, 4, 5, 9, Chen et al. [START_REF] Chen | Optimization of an MRT train schedule: Reducing maximum traction power by using genetic algorithms[END_REF] (P P, dep -dwe, sim) 1, 3, 6, 7, 4, 5, 9 Sanso et al. [START_REF] Sansó | Trains scheduling desynchronization and power peak optimization in a subway system[END_REF] (P P, dwe -int, sim) 2, 3, 6, 7, 4, 5, 9 Albrecht et al. [ of them) and the instant power demand evaluation (by an electrical simulator, a non-linear approximation or a linear approximation). Table 1 classifies different problems studied in the literature using these triples.

In this paper, we shall focus on the problems (G, dwe, lin) and (G, dwe, nonlin), that is, we focus on the problems of modifying solely the dwell times in order to minimize the global energy consumption of a metro line, evaluated using either a linear or a non-linear approximation. In this class of problems, the departure times dep t , the interstation times int t,s , and the braking and acceleration phases brk t,s and acc t,s are given. The modification of the dwell times modifies the arrivals a t,s and departures d t,s in stations, and hence the auxiliary variables.

Complexity

First of all, one can remark that in absence of objective function, the timetable satisfiability problem is polynomial since all the equations of the timetable model are linear. Caprara in et al. showed in [START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF] that minimizing the deviation of a solution timetable comparing to an initial one to satisfy capacity or overtaking constraints is NP-hard. Serafini et al. showed in [START_REF] Serafini | A mathematical model for periodic scheduling problems[END_REF] that the problem of periodically scheduling trains with precedence constraints is NP-complete. Here we show the NP-hardness of several energy optimization timetable problems by polynomial reductions of SAT.

Theorem 1. The problem (G, dep, lin) is NP-hard. Corollary 1. The problem (G, dep, nonlin) is NP-hard. Corollary 2. The problem (G, dep -dwe -int, lin) is NP-hard.
The proofs are given in appendix. The idea of the proof of the first result is to construct a polynomial reduction of SAT to a particular (G, dep, lin) problem. A particular timetable is constructed such that all the acceleration phases can be synchronized with the periodic braking phases of only one special metro x 0 .

In this construction, it is possible to delay the departure time of each metro different from x 0 by one time unit only. This ensures that one unit of energy is saved at this time. A Boolean formula in conjunctive normal form (CNF) with c clauses is then satisfiable if and only if c units of energy can be saved in the timetable.

Example 1. Let us consider the SAT formula (x 2 ∨ x 3 ) ∧ (x 1 ∨ ¬x 2 ∨ x 3 ) ∧ (¬x 1 ∨ y).
The constructed timetable contains four trips x 0 , x 1 , x 2 and x 3 , and is divided in three periods during which metro x 0 has the same behaviour. For each braking phase of x 0 occurring at times i 1 , i 7 and i 13 , either x 1 , x 2 or x 3 can be synchronized to save one energy unit. At each time unit, a metro is either accelerating, consuming 1 energy unit, braking, producing 1 energy unit, coasting or dwelling, producing or consuming nothing. Each trip departure time can be delayed by one time unit. The formula is satisfiable if and only if for each braking phase of x 0 , one other metro can synchronize its acceleration phase with it. The timetable with the energy consumed or produced by each metro at each time is shown in Table 2.

T T i 0 i 1 i 2 i 3 i 4 i 5 i 6 i 7 i 8 i 9 i 10 i 11 i 12 i 13 i 14 i 15 i 16 i 17 x 0 1 -1 0 0 0 0 1 -1 0 0 0 0 1 -1 0 0 0 0 x 1 0 0 1 0 -1 0 1 0 0 0 -1 0 0 1 0 0 -1 0 x 2 1 0 0 0 -1 0 0 1 0 0 -1 0 1 0 0 0 -1 0 x 3 1 0 0 0 -1 0 1 0 0 0 -1 0 0 0 1 0 -1 0 Table 2: Timetable energy problem encoding the satisfiability of (x 2 ∨ x 3 ) ∧ (x 1 ∨ ¬x 2 ∨ x 3 ) ∧ (¬x 1 ∨ x 2 )
. Each cell represents the power produced (-1) or consumed (1) by the trips x i at each time. In this example, to save the 3 energy units produced by x 0 , and satisfy the SAT formula in this encoding, it suffices to delay x 3 by one time unit since x 1 and x 2 already save the 2 other energy units.

Instant Power Demand Model

The energy consumption of a metro line can be computed with an electrical simulator. The electrical circuit is composed of the N stations, the electrical substations (ESS) and the accelerating and braking trains at a given time.

The ESSs are electrical devices that convert the AC power provided by the electricity provider to DC power, directly usable by metros. We assume that they are not revertible, i.e. the regenerative braking can only be used by accelerating metros, and that they are directly connected to the stations. Part of the power provided by the ESSs is lost by Joule effects on the DC network.

Each station is connected to the next one by a resistive cable. The braking metros are connected to the station to which they arrive and are modelled as an ideal power source. The accelerating metros are connected to the station from which they depart and are modelled as an ideal power sink. Figure 1 depicts the circuit associated to a small network with five stations, three ESSs, one braking metro and one accelerating metro.

ESS1 + -V ESS R ESS ESS3 + -V ESS R ESS ESS5 + -V ESS R ESS R 1 R 2 R 3 R 4 P 2,i < 0 P 5,i > 0 S 1 S 2 S 3 S 4 S 5
Figure 1: Electrical circuit associated to a metro line with five stations at time i. It is composed of three electric sub-stations in S 1 , S 3 and S 5 , a braking metro arriving in S 2 and producing P 2,i , and an accelerating metro departing from S 5 and consuming P 5,i . The points in the network are linked by resistive cables.

Network Parameters

The electrical properties of the metro line are fixed and given by the following set of parameters valid at any time:

• V ESS ∈ R + is
the fixed voltage supplied by the ESSs to the line. Typical values are 750V [START_REF] Chen | Optimization of an MRT train schedule: Reducing maximum traction power by using genetic algorithms[END_REF] or 1500V [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF].

• R ESS ∈ R + is the value of the internal resistance of the ESSs.

• R s ∈ R + is the resistance of the electric cable of the metro network between the stations S s and S s+1 , for 1 ≤ s ≤ N -1. This is computed using the linear resistance equation R = ρ.l a , ρ being the resistivity of the third rail, a its section and l its length.

Timetable Parameters

P s,i ∈ R is the net power demand or production, involved by a metro, at time i for every station S s . This value is different from 0 only if there is effectively a metro either braking or accelerating near the station at this particular time. These values are modified according to the acceleration and braking phases of the metros at each time point. We have

P s,i      > 0 if metro accelerating near S s < 0 if metro braking near S s = 0 otherwise 1 ≤ s ≤ N, i ∈ I
The precise power values are supposed to be known beforehand by direct measurement on the metro motors. Figure 2 illustrates the net power demand and production of a metro during an interstation run. In a typical example depicted in Figure 2, the parameters P s,i are given by the power curve for each time point i.

Electric Variables

The energy transfers are modelled by the following variables. By convention, all voltages are positive and the currents can be negative.

• v s,i ∈ R + is the electric potential at station S s , for 1 ≤ s ≤ N , at time i ∈ I.

• i s,i ∈ R is the current flowing through the cable between stations S s and S s+1 , for 1 ≤ s ≤ N -1, at time i ∈ I.

• i ESS s,i
∈ R is the current flowing from the electrical substation connected to the station S s , for 1 ≤ s ≤ N , at time i ∈ I.

• i M ET s,i
∈ R is the current flowing between the network and the metro located in S s , for 1 ≤ s ≤ N , at time i ∈ I.

Constraints and Instant Power Demand Value

The following equations constrain the current and voltage at each metro station. Ohm's law gives

v s,i -v s+1,i = R s .i s,i 1 ≤ s ≤ N -1, i ∈ I (11) V ESS -v s,i = R ESS .i ESS s,i 1 ≤ s ≤ N, i ∈ I (12)
Kirchhoff's current law gives

i s,i + i s+1,i + i ESS s,i + i M ET s,i = 0 1 ≤ s ≤ N, i ∈ I ( 13
)
The satisfaction of the instant power gives rise to a non-linear equation:

P s,i = v s,i .i M ET s,i 1 ≤ s ≤ N, i ∈ I (14)
It is worth remarking that the instant power demand of a metro line P i is not equal to the sum of the net instant power demands of the metros N s=1 P s,i , but to the power supplied by the electrical substations over the line to fulfil the metro power demands:

P i = N s=1 V ESS . max(0, i ESS s,i ) i ∈ I
This value represents the net power consumption of the metro line. The currents i ESS s flowing through each ESS can be negative, i.e. can flow backwards from the line to the grid, but this negative power is not counted in the instant power consumption. Indeed, ESSs possess a rectifier that works as a diode and forces the current to flow only in one direction. In reality, if an electrical substation receives energy, typically when too many metros are braking and none is accelerating, the energy is absorbed by resistors that are placed on the line or on the metros brakes.

Instant Power Demand Approximations

To simplify the evaluation of the power demand, some contributions in the literature have made the choice of directly computing the power transfers between braking and accelerating metros instead of calculating voltages and intensities, and deducing the power demand from it [START_REF] Fournier | A greedy heuristic for optimizing metro regenerative energy usage[END_REF][START_REF] Kim | A mathematical approach for reducing the maximum traction energy: The case of Korean MRT trains[END_REF][START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF][START_REF] Kim | A model and approaches for synchronized energy saving in timetabling[END_REF][START_REF] Fournier | Method and system for timetable optimization utilizing consumption factors[END_REF]. We introduce in this section the notion of power flow network, which is a particular case of a generalized flow network, to model these power transfers. The idea is that setting the flow along the paths of the power flow network is an approximation of the power transfers between braking and accelerating metros, the flow arriving in the sink of the graph representing the power saved by regenerative braking reuse.

Both the electrical simulator and the power flow, model the fact that one accelerating metro can benefit from the regenerative energy of several braking metros and that one braking metro can feed several accelerating metros. The power flow approximation does not intend to reproduce exactly the electrical behaviour of the metro line but rather to give a fast evaluation of it for the optimization algorithm.

The distribution matrix ∆ is the matrix of distribution ratios ∆ s,s = (P -P s )/(P s ) ∈ [0, 1] between each stations S s and S s , with 1 ≤ s < s ≤ N , where P s is the positive net power demand of a metro accelerating in station S s , P s the negative net power production of a metro braking in station S s , and P the resulting power demand of the electrical network as computed by the electrical simulator. The distribution ratio represents the ratio of power a metro braking in station S s effectively transfers to a metro accelerating in station S s . As P s ≥ P , a value of 1 means that the energy is fully transferred and a value of 0 means that the energy is completely lost in resistors as a consequence of Joule effects.

Power Flow Approximation

A generalized flow network is a finite directed graph G(V, E) given with capacities c(u, v) on edges in E and a flow f (u, v) ≤ c(u, v). The graph is given in addition with positive gains γ(u, v) such that, if a a flow f (u, v) is entering at vertex v, then γ(u, v).f (u, v) is going out from v:

v∈V |(u,v)∈E γ(u, v)f (u, v) = v∈V |(v,u)∈E f (v, u), (15)
Two vertices in V are distinguished, the source t which can produce flow and the sink t which can absorb flow.

For each time i ∈ I, the power flow network is the generalized flow network defined by

V = {t, t } ∪ V brk i ∪ V acc i , V brk i = {s, 1 ≤ s ≤ N |P s,i < 0}, i ∈ I, V acc i = {s , 1 ≤ s ≤ N |P s ,i > 0}, i ∈ I, E = {(t, s)} ∪ {(s, s )} ∪ {(s , t )} s ∈ V brk i , s ∈ V acc i and c(t, s) i = P s,i γ(t, s) i = 1 c(s, s ) i = +∞ γ(s, s ) i = ∆ s,s c(s , t ) i = P s ,i γ(s , t ) i = 1
The instant power demand approximation P i can then be defined as:

P i = s ∈V acc i [c(s , t ) i -f (s , t ) i ] (16)
Figure 3 shows an example of a power flow network modelling power transfers between 3 metros braking in stations S 1 , S 3 and S 5 , and 3 metros accelerating in stations S 2 , S 6 and S 7 . The flows are attenuated between braking and accelerating metros according to their distribution ratio ∆ s,s . The flows are not bound and can be set freely between braking and accelerating metros. However, they are bound by the powers produced or demanded at the source or to the sink. The flows effectively arriving to the sink represent the regenerative power that has been saved, while the sum of the capacities' edges arriving to the sink represent the total power demanded by the accelerating metros. Subtracting these two values gives an approximation of the instant power demand.

V brk i V acc i - P 1 , i , 1 -P 3,i ,1 - P 5 , i , 1 +∞, ∆ 1,2 + ∞ , ∆ 1 ,6 + ∞ , ∆ 1 , 7 + ∞ , ∆ 3 ,2 +∞, ∆ 3,6 + ∞ , ∆ 3 ,7 + ∞ , ∆ 5 , 2 + ∞ , ∆ 5 ,6
+∞, ∆ 5,7

P 2 , i , 1

P 6,i ,1 P 7 , i , 1

Power Flow Algorithm

According to the equation ( 16), the power demand equation of the power flow is equal to the capacities minus the flows directed to the sink. We have: c(s , t ) i = P s ,i , s ∈ V acc i and, according to the flow conservation equation ( 15)

f (s , t ) i = s∈V brk i γ(s, s ) i .f (s, s ) i with γ(s, s ) i = ∆ s,s , s ∈ V brk i , s ∈ V acc i
We can also reformulate the flow f (s, s ) i as the ratio of power transferred by the metro braking at station S s multiplied by the power P s,i :

f (s, s ) i = -x s,s ,i .P s,i , s ∈ V brk i , s ∈ V acc i
The power transfer ratio x s,s ,i ∈ [0, 1] is the ratio of the power P s,i , with s ∈ V brk i , transferred from the metro braking at station S s to the metro accelerating at station S s such that x s,s ,i = -f (s, s ) i /P s,i .

The following power flow algorithm computes the power transfer ratios by transferring the produced power of each braking metro in priority to the accelerating metro whose distribution ratio is maximum, until all the produced power is transferred. The algorithm returns the power transfer ratios x s,s ,i , either when all braking metros have transferred their power or when all accelerating metros have their demand fulfilled: Algorithm 1 Power transfer ratios computation at time i Require:

V acc i , V brk i , ∆ s,s , P s,i 1: Initialize vector x s,s ,i ← 0 2: while V brk i = ∅ do 3:
Choose randomly s ∈ V brk i 4: The instant power demand P i can now be reformulated with the non-linear equation:

P IN IT ← P s,i 5: while P s,i < 0 do 6: if V acc i = ∅ then 7: Choose s ∈ V acc i s.t
P i = s ∈V acc i P s ,i + s ∈V acc i s∈V brk i (P s,i .x s,s ,i .∆ s,s ), (17)

Linear Approximation for MILP

To avoid introducing too many variables in MILP, Peña proposed in [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF] a linear approximation of the power flow, which tends to maximize the overlapping times between acceleration and braking phases. The formulation of the power flow is simplified by authorizing only one single transfer between one braking metro and one accelerating metro. The objective function is then the sum, weighted by the distribution matrix, of the overlapping times of these transfers.

The MILP model introduces two new variables to describe the overlaps between metros :

• γ t,s,t ,s ∈ {0, 1} is a boolean variable equal to one if the trip T t , 1 ≤ t ≤ M braking in station S t (s), 2 ≤ s ≤ N , transfers its power to the trip

T t , 1 ≤ t ≤ M accelerating in S t (s ), 1 ≤ s ≤ N -1.
Otherwise it is equal to 0.

• O t,s,t ,s ∈ R is the overlapping time between the braking phase of the trip T t , 1 ≤ t ≤ M at station S t (s), 2 ≤ s ≤ N , and the acceleration phase of the trip

T t , 1 ≤ t ≤ M at station S t (s ), 1 ≤ s ≤ N -1. Constraint M t=1 N s=1 γ t,s,t ,s ≤ 1 1 ≤ t ≤ M, 1 ≤ s ≤ N (18)
ensures that each accelerating metro receives the totality of the power of at most one braking metro, and constraint

M t =1 N s =1 γ t,s,t ,s ≤ 1 1 ≤ t ≤ M, 1 ≤ s ≤ N (19)
ensures that each braking metro transfers the totality of its power to at most one accelerating metro, effectively modelling the fact that metros are only authorized to do one single pairing to transfer their power. Constraint

O t,s,t ,s ≤ brk t,s .γ t,s,t ,s 1 ≤ t < t ≤ M, 2 ≤ s ≤ N, 1 ≤ s ≤ N -1 (20)
ensures that the overlapping time of a braking phase with an acceleration phase cannot be bigger than the braking phase time brk

1 ≤ t < t ≤ M, 2 ≤ s ≤ N, 1 ≤ s ≤ N -1 (22)
ensure that the overlapping time is the minimum value of d acc v,u -a brk t,s and a t,sd v,u when γ t,s,v,u = 1. The member m(1-γ t,s,v,u ), with m a big enough number, ensures that O is never negative. The MILP objective function used in [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF] maximize

M t=1 N s=1 M t =1 N s =1 (O t,s,t ,s .∆ s,s ) (23)
is the sum of all overlapping times of the timetables, weighted by the distribution matrix. The weights tend to synchronize braking and accelerations that are close to each other. Concerning the size of the generated MILP instances, the model contains: The MILP model thus introduces a quadratic number of variables and constraints in the number of trips and stations. A pre-processing have been proposed in [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF] to remove irrelevant constraints and variables. In particular, the number of overlapping variables O t,s,t ,s to consider can be reduced by noticing that two trips T t and T t such that d t,1 > d t ,1 + trt t cannot overlap in any case. However, this pre-processing is not sufficient to handle real data timetables as we will show in Section 8.

• M 2 .N 2 boolean

Greedy Heuristics Algorithm for (G, dwe, nonlin)

The idea of our heuristic algorithm is to shift the acceleration phases to synchronize them with braking phases and to recompute the non-linear objective function to check for improvement. These local moves increase the power transfers between braking and accelerating metros, and possibly decrease the solution timetable global energy consumption.

Braking Phase Neighbourhood

The time neighbourhood of a braking phase N (a t,s ) is defined as the set of acceleration phases that can overlap this braking phase within the given timetable tolerances. This means that every acceleration phase that may start before the end and finish after the beginning of a given braking phase, belongs to the neighbourhood of the latter: N (a t,s ) = {d t ,s |(t = t ) ∧ (d acc t ,s + dwe t ,s > a brk t,s ) ∧ (d t ,s + dwe t ,s < a t,s )}

Acceleration Phase Shift Function

The acceleration phase shift function consists in modifying the departure time d t,s of a trip T t at a station S t (s) to make it correspond to the beginning of the neighbour braking phase a brk t ,s of a trip T t at station S t (s ), and by fixing acc t,s . The function considers the timetable constraints given by the bounds on dwell times, trip times and headways. If the function cannot shift the departure time of the acceleration to the beginning of the braking phase due to tolerance constraints, it will shift it to the closest time which respects these constraints: Shift(d Modifying d t,s consists thus in modifying the duration of dwe t,s , which will shift the future arrivals and departures a t,s and d t,s such that s > s.

Greedy Heuristics Optimization Algorithm

The algorithm comprehensively searches the acceleration phase shift that minimizes the objective function in the time neighbourhood of each braking phase,. All the braking phases are first sorted in chronological order and the acceleration phases are shifted for the earliest braking phases.

For each braking phase, the algorithm computes its time neighbourhood N (a t,s ), modifies each acceleration phase d t ,s by applying the shift function Shift(d t ,s , a brk t,s ) and checks whether this shift is decreasing the objective function. If it does, the best current objective function and the best shift are updated.

When all the acceleration phases have been shifted and evaluated, the algorithm shifts the acceleration phase that minimizes the most the objective function, or does nothing if none of them decreases the objective function. This monotonic behaviour ensures that the algorithm will converge after some iterations and does not worsen the initial timetable. Once an acceleration phase is shifted, it is removed from the pool of neighbour phases and cannot be shifted any more for another braking phase. The following pseudo-code summarizes the greedy heuristics optimization algorithm:

Algorithm 2 Greedy heuristics optimization algorithm Require: T T 1: Sort a t,s , 1 ≤ t ≤ M, 2 ≤ s ≤ N in chronological order 2:
for all a t,s do 3:

Compute initial objective function G IN IT

T T

4:

Initialize best objective function

G BEST T T = G IN IT T T

5:

Initialize best shift d BEST = 0 6:

Compute N (a t,s )

7:

for all d t ,s ∈ N (a t,s ) do 8:

d IN IT ← d t ,s 9:
Shift(d t ,s , a t,s )

10: end if 20: end for 21: return T T , G BEST T T

Compute G T T 11: if G T T < G BEST T T then 12: G BEST T T ← G T T 13: d BEST ← d t ,

Incremental Computation of the Objective Function

Equation [START_REF] Xun | Coordinated train control in a fully automatic operation system for reducing energy consumption[END_REF] shows that the instant power demand P i is a function of the net power demands and productions P s,i , the power transfer ratios x s,s ,i and the distribution matrix ∆ s,s . Since the distribution matrix is precomputed, only the power transfers are modified by the optimization process.

In our dedicated heuristics, the objective function is re-evaluated every time an acceleration phase is shifted. Since a shift consists in modifying the dwell time of a specific trip at a specific station, many time instants remain with the exact same net power demand and production along the timetable. Let P IN IT i be the instant power demand at time i before a dwell time shift and P IN IT s,i be the net power demands and productions at each station S s at time i. After the dwell time shift, only the instant power demands where the power demands and productions have changed need recomputing, as follows:

P i = P IN IT i if P s,i = P IN IT s,i ∀1 ≤ s ≤ N, needs recomputation otherwise
The incremental computation avoids recomputing known values, which increases the computation time of the algorithm by an order of magnitude as shown in Table 3 on six benchmark instances detailed in the following section. Table 3: Computation time in seconds of one run of the greedy heuristics without and with incremental computation of the objective function. The two implementations are compared on the six benchmark instances described in Section 7 representing typical off-peak (op) or peak (p) hours timetables. The lengths of the timetable instances are either 15 or 60 minutes and the number of decision variables #dt,s are given.

Iterative Optimization

After one run of the algorithm, the optimized timetable can be utilized as input for a second run starting from the new solution. The optimization algorithm can thus be executed either once or iteratively until the iterative algorithm stops improving the objective function.

Figure 4 shows the minimization of the objective function using iterative optimization. It shows that the first run largely improves the timetable and that the following iterations lead to further improvements.

Performance Results Compared to MILP and CMA-ES

In this section, the greedy heuristics algorithm is compared to MILP and CMA-ES on six small size benchmark instances on which the three methods give solutions. Both the greedy heuristics and CMA-ES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] are implemented in C++, and the MILP model is solved using CPLEX 12. The machine used for the experiments is a PC with an Intel Core i5 with 3GB of RAM.

For all the following results, a timeout of 1500 seconds was set. If not specified otherwise, the objective function is computed using the electrical simulator described in section 4. The greedy heuristics used the incremental computation of the objective function and iterative optimization.

Benchmark Instances

The six timetables have been drawn from real data and represent relevant portions of the timetable, i.e. peak (p) and off-peak (op) parts of size of 15 minutes and one hour. These instances contain the initial parameters of the timetable (d

IN IT t,1
, dwe IN IT , int IN IT and so on) as well as the tolerances, given by the customer, on which variables can be modified and by how much. The tolerances on dwell times, trip times and headways are equal for all six instances. For the sake of simplicity, the tolerance values are given relatively to the initial timetable instance and dwe t,s = -3 shall be read dwe t,s -dwe IN IT t,s = -3. The departure times, the interstation times, the braking and acceleration phases lengths are fixed :

dwe t,s = -3 1 ≤ t ≤ M, 1 ≤ s ≤ N dwe t,s = 9 1 ≤ t ≤ M, 1 ≤ s ≤ N trt t = -30 1 ≤ t ≤ M trt t = 30 1 ≤ t ≤ M hdw t,s = -30 t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N hdw t,s = 30 t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N 7.2.
Comparison with CMA-ES Evolution strategies are stochastic search algorithms that try to minimize an arbitrary objective function called fitness function. The covariance matrix adaptation evolution strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] applies to vectors of real-valued variables and arbitrary real-valued fitness functions. This algorithm is a multipoint method which at each iteration, samples the search space according to multivariate normal distributions, estimates its covariance matrix, determines a move to make in the most promising direction and updates the multivariate normal distributions for the variables. One important characteristic of CMA-ES compared to other meta-heuristics, is the limited number of parameters that need to be set, namely the initial standard deviation and the termination criteria. The other parameters are automatically adapted during the execution. The CMA-ES

In our experiments, we use the default value for the population size 4 + 3 log(#d t,s ). The optimization is stopped after 10 iterations without improvement of the objective function. The initial distribution has a default variance of (dwe t,s -dwe t,s )/7 for each trip at each station. A quadratic penalty function is added to the objective function for each variable out of its domain, to enforce the algorithm to search solutions within the given tolerances.

Table 4 shows the results of CMA-ES against our heuristics. Due to its stochastic behaviour, CMA-ES has been run 100 times for each instance. The table compiles the average computation time, and both the average and best value found for the objective function over the 100 runs. The results show that the greedy heuristics performs better than the best run of CMA-ES on four of the six benchmark instances. On op1, our heuristics is better than the average result of CMA-ES but not than its best result. Finally, CMA-ES is slightly better than the greedy heuristics in average only for p3. The better performance of the heuristic algorithm is partly due to the incremental computation of the objective function, which cannot be implemented in CMA-ES since all solutions are sampled randomly.

Comparison with MILP

As shown in Table 5, CPLEX is able to prove the optimality on the four smallest instances and outperforms the greedy heuristics on the linear objective function (Equation 23). However, when comparing both optimization methods on the objective function computed with the electrical simulator (Table 6), our greedy heuristics performs better on five of the six instances. This is due to the fact that the linear objective function is less accurate that the one used in the greedy heuristics (Equation 17). The main differences between these two objective functions are that the MILP model is only able to pair one braking with one acceleration (Equations 18, 19), when the power flow objective function is able to dispatch dynamically to different braking and accelerations over time as described in Section 5.2. Thus eventhough our algorithm does not prove optimality, it better approximates the real behaviour of the electricity flows and leads to better solutions.

Performance Results on Real Data

Our greedy heuristics has also been applied on a major city metro line comprising 16 stations for optimizing one full day timetable in two typical situations:

• a weekday timetable comprising 694 trips and 9585 dwell times,

• a Sunday timetable comprising 556 trips and 7679 dwell times.

Both the C++ implementation of CMA-ES and the MILP resolution by CPLEX have failed to tame problems of this size. The MILP model contains, after its pre-processing, 230908 constraints and 165760 variables, whose 52872 are binary, and runs out of memory on CPLEX on a PC with an Intel Core i5 with 3GB of RAM.. The size of this instance is to relate with the size of the problem handled in [START_REF] Peña | Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy[END_REF] which was containing only 17850 constraints and 13860 variables, whose 4780 were binary. For CMA-ES, it fails at computing the global energy consumption of the initial population within 30 minutes. On the other hand, our greedy heuristics is able to compute a solution in 20 minutes. The tolerances on the dwell times, trip times and headways have been first set such that there is no visible change in the quality of service for the passengers. Their relatively small values are as follows:

dwe t,s = -3 1 ≤ t ≤ M, 1 ≤ s ≤ N dwe t,s = 3 1 ≤ t ≤ M, 1 ≤ s ≤ N trt t = -15 1 ≤ t ≤ M trt t = 15 1 ≤ t ≤ M hdw t,s = -15 t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N hdw t,s = 15 t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N
For the Sunday timetable, the trip times and headways tolerances have been enlarged to 20 seconds as follows:

trt t = -20 1 ≤ t ≤ M trt t = 20 1 ≤ t ≤ M hdw t,s = -20 t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N hdw t,s = 20 t ∈ [2, M 2 ] ∪ [ M 2 + 2, M ], 1 ≤ s ≤ N
While the optimized timetable with regular tolerances is saving energy by 7.54%, the solution with increased tolerances can save up to 8.91%, increasing possibilities to synchronize phases better. CPLEX and the greedy heuristics on three full size timetables. CMA-ES and CPLEX did not manage to output a solution. The ratios represent the energy savings compared to the initial timetable energy consumption.

Figure 5: Weekday timetable between 6am and 1am: energy consumption by intervals of 30 minutes compared between the initial the timetable, in red, and timetable computed by the greedy heuristics in green.

Figures 5, 6 and 7 compare the initial and optimized timetable energy consumptions on three real data instances. For the weekday timetable, the two peak hour periods are clearly visible, from 8am to 11am and from 5pm to 9pm. It appears that more energy is saved during these hours. This is due to the fact that, according to the computed distribution matrix ∆ s,s , the energy transfers can be done only between metros that are very close from each other. Since the density of metros on the line is higher during peak hours, it is thus possible to save more energy during these times. The extrapolation of these savings shows that the metro company could save 3.65 GW.h of electrical energy per year.

Conclusion

In this paper, we have proposed a generic mathematical model for metro energy optimization problems and a dedicated heuristics for solving the global energy consumption optimization problem by the sole modification of dwell times. This model has led us to a simple classification by triples of several similar problems of the literature, and to prove their NP-hardness.

We have shown that our heuristic algorithm for the problem (G, dwe, nonlin) performs better than the classical optimization methods used in the literature. On six small size benchmark instances on which MILP and CMA-ES could be run, we have shown that our heuristic algorithm computes better solutions. In particular for the MILP formulation, the results computed by CPLEX were of lesser quality due to the linear approximation of the objective function whereas our heuristics uses a non-linear power flow approximation. It was also shown to perform better than the state-of-the-art metaheuristic CMA-ES on these instances thanks to the possibility of incrementally computing the objective function over iterations.

Furthermore, our dedicated heuristics was the only method able to solve a full day timetable of 7679 variables for the Sunday configuration and 9585 variables for the weekday configuration, decreasing the global energy consumption by 5.15% and 7.54%, and up to 8.91% by increasing the tolerances on the variables. These results show the applicability of this method in an industrial context.

Appendix A. NP-completeness proofs Theorem 1. The problem (G, dep, lin) is NP-hard. Proof 1. We show that there is a polynomial reduction of SAT to (G, dep, lin). Let X = {x 1 , ..., x m } be a set of m variables and ¬X = {¬x 1 , ..., ¬x m } be the set of their negations. Let = φ be a Boolean formula in conjunctive normal form:

φ = n i=1 c i
where c i are clauses of the form mi j=1 l i,j with l i,j ∈ X ∪ ¬X. Let us consider the discrete time domain I = {0, ..., 6n -1}. Let T T be the metro timetable composed of a sequence of m + 1 trips T = (x 0 , x 1 , ..., x m ) running in the same direction, such that each trip x t ∈ T crosses a sequence of unique stations of length n + 1, and that the trips are :

• The trip x 0 cannot be shifted and is crossing stations every 6 time units, departing from its first station at time i = 0:

d 0,0 = 0

-int 0,s = 1 1 ≤ s ≤ n -dwe 0,s = 5 2 ≤ s ≤ n
Reminding the auxiliary variables equations:

the interstation time for a given trip, int t,s = a t,s+1 -d t,s 0 ≤ t ≤ m, 1 ≤ s ≤ n, (A.1)

the dwell time, or stopping time, in a station for a given trip, dwe t,s = d t,s -a t,s 0 ≤ t ≤ m, 2 ≤ s ≤ n, (A.2) let us prove by induction that d 0,s = 6(s -1), 1 ≤ s ≤ n:

Basis: the statement holds for s = 1,

d 0,1 = 0 = 6(1 -1)
Inductive step: if d 0,s = 6(s -1), then d 0,s+1 = 6s.

We can write d 0,s+1 = d 0,s + a 0,s+1 -d 0,s + d 0,s+1 -a 0,s+1

According to (A.1) and (A.2) we have:

d 0,
) if 6l + 4 ≤ i ≤ 6l + 5 ∀ 0 ≤ l ≤ n -1
with k ∈ N, does not change the result of the objective function. Indeed, the timetable is constructed in such a way that m t=0 P t,i ≤ 0, for all 6l + 4 ≤ i ≤ 6l+5 and 0 ≤ l ≤ n-1. Thus, like for the previous timetable, P i = 0, for all 6l+ 4 ≤ i ≤ 6l + 5 and 0 ≤ l ≤ n -1. Finally the problem becomes (G, dep, nonlin) without changing its resolution, thus the problem (G, dep, nonlin) is NP-hard.

Corollary 2. The problem (G, dep -dwe -int, lin) is NP-hard. Proof 3. A problem is of class dwe if at least one dwell time can be modified. Likewise, a problem is of class int at least one interstation time can be modified. To get a problem (G, dep -dwe -int, lin), it suffices to add at the end of the timetable constructed in the proof of Theorem 1, one period of time where only x 0 is running, and where its last dwell time and interstation time can be modified.

Let us thus consider the metro timetable T T crossed by a sequence of m + 1 trips T = (x 0 , x 1 , ..., x m ) crossing n + 1 stations each which encodes a Boolean formula in CNF containing n clauses. To correctly encode the formula, the timetable should have a length of 6n -1 time units. Let us add at the end of the timetable 4 additional time units where only x 0 is crossing a new station S 0 (n + 2) such that:

Figure 2 :

 2 Figure 2: Net power demand and production curve as a function of real time for a metro accelerating and braking between two stations. The points on the curve represent the sampling over discrete time.

Figure 3 :

 3 Figure 3: Example of a generalized flow network. Flows are created by the source t and absorbed by the sink t . The left hand side of the graph contains the vertices corresponding to the braking metros, which are connected to the vertices corresponding to the accelerating metros on the right hand side. Each edge of the graph is characterized by a (capacity,gain) couple.

Figure 4 :

 4 Figure 4: Evolution of the global energy consumption over computation time in seconds on a sample timetable during the iterative optimization process. The first cross represents the global energy consumption of the original timetable and the following ones the global energy consumption by iterating the greedy heuristics.

Figure 6 :

 6 Figure6: Sunday timetable between 6am and 1am with a tolerance of 15 seconds on trips and headways: energy consumption by intervals of 30 minutes compared between the initial timetable, in red, and the timetable computed by the greedy heuristics in green.

Figure 7 :

 7 Figure7: Sunday timetable between 6am and 1am with a tolerance of 20 seconds on trips and headways: energy consumption by intervals of 30 minutes compared between the initial timetable, in red, and the timetable computed by the greedy heuristics in green.

  d 0,n+1 = 6n a 0,n+2 = 6n + 1 Now, the last dwell time of x 0 can be extended by one time unit, encoded in δ dwe = {0, 1}, dwe 0,n+1 = 5 + δ dwe , and the last interstation time of x 0 can be extended by one time unit, encoded in δ int = {0, 1}, int 0,n+1 = 1 + δ int By adding these two tolerances, the problem (G, dep, lin), which encodes the SAT formula and which is NP-hard, has become a problem (G, dep -dweint, lin). Thus the problem (G, dep -dwe -int, lin) is NP-hard.

Table 1 :

 1 Some metro timetabling energy optimization problems from the literature, classified by the problem triple they solve and the corresponding timetable equations.

	8]

  . s = arg max s ∈V acc

					i	(∆ s,s )
	8:		if -P s,i .∆ s,s > P s ,i then
	9:		x s,s ,i ← (P s ,i /∆ s,s )/P IN IT
	10:		P s,i ← P s,i + P s ,i /∆ s,s
	11:		V acc i	← V acc i	\{s }
	12:		else
	13:		x s,s ,i ← P s,i /P IN IT
	14:		P s ,i ← P s ,i + P s,i .∆ s,s
	15:		P s,i ← 0
	16:		end if
	17:	else
	18:		P s,i ← 0
	19:	end if
	20:	end while
	21:	V brk i	← V brk i	\{s}
	22: end while
	23: return vector x s,s ,i

  variables γ, M 2 .N 2 overlapping times variables O and M.N dwell times variables dwe,

	• 3M 2 .N 2 + 2M.N MILP constraints (Equations 18, 19, 20, 21, 22) and
	5M.N + 2M timetable constraints (Equations 1, 2, 3, 4, 5, 6, 7).

Table 4 :

 4 Compared performance in computation time (Time in seconds) and energy consumption (Value in kW.h) between the average and best values found over 100 runs of CMA-ES and the greedy heuristics on six benchmark instances. The instances op represent an off-peak hour timetable and the instances p represent a peak hour timetable, both of either 15 minutes or 60 minutes long.

					Initial		CMA-ES		Greedy Heuristics
	Inst. Length #dt,s	Value	Value Average	Best	Time Value	Time
	op1	15 min	127	2514	2401		2381	256	2394	45.6
	op2	15 min	129	2516	2402		2388	223	2376	38.0
	op3	60 min	449	9956	9724		9716	761	9556	648
	p1	15 min	173	3433	3300		3285	503	3262	178
	p2	15 min	186	3651	3516		3483	669	3442	291
	p3	60 min	670	13067	12696	12675 1030 12713	1500
	Instance Length #dt,s	Initial Greedy Heuristics Value Value	MILP Value Integrality gap
	op1	15 min 127 12.48 s	101.8 s	318.1 s	optimal
	op2	15 min 129 11.48 s	159.1 s	351.6 s	optimal
	op3	60 min 449 45.97 s	817.9 s	1637 s	10.62%
	p1		15 min 173 250.5 s	414.8 s	772.5 s	optimal
	p2		15 min 186 279.1 s	533.8 s	835.8 s	optimal
	p3		60 min 670 1019 s	1576 s	3003 s	20.44%

Table 5 :

 5 Compared performances over the MILP objective function (Value in s) between CPLEX and the greedy heuristics on six benchmark instances. The MILP solutions are given with their integrality gap, optimal standing for 0%. The instances op represent an off-peak hour timetable and the instances p represent a peak hour timetable, both of either 15 minutes or 60 minutes long.

Table 6 :

 6 Compared performances in computation time (Time in seconds) and energy consumption (Value in kW.h) between MILP and the greedy heuristics on six benchmark instances. The instances op represent an off-peak hour timetable and the instances p represent a peak hour timetable, both of either 15 minutes or 60 minutes long.

	Instance Length #dt,s	Initial Value	MILP Value Time Value Greedy Heuristics Time
	op1	15 min	127	2514	2427	0.72	2394	45.6
	op2	15 min	129	2516	2419	1.00	2376	38.0
	op3	60 min	449	9956	9579	1500	9556	648
	p1	15 min	173	3433	3281	460	3262	178
	p2	15 min	186	3651	3494	82.5	3442	291
	p3	60 min	670	13067 12711 1500 12713	1500

  Table 7 summarizes these results.

	Instance	Length #dt,s	Initial	CMA-ES MILP Greedy heuristics
	weekday	1 day	9585	218294	-	-	207052 (-5.15%)
	sunday15	1 day	7679	189953	-	-	175638 (-7.54%)
	sunday20	1 day	7679	189953	-	-	173036 (-8.91%)

Table 7 :

 7 Compared performances in terms of energy consumption, given in kW.h, of CMA-ES,

  s+1 = d 0,s + int t,s + dwe t,s+1 ⇔ d 0,s+1 = 6(s -1) + 1 + 5 ⇔ d 0,s+1 = 6s Corollary 1. The problem (G, dep, nonlin) is NP-hard. Proof 2. A problem is of class nonlin if it exists at least one equation which is non-linear. Modifying the equation (A.7) of the previous timetable into the non-linear equation

	P i = max(0,	m t=0 P t,i ) if 6l ≤ i ≤ 6l + 3
	P i = max(0,	m t=0 P t,i (2+1)k

We thus have:

3) a 0,s = d 0,s-1 + int 0,s-1 = 6(s -2) + 1, 2 ≤ s ≤ n + 1 (A.4)

• Furthermore, the only possible shift applicable on trips x t in this timetable is a delay of their departure time by 1 time unit, denoted δ t ∈ {0, 1}. The trips x t with 1 ≤ t ≤ m are constructed according to the clauses c i of φ as follows:

Like for x 0 , we prove by induction that

Basis: the statement holds for s = 1,

We can write

According to (A.1) and (A.2) we have:

We thus have:

Let the instant power demand function be

where P t,i ∈ R is the power demand or production of the trip x t at time i. The objective function is G T T . Let the trip acceleration and braking phases last one time unit only. Each trip will demand one unit of power when departing from a station, and will produce one unit of power when arriving to a station. The rest of the time, each trip will not demand or produce any power. For all trips x t , the instant power demand or production P t,i at time i can be written as follows:

Now let us prove that the global energy consumption is not modified by the powers produced by trips {x 1 , ..., x m } since they brake it when no trip is accelerating. According to equations (A.6) and (A.8) we have:

Conversely and according to equations (A.3), (A.5) and (A.8), we have for all trips x t ∈ T : P t,i = 1 if i = 6(s -1) + (0 or 1 or 2 or 3), 0 ≤ t ≤ m, 1 ≤ s ≤ n (A.9) Thus, there is no time where the braking of any trip in {x 1 , ..., x m } can be synchronized with the acceleration of any trip in T :

On the other hand, the braking phases of the trip x 0 can be absorbed by the acceleration phases of the other trips, optimizing the objective function. According to equations (A.4) and (A.8) we have:

Also according to equation (A.9), P t,s = 1 if i = d t,s . To synchronize the acceleration of the trip x t at station s with one braking of the trip x 0 we need:

In other terms, the timetable is constructed such that for each trip x t and for each station s we have:

• If the variable x t is in the clause c s , then the acceleration phase of x t at station S t (s) is synchronized with the braking phase of x 0 at station S 0 (s + 1) if and only if δ t = 1. Thus setting δ t = 1 is equivalent to say that the variable x t is true, satisfying all clauses c s containing it.

• If the variable ¬x t is in the clause c s , then the acceleration phase of x t at station S t (s) is synchronized with the braking phase of x 0 at station S 0 (s + 1) if and only if δ t = 0. Thus setting δ t = 0 is equivalent to say that the variable ¬x t is true, satisfying all clauses c s containing it.

• If neither the variable x t nor the variable ¬x t are in the clause c s , then the acceleration phase of x t at station S t (s) cannot be synchronized with any of the braking phases of x 0 .

Every time one braking phase of x 0 is synchronized with the acceleration phase of one trip x t , one unit of power is saved. Given the period of time I and the structure of the timetable, there is n power units that an optimal synchronization could save. Consequently, considering the decision problem of saving n power units, there is a timetable and a set of δ t with 1 ≤ t ≤ m that save n units of power if and only if the set of δ corresponds to a valuation that satisfies φ ; the n synchronizations corresponding to the n satisfied clauses.


Xiaorong Yang

Dan Wang

Lihua Jiang

Jian Ma

Hui Xie

Study on Cloud Service Mode of Digital Libraries based on Sharing Alliance

Keywords: Cloud services, Digital library, Service mode, Sharing alliance

As a new service model, cloud services have become the growth point of the service innovation of digital libraries. By analyzing and summarizing the cloud service mode of the digital libraries at home and abroad, this paper presents three-layer cloud service model of digital libraries which includes data layer, automated management layer and service layer from bottom to top. Based on this model, this paper puts forward the union mechanism established by some libraries which joint provide intelligence research service, knowledge push service, reference service and so on.

Introduction

With the development of cloud computing technology, the cloud service application research becomes a hot issue. The operating mode and technical architecture of cloud service are closely related to Google and Amazon. Google put forward a complete set of distributed parallel cluster infrastructure according to the characteristics of the large scale network data. And Google provides a series of SaaS for individual users and enterprises including the Google search engine, Google maps, photos and videos sharing, social networks, Gmail, Google calendar, Google Apps Market place and so on. In order to make full use of idle IT infrastructure, Amazon provided storage servers, bandwidth and CPU resource to the third party by Amazon Web Services in 2002. In 2006, the network service platform based on cloud computing were provided to enterprises so that they do not need to undertake the underlying work such as hardware maintenance and can focus on their own business.

Information industry is always an important field of information technology application. To apply cloud computing technology, a library can make full use of its resources and technical strength. The service efficiency can be enhanced obviously by innovating service mode.

Cloud Service Status of Domestic and Abroad Libraries

Cloud Service Demand Analysis of Digital Libraries

Service Situation and Problems of Digital Libraries

Due to the monopoly service of domestic companies such as CNKI, WANFANG and VIP in the Chinese literature resources, most of the libraries are not engaged in the electronic processing and network services of Chinese journals, dissertations, conference papers and so on. For readers, only buyers can access to these commercial database of Chinese literature resources such as CNKI. Many researchers have no effective way to obtain the required Chinese document information resources. In the processing and services of foreign literature resources, only CALIS and NSTL established the co-construction and sharing mechanism. It is common in provincial libraries and county libraries that resources are not complete and the service effects are unsatisfactory.

Besides, the functions of cataloguing system and data processing system are unsatisfactory in some libraries because of limited capital and technical strength.

When processing literature data, data format conversion between different software results to inefficient data processing work.

Cloud Service Demand of Digital Libraries

Cloud service pattern will bring the library more flexible, rich and diversified service So the libraries at the same field can jointly build sharing alliance. Member libraries can carry out digital construction and service of Chinese-foreign resources to provide integrated retrieval and services for readers. The scale advantage of the library service will be helpful to ascend influence of the member libraries. Besides the interlibrary information sharing based on the same platform will be conducive to develop knowledge services and applications.

Cloud Service Mode of Digital Libraries based on Sharing Alliance

In the future libraries services will focus on SaaS based on data resources. Several libraries establish the co-construction and sharing system to expand the scope and service mode of information resources. Cloud computing technology is adopted to establish SaaS cloud platform to which member libraries can submit the metadata of their Chinese-foreign resources. Large joint library can provide single sign-on and one-stop service with huge amounts of resources for readers.

Cloud Service Model of Digital Library

The cloud service hierarchy should be built based on the standards system and supported by the operation maintenance and security system. The hierarchy includes data layer, automated management system layer and service layer from bottom to top.

The relationship and structure is shown in figure 1.

Figure 1 Cloud Service Model of Digital Library

The data layer is the underlying structure of the three layer service model of cloud services. This layer consists of data resources built by member libraries of shared alliance such as catalog database, abstract database, citation database, institutional repository and so on. The automated management system layer takes the system management based on data layer . The layer includes a set of middleware service such as user registration, service scheduling, unified authentication, search engine, load management, unified authorization, log statistics and service register.

The service layer is a set of services based on data layer and the automated management system layer. This layer is a software application model which provides cloud application service through the Internet transmission. The service can be divided into three categories: information service, information retrieval and data processing and management. The cloud services of data processing and management provide a cloud platform based on the unified metadata standard for member libraries.

The information cloud service based on co-build data helps member libraries to jointly carry out intelligence research, knowledge push, the original offer, check by services and so on. The information retrieval cloud service provides readers with the retrieval tools to quickly find the resources needed(Chen Chen et .al, 2012).

Cloud Service Mode of Digital Libraries based on Sharing Alliance

The cloud services based on sharing alliance consist of information service, information retrieval and data processing and management. Member libraries can use these cloud services in whole or in part according to need.

Information Service

(1)Intelligence Research Service. Intelligence research software based on co-built basic database can be provided to member libraries. By using the software, member libraries can understand their own competitiveness and analysis hot research field so as to provide the basis for future decision-making.

(2)Knowledge Push Service. Knowledge service software can be provided to member libraries. Thus the readers' access historical data can be analyzed in order to get the interest preference. By being selected, optimized and integrated, information resources can be further developed and pushed to readers. By intelligently filtering and recommending the information, service mode transformation can be promoted from passive service to active service, which the quality of information service can be significantly improved.

(3) Reference Consulting Service. Reference consulting software can be used online.

By using the software, librarians can query advisory forms and provide joint consultation among multiple libraries.

(4)Original Offer Service. Original offer software can be used online. Librarians can manage user account, process orders, transfer documents and manage settlement by using the software. Thus querying, acquisition and transmission services of all kinds of literature resources can be provided for researchers and graduate students.

Information Retrieval

Unified retrieval function based on sharing alliance resources can be provided. By using the resources search engine, readers can retrieve all abstract information of books, journals, conference proceedings, reports, dissertations, standards and patent documents in cloud service platform to complete one-stop retrieval of printed and digital resources. The retrieval system can be used remotely or be embedded into the web systems of member libraries.

Data Processing and Management

(1)Joint Cataloging Services. Joint cataloging software can be provided for member libraries. By using the software, librarians can process metadata cataloging database of Chinese-foreign journals. Readers can log in the joint cataloging system of cloud services platform to search the resources of all member libraries.

(2)Joint Data Processing Services. Joint data processing software can be provided for member libraries. Librarians can adopt the unified standard to process data so as to provide data support for the retrieval service of Chinese-foreign information resources.

Conclusion

The wide application of cloud computing technology has brought a revolution to library information service mode. Cloud service has become a growth point of the digital library service innovation. As SaaS, cloud service based on library sharing alliance can significantly improve the efficiency of access to information and the utilization of information resources. In the future libraries should focus on the cloud service mode and play advantages guided by demand.

1. 1

 1 Cloud Service Application Research Status of Abroad Libraries Online Computer Library Center (OCLC) is a profitable and the largest library cooperation organization all over the world. In 2009, OCLC launched the WorldCat Local based on cloud computing technology which provides a full set of cloud computing information management service for libraries such as WorldCat.org, WorldCatLocal, Questionpoint, CONTENTdm and so on. As a library OCLC firstly applied the cloud computing technology to provide services for users. National Digital Information Infrastructure and Preservation Program (NDIIPP) is a joint action of digital resources preservation initiated by the American Congress. In 2009, NDIIPP and Duraspaces company started the experimental project-DuraCloud together. DuraCloud provided permanent access for digital resources of back issues in the New York Public Library so as to detect the effect of cloud computing technology applied in permanent access of digital resources. This research played an important role in testing cloud computing technology applied in digital Library(Yingjun Lu et .al, 2012).

1. 2

 2 Cloud Service Application Research Status of Domestic Libraries China Academic Library&Information System (CALIS) is a literature guarantee system for colleges and universities launched by The Ministry of Education in 1998. CALIS joined nearly 800 university libraries to provide cloud service including E read, current contents of western journals(CCC), unified authentication system (UAS), interlibrary loan (ILL), distributed collaborative virtual reference system (CVRS) and the unified data exchange system (UES) (Wenqing Wang et .al, 2009). Jilin Province set up the union library which Included 50 public libraries, University libraries and scientific libraries. The union library started to plan and construct a cloud computing service platform in 2010. The platform provided IaaS for member libraries, PaaS for social users and integrated the resources of member libraries to provide unified services as SaaS. National Science and Technology Digital Library (NSTL) is a virtual institution of science& technology literature information service. The institution involves in science, engineering, agriculture and medicine. And it includes 9 member libraries and 50 service stations all over the country. At present the co-construction and sharing among several member libraries is SaaS. And the standard open interface service oriented group users is DaaS. In 2012, the project was started to study the cloud service mode and build the cloud service platform(Xiaodong Qiao et .al, 2010).

  changes. For libraries, to innovating the service mode based on a cloud service alliance has become inevitable. Rich data resources are the basic guarantee of libraries' literature service. Because of their expensive price, foreign electronic resources of every library are not complete. Many libraries have their own distinctive Chinese resources. With consciousness ascend of independent intellectual property rights, some organizations have begun to set up their own institutional repository (IR) not to sell their knowledge assets such as Ph.D. thesis to commercial data companies cheaply. As a representative, Chinese Academy of Sciences has established 77 institutional repository of subordinate institute. Island services are not only unfavorable for the readers but also not conducive to the development of libraries themselves. The service system of digital library based on cloud computing technology can provide efficient platform to process data and service for libraries which lack technical strength.

( 5 )

 5 Check by Services. Check by services software can be used online. Librarians can receive and handle the entrusted orders. Query and delivery services facing all kinds of literature resources can be provided for domestic and abroad researchers (Xiaobo Xiao et .al, 2012).

Acknowledgements

The work is supported by project of national science and technology library "Research on cloud service model and construction of cloud computing services platform for NSTL (2012XM02)", and the special fund project for Basic Science Research Business Fee, AII "Research on user behavior analysis of website of CAAS"(No. 2014-J-007).