N Gionfra

G Sandou

H Siguerdidjane
email: houria.siguerdidjane@centralesupelec.fr

D Faille

P Loevenbruck

A Distributed Consensus Control Under Disturbances for Wind Farm Power Maximization

In this paper we address the problem of power sharing among the wind turbines (WTs) belonging to a wind farm. The objective is to maximize the power extraction under the wake effect, and in the presence of wind disturbances. Because of the latter, WTs may fail in respecting the optimal power sharing gains. These are restored by employing a consensus control among the WTs. In particular, under the assumption of discrete-time communication among the WTs, we propose a distributed PID-like consensus approach that enhances the rejection of the wind disturbances by providing the power references to the local WT controllers. The latter are designed by employing a novel feedback linearization control that, acting simultaneously on the WT rotor speed and the pitch angle, guarantees the tracking of general deloaded power references. The obtained results are validated on a 6-WT wind farm example.

I. INTRODUCTION

Nowadays WTs, and wind farms still capture great research attention, as their role in the electric power supply is fairly changing. New grid requirements that has to be met, as well as an increasing know-how regarding the mentioned power systems lead to the interest for different ways to exploit the wind source. Far from being a classic modus operandi of wind farms, the power maximization problem falls within the latter. This is true when considering wind farms composed of several WTs, as they are very likely to experience the so-called wake effect. Thus considering the aerodynamic coupling among the turbines, and in turns the wind farm as a whole, proves potential gain when maximizing the power production (see e.g. [START_REF] Park | Cooperative wind turbine control for maximizing wind farm power using sequential convex programming[END_REF]), and justifies a growing interest in cooperative methods to control them. As a result, distributed control approaches are preferable to centralized ones if improvement of performance is seek [START_REF] Gebraad | A modelfree distributed approach for wind plant control[END_REF].

Typically the problem of power maximization under wake interaction is handled via a first step of optimization under the assumption of a static system. This approximation is mainly due to the high nonconvexity of the wake model that makes the problem hard to be treated directly under a control perspective. The available approaches mainly deal with either model-free decentralized methods, as in [START_REF] Marden | A model-free approach to wind farm control using game theoretic methods[END_REF], [START_REF] Gebraad | A modelfree distributed approach for wind plant control[END_REF],

philippe.loevenbruck@edf.fr or model-based ones as in [START_REF] Tian | Active power dispatch method for a wind farm central controller considering wake effect[END_REF], and [START_REF] Park | Cooperative wind turbine control for maximizing wind farm power using sequential convex programming[END_REF]. The aforementioned optimization algorithms rely on the existence of local control strategies for individual WTs that can stabilize around the obtained optimal set points [START_REF] Marden | A model-free approach to wind farm control using game theoretic methods[END_REF]. However, as shown in [START_REF] Heer | Model based power optimisation of wind farms[END_REF], and [START_REF] Gionfra | Hierarchical control of a wind farm for wake interaction minimization[END_REF], the actual attainable power gain can be highly affected by the system dynamics, and it validates the need for the design of efficient controllers to support the optimization step. It is important to point out that wind farm power maximization can be alternatively seen as the problem of finding the optimal power sharing of the available wind source among the WTs. Similar power sharing problems for wind farms have been treated in [START_REF] Zhang | Fully distributed coordination of multiple dfigs in a microgrid for load sharing[END_REF], [START_REF] Biegel | Distributed low-complexity controller for wind power plant in derated operation[END_REF], [START_REF] Baros | Distributed torque control of deloaded wind dfigs for wind farm power output regulation[END_REF]. The latter, based on the common assumption that the available wind power is higher than the demanded one, employ different distributed control approaches to deal with the problem. In particular, in [START_REF] Zhang | Fully distributed coordination of multiple dfigs in a microgrid for load sharing[END_REF], classic consensus algorithm is employed to let knowledge of the total available and demanded powers in order to compute the local WT set points. [START_REF] Biegel | Distributed low-complexity controller for wind power plant in derated operation[END_REF] proposes a static optimal consensus law to enable power sharing, and [START_REF] Baros | Distributed torque control of deloaded wind dfigs for wind farm power output regulation[END_REF] capitalizes on a leader-follower approach.

In this paper we propose a distributed control approach to let proper power sharing among the WTs in order to maximize the power generation. To the authors' knowledge such distributed control framework was never applied to the problem addressed in this paper, and, despite having some common ideas, it substantially differs from the mentioned references. Our contribution is two-fold. Firstly, a feedback linearization (FL) control is applied at the WT level to let the distributed problem be treated in the linear systems framework. Moreover, this control step has to allow the WT to track a general deloaded power reference, which is a necessary condition for wind farm power maximization [START_REF] Park | Cooperative wind turbine control for maximizing wind farm power using sequential convex programming[END_REF]. The overall approach can be seen as a novel wind farm level maximum power point tracking (MPPT) technique. [START_REF] Boukhezzar | Nonlinear control of variable speed wind turbines for power regulation[END_REF] proposes a FL controller for variable-speed fixed-pitch WTs to track a desired power reference, but it is only applicable at high wind speed. [START_REF] Gionfra | Combined feedback linearization and mpc for wind turbine power tracking[END_REF] capitalizes on a combined FL, and model predictive control technique enabling general power tracking in the whole WT operating envelope, but it does not particularly simplify the design of a distributed control law. Thus, based on [START_REF] Boukhezzar | Nonlinear control of variable speed wind turbines for power regulation[END_REF] we make use of a new FL application to WTs to let power tracking below the rated wind speed. Secondly, under the assumption of discrete-time communication among the WTs, we employ a distributed PID-like control architecture to force the system to stabilize around the optimal power sharing set points under system disturbances. The PID structure is justified for a simpler P-like protocol would not allow a satisfactory disturbance rejection, if the dynamics of the agents, i.e. the WTs, are general. In the literature, generic distributed controllers are proposed for both continuous and discrete-time in e.g. [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF], and in continuous time in [START_REF] Xi | Output consensus analysis and design for high-order linear swarm systems: partial stability method[END_REF], and [START_REF] Liu | Dynamic output feedback control for consensus of multi-agent systems: an h ∞ approach[END_REF] to meet respectively limited energy, and H ∞ performance constraints. As far as dynamic controllers in continuous-time with a prescribed structure are concerned, one can cite for instance [START_REF] Carli | A pi consensus controller for networked clocks synchronization[END_REF], and [START_REF] Ou | Distributed h ∞ pid feedback for improving consensus performance of arbitrary-delayed multi-agent system[END_REF], where the former proposes a PI-like algorithm for single integrator dynamic agents, and the latter provides a PID-like one for general high-order SISO systems. Our proposed PID-like distributed algorithm applies to general linear discrete-time MIMO systems, and it aims at solving the problem of leaderless consensus under the presence of disturbances. Moreover we provide a possible way of tuning the control parameters based on the results of [START_REF] Wu | Lmi-based multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF] for LMI-based tuning of MIMO centralized PID controllers.

The remainder of the paper is organized as follows. First, some graph theory preliminaries are provided in Section II. The turbine model is presented in Section III. The main control problems, and their objectives are stated in Section IV. We present, and prove our main results concerning the control architecture in Section V, and test its effectiveness on a 6-WT wind farm in Section VI. The paper ends with conclusions, and future perspectives in Section VII.

II. GRAPH THEORY PRELIMINARIES

An undirected graph G is a pair (V , E ), where V = {1, . . . , N} is the set of nodes, and E ⊆ V × V is the set of unordered pairs of nodes, named edges. Two nodes i, j are said to be adjacent if (i, j) ∈ E . Under the assumption of undirected graph, the latter implies that ( j, i) ∈ E too. An undirected graph is connected if there exists a path between every pair of distinct nodes, otherwise is disconnected. The adjacency matrix A = [a i j ] ∈ R N×N associated with the undirected graph G , considered in this paper, is defined by a ii = 0, i.e. self-loops are not allowed, and a i j = 1 if (i, j) ∈ E . The Laplacian matrix L ∈ R N×N is defined as L ii = ∑ j =i a i j and L i j = -a i j , i = j. Considering an undirected graph we make use of the following Lemma 1: [START_REF] Ren | Consensus seeking in multiagent systems under dynamically changing interaction topologies[END_REF] The Laplacian matrix has the following properties: (i) L is symmetric and all its eigenvalues are either strictly positive or equal to 0, and 1 1 1 is the corresponding eigenvector to 0; (ii) 0 is a simple eigenvalue of L if and only if the graph is connected. We will also make use of another Laplacian matrix, according to the following Lemma 2: [START_REF] Lin | Distributed leadless coordination for networks of second-order agents with time-delay on switching topology[END_REF] Let L = li j ∈ R N×N be a Laplacian matrix such that li j = N-1 /N if i = j, and li j = -1 /N otherwise, then the following hold: (i) the eigenvalues of L are 1 with multiplicity N -1, and 0 with multiplicity 1. 1 1 1 and 1 1 1 are respectively the left and right eigenvector associated to eigenvalue 0; (ii) there exists an orthogonal matrix U ∈ R N×N , i.e. U : U U = UU = I, and whose last column is equal to 1 1 1 / √ N, such that for any Laplacian matrix L associated to any undirected graph we have

