Efficient Train Operation via Shrinking Horizon Parametrized Predictive Control

6th 6th IFAC IFAC Conference Conference on on Nonlinear Nonlinear Model Model Predictive Predictive Control Control Madison, WI, USA, USA, August August 19-22, 2018 2018 6th IFAC Conference on Nonlinear Model Predictive Control Madison, WI, 19-22, Available online at www.sciencedirect.com 6th IFAC Conference on Nonlinear Model Predictive Control Madison, WI, USA, August 19-22, 2018 Madison, WI, USA, August 19-22, 2018

ScienceDirect

IFAC PapersOnLine 51-20 (2018) 203–208

Efficient Efficient Train Train Operation Operation via via Shrinking Shrinking Horizon Horizon EfficientParametrized Train Operation via Shrinking Horizon Predictive Control EfficientParametrized Train Operation via Shrinking Predictive ControlHorizon Parametrized Predictive Control Parametrized Predictive Control ∗,∗∗ Hafsa Farooqi ∗∗ , Lorenzo Fagiano ∗∗ , Patrizio Colaneri ∗,∗∗

Hafsa Farooqi ∗ , Lorenzo Fagiano ∗ , Patrizio Colaneri ∗,∗∗ Hafsa Farooqi , Lorenzo Fagiano , Patrizio Colaneri ∗ , Lorenzo Fagiano ∗ , Patrizio Colaneri ∗,∗∗ ∗ Dipartimento Hafsa Farooqi Elettronica, ∗ Dipartimento di Elettronica, Informazione Informazione ee Bioingegneria, Bioingegneria, Politecnico Politecnico di di ∗ Dipartimento di di Elettronica, Informazione e [email protected], Bioingegneria, Politecnico di Milano, Milano, Italy. E-mail addresses: Milano, Milano, Italy. E-mail addresses: [email protected], ∗ Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Milano, Italy. E-mail addresses: [email protected], [email protected], [email protected] [email protected], [email protected] ∗∗ Milano, Italy. E-mail addresses: [email protected], [email protected], [email protected] ∗∗Milano, IEIIT CNR, Politecnico di Milano, - CNR, Politecnico di Milano, Milano, Milano, Italy Italy ∗∗ IEIIT [email protected], [email protected] IEIIT - CNR, Politecnico di Milano, Milano, Italy ∗∗ IEIIT - CNR, Politecnico di Milano, Milano, Italy Abstract: trains is is considered. considered. Abstract: The The problem problem of of driver driver assistance assistance for for the the energy-efficient energy-efficient operation operation of of trains Abstract: of driver assistance forces for theapplied energy-efficient operation trains isspeed considered. The to control traction/braking to while satisfying limits The goal goal is isThe to problem control the the traction/braking forces applied to the the train, train, while of satisfying speed limits Abstract: The problem of driver assistance forarrival theapplied energy-efficient operation of trains is considered. The goal is to control the traction/braking forces to the train, while satisfying speed limits and reaching the next station at the prescribed time. Moreover, the control input has to belong and reaching the next station at the prescribed arrival time. Moreover, the control input has to belong The goal is tothe control the traction/braking forces applied to the train,the while satisfying speed limits and reaching next station at the prescribed arrival time. Moreover, control input has to belong to a discrete set of values and/or operating modes, which a human driver has to implement. to a discrete set of values and/or operating modes, which a human driver has to implement. A A and reaching the next station at the prescribed arrival time. Moreover, the control input has to belong to a discrete setpredictive of valuescontrol and/or(MPC) operating modes, which a human driver has to horizon implement. A nonlinear model approach is proposed, featuring aa shrinking and an nonlinear model predictive control (MPC) approach is proposed, featuring shrinking horizon and an to a discrete setpredictive of values and/or operating modes, which a human driver has to horizon implement. A nonlinear model control (MPC) approach isoptimization proposed, featuring a shrinking and an input-parametrization strategy to aa continuous problem. Theoretical convergence input-parametrization strategy to retain retain continuous optimization problem. Theoretical convergence nonlinear model predictive control (MPC) approach isoptimization proposed, featuring a shrinking horizon and an input-parametrization strategy to retain a continuous problem. Theoretical convergence guarantees are derived, and the approach is tested in realistic simulations. guarantees are derived, and the approach is tested in realistic simulations. input-parametrization retain ais continuous optimization problem. Theoretical convergence guarantees are derived,strategy and the to approach tested in realistic simulations. © 2018, IFAC Federation Control) simulations. Hosting by Elsevier Ltd. All rights reserved. guarantees are (International derived, and the approachof isAutomatic tested in realistic Keywords: Keywords: Model Model Predictive Predictive Control, Control, Input Input Parametrization, Parametrization, Nonlinear Nonlinear control control systems, systems, Train Train control control Keywords: Model Predictive Control, Input Parametrization, Nonlinear control systems, Train control Keywords: Model Predictive Control, Input Parametrization, Nonlinear control systems, Train control 1. cursive 1. INTRODUCTION INTRODUCTION cursive feasibility. feasibility. To To the the best best of of our our knowledge, knowledge, the the combined combined 1. INTRODUCTION cursive To the best of our knowledge, the and combined presence of dynamics, shrinking horizon, inputpresencefeasibility. of nonlinear nonlinear dynamics, shrinking horizon, and input1. most INTRODUCTION cursive feasibility. To the best ofinour knowledge, the and combined Railway is by far the efficient means of transportation presence of nonlinear dynamics, shrinking horizon, inputparametrization strategy is new the literature, and represents strategy is new in the literature, and represents Railway is by far the most efficient means of transportation parametrization presence of nonlinear dynamics, shrinking horizon, and inputRailway ispoint by far the most efficient means of transportation from the of view of energy consumption and therefore parametrization strategy is new in the literature, and represents the main theoretical contribution of this paper, in addition main theoretical contribution of this paper, in addition to to the the from the ispoint of view of energy consumption and therefore the Railway by far the most efficient means of transportation parametrization strategy is new in the literature, and represents from the point of view of energy consumption and therefore aa strategic the main theoretical contribution of this paper, in addition to the content. Realistic simulation results show strategic sector sector in in today’s today’s society. society. The The earliest earliest existing existing works works application-specific application-specific content. Realistic simulation results show the point view of energy consumption and therefore the main theoretical contribution of this paper, in addition to the afrom strategic sector today’s society. The earliest existing on efficient trainofin control include developing strategies byworks con- the application-specific content. Realistic simulation results show effectiveness of the approach in the considered application. the effectiveness of the approach in the considered application. on efficient train control include developing strategies by conaonstrategic sector in today’s society. Thestate earliest existing application-specific content. Realistic simulation results show efficient train control include developing strategies byworks considering the problem as a bounded variable one (see the effectiveness of the approach in the considered application. sidering thetrain problem asinclude a bounded state variable one (see the effectiveness of the approach in the considered application. on efficient control developing strategies by considering the problem as a bounded variable one (see Ichikawa (1968)) and development development of state numerical optimization 2. Ichikawa (1968)) and of numerical optimization 2. PROBLEM PROBLEM STATEMENT STATEMENT AND AND ABSTRACTION ABSTRACTION sidering thesuch problem asones a boundedof state variable(1981) one (see Ichikawa techniques the by 2. PROBLEM STATEMENT AND ABSTRACTION techniques(1968)) such as asand thedevelopment ones presented presentednumerical by Milroy Milroyoptimization (1981) and and Ichikawa (1968)) of numerical optimization techniques such thedevelopment ones by Milroy (1981) and Consider Strobel and Hornasand (1973). In presented these approaches, approaches, linear expres2. PROBLEM STATEMENT ANDa ABSTRACTION Strobel and Horn (1973). In these linear expresConsider an an electric electric train train controlled controlled by by a digital digital control control unit. unit. In In techniques such the ones by non Milroy (1981) and Strobelofand Hornas(1973). In presented these approaches, linear expressions resistance, constant constraints, dependency of Consider an electric train controlled by a digital control unit. In this work, we consider space as the independent variable, while sions of resistance, constant constraints, non dependency of this work, we consider space as the independent variable, while Strobel and Horn (1973). In these approaches, linear expressions of forces resistance, constant constraints, non dependency of time external on the the train position, position, and constant constant slopes and and Consider an train controlled by a Thus, digitalwe control unit. In this work, weelectric consider space as thestates. independent variable, while will be one of the system’s denote with external forces on train and slopes time will be one of the system’s states. Thus, we denote with sions of forces resistance, constant constraints, nonalldependency of this work, we consider space as the independent variable, while external onwere the train position, constant slopes and track considered. In reality, factors will be one of the system’s states. Thus, we denote with ktime ∈ Z the discrete space variable, and with D the sampling track curvatures curvatures were considered. In and reality, all these these factors ∈ Z the discrete space variable, and with Dss the sampling external forces onwere the train position, constant slopes and ktime track curvatures In and reality, all factors significantly affect theconsidered. energy consumption. A these method that be one of space the system’s states. Thus,the denote with the sampling kdistance, ∈ Zwill theso variable, andalong with Dwe that distance at s track significantly affect the energy consumption. A method that distance, sodiscrete that the the actual actual distance along the track at each each track curvatures were considered. In reality, all these factors significantlyvariable affect slopes the energy consumption. A underground method that kdistance, considered has been developed for the sampling ∈ Z the discrete space variable, and with D so that isthe actual distance along makes thes track at each sampling instant equal to the control ss .. This considered variable slopes has been developed for underground instant equal to kD kD This choice choice the at control significantly affect the energy consumption. A (1971). methodMore that sampling considered slopes has been for underground trains with variable short station distance bydeveloped Maksimov distance, so that is actual distance along makes the track each sampling instant isthe equal to kD . This respect. choice makes the control problem formulation easier in We denote with s some trains with short station distance by Maksimov (1971). More problem formulation easier in some respect. We denote with considered variable slopesdistance hashave beenbeen developed for underground trains with short techniques station Maksimov More sampling T recently, several developed considering instant is equal to kD . This choice makes the control s problem formulation easier in some respect. We denote with T (k), x (k)] the state of the train, where x is its travel x(k) = [x recently, several techniques have by been developed(1971). considering 1 (k), x22 (k)] the state of the train, where x11 is its travel x(k) = [x 1 trains with short techniques station Maksimov More T easier recently, several have by been developed considering T all these factors. A recentdistance review of the existing (1971). techniques is time problem formulation in some respect. Wex1denote with (k), x (k)] the state of the train, where is its travel x(k) = [x T and x the train speed (· denotes the matrix transpose 1 2 2 all these factors. A recent review of the existing techniques is and x2 the train speed (·T denotes the matrix transpose T the recently, several techniques have been considering all these factors. A recent et review of thedeveloped existing techniques is time provided by Scheepmaker al. (2017). (k), x (k)] state of the train, where x is its travel x(k) = [x 2 train 1 transpose time and1 xand speed (· 1] denotes the matrix operator), with u(k) ∈ aa normalized traction force, 2 the provided by Scheepmaker al. (2017). operator), and with u(k) ∈ [−1, [−1, normalized traction force, T 1]denotes all theseenergy factors. A recent et review of the existing techniques is time provided by Scheepmaker et al. (2017). For efficient of and xand train speed (·the theapplicable matrix transpose 2=the operator), u(k) ∈ [−1, 1]maximum a normalized traction force, u(k) 11with corresponds to traction For the the energy efficient operation operation of railways, railways, Model Model Predictive Predictive where where u(k) = corresponds to the maximum applicable traction provided by Scheepmaker etapproach, al. (2017). For the energy efficient operation of railways, Model Predictive Control (MPC) is a suitable thanks to its capability to operator), and with u(k)maximum ∈ [−1, 1]maximum a normalized traction force, u(k) = 1 corresponds to the applicable traction and u(k) = −1 to the braking. The input u Control (MPC) is a suitable approach, thanks to its capability to where u(k) = −1 to the maximum braking. The input u is is the the For the energy operation of railways, Model Predictive Control (MPC) is a suitable approach, thanks to its capability to and deal with state efficient and input constraints and economic objectives. where u(k)= =−1 1 corresponds to thetrain maximum applicable traction and u(k) to the maximum braking. The input u is one the available control variable. The has to move from deal with state and input constraints and economic objectives. available control variable. The train has to move from one Control (MPC) is a suitable approach, thanks to its capability to deal with state and input constraints economic MPC has been been already applied in the theand context of this thisobjectives. industrial and u(k) = −1x1to= the maximum braking. Themove input u1 is the available control variable. The train has to from one station at time 0 and reach the next one at time x = MPC has already applied in context of industrial at time x1 = 0 and reach the next one at time x1 = xx ff ,, deal with state and input constraints economic MPC has been applied theand context of thisobjectives. industrial station application, see e.g. et (2014). available control variable. The train has to move from one station at time x = 0 and reach the next one at time x = covering the corresponding distance s . For a given pair application, seealready e.g. Aradi Aradi et al. al. in (2014). 1 f, covering the corresponding distance s ff . For a given 1pair xof of MPC has been applied the context of this industrial station application, seealready e.g. Aradi et al. in (2014). at time x = 0 and reach the next one at time x = x 1 1 f, covering the corresponding distance s . For a given pair of initial and final stations, the track features (slopes, curvature) f In this paper, paper,seewe wee.g. present study on the the application application of of aa new new initial and final stations, the track features (slopes, curvature) application, Aradi aaetstudy al. (2014). In this present on covering the corresponding distance s . For a given pair of f initial and final stations, the track features (slopes, curvature) are known in advance. Thus, in nominal conditions (i.e. with In this paper, we present a study on the application of a new MPC known in advance. Thus, in nominal conditions (i.e. with MPC approach approach to to the the efficient efficient operation operation of of trains. trains. In In our our rere- are and finaladvance. stations,Thus, the track features (slopes, curvature) are in nominal with rated values like mass the In this paper, wetopresent a study onwith the application ofour a new MPC efficient operation of trains. re- initial search, carried in rail transport, ratedknown valuesinof of the the train train parameters, parameters, like its itsconditions mass and and (i.e. the specspecsearch,approach carried out out the in collaboration collaboration with Alstom Alstom rail In transport, are known in advance. Thus, in nominal conditions (i.e. with rated values of the train parameters, like its mass andaccording the specifications of the powertrain and braking systems), MPC approach to(usually the efficient operation of trains. our research, carried out in collaboration with Alstom rail In transport, the control input the traction/braking force applied to ifications of the powertrain and braking systems), according the control input (usually the traction/braking force applied to rated values of the train parameters, like its mass and the specifications of the powertrain and braking systems), according to Newton’s laws and using the forward Euler discretization search, carried out in collaboration with Alstom rail transport, (usuallytothebelong traction/braking applied to to Newton’s laws and using the forward Euler discretization the control train) is isinput constrained to aa set set of offorce discrete values the train) constrained belong to discrete values of the and braking systems), according to Newton’s lawspowertrain andof theof Euleraccurate discretization method, the motion aa reasonably model theoperating control (usuallyto traction/braking applied to ifications train) isinput constrained tothe belong to a set result offorce discrete values or modes, would naturally in method, the equations equations ofusing motion offorward reasonably accurate model or operating modes, which which would naturally result in aa mixedmixedto Newton’s laws and using the forward Euler discretization method, the equations of motion of a reasonably accurate model of this this system system read: read: the train) is constrained totowould belong to at a each set result ofsampling discrete values of or operating modes, which naturally in a instant. mixedinteger nonlinear program be solved integer nonlinear program to be solved at each sampling instant. method, the equations of motion of a reasonably accurate model of this system read: D or operating modes, which would naturally result in a mixedinteger beformulation solved at each instant. To copenonlinear with this this program problem,tothe the wesampling propose features features + xx1 (k) + Dsss xxthis 11 (k To cope with problem, formulation we propose (k system + 1) 1) = = read: of 1 (k) + xx22D(k)  (k) integer nonlinear program to be solved at each sampling instant. s To cope this problem,strategy the formulation weapropose features an input with parametrization that yields yields continuous op F (x(k),u(k))−FB (x(k),u(k))−FR (k,x(k))   x (k + 1) = x1 (k) + x2D(k) an input parametrization strategy that apropose continuous op(x(k),u(k))−FRR (k,x(k))  2s s FT T xxx112 (k (k + 1) = x (k) + D To cope with this problem, the formulation we features T (x(k),u(k))−FB B 2 + 1) = x (k) + an input parametrization strategy that yields a continuous optimization program. Moreover, we adopt a shrinking horizon, Mx (k) 2 (k + 1) = x12 (k) + D 2 (k) F (x(k),u(k))−F (x(k),u(k))−F (k,x(k)) x2s(k) timization program. Moreover, we adopt a shrinking horizon, T B R Mx 2  2 x2 (k + 1) = x2 (k) + Ds  F (x(k),u(k))−FMx an input parametrization strategy that yields a We continuous optimization program. we adopt a shrinking horizon, rather than the most mostMoreover, common receding one. named the 2 (k) B (x(k),u(k))−F R (k,x(k)) (1) rather than the common receding one. We named the (1) x2 (k + 1) = x2 (k) + Ds T timization program. Moreover, we adopt a shrinking horizon, Mx (k)the traction force, (1) rather than the most commonhorizon receding one. We named the where resulting approach Shrinking Parametrized Predictive where M M is is the the total total mass mass of of the the train, train, F FTT 2is is the traction force, resulting approach Shrinking horizon Parametrized Predictive (1) rather than the most common receding one. We named the F where M is the total mass of the train, F is the traction force, resulting approach Shrinking horizon Parametrized Predictive is the braking force, and F the resistive force. Functions Control (SPPC). We present both a nominal approach and two T B R Control (SPPC). WeShrinking present both a nominal approach and two F B is the braking force, and FR the resistive force. Functions where M is the total mass of the train, F is the traction force, resultingones, approach horizon Parametrized Predictive T F is the braking force, and F the resistive force. Functions Control (SPPC). We present both a nominal approach and two (x, u), F (x, u) are nonlinear and they depend on the specific relaxed where state constraints are softened to retain rerelaxed ones, where state constraints are softened to retain re- FBTT (x, u), FBB (x, u) are nonlinear Rand they depend on the specific is u), the Fbraking force, and FRand thethey resistive force. Functions Control (SPPC). We state present both a nominal approach and two nonlinear depend on the specific relaxed ones, where constraints are softened to retain re- FBT (x, B (x, u) are and they depend on the specific relaxed ones, where state constraintsFederation are softened to retainControl) re- F T (x, u), B (x, u) are 2405-8963 © 2018, IFAC (International of Automatic Hosting byFElsevier Ltd.nonlinear All rights reserved.

Copyright © 2018 235 Copyright © under 2018 IFAC IFAC 235 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 235 10.1016/j.ifacol.2018.11.014 Copyright © 2018 IFAC 235

2018 IFAC NMPC 204 Madison, WI, USA, August 19-22, 2018

Hafsa Farooqi et al. / IFAC PapersOnLine 51-20 (2018) 203–208

train and track profile. They include, for example, look-up tables that link the traction and braking forces to the train speed and to the control input value. These functions are derived either experimentally or from complex models of the train and its traction and braking systems. In our research, these are provided by the business unit at Alstom. More details on these functions are omitted for confidentiality reasons. The resistive force FR (k, x) is also nonlinear, and it is the sum of a first term Rv (x2 ), accounting for resistance due to the velocity, and a second term Rg (k), accounting for the effects of slopes and track curvature: FR (k, x) = Rv (x2 ) + Rg (k) Rv (x2 ) = A + Bx2 +Cx22 (2)   D Rg (k) = Ms g tan(α(k)) + rc (k) where the parameters A, B,C, D are specific to the considered train, Ms is the static mass of the train, i.e. the mass calculated without taking into account the effective inertia of the rotating components, rc (k) and α(k) are, respectively, the track curvature and slope at position k, and g is the gravity acceleration. Besides the prescribed arrival time x f and position s f , there are additional state constraints that must be satisfied. These pertain to the limit on the maximum allowed velocity, x2 (k), which depends on the position k, since a different velocity limit is imposed for safety by the regulating authority according to the track features at eachposition.  Overall, by defining the terminal . space step as k f = s f /Ds (where · denotes the flooring operation to the closest integer), the state constraints read: x(0) x(k f ) x2 (k) x2 (k)

= [0, 0]T = [x f , 0]T ≥ 0, k = 0, . . . , k f ≤ x2 (k), k = 0, . . . , k f

(3)

The control objective is to maximize the energy efficiency of the train while satisfying the constraints above. To translate this goal in mathematical terms, different possible cost functions can be considered. In our case, we consider the discretized integral of the absolute value of the traction power over space (with a constant scaling factor D−1 s ): kf

J=

∑ |FT (x(k), u(k))| .

(4)

k=0

This choice tends to produce controllers that minimize the traction energy injected into the system. The braking energy is not penalized, since in our case there is no restriction to the use of the braking system. As already pointed out, the input variable is also constrained in the interval u ∈ [−1, 1]. However, in a driver assistance scenario which is the main focus of this work, the control algorithm is developed to assist a human driver with a suggested value of the input handle, in which case only a smaller set of possible values can be delivered by the controller, in order to facilitate the human-machine interaction. In particular, in this scenario the input constraints are further tightened according to four possible operating modes prescribed by our industrial partner: • Acceleration: in this mode, the input takes the maximum value, i.e. u = 1. • Coasting: this mode implies that the traction is zero, i.e u = 0. • Cruising: in this mode, the train engages a cruise control system that keeps a constant speed, i.e. u is computed 236

by an inner control loop in such a way that FT = FR for positive slopes and FB = FR for negative slopes. • Braking: in this mode the maximum braking force is used, i.e. u = −1.

Finally, a further feature of this application is a relatively small sampling space Ds with respect to the imposed overall space horizon s f , resulting in a rather large number of sampling periodsin the interval [0, s f ]. 2.1 Problem Abstraction The control problem described above can be cast in a rather standard form: kf

min u

∑ (x(k), u(k))

(5a)

k=0

subject to x(k + 1) = f (x(k), u(k))

(5b)

u(k) ∈ U, k = 0, . . . , k f − 1

(5c)

x(k) ∈ X, k = 1, . . . , k f x(0) = x0

(5d) (5e)

(5f) x(k f ) ∈ X f where x ∈ X ⊂ Rn is the system state, x0 is the initial condition, u ∈ U ⊂ Rm is the input, f (x, u) : X × U → X is a known nonlinear mapping representing the system dynamics, and l(x, u) : X × U → R is a stage cost function defined by the designer according to the control objective. The symbol u = {u(0), . . . , u(k f − 1)} ∈ Rm k f represents the sequence of current and future control moves to be applied to the plant. The sets X ⊂ X and U ⊂ U represent the state and input constraints (including the discrete set of allowed inputs or input modes as described above), and the set X f ⊂ X the terminal state constraints, which include a terminal equality constraint as a special case. We recall that a continuous function a : R+ → R+ is a K function (a ∈ K ) if it is strictly increasing and a(0) = 0. Throughout this paper, we consider the following continuity assumption on the system model f . Assumption 1. The function f enjoys the following continuity properties:    f (x1 , u) − f (x2 , u) ≤ ax x1 − x2  , ∀x1 , x2 ∈ X, u ∈ U  f (x, u1 ) − f (x, u2 ) ≤ au u1 − u2  , ∀u1 , u2 ∈ U, x ∈ X (6) where ax , au ∈ K .  In (6) and in the remainder of this paper, any vector norm  ·  can be considered. Assumption (1) is reasonable in most real-world applications, and it holds in the railway application considered here. 3. NOMINAL SPPC APPROACH To solve problem (5), we resort to Nonlinear Model Predictive Control (NMPC) with shrinking horizon and inputparametrization strategy to reduce the computational burden required by the feedback controller as well as to enforce the discrete constraints on the input while retaining a continuous optimization program, instead of the mixed-integer one that would result from direct optimization of the input.

2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018

Hafsa Farooqi et al. / IFAC PapersOnLine 51-20 (2018) 203–208

3.1 Parametrization Setup We adopt a parametrization of the control input that allows us to naturally incorporate the presence of the discrete (switched) driving modes described in Section 2. The main idea is to first split the track (i.e. the whole prediction horizon from 0 to k f ) into sectors. Then, in each sector, we pre-define a switching sequence of the four driving modes, and we optimize over the switching (space) instants of the sequence. In this way, we retain a continuous vector of optimization variables, thus improving the computational efficiency, while still providing the controller with enough degrees-of-freedom to optimize the predicted system behavior. Specifically, we consider a number sn ∈ N of sectors, each one with length Γi , such that sn

∑ Γi = s f .

(7)

i=1

The choice of sectors is carried out by considering characteristics such as the presence of constant velocity limits and the resistance force due to the slopes and track curvature Rg (k), which are known in advance, see Fig. 1 for an example. Regard15

As an example, values of (s(i,1) , s(i,2) , s(i,3) , s(i,4) ) equal to (0, 0, Γi , 0) correspond to the train coasting throughout sector i, and so on. Remark 1. The presented parametrization provides the controller with enough degrees of freedom to choose a single mode or multiple modes of operation in each sector, in order to compensate the presence of uncertainty and to adapt to the track features and constraints. Additionally, the constraint on prescribed arrival position needs to be satisfied: this is done implicitly by imposing (9) for each sector, which implies (7). We are now in position to introduce the optimal control problem to be solved at each time step in our shrinking horizon approach. We named the resulting control approach Shrinking horizon Parametrized Predictive Control (SPPC), presented next. 3.2 Shrinking horizon Parametrized Predictive Control (SPPC) We start by defining the optimization variables available at each step k. The set of indexes identifying the current and future sectors, from the current position kDs until the end of the track s f , is given by: {i : i(k) + 1 ≤ i ≤ sn },

Velocity Limits

where 10

. i(k) = 5 Rg

0

Γ

-5 0

1

Γ

2

Γ4

Γ3

Γ5

5

Γ6

Γ7

10

15

kD s(m)

Fig. 1. Example of sector choice for a track with s f = 15 and sn = 7. Possible sectors based on similar characteristics such as velocity limits and Rg values are depicted. ing the switching sequence (“driving style”) adopted in each sector, denoted with udi , we choose the following: udi = { u(i,1) , u(i,2) , u(i,3) , u(i,4) } 0, −1 }, ∀i ∈ [1, . . . , sn ] ={ 1, uCR ,

4

  

i

(10)

max i s.t. ∑ Γi < kDs ,

if Γ1 < kDs

0,

otherwise

i≥1

i=1

Then, the number N(k) of free variables to be computed corresponds to the number of remaining sectors, equal to (sn − i(k)), times the number of modes in each sector, i.e. 4 in our case. Therefore, we have N(k) = 4(sn − i(k)). We denote the vector of optimization variables with . vN(k) = {δ s(i(k)+1,1) , . . . , δ s(i(k)+1,4) , . . . , δ s(sn ,4) }T ∈ RN(k) . (11) Let us indicate with u( j|k) the input vector at each space sample k + j predicted at step k. Considering the parametrization described in Section 3.1, we can define the function g(vN(k) , j, k) that links, at each step k and for each j = 0, . . . , k f − k − 1, the optimization variables vN(k) with the predicted input u( j|k): . (12) u( j|k) = g(vN(k) , j, k) = u(iˆ( j,k),(ˆ v , j,k)) , N(k)

(8)

where u(i,) is the input issued in the −th phase of the i−th sector (with  ∈ {1, 2, 3, 4}) and uCR identifies the cruising mode (see Section 2), i.e. where the actual input u(k) is computed by a cruise control system in order to maintain a constant speed. The chosen switching sequence is such that in each sector, theoretically, the controller can choose to have a traction phase, a cruising one, a coasting one, and finally a braking phase. Thus, the maximum number of operating modes that can be set in the whole optimization horizon is equal to 4 sn . As anticipated, our optimization variables will be the switching instants, or more precisely the duration of each phase within each switching sequence. We denote these quantities with δ s(i,) , such that: 0 ≤ δ s(i,) ≤ Γi

∑ δ s(i,) = Γi .

205

(9)

=1

237

where (compare with (8)): . iˆ( j, k) =

min i

i=1,...,sn

s.t. i

∑ Γi ≥ (k + j)Ds

i=1

and

. ˆ N(k) , j, k) = (v min  =1,...,4

s.t.

iˆ( j,k)



i=1

=1

∑ Γi + ∑ δ s(iˆ( j,k),) ≥ (k + j)Ds .

Note that the evaluation of (12) is very efficient, since it just amounts to finding, for each space sample k + j, the indexes of the corresponding sector and phase and then to apply the corresponding pre-defined driving mode from (8). Finally, let us denote with x( j|k), j = 0, ..., k f − k the state vectors predicted

2018 IFAC NMPC 206 Madison, WI, USA, August 19-22, 2018

Hafsa Farooqi et al. / IFAC PapersOnLine 51-20 (2018) 203–208

at space sample k + j starting from the one at step k. At each step k ∈ [0, k f − 1], we formulate the following FHOCP: min

vN(k)

k f −k

∑ (x( j|k), u( j|k))

(13a)

j=0

To model uncertainty and disturbances we consider an additive term d(k) acting on the input vector, i.e.: u(k) ˜ = u(k) + d(k)

subject to u( j|k) = g(v, vN(k) , j, k), j = 0, . . . , k f − k − 1

(13b)

u( j|k) ∈ U, j = 0, . . . , k f − k − 1

(13d)

x( j + 1|k) = f (x( j|k), u( j|k)), j = 0, . . . , k f − k − 1 (13c) x( j|k) ∈ X, j = 1, . . . , k f − k

(13e)

x(0|k) = x(k)

(13f)

A(k)vN(k) ≤ b(k)

(13g)

x(k f − k|k) ∈ X f

4. RELAXED SPPC APPROACH: ALGORITHM AND PROPERTIES

(13h)

(17)

where u(k) ˜ is the disturbance-corrupted input provided to the plant. This model represents well all cases where plant uncertainty and exogenous disturbances can be translated into an effect similar to the control input (the so-called matched uncertainty). For example, in our application, equation (17) can describe uncertainty in the train mass, drivetrain specs, track slope and misapplication of the computed input by the human operator in a driver assistance scenario (see Section 2). We consider the following assumption on d: Assumption 2. The disturbance term d belongs to a compact set D ⊂ Rm such that:

Where the matrices A(k), b(k) in (13g) are built to enforce the ∗ (18) d ≤ d, ∀d ∈ D constraints (9). We denote with vN(k) = {δ s∗(i(k)+1,1) , . . . , δ s∗(sn ,4) }T a solution (in general only locally optimal) of (13). Moreover,  we denote with x∗ (k) and u∗ (k) the corresponding predicted where d ∈ (0, +∞). sequences of state and input vectors: This assumption holds in many practical cases and in the considered train application as well. We indicate the perturbed x∗ (k) = {x∗ (0|k), . . . , x∗ (k f − k|k)} (14a) state trajectory due to the presence of d as: u∗ (k) = {u∗ (0|k), . . . , u∗ (k f − 1 − k|k)}

(14b)

where

x∗ (0|k) = x(k) x∗ ( j + 1|k) = f (x∗ ( j|k), u∗ ( j|k)) u

∗

∗ ( j|k) = g(vN(k) ,

j, k)

(14c) (14d)

The SPPC strategy is obtained by recursively solving (13), as described by the following pseudo-algorithm. Algorithm 1. Nominal SPPC strategy (1) At sampling instant k, measure or estimate the state x(k) ∗ and solve the FHOCP (13). Let vN(k) be the computed solution; (2) Apply to the plant the first element of the control sequence ∗ , 0, k); u∗ (k), i.e. u(k) = u∗ (0|k) = g(vN(k) (3) Repeat the procedure from 1) at the next sampling period.  Algorithm 1 defines the following feedback control law: u(k) = µ(x(k)) = u∗ (0|k),

