Zohra Kader

Adnane Saoud

Antoine Girard

Safety controller design for incrementally stable switched systems using event-based symbolic models

In this paper, we investigate the problem of lazy safety controllers synthesis for event-based symbolic models of incrementally stable switched systems with aperiodic time sampling. First of all, we provide a novel event-based scheme for symbolic models design. The obtained symbolic models are computed while considering all transitions of different durations satisfying a triggering condition. In addition, they are related to the original switched system by a feedback refinement relation and thus useful for control applications. Then, using the particular structure of the obtained event-based symbolic model, a lazy safety controller is designed while choosing transitions of longest durations. Secondly, for the same state sampling parameter and desired precision, we show that the obtained event-based symbolic model is related by a feedback refinement relation to the classical symbolic model designed for incrementally stable switched systems with periodic time sampling. Based on this relationship, we prove analytically that the size of the set of controllable states obtained with the lazy safety controller designed for an event-based symbolic model is larger than the one obtained with a safety controller designed for the classical symbolic model. Finally, an illustrative example is proposed in order to show the efficiency of the proposed method and simulations are performed for a Boost DC-DC converter structure.

I. INTRODUCTION

Switched systems have attracted a wide interest of the control community during the last decades [START_REF] Liberzon | Switching in systems and control[END_REF]. The heterogeneous nature of this class of hybrid systems renders their study more complex and thus limits the investigated problems to stability and stabilization [START_REF] Shorten | Stability criteria for switched and hybrid systems[END_REF], [START_REF] Hetel | Local stabilization of switched affine systems[END_REF]. Recently, several approaches based on the use of symbolic models, also called discrete abstractions, for controller design have been proposed [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]. These studies have been motivated by technology advances which demand that more complex control objectives like safety properties, language and logic specifications be considered. Moreover, the use of discrete abstraction for control design becomes more interesting when the obtained symbolic model is finite. Indeed, in this case, the problem of controller design can be efficiently solved *This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 725144). This research has been partially supported by Labex DigiCosme (project ANR-11-LABEX-0045-DIGICOSME) operated by ANR as part of the program "Investissement d'Avenir" Idex Paris Saclay (ANR-11-IDEX-0003-02). 1 {adnane.saoud}@l2s.centralesupelec.fr using the mature methods obtained for supervisory control design for discrete-event systems.

Based on the Lyapunov theory, several constructive approaches of symbolic models for incrementally stable switched systems have been proposed last years. For instance, we can cite the work proposed in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] where a symbolic model has been designed using both state and time discretization. In that paper, the obtained symbolic model is related to the original system by an approximate bisimulation relation.

Here, we are interested in lazy safety controllers synthesis for symbolic models of switched systems with aperiodic time sampling. First of all, we provide a novel event-based scheme for symbolic models design. Contrarily to the event-based symbolic approach proposed in [START_REF] Kader | Symbolic models for incrementally stable switched systems with aperiodic time sampling[END_REF] where only transitions of shorter durations are considered, here the symbolic model is computed while considering all transitions of different durations satisfying a triggering condition. In this paper, we show that the obtained symbolic abstraction is related to the original switched system by a feedback refinement relation [START_REF] Reissig | Feedback refinement relations for the synthesis of symbolic controllers[END_REF] and is thus useful for control applications. Then, using the particular structure of the obtained event-based symbolic model a lazy safety controller that keeps the trajectory of the closed-loop system in the safe set while choosing transitions of longest durations is designed. Secondly, for the same sampling state parameter and desired precision, we show that the obtained event-based symbolic model is related by a feedback refinement relation to the classical symbolic model designed for switched systems with periodic time sampling. Based on this relationship, we provide an analytic proof of the fact that the set of controllable states obtained with a safety controller designed for the classical symbolic model with a periodic time sampling is included in the one obtained with the lazy safety controller designed for the event-based symbolic model. This result has not been shown in the literature.

This paper is structured as follows. In Section II, the class of switched systems under study is described and some required preliminaries and definitions are provided. A novel event-based scheme for symbolic models design for incrementally stable switched systems is proposed in Section III. In addition, a feedback refinement relation from the obtained event-based symbolic model to the classical symbolic model designed for incrementally stable switched systems with periodic time sampling is provided in the same section. In Section IV, a lazy safety controller is designed for the event-based symbolic model. Moreover, using the result obtained in Section III, we provide analytic proof of the fact that the set of the controllable states obtained with the lazy safety controller is larger than the one obtained with a safety controller designed for the classical symbolic model. In Section V, an illustrative example is proposed in order to show the efficiency of the proposed method and simulations are performed for a Boost DC-DC converter structure. Section VI ends the paper with concluding remarks. Proofs of all theorems and lemmas can be find in the Appendices.

Notations.: In this paper we use the notations R, R + 0 and R + to refer to the set of real, non-negative real, and positive real numbers, respectively. Z, N, and N + refer to the sets of integers, of non-negative integers and of positive integers, respectively. card(S ) refers to the cardinal of a set S . x denotes the Euclidean norm of a vector x ∈ R n and x (i) refers to its i-th row. A continuous function γ is said to belong to class K if it is strictly increasing and γ(0

) = 0. It is said to belong to class K ∞ if γ is K and γ(r)
goes to infinity as r tends to infinity. A continuous function β : R + 0 × R + 0 → R + 0 is said to belong to class K L if : for any fixed r, the map β (., s) belongs to the class K , and for each fixed s the map β (r, .) is strictly decreasing and β (r, .) goes to zero as s tends to infinity.

II. PRELIMINARIES AND PROBLEM STATEMENT

A. System description

In this paper we consider the class of switched systems defined as follows:

Definition 1: A switched system is a quadruple Σ = (R n , P, P, F), where:

• R n is the state space;

• P is the finite set of modes P = {1, . . . , m};

• P is a subset of S(R + 0 , P) which denotes the set of piecewise constant and right continuous functions σ from R + 0 to the finite set of modes P, with a finite number of discontinuities on every bounded interval of R + 0 . This guarantees the absence of Zeno behaviors. • F = { f 1 , . . . , f m } is a collection of vector fields indexed by P. The continuous subsystems Σ p of the switched system Σ are defined by the following differential equation:

ẋ(t) = f p (x(t)), ∀p ∈ P. ( 1
)
We assume that for all p ∈ P the vector field f p : R n → R n is locally Lipschitz continuous map and forward complete. Under this assumption, solutions of (1) are unique and defined for all t ∈ R + 0 . Necessary and sufficient conditions for forward completeness of a system have been provided in [START_REF] Angeli | Forward completeness, unboundedness observability, and their lyapunov characterizations[END_REF]. From now on, x(t, x, σ ) will denote the point reached by the trajectory of Σ at time t ∈ R + 0 from the initial state x under the switching signal σ .

B. Incremental stability of switched systems

In order to construct symbolic models for switched systems, we use the notion of incremental stability [START_REF] Angeli | A Lyapunov approach to incremental stability properties[END_REF], [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. Loosely speaking, incremental stability means that independently of their initial states, all the trajectories induced by the same switching signal converge to the same reference trajectory. In [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], the notion of incremental stability of switched systems has been characterized using Lyapunov functions as follows:

Definition 2: A smooth function V : R n × R n → R + 0 is a common δ -GUAS Lyapunov function for system Σ if there exist K ∞ functions α, α and κ ∈ R + such that for all x, y ∈ R n , for all p ∈ P α( xy ) ≤ V (x, y) ≤ α( xy );

(2)

and ∂V ∂ x (x, y) f p (x) + ∂V ∂ y (x, y) f p (y) ≤ -κV (x, y). ( 3
)
In the sequel we assume that: A-1 There exists a common δ -GUAS Lyapunov function V :

R n × R n → R + 0 for subsystems Σ p ; A-2 There exists a K ∞ function γ such that ∀x, y, z ∈ R n , |V (x, y) -V (x, z)| ≤ γ( x -y ). (4)
Assumption A-2 has been already used in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. It has been shown that this assumption is satisfied provided that the dynamics of the switched system are considered on a compact set S ⊂ R n and the Lyapunov function V is of class C 1 on S . In this case, we have

∀x, y, z ∈ S , |V (x, y) -V (x, z)| ≤ c x -y , (5)
with c = max

x,y∈S ∂V ∂ y (x, y) . Thus, ( 4) is verified for linear K ∞ function given by γ(s) = cs. Moreover, note that for all x ∈ R n we have V (x, x) = 0, then

V (x, y) ≤ |V (x, y) -V (x, x)| ≤ γ( x -y ), (6)
for all x, y ∈ R n . Thus, considering that the right inequality in (2) holds with α = γ does not induce any loss of generality.

C. Transition systems

In what follows, we recall the concept of transition systems that allows to describe both switched systems and symbolic models in the same framework: Definition 3:

A transition system is a tuple T = (Q,U, O, ∆) where Q is a set of states, U is a set of inputs, O is a set of outputs, and ∆ ⊆ Q ×U × Q × O is a transition relation. T is said to be metric if the set of outputs O is equipped with a metric d such that d(o 1 , o 2 ) = o 1 -o 2 ,
symbolic if Q and U are finite or countable sets. In this paper, (q ′ , o) ∈ ∆(q, u) will denote the transition

(q, u, q ′ , o) ∈ ∆.
This means that under the input u the trajectory of the transition system starting from the state q will evolve to the state q ′ while providing the output o. Here, we consider that the set of initial states for the transition system coincide with the set of states Q. Given a state q ∈ Q, an input u ∈ U is said to belong to the set of enabled inputs, denoted by Enab ∆ (q), if ∆(q, u) = / 0. A state q ∈ Q is said to be blocking if Enab ∆ (q) = / 0, it is said non-blocking otherwise. T is said to be deterministic if for all q ∈ Q and for all u ∈ Enab ∆ (q), card(∆(q, u)) = 1, nonblocking if all its states are non-blocking. In this paper, only deterministic transition systems are considered.

The dynamics of the switched system (1) can be described by a transition system. Indeed, to the switched system Σ = (R n , P, P, F) we can associate the transition system

T (Σ) = (Q,U, O, ∆) where Q = R n is the set of states, U = P × R + is the set of labels, O is the set of outputs, ∆ is the transition relation defined as follows: ∀q, q ′ ∈ Q, ∀(p, τ) ∈ U, ∀o ∈ O, (q ′ , o) ∈ ∆(q, (p, τ)) if and only if x(τ, q, p) = q ′ and o = q,
i.e., when the mode p is active for a duration τ, the trajectory of the switched system Σ starting from q evolves to the state q ′ . Note that the transition system T (Σ) is not symbolic (U and Q are not finite sets).

In order to design symbolic models for switched systems, we consider the ε-approximate feedback refinement relation. This relation is a simulation relation and it is suitable for controllers design.

Definition 4:

Let T i = (Q i ,U, O, ∆ i ), with i = 1
, 2 be two metric transition systems with the same input set U and the same output set O equipped with the metric d.

Let ε > 0 be a given precision. A relation R ⊆ Q 1 × Q 2 is said to be an ε-approximate feedback refinement relation from T 1 to T 2 if for all (q 1 , q 2 ) ∈ R Enab(q 2 ) ⊆ Enab(q 1 ); for all u ∈ Enab(q 2 ) ∀(q ′ 1 , o 1 ) ∈ ∆ 1 (q 1 , u), ∃(q ′ 2 , o 2 ) ∈ ∆ 2 (q 2 , u) such that d(o 1 , o 2 ) ≤ ε and (q ′ 1 , q ′ 2 ) ∈ R; and for all q 1 ∈ Q 1 ∃q 2 ∈ Q 2 , such that (q 1 , q 2 ) ∈ R.
When ε = 0, R is said to be a feedback refinement relation.

D. Safety controller synthesis

Let the transition system T = (Q,U, O, ∆) and Q s ⊆ Q be a safe set. We consider the synthesis problem that consists in determining a controller that keeps the states of the system inside the set of safety specification Q s . For the system T and the set Q s , a state q ∈ Q s is controllable if there exists an infinite sequence of transitions of T initialized in q and remaining in Q s for all time. This set can be formally defined as follows.

Definition 5: A state q of a deterministic transition system T is controllable with respect to the safe set Q s if q ∈ Q s and for all r ∈ N + , there exists a sequence of inputs u 0 , . . . , u r with (q i+1 , o i+1 ) ∈ ∆(q i , u i ) for all i ∈ {1, . . . , r}, where q 0 = q and q i ∈ Q s for all i ∈ {1, . . . , r}. The set of controllable states is denoted Cont(Q s ). Consider a transition system T = (Q,U, O, ∆), a controller for T is a map C : Q ⇒ U such that for all q ∈ Q, C (q) ⊆ Enab ∆ (q). We define the domain of the controller as dom(C ) = {q ∈ Q | C (q) = / 0}. The controlled transition system T /C is defined by the tuple T /C = (Q ∩ dom(C ),U, O, ∆ C ), where the sets U and O are inherited from the transition system T , and the transition relation is given by:

(q ′ , o) ∈ ∆ C (q, u) iff (q ′ , o) ∈ ∆(q, u) and u ∈ C (q). (7)
Definition 6: A safety controller for the transition system T and the safe set Q s satisfies:

(i) dom(C ) ⊆ Cont(Q s );

(ii) for all q ∈ dom(C ) and for all u ∈ C (q), if (q ′ , o) ∈ ∆ C (q, u), then q ′ ∈ dom(C ). There are in general several controllers that solve the safety problem. A suitable solution to the safety problem is a controller that enables as many actions as possible. This controller C * is said to be a maximal safety controller, in the sense that for any other controller C and for all q ∈ Q, we have C (q) ⊆ C * (q). Given the set of controllable states Cont(Q s ), the maximal safety controller can be defined as follows:

• for all q / ∈ Cont(Q s ), C * (q) = / 0; • for all q ∈ Cont(Q s ), C * (q) = {u ∈ Enab(q) | for (q ′ , o) ∈ ∆(q, u), q ′ ∈ Cont(Q s )}.
Let us remark that for any safety controller C we have that dom(C ) ⊆ Cont(Q s ), while for the maximal safety controller C * , we have dom(C * ) = Cont(Q s ).

III. EVENT-BASED SYMBOLIC MODELS

Let Σ = (R n , P, P, F) be a switched system for which the switching is periodically controlled with a period τ * ∈ R + . Then, a transition system T τ * can be associated to Σ by selecting all its transitions of duration τ * > 0. The transition system

T τ * = (Q τ * ,U, O, ∆ τ * ) is defined by: • Q τ * = R n is the set of states; • U = P × {τ * } is the set of inputs; • O = R n is the set of outputs; • the transition relation ∆ τ * ⊆ Q τ * ×U × Q τ * × O is given as follows: ∀x, x ′ ∈ Q τ * , ∀u ∈ U, ∀o 1 ∈ O, (x ′ , o 1 ) ∈ ∆ τ * (x, u) if and only if x ′ = x(τ * , x, p) and o 1 = x.
In this context, an approach for designing approximately bisimilar symbolic models with T τ * (Σ) has been presented in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. This method is based on the approximation of the state space by the lattice:

[R n ] η = q ∈ R n |q i = k i 2η √ n , k i ∈ Z, i = 1, . . . , n ,
where η ∈ R + is the state space sampling parameter. The quantizer

Q η : R n → [R n ] η is defined by Q η (x) = q if and only if ∀i = 1, . . . , n, q (i) - η √ n ≤ x (i) < q (i) + η √ n . (8)
It can be easily shown that for all

x ∈ R n , Q η (x) -x ≤ η. The symbolic abstraction T τ * η = (Q τ * η ,U, O, ∆ τ * η )
have been constructed as follows:

• Q τ * η = [R n ] η is the set of states; • U = P × {τ * } is the set of inputs; • O = R n is the set of outputs; • the transition relation ∆ τ * η ⊆ Q τ * η ×U × Q τ * η × O is given as follows: ∀q, q ′ ∈ Q τ * η , ∀u ∈ U, ∀o 2 ∈ O, (q ′ , o 2 ) ∈ ∆ τ * η (q, u) if and only if q ′ = Q η (x(τ * , q, p)) and o 2 = q.
It has been shown in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] that if Assumptions A-1 and A-2 are satisfied and if for a desired precision ε ≥ 0 the state space sampling parameter η is such that η ≤ γ -1 ((1e -κτ )α(ε)), then T τ * (Σ) is ε-approximately bisimilar to T τ * η (Σ). Here, we are interested in the symbolic models construction for switched systems for which the switching does not occur periodically. This can be the case when fast switching is needed. In this case we assume that the transition duration can be chosen from a finite set of durations T N τ * = { τ * N , 2τ * N , . . . , τ * } where N ∈ N + is a subsampling parameter. To the switched system Σ = (R n , P, P, F), we associate the transition system T e = (Q e ,U e , O e , ∆ e ) where:

• Q e = R n is the set of states;

• U e = P × T N τ * is the set of inputs;

• O e = R n is the set of outputs; • ∆ e ⊆ Q e ×U e × Q e × O e is the transition relation defined as follows: ∀x, x ′ ∈ Q e , ∀(p, τ) ∈ U e , ∀o ∈ O e , (x ′ , o) ∈ ∆ e (x, u) if and only if x(τ, x, p) = x ′ , o = x.
Let ε ∈ R + be the desired precision of the symbolic model. We first approximate the state space by the lattice [R n ] η where η ∈ R + , and we define the transition system T e η (Σ) = (Q e η ,U e , O e , ∆ e η ) where:

• Q e η = [R n ] η is the set of states; • U e = P × T N τ * is the set of inputs; • O e = R n is the set of outputs; • ∆ e η ⊆ Q e η × U e × Q e η × O e is the transition relation defined by: ∀q, q ′ ∈ Q e η , ∀u = (p, τ) ∈ U e , ∀o ∈ O e , (q ′ , o) ∈ ∆ e η (q, u) if and only if g(τ, q, p) ≤ 0 (9)
where g(τ, q, p) := γ( x(τ, q, p)q ′ ) ) -(1e -κτ )α(ε), (10) and q ′ = Q η (x(τ, q, p)), o = q.

One can easily check that the obtained symbolic model T e η (Σ) is deterministic.

One may remark that if the subsampling parameter is fixed N = 1, then by computing T e η one retrieve the symbolic model T τ * η . Moreover, in this case, one can provide an εapproximate bisimulation relation between T e η and T e [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. Here, we are interested in the case where N ≥ 1 such that the symbolic model allows all the transitions of durations τ ∈ T N τ * satisfying [START_REF] Pola | Integrated design of symbolic controllers for nonlinear systems[END_REF]. In this case, we are able to provide an ε-approximate feedback refinement relation from T e to T e η which is useful for control design. This is shown in the following Theorem.

Theorem 1: Consider a switched system Σ and assume that A-1 and A-2 hold. Let us consider a desired precision ε > 0 and a state sampling parameter η > 0 such that

η ≤ γ -1 ((1 -e -κτ * )α(ε)). (11)
Then, the relation

R = {(x, q) ∈ Q e × Q e η |V (x, q) ≤ α(ε)} (12)
is an ε-approximate feedback refinement relation from T e to T e η .

Note that the result of Theorem 1 is constructive. For a desired precision ε > 0 and a chosen state sampling parameter η satisfying [START_REF] Reissig | Feedback refinement relations for the synthesis of symbolic controllers[END_REF], the transitions durations can be computed numerically while computing the symbolic model since they correspond to the values of τ ∈ T N τ * for which the function g changes sign. Contrarily to the event-based scheme for symbolic models design proposed in [START_REF] Kader | Symbolic models for incrementally stable switched systems with aperiodic time sampling[END_REF] where the symbolic model is designed while choosing only transitions of shorter durations, the proposed symbolic model proposed above provides all the transitions with durations τ ∈ T N τ * satisfying ( 9)- [START_REF] Pola | Integrated design of symbolic controllers for nonlinear systems[END_REF].

The choice of the state sampling parameter η in ( 11), provide us with a useful property relating the event-based symbolic model proposed above and the symbolic model obtained with a fixed time sampling period proposed in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. This property is shown in the following Lemma.

Lemma 1: Consider a switched system Σ and assume that A-1 and A-2 hold. Let ε > 0 be a desired precision and η > 0 a state sampling parameter such that

η ≤ γ -1 ((1 -e -κτ * )α(ε)). (13)
Then, the relation

R ′ = {(q 1 , q 2 ) ∈ Q e η × Q η |q 1 = q 2 } (14)
is a feedback refinement relation from T e η to T τ * η . One may remark that a direct consequence of Lemma 1 is the fact that the event-based symbolic model T e η is nonblocking. Indeed, any transition of the transition system T τ * η is a transition of the symbolic model T e η . The result of Lemma 1 is very interesting in the sense that any controller C τ * designed for T τ * η is a controller for T e η . This is utilized to prove the main result of the next section.

IV. LAZY COMPUTATION OF SYMBOLIC SAFETY

CONTROLLERS

Motivated by the properties of the event-based symbolic model proposed above and inspired from the self-trigger control strategy where the controller determines the mode of the switched system and the duration during which the mode is active [START_REF] Anta | To sample or not to sample: Self-triggered control for nonlinear systems[END_REF], this section is dedicated to the synthesis of lazy safety controllers for event-based symbolic models. Here, using Lemma 1, we show that the size of the set of controllable states obtained with a safety controller designed for a symbolic model with a periodic time sampling is included in the set of controllable states obtained with a lazy safety controller designed for the event-based symbolic model.

The classical approach to compute the maximal safety controller C * is based on a fixed point algorithm [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]. However, the computational complexity grows exponentially with state and input spaces dimension. A lazy safety controller is a controller that keeps all trajectories of the transition system within the safe set, while applying for each state a transition of longest possible duration. For this reason, we introduce a priority relation over the set of inputs, for which we give priority to transitions of longer duration. For u 1 = (p 1 , τ 1 ), u 2 = (p 2 , τ 2 ) ∈ U, priority is given to transitions of longer durations where: u 1 u 2 if and only if τ 1 ≤ τ 2 , u 1 ≺ u 2 if and only if τ 1 < τ 2 and u 1 ≈ u 2 if and only if τ 1 = τ 2 . Using the fact that U is finite and for a given subset U ′ ⊆ U, we define:

max (U ′ ) = {u ′ ∈ U ′ | ∀u ∈ U ′ , u u ′ }. ( 15
)
First, we define a lazy safety controller. Definition 7: A lazy safety controller for the transition system T and the safe set Q s is a safety controller such that (i) for all q ∈ dom(C ), if u ∈ C (q), then for any u ′ ∈ Enab(q) with u ≺ u ′ , (q ′ , o) = ∆(q, u ′ ), it holds that q ′ / ∈ dom(C ). Secondly, let us recall the notion of maximal lazy safety controller introduced in [START_REF] Girard | Safety controller synthesis for incrementally stable switched systems using multiscale symbolic models[END_REF]:

Definition 8: A maximal lazy safety controller for the transition system T = (Q,U, O, ∆) and safety specification Q s is a safety controller C l : X ⇒ U such that:

• all controllable states are in dom(C l ):

Cont(Q s ) = dom(C l );
• for all states q ∈ dom(C l ):

1) if u ∈ C l (q), then for any u ′ ∈ Enab(q) with u ≺ u ′ , (q ′ , o) = ∆(q, u ′ ), it holds that q ′ / ∈ Cont(Q s ); 2) if u ∈ C l (q), then for any u ′ ∈ Enab(q) with u ≈ u ′ , (q ′ , o) = ∆(q, u ′ ), it holds that u ′ ∈ C l (q) if and only if q ′ ∈ Cont(Q s ).
It was shown in [START_REF] Girard | Safety controller synthesis for incrementally stable switched systems using multiscale symbolic models[END_REF], that if the set of inputs is finite and equipped with a priority relation, then there exists a unique maximal lazy safety controller. Interestingly, the domain of the maximal lazy safety controller satisfies dom(C l ) = dom(C * ) = Cont(Q s ). An algorithm for synthesizing the maximal lazy safety controller was given in [START_REF] Girard | Safety controller synthesis for incrementally stable switched systems using multiscale symbolic models[END_REF], it is based on depth first search, where transitions of higher priority are explored first. While in classical safety fixed points algorithms [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF] the abstraction needs to be pre-computed, in the lazy algorithm the abstraction is computed on-the-fly [START_REF] Pola | Integrated design of symbolic controllers for nonlinear systems[END_REF]. The maximal lazy safety controller is a compromise between permissiveness and computational complexity, and represents a suitable solution when computational resources are not sufficient to use classical safety algorithms. Given a switched system Σ and its periodic and eventbased abstractions T τ * η (Σ) and T e η (Σ), and given a safety specification Q s . Our objective is to provide a theoretical comparison between the maximal safety controller for T τ * η (Σ) and Q s , and the maximal lazy safety controller for T e η (Σ) and Q s . Interestingly, we show that the size of the set of controllable states obtained with the lazy safety controller of the event-based symbolic model is much larger than the one of the safety controller designed for the symbolic model with periodic time sampling.