U L U = I N-1 0 (N-1)×1 0 1×(N-1) 0 Λ, U L U = L 1 0 (N-1)×1 0 1×(N-1) 0
where L 1 ∈ R (N-1)×(N-1) is symmetric and positive definite if the graph is connected. Moreover we deduce the following extension of Lemma 2. Lemma 3: Let L ∈ R N×N be the Laplacian matrix associated to an undirected connected graph, and let D ∈ R N×N 0, and symmetric, then the following hold: (i) L DL is a Laplacian generally nonsymmetric matrix, L 0, all its eigenvalues are real, and 0 is a simple eigenvalue with associated eigenvector 1 1 1; (ii) consider the orthogonal matrix U ∈ R N×N defined in Lemma 2, then

U L U = L1 0 (N-1)×1 * 0
where L1 ∈ R (N-1)×(N-1) 0, and its eigenvalues are real.

Proof: We have that DL = D 1/2 (D 1/2 L D 1/2 )D -1/2
, thus DL is similar to a symmetric semi-definite positive matrix, so its eigenvalues are positive real. L preserves the 0 eigenvalue, and its associated eigenvector 1 1 1, as DL 1 1 1 = 0 0 0. 0 is a simple eigenvalue for D is nonsingular, and L has one simple 0 eigenvalue by hypothesis. The last column of U L U has all its entries equal to 0 because the last column of U is 1 1 1 / √ N. Being U L U block triangular, and similar to L , L1 has all real strictly positive eigenvalues. We employ the Kronecker product ⊗, for which we have Lemma 4: [START_REF] Graham | Kronecker products and matrix calculus with applications[END_REF] Suppose that U ∈ R p×p , V ∈ R q×q , X ∈ R p×p , and Y ∈ R q×q . The following hold: (i) (U ⊗V ) (X ⊗Y ) = UX ⊗VY ; (ii) suppose U, and V invertible, then (U ⊗V ) -1 = U -1 ⊗V -1 .

III. WIND TURBINE MODELING

The wind turbine model describes the conversion from wind power to electric power. The wind kinetic energy captured by the turbine is turned into mechanical energy of the turbine rotor, turning at an angular speed ω r and subject to a torque T r . In terms of extracted power, it can be described by the nonlinear function P r = ω r T r = 1 /2ρπR 2 v 3 C p (λ , ϑ ), where ρ is the air density, R is the radius of the rotor blades, ϑ is the pitch angle, v is the effective wind speed representing the wind field impact on the turbine, obtained by filtering the time series of wind data as described by [START_REF] Petru | Modeling of wind turbines for power system studies[END_REF], λ is the tip speed ratio given by λ = ω r R v . C p , nonlinear function of the tip speed ratio and pitch angle, is the power coefficient. This is typically provided in turbine specifications as a lookup table. As far as the turbine parameters are concerned, in this work we make use of the CART (Controls Advanced Research Turbine) power coefficient. This turbine is located at NREL's National Wind Technology Center. Nonetheless, we employ a polynomial approximation of the latter for the purpose of the synthesis of the controller. Referring to a two-mass model as in [START_REF] Boukhezzar | Nonlinear control of variable speed wind turbines for power regulation[END_REF], and as shown in Fig. 1, then, the low speed shaft torque T ls acts as a braking torque on the rotor, the generator is driven by the high speed torque T hs , and braked by the generator electromagnetic torque T g . The drive train turns the slow rotor speed into high speed on the generator side, ω g . Finally J r is the rotor inertia, K r , and K g damping coefficients, n g the gear ratio, and J g the generator inertia. The dynamics of the WT is thus described by J r ωr = T r -K r ω r -T ls , and J g ωg = T hs -K g ω g -T g . In this paper we also consider a first order system to model the pitch actuator, endowed with a sigmoid function σ : R → [ϑ min , ϑ max ] to model the pitch constraints. In addition, for ease of further development we can bring the system equations back on the low speed side, obtaining the simplified overall model

       τ ϑ θs = -ϑ s + ϑ r ϑ = σ (ϑ s ) J t ωr = P r (ω r , ϑ , v) ω r -K t ω r -T g (1)
where J t J r + n 2 g J g , K t K r + n 2 g K g , and where we used the relation n g = ω g/ω r = T ls/T hs . Eventually, neglecting the generator loss, the electric power delivered to the grid is P = T g ω r .

The system inputs are T g , and ϑ r , while the wind speed v acts as a disturbance. The feasible domain of the state variable is

X (ω r , ϑ ) ∈ R 2 : ω r ∈ [ω r,min , ω r,max ], ϑ ∈ [ϑ min , ϑ max ] . IV. PROBLEM STATEMENT At low wind speed, i.e. v ∈ [v min , v n ],
WTs are usually operated according to the well-known MPPT algorithm. The maximum power that a WT can extract from the wind is thus attained for a constant value of ϑ , named here ϑ o , depending on the turbine C p , and by controlling the WT to track the optimal tip speed ratio value

λ o arg max λ P r (v, ϑ o , λ ) = arg max λ C p (λ , ϑ o ). We name C o p C p (λ o , ϑ o ), and P o (v) P r (v, ϑ o , λ o ). For the considered CART turbine λ o ∼ = 8.
Nonetheless, as mention in Section I, when considering the wake effect in the optimization step of a farm of N WTs, the optimal value of C p related to the generic turbine i is such that C * p,i ≤ C o p . As a matter of fact, this implies that a turbine i should track an optimal power reference P * i (v) that satisfies

P * i (v) ≤ P o i (v), i.e
. it has to be deloaded if maximum wind farm power is seek. The reader may refer to the cited works in Section I to see how values C * p,i can be computed. Note that according to the usually employed wake models, the static optimization step needs to be run only when the wind direction changes, as optimal values C * p,i do not depend on the wind speed value [START_REF] Park | Cooperative wind turbine control for maximizing wind farm power using sequential convex programming[END_REF].

Assumption 1: The average wind direction is considered to be slowly varying with respect to the system dynamics. Thus, it is considered to be constant. In the sequel, for consistency of notation we add the index i to the WT variables described in Section III when referring to turbine i variables, and we drop it when the results hold for any WT. We can formulate the control problem in two subproblems, the first of which being Problem 1: Consider the system described by [START_REF] Park | Cooperative wind turbine control for maximizing wind farm power using sequential convex programming[END_REF]. Given an effective wind speed signal v(t), and a time-varying reference trajectory P re f (t), verifying P re f (t) ≤ P o (t) ∀t ≥ 0, find the signals (ϑ r (t), T g (t)) ∀t ≥ 0 such that lim t→∞ |P re f (t) -P(t)| = 0 for every initial condition (ω r (0), ϑ (0)) ∈ X : P(0) ≤ P o (0). Let us now assume that each local WT controller can measure, or estimate, the effective wind speed v m,i (t) such that v i (t) = v m,i (t) + v d,i (t) ∀t ≥ 0. Thus v d,i represents a nonmeasurable time varying disturbance for turbine i.

Assumption 2: We consider small disturbances v d,i with respect to v m,i , and slowly-varying with respect to the dynamics of (1). In nominal condition, i.e. if v d,i ≡ 0, each WT can compute the optimal power reference, as described in [START_REF] Gionfra | Hierarchical control of a wind farm for wake interaction minimization[END_REF], from its maximum available power P o i , according to

P i = C * p,i C o p P o i (2)
We can additionally require the WTs to meet an optimal power sharing condition given by P i

χ i = P k χ k i, k = 1, • • • , N (3)
Indeed, by naming P o ∞ the maximum power that a WT could extract from the wind if there was no wake effect, from [START_REF] Marden | A model-free approach to wind farm control using game theoretic methods[END_REF] we have that

P i/C * p,i = P o i /C o p = γ i P o ∞/C o p , thus C o p/γ i C * p,i P i = P o ∞ , i = 1, • • • , N, yielding P i/γ i C * p,i = P k/γ k C * p,k , i, k = 1, • • • , N. We name χ i γ i C * p,i ∈ R + ,
and where 3 are constant values for any value of v ∞ according to Assumption 1, being v ∞ the free stream wind speed. Despite being redundant information with respect to (2) in nominal conditions, (3) provides additional signals that can be exploited when the system is subject to disturbances. We can now state the second subproblem.

γ i = P o i /P o ∞ = ( v i/v ∞ )
Problem 2: Given N identical WTs, allowed to communicate on an undirected connected graph G ; given optimal values C * p,i , and