(15)

and the resulting model of the closed-loop system is: x(k + 1) = f (x(k), µ(x(k))

x(k ˜ + 1) = f (x(k), ˜ u(k)), ˜ k = 0, . . . , k f

where x(0) ˜ = x(0). To retain recursive feasibility in presence of the disturbance, we soften the constraints in the FHOCP. Since, in our railway application the terminal state constraint is the most important one from the viewpoint of system performance, to be more specific we restrict our analysis to the terminal state constraint in (5f), i.e. x(k f ) ∈ X f without loss of generality. The other constraints (velocity limits) are always enforced for safety by modulating traction or by braking. Let us denote the distance between a point x and a set X as: ∆(x, X) = min x − y. y∈X

Then, we want to derive a modified SPPC strategy with softened terminal state constraint (to ensure recursive feasibility) that guarantees a property of the following form in closed loop: ∆(x(k ˜ f ), X f ) ≤ β (d), β ∈ K .

(19)

That is, the distance between the terminal state and the terminal constraint is bounded by a value that decreases strictly to zero as d → 0. In order to obtain this property, we propose a relaxed SPPC approach using a two-step constraint softening procedure, described next.

(16)

The recursive feasibility of (13) and convergence of the state of (16) to the terminal set is established by construction of the nominal SPPC strategy. In the next section, we introduce a model of uncertainty, whose form is motivated again by the application described in Section 2, and a possible variation of Algorithm 1 to deal with it, along with its guaranteed convergence properties. We term this variation the “relaxed” approach, since it involves the use of suitable soft (i.e. relaxed) constraints to guarantee recursive feasibility. 238

4.1 Two-step relaxed SPPC strategy At each step k we consider a strategy consisting of two optimization problems to be solved in sequence: a) we compute the best (i.e. smallest) achievable distance between the terminal state and the terminal set, starting from the current perturbed state x(k): ˜