Theorem 2: Let the transition systems T τ * η (Σ) and T e η (Σ) for which [START_REF] Reissig | Feedback refinement relations for the synthesis of symbolic controllers[END_REF] holds. Consider the safety specification

Q s ⊆ Q τ * η = Q e η . Let C * e : Q e η ⇒ U e
η be a maximal safety controller

for T e η (Σ) and safety specification Q s , C l e : Q e η ⇒ U e η be a maximal lazy safety controller for T e η (Σ) and safety specification Q s , and C * τ * : Q e η ⇒ U e η be a maximal safety controller for T τ * η (Σ) and safety specification Q s . Then, for all q ∈ Q e η , C * τ * (q) ⊆ C l e (q) ⊆ C * e (q). (16) A direct implication of Theorem 2, is that any transition allowed by the maximal safety controller designed for a symbolic model with periodic time sampling is also enabled by the lazy safety controller designed for an event-based symbolic model. Thus, the size of the set of controllable states obtained with a lazy safety controller designed for an event-based symbolic model is much larger compared to the one obtained with a safety controller designed for a classical symbolic model. One may remark also that we can not compare with the symbolic model with a periodic time sampling τ * N , since the systems T e η (Σ) and T

τ * N η ′ (Σ) did not have the same state space (T e η (Σ) is constructed using a discretization parameter η = γ -1 ((1 -e -κτ * )α(ε)) and T τ * N η ′ (Σ) is contructed using the discretization parameter η ′ = γ -1 ((1 -e -κτ * N )α(ε))).

V. ILLUSTRATIVE EXAMPLE

We consider the boost DC-DC converter which is mathematically modeled as a switched affine system given by:

ẋ(t) = A p(t) x(t) + b (17)
where

x(t) = i L (t) v c (t) T is the state vector with i l (t)
is the current in the inductor and v c is the voltage in the capacitor, b and A p , p ∈ {1, 2} are matrices of appropriate dimension and the numerical values of their parameters have been provided in [START_REF] Beccuti | Optimal control of the boost dc-dc converter[END_REF]. It has been shown in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] the boost DC-DC converter is an incrementally stable system. In order to design our symbolic model we consider the Lyapunov function and the parameters provided [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. We consider the time sampling parameter τ * = 0.5sec and the subsampling parameter N = 50. Let us consider a desired precision for the symbolic model as ε = 0.1 and the state sampling parameter is such that η = γ -1 ((1e -κτ * )α(ε)) = 9.7 × 10 -4 . Let the safe set of the symbolic abstraction

Q s = [R 2 ] η ∩ [1.3 1.6] × [5.6 5.8].
For the obtained symbolic model we apply a lazy safety control strategy. Simulations results are presented in Figures 2,3. Moreover, in order to illustrate the result of Theorem 2, we have designed a symbolic model with a fixed time sampling period τ * = 0.5sec and the same desired precision ε = 0.1. To the obtained symbolic model we have designed a safety controller to ensure the same control objective as for the lazy safety controller. The obtained symbolic controller is shown in Figure 1. From Figures 23, we can observe that the control objective is satisfied and the trajectories of the closed-loop system remain inside the safe set. We can see from Figures 23that when the trajectory of the closed-loop system is far from the boundary of the safe set the time sampling parameter is equal to τ * and as the trajectory get closer to the unsafe set the sampling parameter becomes smaller (τ < τ * ) in order to allow fast switching and keep the trajectory in the safe set. Comparing Figures 2 and1, we can observe that all the transitions allowed by the safety controller designed for the symbolic model with a periodic time sampling are enabled by the lazy safety controller designed for the event-based symbolic model. It can be seen clearly that the size of the set of the controllable states obtained with the event-based approach is much larger than the one obtained with the classical abstraction method. This observations are consistent with the theoretical result.

VI. CONCLUSION

This paper has provided a novel event-based scheme for symbolic models design. The obtained symbolic models have been shown to be related to the original switched system by a feedback refinement relation. Then, using the particular structure of the obtained event-based symbolic model, a lazy safety controller has been designed while choosing transitions of longest durations. We then prove analytically that the size of the set of controllable states obtained with the lazy safety controller designed for the event-based symbolic model is larger than the one obtained with a safety controller designed for the classical symbolic model has been provided. Finally, simulations have been performed for a Boost DC-DC converter structure in order to show the efficiency of the proposed method.

APPENDIX I PROOF OF THEOREM 1

Proof: In order to prove Theorem 1, we will follow the statements of Definition 4. First let us remark that for all (x, q) ∈ R, we have Enab ∆ e η (q) ⊆ Enab ∆ e (x) = U e . Thus, the first condition of Definition 4 holds.

Let us consider (x, q) ∈ R, from the left inequality in (2) we obtain α( xq ) ≤ V (x, q) ≤ α(ε).

This leads to

x -q ≤ α -1 (V (x, q)) ≤ ε. (19)
Therefore, d(x, q) ≤ ε. Now, let u ∈ Enab ∆ e η (q) and (x ′ , o) ∈ ∆ e (x, u). There exists q ′ ∈ [R n ] η such that x(τ, q, p)q ′ ≤ η. In order to demonstrate that (x ′ , q ′ ) ∈ R, it is sufficient to prove that

V (x ′ , q ′ ) ≤ α(ε). (20)
From (4), we obtain

V (x ′ , q ′ ) ≤ V (x ′ , x(τ, q, p)) + γ( x(τ, q, p) -q ′ ). (21)
Using the fact that V is a δ -GUAS Lyapunov function for the switched system Σ that satisfies (3), we obtain

V (x ′ , q ′ ) ≤ V (x ′ , x(τ, q, p)) + γ( x(τ, q, p) -q ′ ) ≤ V (x(τ, x, p), x(τ, q, p)) + γ( x(τ, q, p) -q ′ ) ≤ V (x, q)e -κτ + γ( x(τ, q, p) -q ′ ). (22) Since (x, q) ∈ R, the last inequality leads to V (x ′ , q ′ ) ≤ V (x, q)e -κτ + γ( x(τ, q, p) -q ′ ) ≤ α(ε)e -κτ + γ( x(τ, q, p) -q ′ ). (23)
From ( 9) and ( 10), (23) becomes

V (x ′ , q ′ ) ≤ α(ε)e -κτ + γ( x(τ, q, p) -q ′ ) ≤ α(ε). (24)
Therefore, (20) is verified. Thus, (x ′ , q ′ ) ∈ R. Finally, the last requirement in Definition 4 can be shown directly from the construction of the state space approximation [R n ] η . Indeed, we have that for all x ∈ Q e there exists q ∈ Q e η given by q = Q(x) such that xq ≤ η. Assuming that the right inequality of (2) holds for α = γ, we obtain V (x, q) ≤ γ( xq ) ≤ γ(η).

(25)

Using [START_REF] Reissig | Feedback refinement relations for the synthesis of symbolic controllers[END_REF], we obtain V (x, q) ≤ γ( xq ) ≤ γ(η) ≤ (1e -κτ * )α(ε) ≤ α(ε).

(26) Thus, (x, q) ∈ R, which ends the proof.

APPENDIX II PROOF OF LEMMA 1

Proof: First let us show that for all (q 1 , q 2 ) ∈ R ′ , Enab ∆ τ * η (q 2 ) ⊆ Enab ∆ e η (q 1 ). We have (q 1 , q 2 ) ∈ R ′ . Thus q 1 = q 2 = q. Let u = (p, τ * ) ∈ Enab ∆ τ * η (q). Therefore, ∆ τ * η (q, u) = / 0. Then, there exists q ′ ∈ [R n ] η such that q ′ = Q η (x(τ * , q, p)) and

x(τ * , q, p)q ′ ≤ η.

Now let us verify that g(τ * , q, p) ≤ 0 with the function g defined in [START_REF] Pola | Integrated design of symbolic controllers for nonlinear systems[END_REF] and g(τ * , q, p) = x(τ * , q, p)q ′ -(1e κτ * )α(ε).

Using (27), we obtain g(τ * , q, p) ≤ γ(η) -(1e κτ * )α(ε)

From ( 13), the last inequality leads to g(τ * , q, p) ≤ γ(η) -(1e κτ * )α(ε) ≤ 0. (28)

Therefore, (q ′ , o) ∈ ∆ e η (q, u). Thus, u = (p, τ * ) ∈ Enab ∆ e η (q), which proves the first requirements of Definition 4. Now let u = (p, τ * ) ∈ Enab ∆ τ * η (q). Consider (q ′ 1 , o 1 ) ∈ ∆ e η (q, u). Then, there exists q ′ 1 = Q η (x(τ * , q, p)) and o 1 = q. On the other hand, since u = (p, τ * ) ∈ Enab ∆ τ * η (q) there exists q ′ 2 = Q η (x(τ * , q, p)) = q ′ 1 and o 2 = q = o 1 . Therefore, (q

′ 1 , q ′ 2 ) ∈ R ′ . Moreover, d(o 1 , o 2 ) = 0.
Finally, thanks to the fact that Q τ * η = Q e η the third condition in Definition 4 is satisfied for all q ∈ Q e η .

APPENDIX III PROOF OF THEOREM 2

Proof: Let us remark that the inclusion of C l e in C * e follows directly from the fact that the maximal lazy safety controller is a safety controller.

For the first inclusion, we first prove that dom(C * τ * ) ⊆ dom(C l e ). We have from Lemma 1 that the relation R ′ defined in ( 14) is a feedback refinement relation from T e η to T τ * η . Since C * τ * is the maximal safety controller for T τ * η (Σ) and safety specification Q s and using the fact that R ′ is a feedback refinement relation [START_REF] Reissig | Feedback refinement relations for the synthesis of symbolic controllers[END_REF], we have that C * τ * is a safety controller for T e η (Σ) and safety specification Q s . Hence, for all q ∈ Q e η , C * τ * (q) ⊆ C * e (q), which implies that dom(C * τ * ) ⊆ dom(C * e ) = dom(C l e ).

We have that C * τ * is a safety controller, then conditions (i) and (ii) of Definition 6 are directly satisfied. Now let q ∈ Q τ * η , and u = (p, τ * ) ∈ Enab ∆ τ * η (q). In addition, since there are no u ′ ∈ Enab ∆ τ * η (q) such that u ≺ u ′ , the condition (i) of Definition 7 is immediately satisfied. Then, C * τ * is a lazy safety controller for T e η and Q s . Therefore, for all q ∈ Q e η , C * τ * (q) ⊆ C l e (q). Which ends the proof.

Fig. 1 .

 1 Fig. 1. Safety controller designed for a symbolic model of the boost DC-DC converter with a fixed time sampling period τ * = 0.5sec (dark gray: mode 1, light gray: mode 2, medium gray: both modes are acceptable, white: uncontrollable states).

Fig. 2 .

 2 Fig. 2. Lazy safety controller for the event-based symbolic model of the boost DC-DC converter (dark gray: mode 1, light gray: mode 2, medium gray: both modes are acceptable, white: uncontrollable states); Symbolic states of the closed-loop boost DC-DC converter (blue stars).

Fig. 3 .

 3 Fig. 3. Evolution of the state variables of the boost DC-DC converter with the lazy safety controller starting at x = [1.45,5.77] T ; The sampling instants generated while computing the lazy safety controller; Switching signal generated by the lazy safety controller


Ji Ma
email: ji.ma@recherche.enac.fr

Mohammed Sbihi
email: mohammed.sbihi@enac.fr

Daniel Delahaye
email: daniel@recherche.enac.fr

Optimization of departure runway scheduling incorporating arrival crossings

Keywords: Aircraft Sequencing Problem, ILP, Simulated Annealing

The runway is a key airport resource and an efficient runway operation is critical to enhance the airport efficiency and to reduce delays. This paper addresses the problem of scheduling aircraft departures incorporating arrival crossings. Constraints for wake turbulence separations, flight time window restrictions, and holding queue capacity at runway threshold are explicitly considered. We present two Integer Linear Programming (ILP) models and a simulated annealing (SA) algorithm. Comparison tests are conducted for the south side of Paris Charles De Gaulle Airport and show that the three proposed methods can significantly improve the solutions based on the simple first-come-first-served rule. Moreover, SA has sufficiently small computation time for real-time deployment.

Introduction

Continuous air traffic growth has led to congestions in many major airports. Limited capacity is the major cause of congestion. Surface congestion may result in significant delays and operational inefficiencies. Of all airport ground resources, the runway is the main bottleneck, as it is shared by all types of operations (arrivals, departures and crossings) indicated by field observations [START_REF] Anagnostakis | A multi-objective, decomposition-based algorithm design methodology and its application to runaway operations planning[END_REF]. Efficient runway operations planning and control is critical to enhance the airport efficiency and to reduce delays while satisfying all system constraints. One important aspect of runway operations is the departure sequencing and scheduling, termed by departure management (DMAN).

The DMAN is a tactical planning tool supporting the departure schedule. It assists in optimizing take-off time while ensuring safety by taking into account several operational constraints: the minimum wake turbulence separation between successive aircraft, the time window of runway usage for individual flights, the interactions with landing runways with regard to airport layout, etc. The minimum separation requirement is a key factor for deciding the runway capacity. The International Civil Aviation Organization (ICAO) divides aircraft into four different wake turbulence weight categories (Light, Medium, Heavy, and Super Heavy) based on the maximum certificated take-off mass. Due to wake vortex, an aircraft is perturbed by the precedent one and a minimum separation must be respected to ensure safety. The separation time depends on the categories of the leading and trailing aircraft. For example, a Heavy aircraft followed by a Light aircraft generates longer separation than a Light aircraft followed by a Heavy aircraft. The time window constraint represents the earliest and latest time of arrival at the runway threshold. The time window can be applied to both the flights with CFMU (Central flow management Unit) slots and the flights without CFMU slots to ensure fairness. CFMU slot is a take-off slot that can be assigned to a departing flight in order to avoid congestion in some particular airspace. It is defined as the period between 5 minutes before and 10 minutes after the CTOT (Calculated Take-Off Times). It can be transfered to the gate by subtracting the taxi time. For the flights without CFMU slots, the earliest and latest time of the time window represent the EOBT (Estimated Off-Block Time) and a maximum gate holding time with regard to the EOBT, respectively.

The Estimated Off-Block Time (EOBT) is the estimated time when an aircraft will be ready after all doors are closed, boarding bridge removed, and is ready to push back upon reception of clearance from controller. Next, a Target Startup Approval Time (TSAT) is the time provided by controller when an aircraft can expect to receive push back approval taking into account EOBT, and the traffic situation etc. TSAT assignment is important to reduce aircraft waiting time with engine on. If many aircraft arrive at the runway threshold waiting for take-off, a long queue is formed and leads to extra fuel burn. It is more efficient to delay departures before the engines start up so that the traffic arrives at the runway smoothly.

Air traffic controller officers are responsible for sequencing traffic and ensuring minimum separations for departures. In large airports, departure management is handled by three types of controllers: the apron controller (or ramp controller), the ground controller, and the local controller. The apron controller dictates aircraft push back time from the gate and directs aircraft into gates after landing. They also coordinate with airline operators in implementing ground delays if necessary. The ground controller handles all intermediary taxiing routes to avoid conflicting movements of aircraft, and gives the pilot instructions on reaching the runway holding area via the appropriate taxiways. The local controller clears aircraft for take-off or landing ensuring the prescribed runway separation.

Runway can be used in segregated mode (only landings or only take-offs) or in mixed mode (both takeoffs and landings). The mixed mode can achieve higher runway throughput since separation between arrival and departure is shorter. Nevertheless, segregated mode is more often used due to airport layout. Departures and arrivals are highly coupled processes with complex interactions with regard to runway configurations, weather, mix of aircraft, etc. Departing and arrival flights may cross an active runway. Furthermore, controllers delay runway-crossing clearances until a group of aircraft has accumulated at various holding points [START_REF] Anagnostakis | A multi-objective, decomposition-based algorithm design methodology and its application to runaway operations planning[END_REF]. However, there is an upper limit on how many aircraft can hold at the same crossing point. For safety reasons, arrivals always have priority over departures and crossings.

Aircraft Sequencing Problem (ASP) has been studied intensively in the past decades. Research has been conducted both in the US [START_REF] Anagnostakis | A conceptual design of a departure planner decision aid[END_REF][START_REF] Balakrishnan | Efficient and equitable departure scheduling in real-time: new approaches to old problems[END_REF][START_REF] Gupta | Incorporating active runway crossings in airport departure scheduling[END_REF] etc. and in Europe [START_REF] Böhme | Tactical departure management with the eurocontrol/dlr dman[END_REF][START_REF] Atkin | Hybrid metaheuristics to aid runway scheduling at london heathrow airport[END_REF][START_REF] Deau | Airport surface management and runways scheduling[END_REF]) etc. Bennell et al. (2011) presented an extensive overview on landing and take-off scheduling problems. [START_REF] Lieder | Scheduling aircraft take-offs and landings on interdependent and heterogeneous runways[END_REF] recently summarized different models and solution approaches on the ASP that consider heterogeneous or interdependent runways, featuring articles up to 2015.