χ i , i = • • • , N; find P re f i (t) ∀t ≥ 0, i = • • • , N such
that every P i tracks (2), while minimizing the error

| P i/χ i -P k/χ k |, i, k = • • • , N, under the presence of v d,i (t).
Note that a similar idea of constant weighting factors is used in [START_REF] Baros | Distributed torque control of deloaded wind dfigs for wind farm power output regulation[END_REF] to deal with wind farm power regulation, and in [START_REF] Schiffer | Voltage stability and reactive power sharing in inverter-based microgrids with consensusbased distributed voltage control[END_REF] for reactive power control in microgrids.

V. CONTROL DESIGN A. FL step for local WT control

According to the optimization step, it turns out that every WT causing a reduction of available wind power of another one, is very likely to be subject to an optimal C p value such that C * p,i < C o p , i.e. strictly inferior. Thus, WTs whose C * p,i verifies C * p,i = C o p should simply perform classic MPPT regardless the disturbances of the system and the other WTs operating points, and they can be controlled with classic local controllers. In the sequel we only consider WTs that have to be strictly deloaded with respect to their P o i . Following [START_REF] Boukhezzar | Nonlinear control of variable speed wind turbines for power regulation[END_REF], the local control is composed of a first loop to control ω r . We impose a first order dynamics to the rotor speed tracking error ε ω ω re fω r : εω + a 0 ε ω = 0, by choosing a 0 ∈ R + . By naming w a 0 ω re f + ωre f , this is attained via

T g = T r -(K t -a 0 J t )ω r -J t w (4)
Differently from [START_REF] Boukhezzar | Nonlinear control of variable speed wind turbines for power regulation[END_REF], we choose to regulate the power output P by acting on the pitch angle. We impose a first order dynamics to the electric power tracking error ε p P re f -P:

εp + b 0 ε p = 0 ( 5
)
by choosing b 0 ∈ R + . This is attained via FL on (1) by choosing the feedback linearizing input

ϑ r = 1
β (ω r , ϑ s , v) Ṗre f -ω r ∂ T r ∂ v v + ω r τ ϑ ∂ T r ∂ ϑ dσ dϑ s ϑ s + J t ẇω r + 2(K t -a 0 J t )ω r -T r + J t w -ω r ∂ T r ∂ ω r (-a 0 ω r + w) +b 0 ε p ) (6)
where

∂ T r/∂ ω r , ∂ T r/∂ ϑ , and ∂ T r/∂ v are functions of (ω r , ϑ , v),
and where β (ω r , ϑ s , v)

ω r τ ϑ dσ dϑ s ∂ T r ∂ ϑ
. As pointed out in our previous work [START_REF] Gionfra | Combined feedback linearization and mpc for wind turbine power tracking[END_REF], there exist points in which β = 0, called singular points, i.e. points in which [START_REF] Gionfra | Hierarchical control of a wind farm for wake interaction minimization[END_REF], solution of the FL problem with respect to output P, is not defined. These points are determined by the so-

lution of ∂C q ∂ ϑ (ω r , ϑ , v) = ∂C q ∂ ϑ (λ , ϑ ) = 0, being C q C p/λ , for β (ω r , ϑ s , v) = ω r 2τ ϑ ρπR 3 v 2 dσ dϑ s ∂C q ∂ ϑ (ω r , σ (ϑ s ), v) ∼ = ω r 2τ ϑ ρπR 3 v 2 ∂C q ∂ ϑ (ω r , ϑ , v)
in the domain of interest of ϑ , and ω r , v > 0. In Fig. 2 the white area represents Λ = {(λ , ϑ ) : (ω r , ϑ ) ∈ X ∧ β < 0}. If ω re f is chosen to let λ be in a neighborhood of λ o , and ϑ > 0 • in order to let the WT to be deloaded, then it is clear that β is negative valued in the points of functioning of interest. In order to ensure that the trajectories of the closed loop system, defined by ( 1), ( 4), [START_REF] Gionfra | Hierarchical control of a wind farm for wake interaction minimization[END_REF], do not pass through singular points, differently from [START_REF] Gionfra | Combined feedback linearization and mpc for wind turbine power tracking[END_REF], we consider a modified FL function for ϑ r , by substituting the β function appearing in [START_REF] Gionfra | Hierarchical control of a wind farm for wake interaction minimization[END_REF] with

β ω r 2τ ϑ ρπR 3 v 2 dσ dϑ s ∂C q ∂ ϑ (λ , ϑ ) -ε(λ , ϑ ) ε(λ , ϑ )    α max ∂C q ∂ ϑ (λ , ϑ ), 0 if ∂C q ∂ ϑ (λ , ϑ ) = 0 ε 1 otherwise (7)
where ε 1 is a small positive value, and α > 1 is a tunable parameter to let some margin to have β negative valued in the system trajectories. Thus we obtain Λ = (λ , ϑ ) : (ω r , ϑ ) ∈ X ∧ β < 0 , shown in Fig. 2. The idea is to perform an approximated FL only when the system trajectories come close to a singular point. Clearly, in this case, the chosen ϑ r no longer guarantees satisfaction of (5). Nonetheless, under proper choice of ω re f , and deloaded mode of functioning, approximation (7) may occur only during transients. We can summarize the results in this subsection by stating the following Theorem 1: Given system (1), controlled via (4), and ( 6), where the β function is substituted with [START_REF] Zhang | Fully distributed coordination of multiple dfigs in a microgrid for load sharing[END_REF]. For any initial condition (ω r (0), ϑ (0)) ∈ Λ, the system trajectories are bounded if parameters b 0 , ε 1 , and α are chosen such that ε 1 > 0 is sufficiently small, α > 1, and b 0 >α /1-α. In addition, if ∃t ≥ 0 such that

∂C q ∂ ϑ (λ (t), ϑ (t)) < 0 ∀t ≥ t, then lim t→∞ |P re f (t) -P(t)| = 0.
Proof: First of all, initial conditions in Λ imply β (0) < 0, then ε 1 > 0, α > 1 allow β (t) < 0 ∀t ≥ 0. In particular β (t) = 0 ∀t ≥ 0, thus (6) is well-defined. Note that initial conditions required in Problem 1 satisfy β (0) < 0, according to Λ in Fig. 2. The system dynamics in closed loop is given by εp

= -b 0 + 1 - β β + β ε ε p + 1 - β β + β ε ϕ(ς ) (8)
where we named β ε ββ , and ϕ(ς ) the function composed of all the terms appearing in the right factor of (6) deprived of the term b 0 ε p , and being ς (ω r , v, ϑ s , v, w, ẇ, Ṗre f ). The term (1β /β+β ε )ϕ(ς ) is bounded in the trajectories thanks to the choice of β ε , and being ϕ a continuous function on a compact set. The latter is compact because the wind is limited, w, ẇ, Ṗre f are chosen to be so, ω r is bounded thanks to (4), and term dσ dϑ s ϑ s is bounded. Thus it will be considered as a bounded input of (8) to simplify the analysis. Finally the autonomous system of ( 8) is stable if, for instance, we choose b 0 >α /1-α, and ε 1 :

ε -1 1 > 1 -b 0 .
This can be proved by choosing a Lyapunov function for the uncertain system εp = (-b 0 + θ )ε p , where θ ∈ min -b 0 + 1ε -1 1 , -b 0α /1-α , 0 . Eventually, if for some t ≥ 0 : ∂C q ∂ ϑ < 0 ∀t ≥ t, then (8) reduces to εp = -b 0 ε p , thus P → P re f for t → ∞.

Remark 1: Concerning ω re f , we make the choice to use the signal ω o λ o v /R sufficiently filtered of its high frequency components. There are different motivations to support this choice. First of all, if v varies rapidly so it does ω o , then if we consider ω re f = ω o , its variation would directly effect ϑ r via [START_REF] Gionfra | Hierarchical control of a wind farm for wake interaction minimization[END_REF], and in turns ϑ . This fact risks to make ϑ hit the saturation constraints of the sigmoid function, and more in general, to not let the constraints on θ be respected, as in this framework they are only verified a posteriori. Secondly, if ω re f varies too rapidly, by empirical results it turns out that closed-loop system trajectories are more likely to approach singular points, letting the activation of ε(λ , ϑ ) defined in [START_REF] Zhang | Fully distributed coordination of multiple dfigs in a microgrid for load sharing[END_REF], and not allowing satisfaction of [START_REF] Heer | Model based power optimisation of wind farms[END_REF]. On the other hand, filtering ω o let (6) be defined. The physical explanation of this fact is that there exist infinite pairs (ω, ϑ ) ∈ X to deload a WT, i.e. to track P re f , [START_REF] Gionfra | Combined feedback linearization and mpc for wind turbine power tracking[END_REF].