2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018

Hafsa Farooqi et al. / IFAC PapersOnLine 51-20 (2018) 203–208



γ = arg min γ vN(k) ,γ

subject to u( j|k) = g(vN(k) , j, k), j = 0, . . . , k f − k − 1 x( j + 1|k) = f (x( j|k), u( j|k)), j = 0, . . . , k f − k − 1 u( j|k) ∈ U, j = 0, . . . , k f − k − 1 x( j|k) ∈ X, j = 1, . . . , k f − k x(0|k) = x(k) ˜ A(k)vN(k) ≤ b(k) ∆(x(k f − k|k), X f ) ≤ γ (20) b) we optimize the input sequence using the original cost function, and softening the terminal constraint by γ: k f −k

˜ j|k), u( j|k)) min ∑ (x(

vN(k) j=0

subject to u( j|k) = g(vN(k) , j, k), j = 0, . . . , k f − k − 1 x( j + 1|k) = f (x( j|k), u( j|k)), j = 0, . . . , k f − k − 1 u( j|k) ∈ U, j = 0, . . . , k f − k − 1 x( j|k) ∈ X, j = 1, . . . , k f − k x(0|k) = x(k) ˜ A(k)vN(k) ≤ b(k) ∆(x(k f − k|k), X f ) ≤ γ (21) By construction, both problems are always feasible (with the caveat that state constraints are considered to be always feasible, as discussed above, otherwise the softening shall be applied r to these constraints as well). We denote with vN(k) , xr (k) and r u (k) the optimized sequences of decision variables, state and inputs resulting from the solution of (21). The sequences xr (k) r and ur (k) are computed from vN(k) and x(k) ˜ as reported in (14). The resulting relaxed SPPC strategy is implemented by the following pseudo-algorithm. Algorithm 2. Two-stage relaxed SPPC strategy (1) At sampling instant k, measure or estimate the state x(k) ˜ and solve in sequence the optimization problems (20)r (21). Let vN(k) be the computed solution; (2) Apply to the plant the control vector u(k) = ur (0|k) = r g(vN(k) , 0, k); (3) Repeat the procedure from (1) at the next sampling period.  Algorithm 2 defines the following feedback control law: ˜ = ur (0|k), u(k) = µ r (x(k) and the resulting closed-loop dynamics are given by: ˜ + d(k)). x(k ˜ + 1) = f (x(k), ˜ µ r (x(k))

(22) (23)

The next result shows that the closed-loop system (23) enjoys a uniformly bounded accuracy property of the form (19), provided that the nominal SPPC problem (13) is feasible at k = 0. Theorem 1. Let Assumptions 1 and 2 hold and let the FHOCP (13) be feasible at step k = 0. Then, the terminal state x(k ˜ f ) of system (23) enjoys property (19) with ∆(x(k ˜ f ), X f ) ≤ β (d) =

k f −1

∑ βk f −k−1 (d)

(24)

k=0

where β0 (d) = au (d) βk (d) = au (d) + ax (βk−1 (d)), k = 1, . . . , k f − 1

207

(25) 239

The proof is by induction and is very similar to the proof of Theorem 1 in Farooqi et al. (2018) and hence has been removed for brevity. Theorem 1 indicates that the worst-case distance between the terminal state and the terminal set is bounded by a value which is zero for d = 0 and increases strictly with the disturbance bound. In the considered railway application this means that, for example, the worst-case accuracy degradation in reaching the terminal station on time due to the changing mass of the train which is a function of the passenger load, is proportional to the largest difference between the nominal mass and its change due to the varying passenger load. This result provides a theoretical justification to the proposed twostep relaxed SPPC approach. The bound (24) is conservative, since it essentially results from the accumulation of worstcase perturbations induced by the disturbance on the open-loop trajectories computed at each k. As we show in our simulation results, in practice the resulting closed-loop performance are usually very close to those of the nominal case, thanks to recursive optimization in the feedback control loop. 4.2 Multi-objective relaxed SPPC strategy As an alternative to the two-step approach described above, one can also consider a multi-objective minimization: min

vN(k) ,β

k f −k

˜ j|k), u( j|k)) + ωγ ∑ (x(

j=0

subject to u( j|k) = g(vN(k) , j, k), j = 0, . . . , k f − k − 1 x( j + 1|k) = f (x( ˜ j|k), u( j|k)), j = 0, . . . , k f − k − 1 (26) u( j|k) ∈ U, j = 0, . . . , k f − k − 1 x( j|k) ∈ X, j = 1, . . . , k f − k x(0|k) = x(k) ˜ A(k)vN(k) ≤ b(k) ∆(x(k f − k|k), X f ) ≤ γ where ω is a positive weight on the scalar γ. Problem (26) can be solved in Algorithm (2) in place of problems (20)(21). In this case, the advantage is that a trade-off between constraint relaxation and performance can be set by tuning ω. Regarding the guaranteed bounds on constraint violation, with arguments similar to those employed in Fagiano and Teel (2013) one can show that, at each k ∈ [0, k f − 1], for any ε > 0 there exists a finite value of ω such that the distance between the terminal state and the terminal set is smaller than γk f −k−1 (d) + ε. Thus, with large-enough ω, one can recover the behavior obtained with the two-step relaxed SPPC approach. The theoretical derivation is omitted for the sake of brevity, as it is a rather minor extension of the results of Fagiano and Teel (2013). 5. SIMULATION RESULTS We tested the proposed strategies in realistic simulations with real train data provided by Alstom for a section of the Amsterdam metro rail, in particular the track between Rokin and Central Station. The parametric values of the train used in the controller are (see (1)-(2)) M = 142403 kg, Ms = 131403 kg, A = 3975.9 N, B = 24.36 Nsm−1 and C = 4.38 Nsm−2 , while the maximum traction and braking forces allowed for this particular train are in the form of look up tables (see Fig. 2). These forces are of the form FT (x(k), u(k)) = FT max (x2 (k))u(k),

2018 IFAC NMPC 208 Madison, WI, USA, August 19-22, 2018

Hafsa Farooqi et al. / IFAC PapersOnLine 51-20 (2018) 203–208

and FB (x(k), u(k)) = FBmax (x2 (k))u(k). The input variable is constrained in the set [−1, 1].

25 20

2

x2(m/s)

Maximum Traction and Braking Force (N)

10

5

1.5

10

1

5

0.5

0

0

0

10

20

30

40

50

60

70

80

Fig. 2. Maximum traction force FT max (dashed) and braking force FBmax (dotted) allowed for the train considered in the simulation. The considered track has zero curvature, slopes as plotted in Fig. 3, and velocity limits reported in Fig. 4.

200

400

600

800

1000

Fig. 4. Simulation results. Velocity profiles as a function of the train position obtained with different SPPC strategies: nominal (dotted line), two-step relaxed (dashed), multiobjective (dash-dot). The “all-out” solution (x f = 74.2 s) is shown with a thick solid line, and velocity limits with a thin solid line. 6. CONCLUSIONS We considered the problem of energy-efficient train operation, involving a finite terminal position and corresponding state constraint, and a discrete set of allowed input values or operating modes. To address this problem, we proposed a MPC approach with a shrinking horizon and a particular input parametrization, which allows one retain a continuous optimization problem at each step. We derived convergence guarantees in both nominal conditions and under uncertainty, and showcased the technique in realistic simulations.

40 30 20 10 0 -10 -20

ACKNOWLEDGEMENTS

-30 -40

0

kD s(m)

x2(km/h)

α(x 1) (deg)

15

0

200

400

600

800

The authors thank Alstom rail transport for providing the model and track data employed in the simulation study.

1000

kD s(m)

Fig. 3. Slopes of the Amsterdam metro track segment considered in the simulation. The train has to reach the next station at s f = 1106 m in x f = 76 s. The track is divided into sn = 15 sectors. The employed sampling distance is Ds = 18.4 m. We compare a nominal situation, where Algorithm 1 is applied and no model uncertainty is present, with a case where we employ Algorithm 2 (with both the two-stage approach and the multi-objective optimization variants) in presence of random parameter uncertainty (±10% of each model parameter). For the multi-objective approach, we set ω = 160 (see (26)). The obtained velocity profiles are presented in Fig. 4. From the plots, it is evident that in order to save energy, after accelerating to a certain velocity, the train mostly coasts or cruises taking advantage of the slopes. All the velocity constraints are always satisfied. Regarding the obtained final time between the terminal state and the target one at k = k f , in the presence of uncertainty this is delayed by 2 s (i.e. x f = 78 s) for both the 2-stage approach and the multi-objective one. This result is perfectly compatible with the desired performance in this application. Fig. 4 presents also the “all-out” solution, which achieves the shortest arrival time compatible with all the constraints. This corresponds to x f = 74.2 s. 240

REFERENCES Aradi, S., B´ecsi, T., and G´asp´ar, P. (2014). Design of predictive optimization method for energy-efficient operation of trains. In European Control Conference (ECC), 2490–2495. Fagiano, L. and Teel, A. (2013). Generalized terminal state constraint for model predictive control. Automatica, 49(5), 2622–2631. Farooqi, H., Fagiano, L., and Colaneri, P. (2018). On shrinking horizon move-blocking predictive control. Submitted to Automatica, Preliminary version available on arXiv, identifier/1803.09676. URL http://arxiv.org/abs/1803.09676. Ichikawa, K. (1968). Application of optimization theory for bounded state variable problems to the operation of train. Bulletin of JSME, 11(47), 857–865. Maksimov, V. (1971). Optimal control of automatic subway trains. Proceedings of Moscow Railway Engineering Institute (Trudy MIIT), (388), 103–114. Milroy, I.P. (1981). Minimum-energy control of rail vehicles. South Australian Institute of Technology. Scheepmaker, G.M., Goverde, R.M., and Kroon, L.G. (2017). Review of energy-efficient train control and timetabling. European Journal of Operational Research, 257(2), 355– 376. Strobel, H. and Horn, P. (1973). On energy-optimum control of train movement with phase constraints. Electric, Informatics and Energy Technique Journal, 6, 304–308.