ASP can be divided into Aircraft Landing Problem (ALP), Aircraft Take-off Problem (ATP), and a combination of both problems. Numerous physical and operational constraints as well as the layout of runway system must be considered. Various exact and heuristic approaches have been proposed. The solution quality and the computational time of the approach are two important factors. In real situations, controllers can use only algorithms which can quickly find a good solution. Optimal solutions arising from lengthy computation times are of little practical use [START_REF] Bennell | Airport runway scheduling[END_REF].

In the rest of this section, we mainly review the ATP and its interactions with the arrival aircraft. Considering the exact approach, [START_REF] Balakrishnan | Scheduling aircraft landings under constrained position shifting[END_REF], 2007, 2010) introduce a dynamic programming algorithm with Constrained Position Shift (CPS) that set a limit on the number of positions an aircraft can occupy in the sequence for the ATP based on their model for the ALP. They also provide complexity analyses of their algorithms and extend their approach for multiple runways and active runway crossings. [START_REF] Montoya | Multiobjective departure runway scheduling using dynamic programming[END_REF] formulate ATP as a multi-objective optimization problem with respect to total aircraft delay and runway throughput and use dynamic programming to solve. [START_REF] Gupta | A mixed integer linear program for airport departure scheduling[END_REF] present a Mixed Integer Linear Program (MILP) to handle the departure queuing area. A following work by the same authors [START_REF] Gupta | Incorporating active runway crossings in airport departure scheduling[END_REF] presented a MILP formulation for incorporating active runway crossings in departure scheduling. However, the computational performance was not good due to the lack of good lower bounds. Moreover, the MILP formulation only provides optimal solution for each short time frame (e.g. 15 minutes), which is addressed and improved in [START_REF] Bae | Modified mixed integer linear program for airport departure scheduling[END_REF].