Remark 2: One of the reasons why the described approach should not be considered for classic MPPT mode of functioning lies in the fact that it is impossible to simultaneously track ω o , and P o under any control. This is easily seen considering ( 1). In addition, for the purpose of this analysis, let us neglect K t . Then suppose that, thanks to a controller, ω r = ω o , and P r = P o ∀t ≥ 0. If v is not constant, then | ωr | > 0, which implies |P o -P| > 0, and in particular P = P o . This basically means that condition P = P o cannot be forced via (6) because it is not defined. So if ( 6) is employed, trajectories would pass through singular points.

B. Additional Local Control Settings

From now on we carry out the analysis under the following Assumption 3: Trajectories of the closed-loop system described by ( 1), ( 4), [START_REF] Gionfra | Hierarchical control of a wind farm for wake interaction minimization[END_REF] verify ∂C q ∂ ϑ < 0. As mentioned in Section IV, let us assume that turbine i local controller is able to measure v m,i such that v i = v m,i + v d,i . The effect of v d,i on the closed loop dynamics can be approximated as

εp,i = -b 0 ε p,i + µ 1 ( ζi )v d,i + µ 2 ( ζi )v 2 d,i + µ 3 ( ζi ) vd,i (9)
obtained via first order Taylor expansion of the functions depending on v i , in a neighborhood of v m,i , e.g. T r

(v i ) ∼ = T r (v m,i ) + ∂ T r ∂ v (v m,i )v d,i
, and where ζi (ω r,i , ϑ i , v m,i , vm,i ). Functions µ 1 , µ 2 , µ 3 are not reported in this paper for the sake of brevity. According to Assumption 2 we neglect the last term of [START_REF] Baros | Distributed torque control of deloaded wind dfigs for wind farm power output regulation[END_REF]. Moreover the contribution of term µ 2 v 2 d,i can be neglected with respect to µ 1 v d,i . On the compact set on which µ 1 is defined, the function satisfies µ 1,min ≤ µ 1 ≤ 0, thus in the sequel we treat µ 1 as a parametric uncertainty, and we drop its dependency on ζi for ease of notation. Being interested in a discrete-time communication set-up among the WTs we shall consider the discretized system of ( 9)

ξ k+1 i P k+1 i = 0 0 -1 (1 -T s b 0 ) ξ k i P k i + 1 1 + T s b 0 P re f ,k i + 0 µ 1 T s v k d,i (10)
where we used Euler approximation using sampling time T s , we approximated Ṗre

f i ∼ = (P re f ,k i -P re f ,k-1 i ) /T s , we named ξ k i P re f ,k-1 i
, and apex k stands for time kT s . Before providing the distributed controller, we add an additional local PI loop to [START_REF] Boukhezzar | Nonlinear control of variable speed wind turbines for power regulation[END_REF] to be tuned to enhance rejection of v d,i . As it will be clear, the latter also has an important role on the wind farm consensus. Naming K l I , and K l P respectively the integral, and proportional gains of the PI, we can write [START_REF] Boukhezzar | Nonlinear control of variable speed wind turbines for power regulation[END_REF] in closed-loop as

x k+1 i = Ax k+1 i + B 2 u i + B f w P f w i + B w v d,i
, where

A    1 0 -K l I T s 1 0 -K l P (1 + T s b 0 ) -1 (1 -T s b 0 -K l P (1 + T s b 0))    , B w   0 0 µ 1 T s   B 2   K l I T s 1 (1 + T s b 0 )   , B f w   K l I T s K l P (1 + T s b 0 )K l P   (11)
and where we named x i [δ i , ξ i , P i ] , being δ i the state of the integral action, P f w i a forward signal, and u i is left as a degree of freedom to let distributed control. Note that other choices for the placement of u i would have been possible.

Remark 3: Forward signal P f w i is set to be equal to (2), as in steady state it respects the optimal power sharing provided by the optimization step. However, the latter is based on the computation P o i , in turns based on the wind measure v m,i . Thus if v d,i = 0 then P f w i is a set point for [START_REF] Gionfra | Combined feedback linearization and mpc for wind turbine power tracking[END_REF] that may not respect (3). This motivates next subsection analysis, where optimal power sharing is restored via a distributed algorithm by acting on u i .

C. Distributed PID-like Consensus

Let us defined the consensus problem addressed in a general formulation. Consider N identical agents governed by general discrete-time linear dynamics, according to

x + i = Ax i + B 2 u i + B 1 ω i , y i = Cx i i = 1, • • • , N (12)
where

A ∈ R n×n , B 2 ∈ R n×l , B 1 ∈ R n×h , C ∈ R m×n , x i x i (k) ∈ R n and x + i x i (k + 1) ∈ R n
are respectively the agent state at the current step k, and at the next step k + 1,

u i u i (k) ∈ R l is the agent control, ω i ω i (k) ∈ R h its disturbance.
Being y i the measured and the controlled output, we additionally require l ≥ m. Let A be Schur stable, and let the agent communicate on an undirected graph whose Laplacian is L . Thus we address the problem finding a distributed control law for u i such that

y i/χ i -y k/χ k is minimized for i, k = 1, • • • , N with respect to the disturbance ω [ω 1 , • • • , ω N ] .
If the such error is zero, then we say that weighted consensus is achieved. We consider χ i ∈ R + to simplify the analysis. Results can be extended to case of higher dimensional weights for the general case of m > 1. By naming D diag( 1 /χ i ), we additionally define matrix L DL , which satisfies Lemma 3, and whose positive minimum nonzero and maximum eigenvalues are respectively λ L , and λ L . In this work we focus on local controllers of the form

x + c i = A c x c i + B c s i , u i = C c x c i + D c s i i = 1, • • • , N (13)
where x c i x c i (k) ∈ R 2l is the agent controller state, and

A c = I l I l 0 l×l 0 l×l 2l×2l B c = (K i -K d ) K d 2l×m C c = I l 0 l×l l×2l D c = [(K p + K i + K d )] l×m (14)
where K p , K i , K d ∈ R l×m are gain matrices to be tuned, and where s i s i (k) ∈ R m :

s i N ∑ k=1 a ik y i χ i - y k χ k (15)
Thus the closed-loop system for agent i has dimension n n+2l. As shown by [START_REF] Wu | Lmi-based multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF], ( 13) is a state representation of the discrete-time PID MIMO controller, whose z-transform

is u i (z) s i (z) = K p + K i z z -1 + K d z -1 z .
The problem can now be restated as the one of finding the matrices B c , and D c such that the effect of disturbance ω on the weighted consensus is minimized. Before stating the result we introduce the following Definition 1: System ( 12) is said to achieve fast weighted consensus with performance index τ ∈ R + if for any time-constant disturbance ω, and any initial condition, lim k→∞ y i/χ iy k/χ k = 0 for i, k = 1, • • • , N, and (1e -1 )% of consensus is achieved with a time constant inferior to τT s .

Theorem 2: Given the system described by [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF], where N agents can communicate on an undirected connected graph; consider the distributed protocol of equations ( 13),( 14), [START_REF] Carli | A pi consensus controller for networked clocks synchronization[END_REF]; then the systems achieve fast weighted consensus with per-

formance index τ = - 1 log(R)
, where R ∈ R : 0 ≤ R < 1, if there exist two symmetric positive definite matrices P, P ∈ R n× n such that the LMI conditions of Theorem 2 in [START_REF] Wu | Lmi-based multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF] are simultaneously satisfied for two LTI systems whose matrices are respectively (A, B 2 , λ L C), and (A, B 2 , λ L C), and where the real constants (a, b) to be set in Theorem 2 in [START_REF] Wu | Lmi-based multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF] are chosen to be (a, b) = (0, R). Proof: The closed-loop dynamics for the generic agent i, by using ( 12), [START_REF] Xi | Output consensus analysis and design for high-order linear swarm systems: partial stability method[END_REF], and by defining the augmented state ξ i x i , x c i ∈ R n, and matrices C [C 0 m×2l ], is given by ξ Considering the change of coordinates δ = L ⊗ I n ξ , it yields

+ i = Âξ i + B N ∑ k=1 a ik ( ξ i/χ i -ξ k/χ k ) + Bω i , y i = Cξ i , where  = A B 2 C c 0 A c , B = B 2 D c C2 B c C2 , B = B 1 0 By naming ξ ξ 1 , • • • , ξ N , y y 1 , • • • , y N ,
δ + = L ⊗ I n I N ⊗ Â + L ⊗ B ξ + L ⊗ I n I N ⊗ B ω = L ⊗ Â + L L ⊗ B δ + 1 1 1 ⊗ 1 N N ∑ k=1 ξk + L ⊗ B ω = L ⊗ Â + L L ⊗ B δ + L ⊗ B ω
where we used points (i) of Lemma 2, 3, and 4. According to the (ii) point of Lemma 2, we employ the orthogonal matrix U ∈ R N×N to define the change of coordinates: δ

U ⊗ I n δ , ω U ⊗ I h ω, ζ U ⊗ I m ζ
, so that the system equations in the new coordinates are given by

         δ + = U ⊗ I n L ⊗ Â + L L ⊗ B (U ⊗ I n) δ + U ⊗ I n L ⊗ B ω = Λ ⊗ Â + ΛU L U ⊗ B δ + Λ ⊗ B ω ζ = U ⊗ I m I N ⊗ C (U ⊗ I n) δ = I N ⊗ C δ (17)
As shown in Lemma 2, and 3, being the last rows of Λ, and U L U zeros, we can split [START_REF] Wu | Lmi-based multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF] in two by dividing the system variables as

δ = [ δ 1 , δ 2 ] , ω = [ ω 1 , ω 2 ] , and ζ = [ ζ 1 , ζ 2 ]
. It follows that, to conclude on system stability, we can study the reduced order system described by δ + 1 = I N-1 ⊗ Â + L1 ⊗ B δ1 + I N-1 ⊗ B ω1 , and ζ1 = I N-1 ⊗ C δ1 . From Lemma 3, it exists an invertible matrix V ∈ R (N-1)×(N-1) :

V -1 L1 V Λ = diag(λ 1 , • • • , λ N-1 ), where 0 < λ L ≤ λ i ≤ λ L for i = 1, • • • , N -1.
Thus we can define a further change of coordinates, such that δ1 V -1 ⊗ I n δ1 , ω1

V -1 ⊗ I h ω1 , and ζ1 V -1 ⊗ I m ζ1 . The latter yields

δ + 1 = I N-1 ⊗ Â + Λ ⊗ B δ1 + I N-1 ⊗ B ω1 ζ1 = I N-1 ⊗ C δ1 (18)
Note that the transfer function matrix of (18) satisfies

T ζ1 ω1 (z) ∞ = T ζ1 ω1 (z) ∞ = T ζ ω (z) ∞ = T ζ ω (z) ∞ .
We can now separate [START_REF] Ren | Consensus seeking in multiagent systems under dynamically changing interaction topologies[END_REF] in N -1 subsystems, each of them being governed by

           δ + 1 i = x+ 1 i x+ 1,c i = (A + B 2 D c (λ i C)) B 2 C c B c (λ i C) A c x1 i x1,c i + B 1 0 ωi ζ1 i = C x1 i (19)
System [START_REF] Lin | Distributed leadless coordination for networks of second-order agents with time-delay on switching topology[END_REF] can be equivalently seen as the closed-loop form of the two following systems

x+ 1 i = A x1 i + B 2 ũi + B 1 ωi ỹi (λ i C) x1 i x+ 1,c i = A c x1,c i + B c ỹ1 i ũi C c x1,c i + D c ỹ1 i (20)
Thus, we can reformulate the problem as the one finding matrices B c , and D c such that for i = 1, • • • , N -1 the closed-loop system of ( 20) is Schur stable when ω i = 0. We now invoke Theorem 2 in [START_REF] Wu | Lmi-based multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF] where it is shown that given two real constants (a, b), if there exists a symmetric positive definite matrix P i such that a given LMI condition is satisfied, then system (20) is stable with all its eigenvalues λ laying in the complex plane region defined by

F D (ℜ[λ ], ℑ[λ ]) : (ℜ[λ ] + a) 2 + ℑ[λ ] 2 < b 2 .
Such LMI condition happens to be affine in the system matrices, variables and matrix P i . We make use of this fact to provide sufficient conditions for which it exists a controller of the considered form such that the mentioned LMI is simultaneously verified for

i = 1, • • • , N -1. Since the generic eigenvalue of L1 : λ i is such that λ L ≤ λ i ≤ λ L , then it always exists α i ∈ R : 0 ≤ α i ≤ 1 so that λ i = α i λ L + (1 -α i ) λ L .
Notice that the systems to be stabilized, appearing in the first set of equations in [START_REF] Graham | Kronecker products and matrix calculus with applications[END_REF], can be seen as one single system with an uncertain measurement matrix, whose parameter is λ i . In other words, C i λ i C, and ∃α i :

C i = α i C min + (1 -α i )C max ,
where C min λ L C, and C max λ L C, i.e. it can be written as a convex combination of the extreme matrices C min , and C max . Thus, we make use of classic results of robust linear control, and in particular by introducing an affine parameter dependent Lyapunov matrix P(α i ) α i P + (1α i ) P, where P, P are Lyapunov matrices solution of simultaneous LMI of Theorem 2 in [START_REF] Wu | Lmi-based multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF] written for respectively C min , and C max . Eventually, it is easy to show that if P, P exist, then the controller solves the problem ∀λ ∈ R : λ L ≤ λ ≤ λ L , and in particular for λ = λ i , for i = 1, • • • , N. Such a controller is easily found from the solution of the aforementioned LMI condition. Indeed among the LMI variables there are matrices B c , and D c , from which it is easy to calculate the PID gains K p , K i , and K d by employing relations in [START_REF] Liu | Dynamic output feedback control for consensus of multi-agent systems: an h ∞ approach[END_REF]. If such a solution exists then the eigenvalues of system ( 16) are guaranteed to lay in F D . In this framework we are interested in speeding up the system response to ω. For this reason we set a = 0, and b = R, where R : 0 ≤ R < 1. Thus, all system eigenvalues are guaranteed to have a module inferior to R. As a result, the system has the slowest time-constant inferior to -T s/log(R).

Notice that the tools developed in this subsection can be directly applied to solve the wind farm weighted consensus. This is obtained by choosing as A, B 2 the homonym matrices of ( 11), C = [0 0 1], i.e. measure and control of P i , and B 1 = [B f w B w ], i.e. P f w i , and v d,i are both considered as disturbances with respect to the weighted consensus. Eventually, in order not to contrast Assumption 3, we require the additional Assumption 4: All v d,i are such that the corresponding P re f i , depending on (2), on the local PI control, and on the distributed PID, is always lower than the real maximum power that turbine i can extract from the wind.

Remark 4:

The distributed PID minimizes the influence of all P f w i , and v d,i on the power sharing error among P i . This can be seen as forcing given relative distances among the P i variables. However, the absolute values of all the P i are influenced by P f w i , and v d,i . This motivates the internal PI control pushing P i towards P f w i . Note that by construction P f w [P f w 1 , • • • , P f w N ] verifies L DP f w = 0 0 0.

VI. SIMULATIONS

For the following simulations we considered the real CART turbine C p . This represents the only source of modelplant mismatch. Test under more realistic uncertainties goes beyond the scope of this paper, and it will not be addressed. First of all, let us show the FL local controller behavior. During 600 s, the WT is excited by the effective wind speed signal of Fig. 3. Inputs signal shown in Fig. 4, generated by the FL controller, allow the WT to track a deloaded power reference of 50% with respect to P o . |P -P re f | goes to zero with a time constant depending on the chosen b 0 in (5), (see Fig. 3). Small persistent error oscillations are due to the C p mismatch. Concerning the wind farm power sharing simulations, we consider 6 aligned WTs that communicate with their direct neighbor WT. Wind speed signal v ∞ , blowing in front of turbine 1 is chosen as the previous simulation one. Turbine 6, being the last one, is required to operate in classic MPPT mode, thus it does not intervene in the consensus control, and its P signal will not be reported. Wind disturbances v d,i , and the controlled WT powers are shown in Fig. 5. For this problem, weighted consensus is achieved with performance index τ = 24.5. In order to show consensus achievement we provide two additional figures. Naming P [P 1 , • • • , P N ] , the first one shows signals L DP, which in ideal conditions should have all zeros entries in steady state. The second one shows DP, where its entries should ideally reach a common value. These simulations are shown in the bottom of Fig. 5, and in a zoomed window in Fig. 6. Note that C p mismatch, as well as temporary dissatisfaction of Assumption 3 cause persistent small oscillations on the reached weighted consensus.

VII. CONCLUSIONS

We presented a novel distributed approach to control a wind farm for power maximization under wake effect, based on two control layers. First, the proposed FL local controller allows a WT to track a deloaded power reference by acting on both the rotor speed, and the pitch angle. Then, a PIDlike discrete-time weighted consensus control is developed to let the wind farm turbines keep the optimal power sharing under the presence of wind disturbances. Indeed, because of the latter, without this second control layer, WTs would track an absolute power reference which may not respect the optimization step gains. In the near future work we are interested in extending the proposed framework to relax the assumptions on the local controller performance, and to treat other sources of system disturbances.

Fig. 1 :

 1 Fig. 1: Two-mass model of the turbine mechanics.

  (a) Λ: (white area), set of (λ , ϑ ) such that β (λ , ϑ ) < 0.(b) Λ: (white area), set of (λ , ϑ ) such that β (λ , ϑ ) < 0.