Different heuristic approaches have been proposed. Atkin et al. (2008a) applied three different heuristic methods for the ATP taking into account some specific constraints concerning London Heathrow airport. Further work can be found in [START_REF] Atkin | Hybrid metaheuristics to aid runway scheduling at london heathrow airport[END_REF](Atkin et al., , 2008b[START_REF] Atkin | An examination of take-off scheduling constraints at london heathrow airport[END_REF]. [START_REF] Anagnostakis | Runway operations planning: a two-stage solution methodology[END_REF] proposed a two-stage heuristic algorithm to solve the departure and runway crossing problem. [START_REF] Mori | Development of a pushback time assignment algorithm considering uncertainty[END_REF] applied a tabu search technique to develop a pushback time assignment algorithm in order to reduce take-off delay.

Regarding comparative studies between exact and heuristic approaches, [START_REF] Furini | Improved rolling horizon approaches to the aircraft sequencing problem[END_REF] present two MILP models and a tabu search heuristic for finding solutions to the ASP in a short computation time. With a time limit of 15 seconds, the tabu search heuristic outperforms the two MILP models. [START_REF] Malik | Exact and heuristic algorithms for runway scheduling[END_REF] proposed one exact and two heuristic algorithms to solve the ATP, and they conclude that heuristics are more suitable than the exact approach to produce good quality solutions in a relatively short computation time. D' Ariano et al. (2010) compare branch and bound (BB) and heuristic algorithms, setting a time limit of two minutes for the algorithm. Their results show that BB outperforms their heuristics for instances up to 60 minutes of traffic prediction.

In this paper we focus on the departure sequencing problem incorporating arrival crossings. Different from the literature works, besides specific constraints such as wake turbulence separations, flight time window restrictions etc, we consider explicitly the holding queue capacity at runway threshold. Moreover, holding point for arrivals is introduced as a decision variable to allow simultaneous crossings. Gate and runway threshold are two control points to hold departures. Two ILP models and a metaheuristic method are proposed to compare the solution quality and the computational performance. The paper is organized as follows. In Section 2, the problem statement is presented. In Section 3, two ILP formulations are proposed in detail, followed by the metaheuristic method in Section 4. Finally, computational results and conclusions are summarized in Section 5 and Section 6 respectively. Fig. 1: Paris CDG airport runway throughput chart in July 2017. Southern pair of runways (26L/26R) is in higher demand compared to northern side (27L/27R). July 13th is excluded due to abnormal data record.

The Problem Statement

Paris CDG airport is one of the busiest passenger airports in Europe, composed of four parallel runways (two for landings and two for departures) and three terminals. Fig. 1 illustrates one month runway throughput of CDG in July 2017; it handled more than 1,400 flights/day. The southern pair of runways is in higher demand compared to the northern side due to the fact that it is closer to the main terminal (Terminal 2) of CDG. In this paper, we focus on the departure scheduling problem incorporating arrival crossings on the southern pair of runways: 26L for arrival and 26R for departure. As shown in Fig. 2, there are three main landing runway exits: V2, V3, and V4 and three holding points before crossing: S1, S2, and S3. In practice, the landing runway exit is chosen by the pilot with regard to the aircraft type and how fast it can brake after touching down. Controllers prefer to give priority to arrival aircraft to cross the departure runway without stopping when the departure runway is not fully occupied. When the departure traffic demand is high, controllers indicate arrival aircraft to wait at different holding points and to cross simultaneously. For departures, controllers can hold aircraft at the gate with engines off, or have them waiting at the runway threshold before taking off. In order to ensure safety, the minimum wake turbulence separation between successive aircraft must be respected. Moreover, several operational constraints such as time window of runway usage for individual flights, maximum runway queue capacity, and interactions between arrival crossings and departures need to be considered.

To be more precise, let F = A ∪ D be the set of aircraft (arrivals and departures) that will use the runway during the time horizon T , which is discretized into time intervals = {1, . . . T }. Flight taking off/crossing at the time interval t ∈ T refers that it occurs at the beginning of time interval t.

We need to define the following notations for our problem formulation:

i ∈ A, R a i = {1, . . . , R a i }, T ihr possible runway usage time for flight i ∈ A if i uses the holding point h, T ihr = T 0 ih + rρ, r ∈ Ri, sij
minimum separation in runway usage time of two successive flights i and j with regard to their wake turbulence categories, as shown in Table 2, M large positive constant.

Our aircraft runway scheduling consists in:

• Choosing a pushback delay d i , i.e., the difference between the EOBT and the modified pushback time, and a holding duration w i at the runway threshold for each departure i ∈ D. With these decisions, we know the arrival time at the runway threshold t T i = P i + d i + X i , and the runway usage time t i = t T i + w i ; • Choosing a holding point h i and a holding duration w i for each arrival i ∈ A. With these decisions, we In accordance with the operational practices, pushback delay and holding time are given as multiples of time slot ρ. A runway schedule is feasible if and only if for each pair of distinct flights (i, j) ∈ F × F with t i ≤ t j , we have

• t j -t i ≥ s ij and t T j ≥ t T i if (i, j) ∈ D×D, i.e.
, the minimum wake turbulence separation is respected, and the first aircraft arriving at the holding point departs first; • t j -t i ≥ s ij and t T j ≥ t T i if (i, j) ∈ A × A and h i = h j , i.e., two aircraft using the same holding point follow the FCFS order, and the minimum separation is respected;

• t j -t i ≥ s ij if (i, j) ∈ A × D ∪ D × A, i.e.
, the minimum separation is guaranteed between arrival and departure, vice versa.

Our problem consists of finding a feasible schedule while minimizing the total delay i∈D

(d i + w i ) + i∈A w i .
The subproblem of scheduling aircraft take-offs can be seen as a machine scheduling problem with release times, due dates, and sequence-dependent processing times (but zero setup times) [START_REF] Chaudhry | A research survey: review of flexible job shop scheduling techniques[END_REF]. A related problem, the simple total tardiness problem (without sequence dependent setup times), has been shown to be NP-hard [START_REF] Du | Minimizing total tardiness on one machine is np-hard[END_REF]. Consequently, our runway scheduling problem is NP-hard too.

The ILP Formulations

Two ILP formulations are given in this section. First, we present Model A that uses the number of time slots as decision. This model involves many so called "big-M" constraints that makes the LP relaxation rather weak. In order to reduce this number and to strengthen the formulation, we propose Model B termed as delay-indexed formulation.

Model A: time slot based formulation

The Decision Variables. We define the following decision variables:

d i ≡ number of pushback delay slots for departure flight i ∈ D; w i ≡ number of holding time slots for flight i ∈ F; t i ≡ runway usage time for flight i ∈ F; t T i ≡ time arriving at the runway threshold for flight i ∈ F; The Objective Function. Our objective is to minimize the total delay:

δ ij = 1, if i uses runway before j, 0, otherwise; for i ∈ F, j ∈ F, i = j; α is = 1, if i
Min i∈D (d i + w i ) + i∈A w i .
The Constraints.

h∈H

z ih = 1, ∀i ∈ A, (1)
t T i = T 0 i + ρd i , ∀i ∈ D, (2)
t T i = h∈H T 0 ih z ih , ∀i ∈ A, (3)
t i = t T i + ρw i , ∀i ∈ F, (4) δ ij + δ ji = 1, ∀i, j ∈ F, i < j, (5)
t T j -t T i ≥ -M (1 -δ ij ), ∀i, j ∈ D, i = j, (6)
t j -t i ≥ s ij -M (1 -δ ij ), ∀i, j ∈ D, i = j, (7)
t T j -t T i ≥ -M (3 -δ ij -z ih -z jh ), ∀h ∈ H, ∀i, j ∈ A, i = j, (8) t j -t i ≥ s ij -M (3 -δ ij -z ih -z jh ), ∀h ∈ H, ∀i, j ∈ A, i = j, (9) t j -t i ≥ s ij -M (1 -δ ij ), ∀(i, j) ∈ D × A ∪ A × D, (10) α is ≤ α i(s+1) , ∀i ∈ D, ∀s ∈ T \{T }, (11)
β is ≥ β i(s+1) , ∀i ∈ D, ∀s ∈ T \{T }, (12)
β ihs ≥ β ih(s+1) , ∀i ∈ A, ∀h ∈ H, ∀s ∈ T \{T }, (13)
t i - s∈T β is = 0, ∀i ∈ D, (14)
t T i - s∈T (1 -α is ) = 0, ∀i ∈ D, (15)
s∈T

β ihs ≤ M z ih , ∀i ∈ A, ∀h ∈ H, (16)
t i - h∈H s∈T β ihs = 0, ∀i ∈ A, ( 17
) i∈D (β is + α is -1) ≤ C d , ∀s ∈ T , (18)
i∈A|T 0 ih ≤s

β ihs ≤ C a , ∀s ∈ T , ∀h ∈ H, (19)
0 ≤ d i ≤ R p i , integer, ∀i ∈ D, (20) 0 ≤ w i ≤ R w i , integer, ∀i ∈ D, (21) 0 ≤ w i ≤ R a i , integer, ∀i ∈ A, (22)
t T i , integer, ∀i ∈ F, ( 23
) t i , integer, ∀i ∈ F, (24) δ ij ∈ {0, 1}, ∀i, j ∈ F, i = j, (25) α is ∈ {0, 1}, ∀i ∈ D, ∀s ∈ T , (26)
β is ∈ {0, 1}, ∀i ∈ D, ∀s ∈ T , (27) z ih ∈ {0, 1}, ∀i ∈ A, ∀h ∈ H, (28)
β ihs ∈ {0, 1}, ∀i ∈ A, ∀h ∈ H, ∀s ∈ T . (29)
Constraints (1) ensure only one holding point for each arrival. Constraints (2) link the time arriving at runway threshold and the pushback delay for departures. Constraints (3) link the time arriving at holding point and the holding point decision variables for arrivals. Constraints (4) state the runway usage time as a sum of time to reach the holding point and the holding time. Constraints (5) guarantee that given any two aircraft, one leads the other. Constraints ( 6) and ( 7) ensure the first-come-first-served (FCFS) order at the runway threshold and for take-off as well as the minimum separation for departures. There are two cases to consider here: if δ ij = 1, i.e., i uses runway before j , then Constraints (6) and Constraints (7) become t T j ≥ t T i and t j -t i ≥ s ij , ensuring that the FCFS order and the separation are respected. If δ ij = 0, j uses runway before i, then Constraints (6) and Constraints (7) become t T j -t T i ≥ -M and t j -t i ≥ s ij -M , i.e., t T j -t T i and t j -t i are bigger than some large negative values, thereby ensuring that this constraint is effectively inactive. Similarly, Constraints (8) and Constraints ( 9) guarantee that at the same holding point, the crossing order of two landing aircraft keeps the same with the order they reach the holding point, and ensure the minimum separation while crossing. Constraints (10) guarantee the minimum separation between arrival and departure (vice versa). Constraints ( 11), ( 12) and (13) ensure the connectivity in time. For example, if flight i arrives at the runway threshold before time s, then for all the time intervals s ≥ s, α is has to have a value of 1. Constraints ( 14) and Constraints (15) link the two binary variables to the runway usage time and the arrival time at runway threshold for departures. Constraints ( 16) ensure that β ihs = 0 if flight i does not use holding point h. Constraints ( 17) is similar to Constraints ( 14), that link the binary variable β ihs to the runway usage time for arrivals. Constraints (18) ensure that at every time the number of aircraft waiting at the departure runway threshold will not exceed the departure capacity. The term β is + α is -1 is equal to 1 if and only if flight i is waiting at runway threshold during the time interval s. Likewise, Constraints (19) ensure that at every time the number of crossing aircraft waiting at the holding point will not exceed the holding capacity. Notice that, contrary to the departures, we do not need to introduce variables α ihs for arrivals since the arrival time at runway threshold is known as soon as the holding point is given.

Considering the LP relaxation of the problem, we would like M to be as small as possible, and this can be achieved by tailoring M for each constraint. For example in Constraints ( 6), M can be set to max(0, T 0 i + R p i -T 0 j ). Moreover, as we know the time range where each aircraft can use the runway, we can strengthen the formulation by fixing the value of some variables:

α is = 1, for s ∈ {T iR p i , . . . , T }, 0, for s ∈ {1, . . . , T 0 i -1}; ∀i ∈ D;
β is = 1, for s ∈ {1, . . . , T 0 i }, 0, for s ∈ {T iRi + 1, . . . , T }; ∀i ∈ D;

β ihs = 0, for s ∈ {0, . . . ., T 0 ih -1}; ∀i ∈ A, h ∈ H.

Model B: delay-indexed formulation.

The Decision Variables. We define the following binary decision variables:

x ir = 1, if the take-off time of flight i is delayed by r slots, 0, otherwise; The Constraints. In order to introduce the constraints, we define the following sets to represent some infeasible assignments pairs of delays (and crossing point for arrivals). Given two distinct flights i, j ∈ F, we introduce:

for i ∈ D, r ∈ R i x ihr = 1, if i crosses via
C ij = {(r, q) ∈ R i × R j |T jq -T ir < s ij and T ir -T jr < s ji } if (i, j) ∈ D × D
the set of delay assignments which violates the separation requirements between two departure flights;

C ijh = {(r, q) ∈ R i × R a j |T jhq -T ir < s ij and T ir -T jhq < s ji } if (i, j) ∈ D × A {(r, q) ∈ R a i × R a j |(T 0 ih -T 0 jh )(T ihr -T jhq ) < 0 or (T jhq -T ihr < s ij and T ihr -T jhq < s ji )} if (i, j) ∈ A × A
the set of delay assignments which violates the separation requirements between departure and arrival or between two arrivals using the same holding point. The first condition in the set definition expresses the violation of the FCFS order for arrival holding point. Notice that it is not possible to express the FCFS order for departure queue directly in the set definition. It is ensured by Constraints ( 38) and (39).

r∈Ri

x ir = 1, ∀i ∈ D, ( 30
)
h∈H r∈R a i x ihr = 1, ∀i ∈ A, ( 31
) r∈Ri rx ir - r∈R p i ry ir ≥ 0, ∀i ∈ D, ( 32
) r∈Ri rx ir - r∈R p i ry ir ≤ R w i , ∀i ∈ D, (33)
x ir + x jq ≤ 1, ∀(r, q) ∈ C ij , ∀i, j ∈ D, i < j, (34) x ir + x jhq ≤ 1, ∀(r, q) ∈ C ijh , ∀h ∈ H, ∀i ∈ D, ∀j ∈ A, (35) x ihr + x jhq ≤ 1, ∀(r, q) ∈ C ijh , ∀h ∈ H, ∀i, j ∈ A, i < j, ( 36
) δ ij + δ ji = 1, ∀i, j ∈ D, i < j, ( 37
) r∈Rj T jr x jr - r∈Ri T ir x ir ≥ -M δ ij , ∀i, j ∈ D, i = j, ( 38
)
r∈R p j T jr y jr - r∈R p i T ir y ir ≥ -M δ ij , ∀i, j ∈ D, i = j, (39)
α is ≤ α i(s+1) , ∀i ∈ D, ∀s ∈ T \{T }, (40) β is ≥ β i(s+1) , ∀i ∈ D, ∀s ∈ T \{T }, (41) β ihs ≥ β ih(s+1) , ∀i ∈ A, ∀h ∈ H, ∀s ∈ T \{T }, ( 42
) r∈Ri T ir x ir - s∈T β is = 0, ∀i ∈ D, ( 43
)
r∈R p i T ir y ir - s∈T (1 -α is ) = 0, ∀i ∈ D, ( 44
)
r∈R a i T ihr x ihr - s∈T β ihs = 0, ∀i ∈ A, ∀h ∈ H, ( 45
) i∈D (β is + α is -1) ≤ C d , ∀s ∈ T , ( 46
) i∈A β ihs ≤ C a , ∀s ∈ T , ∀h ∈ H, (47)
x ir ∈ {0, 1}, ∀i ∈ D, ∀r ∈ R i , (48) x ihr ∈ {0, 1}, ∀i ∈ A, ∀h ∈ H, ∀r ∈ R a i (49) y ir ∈ {0, 1}, ∀i ∈ D, ∀r ∈ R p i , (50)
δ ij ∈ {0, 1}, ∀i, j ∈ D, i = j, (51) α is ∈ {0, 1}, ∀i ∈ D, ∀s ∈ T , (52)
β is ∈ {0, 1}, ∀i ∈ D, ∀s ∈ T , (53)
β ihs ∈ {0, 1}, ∀i ∈ A, ∀h ∈ H, ∀s ∈ T , (54)
Constraints ( 30) and ( 31) require exactly one time period of runway usage for departure and arrival flights. Constraints (32) guarantee that for departures, the arrival time at the runway threshold is always smaller than the take-off time. Constraints (33) state the maximum holding time for departures. Constraints (34), Constraints (35),and Constraints (36) ensure the separation criteria between two departures, arrival and departure, and two arrival crossings respectively. Constraints (37) guarantee that given any two aircraft, one leads the other. Constraints (38) and Constraints (39) guarantee the same order waiting at the runway threshold and for take-off for departures. are the same as Constraints (11-13) in Model A. Constraints (43) and Constraints (44) link the binary variables for departures. Constraints ( 45) is similar to Constraints (43), that link the binary variables for arrivals. Constraints (46-47) are the same as Constraints (18-19) in Model A.

One disadvantage of this discrete time formulation is a relatively large number of variables (in particular, Constraints (34-36)). In order to reduce the number of variables, we reformulate Constraints (34) and Constraints (35) to Constraints (55),and Constraints (36) to Constraints (56),

x ir + j∈D,j =i q|(r,q)∈Cij x jq + h∈H j∈A q|(r,q)∈Cijh x jhq ≤ 1 + M (1 -x ir ), ∀r ∈ R i , ∀i ∈ D, (55)
x ihr + j∈A,j =i q|(r,q)∈Cijh

x jhq ≤ 1 + M (1 -x ihr ), ∀r ∈ R a i , ∀h ∈ H, ∀i ∈ A. (56)
There are two cases to consider: if x ir = 1, then the take-off time of flight i is delayed by r slots, Constraints (55) become j∈D,j =i q|(r,q)∈Cij x jq + h∈H j∈A q|(r,q)∈Cijh x jhq ≤ 0, ensuring that i is separated with respect to all the other flights. If x ir = 0, then Constraints (55) become j∈D,j =i q|(r,q)∈Cij x jq + h∈H j∈A q|(r,q)∈Cijh x jhq ≤ 1 + M , which inactivates the constraints. The same reasoning applies to Constraints (56).

The decision variables of Model A can be linked to the decision variables of Model B using the following equations (57-60):

w i = r∈Ri ix ir , ∀i ∈ D, (57)
d i = r∈R p i iy ir , ∀i ∈ D, (58)
w i = h∈H r∈R a i ix ihr , ∀i ∈ A, (59)
z ih = r∈R a i x ihr , ∀i ∈ A, ∀h ∈ H, (60)
It is possible to remove all the "big-M" constraints in Model B by introducing the decision variables x irq instead of x ir and y ir , where

x irq = 1,
if the pushback time of flight i is delayed by r slots and flight i holds q slots at the runway threshold, 0, otherwise;

for i ∈ D, r ∈ R p i , q ∈ R w i
However, experimental tests showed that the benefit is counterbalanced by the large increase of the number of variables and constraints.

After describing our two ILP formulations, in the next section we consider also the metaheuristic method to tackle the aircraft runway scheduling problem.

Simulated Annealing

In this section we present a simulated annealing algorithm to solve our aircraft runway scheduling problem.

Simulated Annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]) is a meta-heuristic that simulates the annealing of a metal, in which the metal is heated up and slowly cooled down to move towards an optimal energy state. SA is applied to many optimization problems, in particular the job shop scheduling problem (Cruz- et al., 2017). In SA, the objective function to be minimized is analogous to the energy of the physical problem. Algorithm 1 illustrates the SA procedure, where the parameters are • S 0 : the initial solution;