Fig. 2 :

 2 Fig. 2: Singular points with and without β approximation.

  gathering together the closed-loop agents dynamic, and performing the change of coordinates ξ = (D ⊗ I n)ξ , it yieldsξ + = I N ⊗ Â + DL ⊗ B ξ + I N ⊗ B ω ȳ = I N ⊗ C ξ (16)where we named ω (D ⊗ I h)ω , ȳ (D ⊗ I m)y, and we used point (i) of Lemma 4. Similar to[START_REF] Liu | Dynamic output feedback control for consensus of multi-agent systems: an h ∞ approach[END_REF], we define ζ i ȳi -1 N ∑ N k=1 ȳk , and δ i ξi -1 N ∑ N k=1 ξk , thus ζ i = Cδ i . Note that if ζ i = 0 for i = 1, • • • , N then ȳi = ȳk , i.e. weighted consensus is achieved. If we now name δ δ 1 , • • • , δ N , and ζ ζ 1 , • • • , ζ N , we have that ζ = I N ⊗ C δ , and δ = ξ -1 1 1⊗ 1 N ∑ N k=1 ξk = L ⊗ I n ξ , where L satisfies the conditions of Lemma 2. Thus ζ = I N ⊗ C L ⊗ I n ξ = L ⊗ C ξ .

Fig. 3 :Fig. 4 :

 34 Fig. 3: Wind speed signal.

  (a) Wind disturbances. (b) WTs electric power. (c) Consensus signals: L DP. (d) New coordinate power signals: P = DP.

Fig. 5 :

 5 Fig. 5: Distibuted control simulation.

Fig. 6 :

 6 Fig. 6: Zoom on system variables during distributed control simulation.

ACKNOWLEDGMENT

This study has been carried out in the RISEGrid Institute (www.supelec.fr/342p38091/risegrid-en.html), joint program between CentraleSupélec and EDF ('Electricité de France') on smarter electric grids.


Juliette Pochet

Guillaume Sandou

Sylvain Baro

Automatic Train Supervision for a CBTC Suburban Railway Line Using Multiobjective Optimization

Keywords: CBTC, railway, traffic management, suburban, rescheduling, optimization I

Communication-based train control (CBTC)

systems have been deployed on subway lines to increase capacities on existing infrastructures. For the same purpose, CBTC systems are to be deployed on suburban railway lines where operating principles and constraints are significantly different. In this paper, a regulation method for CBTC trains on a suburban line is presented. This method is designed to combine CBTC functionalities with suburban operating principles. It includes a traffic management method in station, and a rescheduling method in case of disturbances. The proposed regulation method is integrated into the railway system simulation tool SIMONE developed by SNCF. This simulation tool includes models the whole CBTC system, as well as the classic signaling system, train dynamics and railway infrastructures. Models of these different agents are described. The integration of the proposed regulation method into the tool SIMONE allows evaluating performances while taking into account the functional complexity of a CBTC railway system. The approach is illustrated with a realistic case: simulations of a CBTC traffic on the urban part of a railway line in the Paris region network are described. The proposed regulation method shows interesting results in disturbed situations according to the railway operating principles.

INTRODUCTION

In peak times, subway and railway companies operate short headways in order to increase the capacity. As a result, the system reaches its maximal capacity and is very sensitive to disturbances, even small ones. The operational margins are limited or non-existent. More and more, operating companies turn to new intelligent signaling and control systems to increase the line capacity and the traffic robustness without expensive modification of infrastructures.

Communication-based-train-control (CBTC) systems are signaling and control systems that have been developed in the last decades. Currently, they are mainly deployed on subway lines to increase the capacity and improve performances in high-density areas. In such systems, CBTC trains are accurately localized by means of communication devices (regardless of infrastructure detection devices). Besides, CBTC trains are controlled by an automatic pilot system. Through localization and automatic pilot, CBTC systems operate a moving block system and achieve shorter headways than the fixed block signaling system.

In a near future, CBTC systems will be deployed on suburban railway lines which cross big cities. Suburban railway lines are generally characterized by a high-density traffic in the urban part of the line, and a low-density traffic in the suburban parts. The CBTC area can be limited to high density traffic areas where there is a need to gain capacity. Thus, in such a configuration, the driver is in charge of piloting the train in manual mode in the suburbs, while the automatic pilot runs the train in the high density area.

Most CBTC systems include automatic train supervision (ATS) module to manage the traffic and reschedule CBTC trains in real time according to traffic monitoring information.

Suburban railway lines have specific operation principles. In the general case, they differ from subway lines on four mains aspects: (1) on a suburban line, operation management generally seeks to combine two objectives: on one hand, like on subway lines, management aims to maintain a constant headway between trains in the high density area of the line. On the other hand, unlike subway lines, operation management seeks to respect the timetable in the suburbs where the traffic is not dense. In disturbed situations, a compromise can be necessary. (2) On a suburban line, the network is more complex. The ATS module can be charged with choosing itineraries to link remarkable nodes and with choosing platforms stations. These choices are made in real time and depend on the traffic state. (3) On a suburban line partially equipped with a CBTC system, it is unlikely that platform screen doors will be deployed, for technical and economic reasons. Thus the ATS will face more disturbances such as passengers blocking the doors or going to the tracks. (4) A suburban railway line is part of a large regional network. In nominal mode, tracks can be shared with trains which are not controlled by the ATS. They can either be non-equipped trains which use the same tracks or CBTC-trains in manual mode.

Supervision and rescheduling of a CBTC suburban line are challenging topics in the field of traffic management. In this paper, we present a regulation approach based a traffic management method in station, and on a multi-objective optimization method for rescheduling after disturbances. This approach has been integrated in the ATS module of a simulation tool called SIMONE. SIMONE is developed by SNCF, the French national railway operator. This tool models the whole CBTC system, as well as the classic signaling system, train dynamics and railway infrastructures. Section II briefly outlines some related work in literature. Section III introduces the simulation tool SIMONE. Section IV focuses on the traffic model in the ATS module of SIMONE. Section V presents the proposed regulation approach. Section VI presents some experimental optimization results. Section VII provides some discussion and points to future work.

II. RELATED WORK

A. Rescheduling methods

In the literature of real-time traffic management, various methods have been proposed to solve the rescheduling problem in case of disturbances. Many approaches are based on multiobjective optimization and have shown interesting results on simulation tools. They tend to minimize delay as well as additional objectives such as travel time, energy cost, or passenger waiting time, depending on the application. There are four types of control actions: re-routing, modifying running times or dwell times (equivalent to modifying arrival and departure times), reordering at junctions, and re-servicing. Each rescheduling method generally uses one or two control actions, adapted to the kind of infrastructures and applications that are considered.

Most railway traffic studies apply re-ordering (or dispatching) as a traffic management action [START_REF] Ho | Railway Junction Traffic Control by Heuristic Methods[END_REF], [START_REF] Wegele | Automatic Dispatching of Train Operations Using a Hybrid Optimisation Method[END_REF], [START_REF] Fan | Railway "Traffic Rescheduling Approaches to Minimise Delays in Disturbed Conditions[END_REF], [START_REF] Tornquist | Train Traffic Deviation Handling Using Tabu Search and Simulated Annealing[END_REF], [START_REF] Kersbergen | Distributed model predictive control for rescheduling of railway traffic[END_REF]. In complex networks, reordering is combined with rerouting [START_REF] D'ariano | Improving Real-Time Train Dispatching: Models, Algorithms and Applications[END_REF], [START_REF] Corman | Real-time Railway Traffic Management: Dispatching in Complex, Large and Busy Railway Networks[END_REF], [START_REF] Flamini | Job-Shop Scheduling Algorithms with Applications to Railway Traffic Control[END_REF], [START_REF] Tornquist | N-tracked Railway Traffic Re-scheduling during Disturbances[END_REF] to make use of the network alternatives. In [START_REF] Fuchsberger | Algorithms for Railway Traffic Management in Complex Central Station Areas[END_REF], the rescheduling problem is adapted to complex central stations. The proposed approach assigns routes and running times to trains, in order to fully determine the block allocations.

This corresponds to a microscopic regulation method adapted to high-density area.

Usually, subway traffic management method are designed for closed near-linear networks, and so focus on modifying running and dwell times (or departure and arrival times) [START_REF] Breusegem | Traffic Modeling and State Feedback Control for Metro Lines[END_REF], [START_REF] Fernandez | SIRO: An Optimal Regulation System in an Integrated Control Centre for Metro Lines[END_REF], [START_REF] Assis | Generation of Optimal Schedules for Metro Lines Using Model Predictive Control[END_REF]. Re-servicing is only considered to manage urban or multi-modal networks offering alternative routes to passengers [START_REF] Zidi | SARR: Support System for Regulation and Reconfiguration of Multimodal Transportation Networks[END_REF]. Re-ordering and rerouting are not relevant for this kind of application.