Algorithm 1 Simulated Annealing Require: Initial state S 0 , T 0 ; S c ← S 0 , T ← T 0 , i ← 0; while T > T f do for i = 0 to N T do Generate a neighborhood solution, S N ; if E(S N ) ≤ E(S C ) then S C ← S N ; else if random(0,1) < exp -(E(SN )-E(SC )) T then S C ← S N ; end if end for T ← βT ; end while return S c

Chávez

• S C : the current solution;

• S N : the neighborhood solution generated from S C ; • T 0 : the initial temperature; • T : the current temperature; • T f : the final temperature; • E: the SA fitness function;

• N T : the total number of transitions in each temperature step; • β: the geometrical temperature reduction coefficient.

A global parameter T is used to simulate the cooling process. A current solution may be replaced by a random "neighborhood" solution with a probability. We start the cooling process from a high initial temperature T 0 (which can be determined by a heating process or defined by user), the current solution changes almost randomly at a high temperature, thus the algorithm is able to escape from local minima. The decrease of temperature may follow different laws, such as linear law, geometric law, etc. At each temperature step, a number of transitions are executed. The probability to accept a degrading solution become smaller and smaller when T decreases. Therefore, at the final stages of the annealing process, the system will converge to a near-global or global optimum.

In order to apply SA to our problem, the neighborhood generation is critical and needs to be defined precisely. Our neighborhood choice is similar to a Traveling Salesman Problem with a constrained search space. The algorithm starts with an initial sequence sorted by the earliest runway usage time. First, we randomly choose one flight, and we search for a list of sequence positions that can be used by the current chosen flight without exceeding the maximum delay. Then, we choose one flight in the list and apply three exchange strategies with the current one: swap, inversion, and circular shift. As shown in Fig. 3, swap move exchanges the positions of two aircraft, inversion move inverts the order of aircraft between two positions, and circular shift move changes the first aircraft to the final position, while shifting all other aircraft to the previous position. Note that it is necessary to relax the maximum pushback delay and maximum holding time constraints in order to allow SA to search more freely in the state space. Nevertheless, a penalty coefficient is added if one aircraft cannot meet the maximum delay constraints. At the end of the procedure, SA can always find a solution satisfying the maximum delay requirement. A comparison test with one hour traffic is made to find the most efficient neighborhood generation method. As indicated in Fig. 4, the strategy combining circular shift move and swap move can achieve less delay than other strategies. Thus, in the remaining part of the paper, we choose circular shift move and swap move as our neighborhood generation method. Fig. 4: Comparison of different neighborhood methods on solution qualities. We run the proposed SA algorithm 10 times with the same input parameters for each method. The error bar represents the minimum delay, the average delay, and the maximum delay respectively.

After the previous step, we obtained a new sequence, and we now have to update the flight schedules: First, we set the earliest runway usage time of the first flight in the sequence as its actual runway usage time. Then, for each aircraft, its runway usage time is set with respect to its earliest available time and the minimum separation with its precedent aircraft. Moreover, the best holding point for arrivals is determined. More precisely, we calculate for each holding point the required waiting time for the flight and choose the holding point with the shortest waiting time. Table 3 lists some SA parameters related to our problem.

Computational results

In this section we present the results of some computational experiments. We test our methodology on both real data cases at Paris CDG Airport and randomly generated data. Numerical results with different settings of (user-defined) algorithm parameters are presented and discussed. The overall process is run on a 2.50 GHz core i7 CPU, PC under Linux operating system, implemented in the Java programming language. The ILP models were solved with GUROBI solver [START_REF] Optimization | Gurobi optimizer reference manual[END_REF] with the default parameters settings. Table 4 lists some parameters related to our mathematical formulation. The actual pushback time from historical data is considered as the EOBT in this paper. In our test cases, all the flights are without CFMU slots due to data unavailability. In the current system, controllers schedule aircraft take-offs in a FCFS order with regard to their earliest arrival time at runway threshold. In this study, we set a FCFS sequence as the baseline based on the earliest arrival time at runway threshold. We assume that all the arrival flights use the closest holding point after they land, i.e., S3. There is no limit on the maximum holding number and the maximum waiting time at the runway threshold. Moreover, a landing aircraft is assumed to cross the departure runway without any holding time. In a real-world case, arrivals do have priority on departures, as departures can be held on the ground. The arrival holding area is limited in terms of capacity and should be vacated as early as possible.

To illustrate how we calculate the FCFS sequence and compare with the optimized sequence, we give an example by generating a random instance consisting of 15 aircraft (5 arrivals and 10 departures). As shown in Table 5, given the initial landing time or initial pushback time, we calculate the earliest arrival time at runway threshold, t T i . For departures, t T i = P i + X i ; For arrivals, t T i = L i + O i + T e,h . We (b) Optimal sequence with a preference to wait at the holding point.

(c) Optimal sequence with preference to wait at the gate.

Fig. 6: Comparison of the FCFS sequence, optimal sequence with a preference to wait at the holding point, and optimal sequence with preference waiting at the gate. holding time. The departure queue is reduced to 3 aircraft. In Fig. 6c, with a preference of holding at gate, aircraft taxi to the runway threshold and take off smoothly. Next, we generate some random data set to check if our methods can be suitable for dense traffic scenarios. The mix of arrivals and departures, and the fleet mix of Medium and Heavy aircraft are two important factors in departure scheduling. Thus, we generate 15 random high traffic demand scenarios in a time window of 30 minutes. The number of arrivals/departures are set to be 10/20, 15/15, and 15/20. A mix of Heavy aircraft from 0% to 40% is set. Results are shown in Table 9: the solution quality decreases when the number of departures increases; SA keeps the shortest computational time; Model A is more suitable for cases with low departure rate.

In conclusion, Model B is suitable for finding an optimal solution while Model A is suitable for finding a feasible solution within a short computation time. In fact, by imposing a time limit, Model A solves most instances and sometimes reaches optimality. However, when the runway is in high demand, Model A is not able to find a good quality solution. The SA is suitable for finding a near-optimal solution in a short computational time.

Conclusions

The runway system is a major source of delay in the departure process It is critical to achieve efficient scheduling of aircraft taking into account specific operational constraints. In this paper, we have developed two ILP models and one metaheuristic algorithm for departure runway scheduling incorporating

Fig. 2 :

 2 Fig. 2: CDG south side runway layout. Runway 26L is for landing and 26R is for taking off. The arrival aircraft have to cross the departure runway to reach the taxi area. Three arrival holding points before crossing 26R: S1, S2, and S3. Three main runway exits: V2, V3, and V4.

Fig. 3 :

 3 Fig.3: Illustration of three neighborhood generation methods: swap, i.e., exchanging the positions of two aircraft; inversion, i.e., inverting the order of aircraft between positions; circular shift, i.e., moving the first aircraft to the final position, while shifting all other aircraft to the previous position.

  (a) FCFS sequence. Each line represents one flight, the four points from top to bottom for departures indicate respectively initial pushback time, actual pushback time, arriving time at the runway threshold, and take-off time; the three points from top to bottom for arrivals indicate respectively landing time, arriving time at the runway threshold, and crossing time. The departure waiting queue is shown at the bottom of the figure.

Table 1 :

 1 Average time from runway exit of 26L to holding point, T e,h (in seconds). The average values are calculated based on one-month real traffic data.

	S1	S2	S3
	V2 104	70	59
	V3 118 102	75
	V4 158 142 115

Table 2 :

 2 Single-runway separation requirements for departures and arrival crossings (in seconds).

	Category	Trailing Aircraft g Heavy Medium Light Crossing
		Heavy	96	120	120	60
	Leading	Medium	60	60	60	60
	Aircraft f	Light	60	60	60	60
		Crossing	40	40	40	10
	know the arrival time at holding point t T

i = L i + O i + d ei,hi , and the runway usage time t i = t T i + w i .