Rescheduling methods are tested in simulation tools to evaluate their performances and their limits. Most of the times, simulation tools consist of a description of the tracks based on a node-link scheme, a description of train movements based on arrival and departure events, and a description of the signaling system based on a set of traffic constraints. These models allow describing traffic states with a flexible precision level and can be easily used for optimization process. The impact of the railway system components is not directly described; it appears through traffic constraints and delays.

To evaluate performances of a CBTC equipped line, it is not sufficient to model the CBTC system with constraints and delays. More precise models are needed to take into account its functional architecture, as well as some specific technical choices that can vary from a manufacturer to another. These features can have a strong impact on performances.

B. CBTC System

In the literature on CBTC systems, studies tend to focus on modeling the complex architecture of the system [START_REF] Chen | Development of the New CBTC System Simulation and Performance Analysis[END_REF], on analyzing availability and performance of communication processes [START_REF] Sun | Performance improvement in Communication-Based Train Control (CBTC) Systems using cognitive control[END_REF], or on designing speed profiles [START_REF] Fernandez-Rodriguez | Design of Robust and Energy-Efficient ATO Speed Profiles of Metropolitan Lines Considering Train Load Variations and Delays[END_REF], [START_REF] Bu | Predictive Function Control for Communication-Based Train Control (CBTC) Systems[END_REF].

CBTC systems are broadly deployed on closed near-linear networks. They generally operate simple but efficient regulation strategies based on maintaining a constant headway to handle high passenger volume [START_REF] Ning | An Integrated Control Model for Head-way Regulation and Energy Saving in Urban Rail Transit[END_REF]. Although usual CBTC traffic management strategies are efficient, they cannot be transposed to suburban railway context because the operation principles are not the same. Our work aims to develop a multiobjective optimization strategy for a suburban railway line, partially equipped with a CBTC system.

III. SIMULATION TOOL SIMONE

In previous work [START_REF] Pochet | Supervision and rescheduling of a mixed CBTC traffic on a suburban railway line[END_REF], the proposed method was firstly tested in a prototype traffic model to prove the viability of the approach on a reduced size and reduced complexity problem. This academic case study includes a description of the tracks based on a node-link scheme, a description of train movements based on arrival and departure events, and a description of the signaling system based on a set of traffic constraints. As described in section II, this kind of modelling is broadly used in literature to test traffic management strategy.

Our current work is supported by a new simulation tool SIMONE, developed by SNCF Reseau. This tool intends to simulate railway traffic including CBTC trains. In particular, this tool simulates the functional complexity of the different components of the railway system as well as their interactions. In this way, SIMONE takes into account the impact of functional complexity on performances.

The tool SIMONE is based on a multi agent simulation environment. A simulation core activates agents and assigns triggers such as messages, events or cycles. In this section, we describe some main agents and the models they embed. Due to space limitation, the models of tracks and vehicles are not described here.

A. Signaling system

In SIMONE, the classic signaling system is modeled with a set of state machines. The scope of the model consists of all the common components which have an impact on performances, and includes sensors (pedals, track circuits), interlocking, points, stop lights.

Each component is described by a state machine, and each state machine can be individually configured to include delays. On top of the component state machines, the main state machine models the software layer in charge of setting routes. It reproduces the functional sequences that lead to setting a route by activating the different components of the model. The SIMONE signaling system model is currently based on the French main line signaling system. In particular, it described interlocking and sectional route releasing concepts.

B. On-ground ATC module

On-ground ATC module is a wayside equipment. It localizes all trains (CBTC and non CBTC trains) on the map of the track. It also defines the movement authority of CBTC trains.

In the SIMONE framework, the on-ground ATC has a precise knowledge of the tracks, receives information regarding the tracks' state from the signaling system, and receives localization of trains from communication devices. The main functions of the on-ground ATC model include monitoring trains on a map of the tracks, computing the movement authority for each CBTC trains, sending the movement authorities to the on-board ATC modules, sending information regarding the tracks' state.

C. On-board ATC modules

The on-board ATC module is embedded into the CBTC train, it is responsible for respecting the safety area of the train called movement authority and computed by the on-ground ATC module. It drives the train and supervises it (w.r.t. safety). It computes its speed and localization and sends its localization report to the on-ground ATC module.

In the SIMONE tool, the main functions of the on-board ATC model include sending the localization of the train to the on-ground ATC, computing the target speed to respect the movement authority, computing traction or breaking command, operating turnarounds. The safe braking model, in compliance with IEEE Standards [START_REF]IEEE Standard for Communications-Based Train Control (CBTC) Performance and Functional Requirements[END_REF], consists of 5 phases: CBTC response time, traction disable response time, cruise time, emergency brake build-up, and emergency braking at guaranteed emergency brake rate.

D. Automatic Train Supervision (ATS)

ATS module supervises operations, adapts the CBTC trains flow and controls the signaling system, in order to reach the performance objectives.

In the SIMONE framework, the main functions of the ATS model include supervising trains through a traffic model, updating the traffic model with monitoring information, sending route commands to the signaling system, sending mission and schedule commands to the on-board ATC modules, rescheduling missions and schedule in case of disturbances. More details on the ATS module will be given in section IV and section V.

In the next section, we describe the ATS traffic model implemented in the SIMONE tool.

IV. ATS TRAFFIC MODEL

In traffic management systems, regulation decisions are based on monitoring traffic information. Monitoring data are generally used to update a traffic model in case of disturbances. The traffic model helps to evaluate both the impact of disturbances and the impact of regulation decisions. Reference [START_REF] Hansen | Railway Timetable and Traffic[END_REF] compiles the state-of-the-art in infrastructure and operation modeling. This section describes the formulations of infrastructure configuration and railway operations used in the presented work.

The presented traffic model is designed to be part of a supervision module, like an ATS module. It is, to some extent, based on the model previously introduced in [START_REF] Pochet | Supervision and rescheduling of a mixed CBTC traffic on a suburban railway line[END_REF]. It has been adapted to be integrated in the ATS module of the SIMONE simulation tool. In this way, the traffic model is designed for traffic supervision and has been integrated into the ATS module of a modelling of a railway system. It is composed of of a description of the tracks based on a node-link scheme, a description of train movements based on arrival and departure events, and a description of the signaling system based on a set of traffic constraints.

A. Description of tracks

The description of tracks is the base structure of the traffic model; traffic events will be monitored and simulated according to this base structure. Tracks description is based on a node-link scheme, adapted to the data structure of the SIMONE tool. The density of nodes describing a track section depends on the required level of accuracy. Thus, it is possible to find a trade-off between the size of the model and the precision in critical areas.

The presented model considers two types of nodes. Stopping nodes are located on tracks, in a station or on a sidetrack. A stopping node is specific to a train and models its target functional stopping node. It is characterized by a station (or a trackside), a platform, and a position along the platform. As a result, a station hosts several stopping nodes. If stopping nodes are computed in real-time, to take into account the traffic state, the platform and position attributes can be modified in real-time.

On-line nodes are located on tracks, outside of stopping areas. On-line nodes can be used in a timetable to model a milestone where the reference schedule should be respected. On-line nodes can also be used in a timetable to force a route else than the nominal route.

B. Timetables

To supervise traffic in real time, the ATS module uses timetables based on the structure of the link-node scheme. A timetable is composed of a collection of triplets: node identifier, arrival time and departure time, for each train.

The supervision functions use three types of timetable: the reference, the monitoring and the evaluation ones. The reference timetable consists of the objective schedule followed in nominal operations. The monitoring timetable reflects how the traffic actually evolves during operations. The ATS updates the monitoring schedule with arrival and departure data coming from the on-ground ATC for CBTC-trains, and with track circuit states coming from the signaling system for all trains (CBTC and non-CBTC trains). Then the future steps of train missions are shifted in consequence. The evaluation timetable is used to evaluate solutions to the rescheduling problem.

Moreover, the ATS checks, through the evaluation timetable, that solutions do not break traffic constraints.

C. Traffic constraints

The simulation environment comprises precise models of the signaling system and the CBTC equipment that are in charge of ensuring safety. These models reproduce the real architecture and functional sequence. To evaluate solutions to the rescheduling problem, the ATS needs to simulate the traffic and to take into account safety requirements. In order to limit the computational time, the ATS module includes a more abstract, and less complex, model of safety based on a set of constraints. The signaling system and the CBTC functions are not explicitly described in the supervision traffic model (it is modeled in the SIMONE simulation environment). When the ATS needs to evaluate the feasibility of a solution, it checks traffic constraints at the nodes of the node-link scheme, using the evaluation timetable.

1) Order Constraint

In the presented model, the ATS uses a First Come First Served (FCFS) strategy. The order constraint is checked at every node. It confirms that arriving order is compliant with departing order at previous nodes.

2) Capacity Constraint

The capacity constraint is checked at every station. For each station, and for each train, it confirms that the arrival time is compliant with the station capacity and the departure times of previous trains.

3) Headway constraint