Table 3 :

 3 Empirically-set parameter values of SA

	Parameter	Value
	Geometrical temperature reduction coefficient, δ	0.95
	Number of transitions at each temperature step, NT	100
	Initial temperature, T0
	Final temperature, T f	0.0001*T0

Table 4 :

 4 Chosen parameter values specifying the optimization problem.

	Parameter

Table 5 :

 5 Example of FCFS sequence and optimized sequence with a random instance of 5 arrivals and 10 departures. H refers to Heavy, M refers to Medium, C refers to Crossing. All the times are in seconds.

Table 8 :

 8 Comparison of heuristics for solving one-hour real traffic on February 18th, 2016. Best results are in bold.

	Time Window	Number of flights All Dep. Arr.	Gap (%) FCFS Model A Model B	SA	Model B Optimum	Computation time (s) Model A Model B	SA
	6:00-7:00	41	19	22	61.23	11.79	-	0.64	86.73	10.00	10.00	1.51
	7:00-8:00	48	13	35	71.16	0.00	53.67	3.81	14.57	5.35	10.00	2.31
	8:00-9:00	34	17	17	57.49	0.00	0.00	0.00	8.20	1.25	8.20	1.05
	9:00-10:00	58	35	23	85.31	41.09	-	6.55	353.74	10.00	10.00	3.28
	10:00-11:00	48	21	27	79.81	8.46	-	1.81	106.17	10.00	10.00	2.40
	11:00-12:00	42	20	22	26.98	1.57	-	0.00	51.01	10.00	10.00	1.68
	12:00-13:00	42	30	12	33.09	7.31	-	8.48	124.70	10.00	10.00	1.35
	13:00-14:00	37	14	23	89.82	0.00	0.00	1.96	6.22	1.36	6.22	1.33
	14:00-15:00	41	23	18	47.05	0.00	69.09	2.54	24.38	6.80	10.00	1.40
	15:00-16:00	40	20	20	45.77	0.00	-	0.99	16.02	5.39	10.00	1.47
	16:00-17:00	38	10	28	58.86	0.00	4.21	0.00	11.51	4.76	10.00	1.49
	17:00-18:00	46	25	21	84.14	9.09	-	0.00	25.20	10.00	10.00	1.99
	18:00-19:00	45	19	26	69.45	0.47	-	0.61	80.24	10.00	10.00	2.02
	19:00-20:00	46	26	20	73.77	5.26	-	2.17	34.19	10.00	10.00	1.89
	20:00-21:00	41	21	20	74.54	0.00	-	2.36	28.71	10.00	10.00	1.54

Acknowledgments

This work has been partially supported by CAUC, by China Scholarship Council (CSC) and by the Natural Science Foundation of Tianjin grant 17JCYBJC43100. We would like to thank Serge Roux and Andrea Hakala for helpful discussions and feedback on the manuscript. We would like to thank SNA-RP/CDG-LB for providing the traffic data. The authors thank the anonymous reviewers for their insightful comments.

Data

assume that all the arrival flights use the closest holding point after they land, i.e., S3. Then, we build a FCFS sequence by sorting t T i and we obtain the runway usage time t i with regard to the separation requirements. After that, the optimized sequence is presented in Table 5. We observed a decrease of total delay from 4340 s in the baseline case to 2950 s in the optimized case. We chose a heavy traffic summer day, July 11th, 2017 as our test case. An optimization time window of one hour is applied from 8:00 to 20:00; the rest of the periods involve only a few flights. As illustrated in Table 6, real traffic data provide various mixes of arrivals and departures: three fully loaded periods of 8:00-9:00, 11:00-12:00, and 19:00-20:00 with more than 30 departures per hour per runway, and 14:00-15:00 with 30 arrivals per hour per runway are observed.

Table 6 shows the comparison of optimization results for four cases: FCFS, Model A, Model B, and SA. We list the gap, i.e., the percent increase in total delay over the optimal solution. The computation time in seconds is listed in the last four columns. For SA, the gap and the computation time are averaged over 10 times random runs. In order to be fair and practical in the real world application, we set the computational time limit to be 10 seconds for all the algorithms. As shown in Table 6, first we can see a significant reduction of delays from Model A and SA compared to FCFS sequence. Model B cannot find a feasible solution in most of cases within the time limit. Model A can find a near-optimal solution. Nevertheless, it is difficult to prove the optimality due to the poor lower bound when the number of aircraft is large and the traffic scenario is dense. Therefore, we still use Model B to find the optimal solution with a reasonable time from 12 seconds up to 156 seconds, as shown in the first column of computational times in Table 6. In contrast, SA can find good quality solutions with most of the gaps less than 10 % within 5 seconds.

In order to test the performance of algorithms facing large instances, we set our time window to be two hours and the algorithm running time limit to be 20 seconds. The number of flights ranges from 65 to 103, which is challenging to find a good solution in a short period of time. As shown in Table 7, first we observe a large gap between FCFS solution and optimal solution. After optimization, SA can still find a near-optimal solution within the time limit compared to other models. In 18:00-20:00, the computational time is 384 seconds for Model B, even Model A has a gap of 39 % compared to optimality. This proves that SA is more suitable and more practical for dealing with large instances. We test another heavy traffic day on February 18th, 2016. A departure queue of 9 aircraft with 3 arrivals waiting at the holding point at 9:30 is illustrated in Fig. 5. As shown in Table 8, in 9:00-10:00, there are 58 aircraft with 35 departures and 23 arrivals.

The optimization results are listed in Table 8. SA still keeps a similar performance compared to the first test. Model A also shows very good performance to find the optimal solution within 10 seconds, except in 9:00-10:00 with a gap of 41%. Note that 9:00-10:00 is the period in highest demand of the day. Below, we analyze in detail the optimization results of time window 9:00-10:00. Fig. 6 shows the sequencing results in 9:00-10:00 for three cases: the FCFS sequence, the optimal sequence with a preference to wait at the holding point, and the optimal sequence with a preference to wait at the gate. One can observe a large departure delay as well as a long departure queue for FCFS case illustrated at the bottom of the Figure 6a. While in Fig. 6b, much less departure delays are achieved with slight modification of arrival arrival crossings. Specific constraints such as wake turbulence separations, flight time window restrictions, and holding queue capacity at runway threshold are explicitly considered. Different comparison tests on real and random data were launched and demonstrate that: the time slot based formulation can find an optimal or near-optimal solution quickly, but it is hard to prove optimality facing high demand periods; the delay-indexed formulation can prove optimality in a reasonable time; SA is a good candidate for having a near-optimal solution to reduce flight delays, and the run times are sufficiently small for real-time application. The runway sequencing algorithm can be implemented as a decision support tool for controllers in actual operations. However, the benefits of the runway sequencing algorithm may be limited by the uncertainty arising from real operations. For real-time application, flight information needs to be updated frequently and accurately as time passes. For example, the taxi time to the runway threshold can be varied with regard to the current aircraft position at the airport after pushback. A proper prediction of the travel time on the surface can help the tactical runway sequencing algorithms better organize the sequence and decrease flight delays. Future extensions of the model could develop a robust algorithm facing uncertain pushback time and taxi time changes.


Daniel Coggia
email: daniel.coggia@inria.fr

Alain Couvreur
email: alain.couvreur@inria.fr

Dga

On the security of a Loidreau's rank metric code based encryption scheme

Keywords: Rank metric codes, Gabidulin codes, code based cryptography, cryptanalysis

We present a polynomial time attack of a rank metric code based encryption scheme due to Loidreau for some parameters.

Introduction

To instantiate McEliece encryption scheme, one needs a family of codes with random looking generator matrices and an efficient decoding algorithm. If the original proposal due to McEliece himself [START_REF] Mceliece | A Public-Key System Based on Algebraic Coding Theory[END_REF] relies on classical Goppa codes endowed with the Hamming metric, one can actually consider codes endowed with any other metric. The use of F q m -linear rank metric codes, first suggested by Gabidulin et. al. [START_REF] Gabidulin | Ideals over a noncommutative ring and their applications to cryptography[END_REF] is of particular interest, since the F q m -linearity permits a very "compact" representation of the code and hence permits to design a public key encryption scheme with rather short keys compared to the original McEliece proposal.

Compared to the Hamming metric world, only few families of codes with efficient decoding algorithms are known in rank metric. Basically, the McEliece scheme has been instantiated with two general families of rank metric codes, namely Gabidulin codes [START_REF] Delsarte | Bilinear forms over a finite field, with applications to coding theory[END_REF][START_REF] Gabidulin | Theory of codes with maximum rank distance[END_REF] and LRPC codes [START_REF] Gaborit | Low rank parity check codes and their application to cryptography[END_REF].

In [START_REF] Loidreau | A new rank metric codes based encryption scheme[END_REF], Loidreau proposed the use of codes which can somehow be regarded as an intermediary version between Gabdiulin codes and LRPC codes. These codes are obtained by right multiplying a Gabidulin code with an invertible matrix whose entries are in F q m and span an F qsubspace of small dimension λ. This approach can be regarded as a "rank metric" counterpart of BBCRS scheme [START_REF] Baldi | Enhanced public key security for the McEliece cryptosystem[END_REF] in Hamming metric.

In the present article, we explain why the case λ = 2 and dim C pub n/2 is weak and describe a key recovery attack in this situation.

Note. The material of the present article has been communicated to Pierre Loidreau in april 2016. The article [START_REF] Loidreau | A new rank metric codes based encryption scheme[END_REF] is subsequent to this discussion and proposes parameters which avoid the attack described in the present article.

Prerequisites

Rank metric codes

In this article m, n denote positive integers and q a prime power. A code of dimension k is an F q m -subspace of F n q m whose dimension as an F q m -vector space is k. Given a vector x ∈ F n q m , the rank weight or rank of x, denoted as |x| R is the dimension of the F q -vector sub-space of F q m spanned by the entries of x. The support of a vector x ∈ F n q m , denoted supp(x) is the F q -vector space spanned by the entries of x. Hence the rank of x is nothing but the dimension of its support. The rank distance or distance of two vectors x, y ∈ F n q m is defined as

d R (x, y) def = |x -y| R .
Given a code C ⊆ F n q m , the minimum distance of C is defined as

d min (C ) def = min x∈C \{0} {|x| R } .

q-polynomials and Gabidulin codes

A q-polynomial or a linear polynomial is an F q m -linear combination of monomials X, X q , X q 2 , . . . , X q s , . . . Such a polynomial induces a function F q m → F q m which is F q -linear. The q-degree of a q-polynomial P , denoted by deg q (P ) is the integer s such that the degree of P is q s . In short:

P = deg q P i=0 p i X q i , p i ∈ F q m , p deg q P = 0.
The following very classical result is crucial in what follows.

Proposition 1. Let P ∈ F q m [X] be a q-polynomial. Then, the set of roots of P in F q m is an F q -vector space of dimension less than or equal to deg q P .

The space of q-polynomials is denoted by L and, given a positive integer s, the space of q-polynomials of degree less than s is denoted by

L <s def = P ∈ L | deg q P < s .
Given positive integers k, n with k n m and an n-tuple a = (a 1 , . . . , a n ) of F q -linearly independent elements of F q m , the Gabidulin code G k (a) is defined as

G k (a) def = {(f (a 1 ), . . . , f (a n )) | f ∈ L <k }.
Such codes are known to have minimum distance n -k + 1 and to benefit from a decoding algorithm correcting up to half the minimum distance (see [START_REF] Loidreau | A Welch-Berlekamp like algorithm for decoding Gabidulin codes[END_REF]).

The n-tuple a is referred to as the support of the code. Note that the support is not unique as shown by the following lemma which will be useful for our attack.

Lemma 1. Let α ∈ F q m . Then G k (a) = G k (αa).

The component-wise Frobenius map

In what follows, we will frequently apply the component-wise Frobenius map or its iterates to vectors or codes. Hence, we introduce the following notation. Given a vector v ∈ F n q m and a nonnegative integer s, we denote by v [s] the vector:

v [s] def = (v q s 1 , . . . , v q s n ).
Similarly, given a code C ⊆ F n q m and a positive integer s, the code

C [s] denotes the code C [s] def = {c [s] | c ∈ C }.

Overbeck's distinguisher

In [START_REF] Overbeck | Structural attacks for public key cryptosystems based on Gabidulin codes[END_REF], Overbeck proposes a general framework to break cryptosystems based on Gabidulin codes. The core of his attack is that a simple operation permits to distinguish Gabidulin codes from random ones. Indeed, given a random code C ⊆ F n q m of dimension k < n/2, the expected dimension of the code C + C [1] equals 2k and, equivalently C ∩ C [1] is likely to be equal to 0. More generally, we have the following statement.

Proposition 2. If C rand is a code of length n and dimension k chosen uniformly at random, then for a nonnegative integer a and for a positive integer s < k, we have

P dim F q m C rand +C rand [1] +• • •+C rand [s]
min(n, (s+1)k)-a = O(q -ma ).

On the other hand, for a Gabidulin code, the behaviour with respect to such operations is completely different as explained in the following statement.

Proposition 3. Let a ∈ F n q m be a word of rank n, k n and s be an integer. Then, [1] = G k-1 (a [1] );

G k (a) ∩ G k (a)
G k (a) + • • • + G k (a) [s] = G k+s (a).

Loidreau's scheme

In order to mask the structure of Gabidulin codes and to resist to Overbeck's attack, Loidreau suggested in [START_REF] Loidreau | A new rank metric codes based encryption scheme[END_REF] the following construction. Denote by G a random generator matrix of a Gabidulin code G k (a). Fix an integer λ m and an F q -vector subspace V of F q m of dimension λ. Let P ∈ GL(n, F q m ) whose entries are all in V. Then, let

G pub def = GP -1 .
We have the following encryption scheme.

Public key: The pair (G pub , t) where

t def = n-k 2λ .
Secret key: The pair (a, P ). Encryption: Given a plaintext m ∈ F k q m , choose a uniformly random vector e ∈ F n q m of rank weight t. The ciphertext is

c def = mG pub + e.
Decryption: Compute, cP = mG + eP .

Since the entries of P are all in V then, the entries of eP are in the product space supp(e)

• V def = uv | u ∈ supp(e), v ∈ V Fq .
The dimension of this space is bounded from above by tλ n-k 2 . Therefore, using a classical decoding algorithm for Gabidulin codes, one can recover m.

3 A distinguisher when λ = 2

Context

The goal of this section is to establish a distinguisher for Loidreau's cryptosystem instantiated with λ = 2 and a public code C pub of dimension k n 2 . Similarly to Overbeck's attack, this distinguisher reposes on Propositions 2 and 3.

Similarly to the attacks of BBCRS system [START_REF] Couvreur | Distinguisher-based attacks on public-key cryptosystems using Reed-Solomon codes[END_REF][START_REF] Couvreur | A polynomial-time attack on the BBCRS scheme[END_REF], it is more convenient to work on the dual of the public code because of the following lemmas. Lemma 1]). Any full-rank generator matrix H pub of C ⊥ pub can be decomposed as H pub = H sec P T where H sec is a parity-check of the Gabidulin code G k (a).

Lemma 2 ([10, Page 52]). The code G k (a) ⊥ is a Gabidulin code G n-k (b) for some b ∈ F n q m of rank n. Lemma 3 ([4,
The convenient aspect of the previous lemma is that the matrix P has its entries in a small vector space, while its inverse has not.

The case λ = 2

We suppose in this section that the vector space V ⊆ F n q m in which the matrix P has all its entries has dimension 2:

λ = dim Fq V = 2.
Note that, w.l.o.g, one can suppose that 1 ∈ V. Indeed, if V is spanned over F q by α, β ∈ F q m \ {0}, then one can replace H sec by H sec = αH sec which spans the same code and P = α -1 P has entries in V = Span 1, α -1 β and H pub = H sec P T .

Thus, from now on, we suppose that V = Span{1, γ} for some γ ∈ F q m \ F q . Consequently, P T can be decomposed as

P T = P 0 + γP 1 ,
where P 0 , P 1 are square matrices with entries in F q . For convenience, we suppose from now on that P 0 and P 1 are both invertible. Note that this actually holds with a high probability. If one of these matrices was not invertible, then the attack could probably be performed after minor adjustments.

We have seen that C ⊥ sec = G n-k (a) for some a ∈ F n q m with |a| R = n. We define g def = aP 0 and h def = aP 1 .

Lemma 4. The code C ⊥ pub is spanned by g + γh, g [1] + γh [1] , . . . ,

g [n-k-1] + γh [n-k-1] .
We can now state a crucial result.

Theorem 1. The dual of the public code satisfies:

dim F q m C ⊥ pub + C ⊥ pub [1] + C ⊥ pub [2] 2 dim F q m C ⊥ pub + 2.
As a conclusion, thanks to Proposition 2, we deduce that C ⊥ pub is distinguishable in polynomial time from a random code as soon as 2k-2 > n.

The attack

In this section, we derive an attack from the distinguisher defined in Section 3. In what follows, we suppose that λ = 2 and the public code has rate larger than 1/2 so that the distinguisher introduced in Section 3 works on it. Recall that C ⊥ pub = G n-k (a)P for some a ∈ F n q m whose entries are F q -independent and P is of the form P 0 + γP 1 for P 0 , P 1 ∈ M n (F q ) and γ ∈ F q m \ F q . Finally recall that

g def = aP 0 and h def = aP 1 .
In addition, we make the following assumptions:

(1) P 0 , P 1 ∈ GL(n, F q );

(2) γ is not contained in any subfield of

F q m ; (3) dim C ⊥ pub + C ⊥ pub [1] + C ⊥ pub [2] = 2(n -k) + 2.
Assumption (1) has already been discussed in § 3.2. Assumption ( 2) is reasonable in order to prevent against possible attacks based on an exhaustive search of γ. Finally, Assumption ( 3) is what typically happens according to our experiments using Magma [START_REF] Bosma | The Magma algebra system I: The user language[END_REF]. The aim of the attack is to recover the triple (γ, g, h), or more precisely, to recover a triple (γ , g , h ) such that

C ⊥ pub = g [i] + γ h [i] | i = 0, . . . , n -k -1 . (1)
Actually, the triple (γ, g, h) is far from being unique and any other triple satisfying (1) permits to decrypt messages (see further § 4.3). Let us describe an action of PGL(2, F q ) on such triples.

Proposition 4. Let a, b, c, d ∈ F q such that ad -bc = 0 and δ ∈ F q m such that γ = aδ+b cδ+d . Then, the triple (δ, dg + bh, cg + ah) satisfies (1).

Step 1: using the distinguisher to compute some subcodes

As shown in the proof of Theorem 1, C ⊥ pub + C ⊥ pub [1] is spanned by: g + γh and g [1] , h [1] , . . . , g [t] , h [t] and g

[t+1] + γ q h [t+1] ,
where

t def = n -k -1.
Then, by iterating intersections

(C ⊥ pub + C ⊥ pub [1] ) ∩ (C ⊥ pub [1] + C ⊥ pub [2] ) ∩ • • • ∩ (C ⊥ pub [t] + C ⊥ pub [t+1] ),
we obtain the code spanned by

g [t] + γ q t h [t] and g [t+1] + γ q h [t+1] .
Notice that Assumption (3) permits to prove that this intersection has exactly dimension 2. Applying the inverse of the t-th Frobenius, we get the code spanned by g + γh and g [1] + γ q 1-t h [1] .

Note that t = n -k -1 < n m. Hence γ q 1-t = γ q because of Assump- tion (2)
. Next, one can compute

C ⊥ pub ∩ g + γh, g [1] + γ q 1-t h [1] = g + γh (2)
and g + γh, g [1] + γ q 1-t h [1] + g + γh [1] = g + γh, g [1] + γ q 1-t h [1] , g [1] + γ q h [1] = g + γh, g [1] , h [1] .

Similarly, we compute the intersection with 1] and get

C ⊥ pub [-1] def = C ⊥ pub [m-
C ⊥ pub [-1]
∩ g + γh, g [1] , h [1] = g [1] + γ q m-1 h [1] .

Applying the inverse Frobenius to the last code, we get g + γ q m-2 h . Since, from (2), we also know g + γh , one can compute g + γh + g [1] + γ q m-1 h

[1] [-1] = g + γh, g + γ q m-2 h = g, h . (4)
Next, for any i ∈ {1, . . . , t} one can compute

C pub ⊥ ∩ g, h [i] = g [i]
+ γh [i] .

By applying the i-th inverse Frobenius to the previous result, we obtain the space g + γ q -i h for any i ∈ {1, . . . , t}. In summary, we know the spaces g + γh , g + γ q -1 h , . . . , g + γ q -t h .

In addition, from Lemma 1, the vector g is determined up to some multiplicative constant. Therefore, one can choose an arbitrary element of g + γh and suppose that this element is g + γh.

Step 2. Finding γ

In summary, the vector g + γh and the spaces g + γ q i h for any i ∈ {-1, . . . , -t} are known. To compute γ, we will use the following lemma.

Lemma 5. For i, j ∈ {1, . . . , t}, i = j, there exists a unique pair

(u ij , v ij ) ∈ g + γ q -i h × g + γ q -j h such that u ij + v ij = g + γh.
The pairs of vectors (u ij , v ij ) can be easily computed. Thus, from now on, we suppose we know them. In addition, despite γ, g + γ q -i h and g + γ q -j h are unknown, a calculation permits to show that u ij , v ij have the following expressions.

u ij = γ q -j -γ γ q -j -γ q -i •(g+γ q -i h) and v ij = γ -γ q -i γ q -j -γ q -i •(g+γ q -j h). (5)
Consider the vectors u 12 and u 13 . They are collinear since, from (5), they are both multiples of g + γ q -1 h. Therefore, one can compute the scalar α such that u 12 = α • u 13 . From (5) we deduce that γ satisfies the following relation.

γ q -2 -γ γ q -2 -γ q -1 = α • γ q -2 -γ γ q -2 -γ q -1 • (6)
Or equivalently, γ is a root of the polynomial

Q γ (X) def = (X q -X q 3 )(X -X q 2 ) -α q 3 (X -X q 3 )(X q -X q 2 ).
One can easily check that (X q -X) q+1 divides Q γ and we set

P γ (X) def = Q γ (X q -X) q+1 •
The element γ we look for is a root of P γ but actually, the forthcoming Proposition 5 provides the description of the other roots. We first need a technical Lemma. Lemma 6. Let a, b, c, d ∈ F q , let i, j be two nonnegative integers and set A(X) = X q i -X q j ∈ F q [X]. Then,

A aX + b cX + d = ad -bc (cX + d) q i +q j • A(X).
Proposition 5. The set of roots of P γ equals the orbit of γ under the action of PGL(2, F q ). Equivalently, any root of P γ is of the form aγ+b cγ+d for a, b, c, d ∈ F q such that ad -bc = 0.

Thanks to Propositions 4 and 5, we deduce that choosing an arbitrary root γ of P γ provides a candidate for γ and there remains to compute g , h providing our triple. The triple can be deduced from the knowledge of g + γh = g + γ h and the computation of g + γ q -1 h which can be proved to satisfy g + γ q -1 h = γ q -2 -γ q -1 γ -2 -γ u 12 .

End of the attack

Given the pair g , h , compute the matrix Q ∈ GL(n, F q ) such that h = g Q. Then, C ⊥ pub = G n-k (g ) • (I n + γ Q) and this representation of the dual provides all the elements necessary to decode, that is to decrypt any ciphertext.

Conclusion

We provided a distinguisher à la Overbeck for the public keys of Loidreau's scheme when λ = 2 and the public code has rate R pub 1 2 . From this distinguisher, we are able to derive a polynomial time key recovery attack. This attack can probably be extended to other values of λ when the public code rate satisfies R pub 1 -1 λ . Therefore, such parameters should be avoided in Loidreau's scheme.

Acknowledgements

The second author was funded by French ANR projects ANR-15-CE39-0013 Manta and ANR-17-CE39-0007 CBCrypt.