The headway constraint is checked at every on-line node. For each train crossing the considered on-line node, it confirms that the arrival time is compliant with the departure time of the previous train according to the minimum headway.

The minimum headway between two trains depends on their communication mode (communicating mode / noncommunicating mode) and their driving mode (automatic pilot mode / manual driving mode). Non CBTC-trains are considered to be in non-communicating mode and in manualdriving mode, at any time.

4) Running time constraint

The running time and dwell time constraints define the minimum values for running times and dwell times [START_REF] Hansen | Railway Timetable and Traffic[END_REF].

The running time constraint is not explicitly defined in the ATS model. Instead, the functions that are able to modify running times select the corresponding values in a dictionary of possible values. The possible values of running times correspond to the different running types available in the onboard ATC. A running type consists in a coefficient which multiplies the maximum speed of a train on a link of the nodelink scheme. For example, we consider that the nominal running type correspond to a multiplying coefficient of 0.9. The running times dictionary is computed off-line, for each train according to its features, and for each link of the nodelink scheme crossed by the considered train.

Based on the presented traffic model, the ATS can evaluate traffic indicators, simulate solutions to disturbances and take decisions. Regulation will be discussed in section V.

V. ATS REGULATION METHOD

The presented work aims to develop a regulation method for a suburban railway line, partially equipped with a CBTC system. This method is designed to be part of the ATS module of a CBTC system. It is expected to supervise all trains in the CBTC area and to reschedule CBTC trains, while taking into account the specific principles of suburban operations. Thus, the proposed method is based on the state-of-the-art of railway traffic management as well as on the state-of-the-art of CBTC systems. Problem formulation and solution are adapted to the specific features of CBTC mixed traffic management on a suburban railway. The presented regulation method includes two main modules: a train management method in multiplatform stations, and a rescheduling method for CBTC trains using a multi-objective optimization. A first version of the rescheduling method was presented in [START_REF] Pochet | Supervision and rescheduling of a mixed CBTC traffic on a suburban railway line[END_REF].

A. Train management method in multi-platform stations

In order to make full use of platforms in stations, the proposed regulation method includes a train management method in multi-platform stations. This method is composed of three sub-functions that can be independently activated.

1) Platform choice according to expected occupation

When the ATS module defines its timetables and its nodelink scheme based on the reference data, the platform attributes of stopping nodes of trains are not computed. The ATS computes the platform choice only when the considered CBTC train is in the previous station. This computation takes into account a set of weights that are predefined in the reference timetable, as well as the expected occupation of platforms on the ATS monitoring schedule.

2) Delay of route command in stations

After computing a platform choice for a CBTC train, the ATS computes the best route to reach the platform. Then, the command is sent to the signaling system as soon as the trigger area corresponding to the route is occupied. Generally this trigger area consists of the platform area. Once a route is commanded to the signaling system and set, it is quite long and complex to undo this decision in favor of another train with another route. It is due to the approach locking [START_REF] Schoen | Railway Signaling and automation[END_REF]. In order to allow flexible departure order at station, the ATS delays the sending of commands to the signaling system. It sends command just before departure. The time of sending before departure is defined by an ATS parameter.

3) Delay of platform choice

When a CBTC train is stopped in a station, it may be blocked for various reasons (for example change of departure order or disturbances). As a result, its departure time is delayed; it will be shifted in the ATS monitoring timetable. In order to take into account the latest occupation prediction at the next station, the platform choice is computed just before departure. The time of computation before departure is defined by an ATS parameter.

B. Rescheduling through multi-objective optimization

In case of disturbances, the ATS module executes a rescheduling method to adapt train runs in the aim of improving some traffic indicators. In the presented work, the rescheduling method is based on a model predictive control approach (MPC) and its receding horizon principle. It adapts the running types of CBTC trains; Non-CBTC trains are not controlled by the ATS, but are taken into account in the supervision functions, and in the evaluation functions.

The rescheduling method uses a multi-objective centralized genetic algorithm to compute the control inputs. This algorithm is partly based on the one presented and tested in [START_REF] Pochet | Supervision and rescheduling of a mixed CBTC traffic on a suburban railway line[END_REF]. The solutions are evaluated according to 2 objectives: distance to objective headway and punctuality at the end of the CBTC area. We also consider distance to traffic constraints satisfaction as a third traffic indicator. It is used as a strong penalty applied to both objective values. The considered traffic constraints were presented in section IV.

In the genetic algorithm, a chromosome is a list of running types. It is built on the basis of the monitoring timetable. For each train, according to its position on the node-link scheme, the first link on which the running type can be modified is computed. From this first link to the end of the optimization horizon, each running type corresponding to a step of the nodelink scheme is part of the chromosome.

The selection process uses a Pareto ranking process. The random mutation process uses values generated from the running times dictionary mentioned in section IV. The crossover process uses a 2-points-crossover operator.

The proposed regulation method is designed to combine CBTC functionalities with suburban operating principle. It includes a train management method in multi-platform stations to make use of available capacity in a flexible way. It also includes a rescheduling method which is executed in case of disturbances (Fig. 1) and provides a global response to disturbances, for CBTC trains in the CBTC area and according to the suburban operation objectives.

VI. EXPERIMENTS

The regulation method, described in section V, has been tested in the SIMONE simulation tool described in section III.

In the tool SIMONE, a scenario is composed of: a precise description of tracks, the corresponding signaling system model (described in section III), vehicles parameters, CBTC parameters, reference operation data and a set of predefined disturbances. Reference operation data consist of a reference timetable and reference weights for platform choice computation.

In the presented experiments, the scenarios take place on the future extended version of the suburban line "RER E" in Paris region. The description of tracks exactly corresponds to the real track map, from the sidetracks of station "Nanterre-La-Folie" (Fig. 2) to the sidetracks of station "Rosa Parks". This 18km-long area corresponds to the CBTC control area and includes 6 stations. The signaling system model is also based on the real system configuration. The vehicle parameters are based on validated vehicle description data. CBTC parameters are based on the state-of-the-art of existing CBTC systems.

The proposed method is tested on a set of 4 train paths in each direction. The purpose of these tests is to evaluate how much a disturbed situation can be solved or at least mitigated. In order to analyze the results, the method is compared with a common basic strategy (denoted individual delay strategy) which consists in individually adapting the running types of trains according to their own delays. The individual delay strategy was implemented into the ATS module of the CBTC system in the SIMONE tool. Both compared methods include the traffic management method in stations; they only differ on their rescheduling approach.

The disturbances in scenarios represent the second train being blocked in a station, up to 60 seconds over the expected dwell time. Regulation methods are evaluated according to two indicators: the total delay and the total distance to the objective headway. For scenarios with a blocking time less than 40 seconds, the proposed method shows equivalent results to the individual delay strategy. It proposes a compromise between the two objectives. The individual delay strategy is able to partly recover the primary delay. The latter is faintly spread. For scenarios with a blocking time over 40 seconds, the proposed method has the ability to outperform the individual delay strategy. Table I. presents experimental results of two scenarios. It shows the average improvement of each regulation method, for both indicators, in comparison with the same scenario without speed regulation. The proposed rescheduling method provides a global response to disturbances and takes both indicators in consideration.

These experiments points to the usefulness of multiobjective optimization in the traffic management functions for CBTC railway lines. Multi-objective optimization is able to provide good solutions to disturbed situations. In the case of complex situations, where delays spread across the line, it outperforms the individual delay strategy by better taking into account the system features. In future work, it would be interesting to adapt the regulation method to the monitored situation. The ATS could host several regulation methods and use the more relevant one, depending on the extent of disturbances.

VII. CONCLUSION

This paper has presented the design of automatic train supervision functions for a railway line equipped with a CBTC system. Our work is supported by a new simulation tool SIMONE, developed by SNCF Reseau. This tool precisely describes the main agents of the railway system and the CBTC system. By means of the SIMONE simulation tool, we developed an ATS module, including a traffic management method in multi-platform stations and a rescheduling method based on a multi-objective optimization algorithm. The proposed rescheduling method was tested in SIMONE. Experiments showed an interest in designing CBTC regulation method using a multi-objective optimization, according to suburban railway principle.

Fig. 1 .

 1 Fig. 1. Simplified functional sequence of the ATS module.

Fig. 2 .

 2 Fig. 2. Sidetracks and platforms of station "Nanterre-La-Folie" in SIMONE

TABLE I .

 I

		EXPERIMENTAL RESULTS
	Scenario	Regulation method Punctuality Headway
	Blocking time = 40s,	Individual strategy	21%	4%
	in "Saint-Lazare" station	Proposed method	11%	16%
	Blocking time = 60s,	Individual strategy	-2%	4%
	in "Saint-Lazare" station	Proposed method	33%	23%
	Blocking time = 60s,	Individual strategy	1%	-2%
	in "La Defense" station	Proposed method	37%	26%