Olivier Finkel
email: finkel@math.univ-paris-diderot.fr

Stevo Todorčević

Automatic Ordinals

Keywords: ω-tree-automatic structures, ω n -automatic structures, ordinals

We prove that the injectively ω-tree-automatic ordinals are the ordinals smaller than ω ω ω . Then we show that the injectively ω nautomatic ordinals, where n ≥ 1 is an integer, are the ordinals smaller than ω ω n . This strengthens a recent result of Schlicht and Stephan who considered in [SS11] the subclasses of finite word ω n -automatic ordinals. As a by-product we obtain that the hierarchy of injectively ω n -automatic structures, n ≥ 1, which was considered in [FT12], is strict.

Introduction

An automatic structure is a relational structure whose domain and relations are recognizable by finite automata reading finite words. Automatic structures have very nice decidability and definability properties and have been much studied in the last few years, see [START_REF] Blumensath | Automatic structures[END_REF][START_REF] Blumensath | Finite presentations of infinite structures: Automata and interpretations[END_REF][START_REF] Khoussainov | Automatic structures: Richness and limitations[END_REF][START_REF] Rubin | Automatic Structures[END_REF][START_REF] Rubin | Automata presenting structures: A survey of the finite string case[END_REF]. Blumensath considered in [START_REF] Blumensath | Automatic Structures[END_REF] more powerful kinds of automata. If we replace automata by tree automata (respectively, Büchi automata reading infinite words, Muller or Rabin tree automata reading infinite labelled trees) then we get the notion of tree-automatic (respectively, ω-automatic, ωtree-automatic) structures. Notice that an ω-automatic or ω-tree-automatic structure may have uncountable cardinality. All these kinds of automatic structures have the two following fundamental properties. (1) The class of automatic (respectively, tree-automatic, ω-automatic, ω-tree-automatic) structures is closed under first-order interpretations. (2) The first-order theory of an automatic (respectively, tree-automatic, ω-automatic, ω-tree-automatic) structure is decidable.

A natural problem is to classify firstly automatic structures using some invariants. For instance Delhommé proved that the automatic ordinals are the ordinals smaller than ω ω , and that the tree-automatic ordinals are the ordinals smaller than ω ω ω , see [START_REF] Delhommé | Automaticité des ordinaux et des graphes homogènes[END_REF][START_REF] Rubin | Automatic Structures[END_REF][START_REF] Rubin | Automata presenting structures: A survey of the finite string case[END_REF]. Kuske proved in [START_REF] Kuske | Is Ramsey's Theorem omega-automatic? In 27th[END_REF] that the ω-automatic ordinals are the automatic ordinals, i.e. the ordinals smaller than ω ω .

In the first part of this paper we characterize the (injectively) ω-tree automatic ordinals, proving that they are the ordinals smaller than ω ω ω . This seems to be the first complete characterization of a class of ω-tree automatic structures. The proof uses some results of Niwinski [START_REF] Niwinski | On the cardinality of sets of infinite trees recognizable by finite automata[END_REF] on the cardinality of regular tree languages, a recent result of Barany, Kaiser and Rabinovich [START_REF] Bárány | Cardinality quantifiers in MLO over trees[END_REF] on (injectively) ω-tree automatic structures, the result of Delhommé on tree automatic ordinals, and some set theory.

The ω n -automatic structures, presentable by automata reading ordinals words of length ω n , where n ≥ 1 is an integer, have been recently investigated in [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF]. We show here that the injectively ω n -automatic ordinals are the ordinals smaller than ω ω n . This strengthens a recent result of Schlicht and Stephan who considered in [START_REF] Schlicht | Automata on ordinals and linear orders[END_REF] the subclasses of finite word ω n -automatic ordinals. As a by-product we obtain that the hierarchy of injectively ω nautomatic structures, n ≥ 1, which was considered in [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF], is strict.

The paper is organized as follows. In Section 2 we recall basic notions and some properties of automatic structures. We characterize the (injectively) ωtree automatic ordinals in Section 3. We consider the ω n -automatic ordinals in Section 4. Some concluding remarks are given in Section 5.

Automatic structures

When Σ is a finite alphabet, a non-empty finite word over Σ is any sequence x = a 1 a 2 . . . a k , where a i ∈ Σ for i = 1, . . . , k , and k is an integer ≥ 1. The length of x is k. The empty word has no letter and is denoted by ε; its length is 0. For x = a 1 a 2 . . . a k , we write x(i) = a i . Σ ⋆ is the set of finite words (including the empty word) over Σ.

We introduce now languages of infinite binary trees whose nodes are labelled in a finite alphabet Σ. A node of an infinite binary tree is represented by a finite word over the alphabet {l, r} where r means "right" and l means "left". Then an infinite binary tree whose nodes are labelled in Σ is identified with a function t : {l, r} ⋆ → Σ. The set of infinite binary trees labelled in Σ will be denoted T ω Σ . A tree language is a subset of T ω Σ , for some alphabet Σ. (Notice that we shall mainly consider in the sequel infinite trees so we shall often use the term tree instead of infinite tree).

We assume the reader has some knowledge about Muller or Rabin tree automata and regular tree languages. We recall that the classes of tree languages accepted by non-deterministic Muller, Rabin, Street, and parity tree automata are all the same. We refer for instance to [Tho90, [START_REF] Perrin | Infinite words, automata, semigroups, logic and games[END_REF][START_REF] Grädel | Automata, Logics, and Infinite Games: A Guide to Current Research [outcome of a Dagstuhl seminar[END_REF] for the definition of these acceptance conditions.

Notice that one can consider a relation

R ⊆ T ω Σ 1 × T ω Σ 2 × . . . × T ω Σ k
, where Σ 1 , Σ 2 , . . . Σ k , are finite alphabets, as a tree language over the product alphabet

Σ 1 × Σ 2 × . . . × Σ k . Let now M = (M, (R M i ) 1≤i≤k
) be a relational structure, where M is the domain, and for each i ∈ [1, k] R M i is a relation of finite arity n i on the domain M. The structure is said to be ω-tree-automatic if there is a presentation of the structure where the domain and the relations on the domain are accepted by Muller tree automata, in the following sense.

Definition 2.1 (see [START_REF] Blumensath | Automatic Structures[END_REF]) Let M = (M, (R M i ) 1≤i≤k ) be a relational structure, where n ≥ 1 is an integer, and each relation R i is of finite arity n i . An ω-tree-automatic presentation of the structure M is formed by a tuple of Muller tree automata (A, A = , (A i ) 1≤i≤k ), and a mapping h from L(A) onto M, such that:

1. The automaton A = accepts an equivalence relation E ≡ on L(A), and

2. For each i ∈ [1, k], the automaton A i accepts an n i -ary relation R ′ i on L(A) such that E ≡ is compatible with R ′ i , and
3. The mapping h is an isomorphism from the quotient structure

(L(A), (R ′ i ) 1≤i≤k )/E ≡ onto M.
The ω-tree-automatic presentation is said to be injective if the equivalence relation E ≡ is just the equality relation on L(A). In this case A = and E ≡ can be omitted and h is simply an isomorphism from (L(A), (R ′ i ) 1≤i≤k ) onto M. A relational structure is said to be (injectively) ω-tree-automatic if it has an (injective) ω-tree-automatic presentation.

We now recall a recent decidability result of [START_REF] Bárány | Cardinality quantifiers in MLO over trees[END_REF] about injectively ω-tree-automatic structures.

The quantifier ∃ is usual in first-order logic. In addition we can consider the infinity quantifier ∃ ∞ , and the cardinality quantifiers ∃ κ for any cardinal κ. The set of first-order formulas is denoted by FO. The set of formulas using the additional quantifiers ∃ ∞ , ∃ ℵ 0 and ∃ 2 ℵ 0 is denoted by FO(∃ ∞ , ∃ ℵ 0 , ∃ 2 ℵ 0 ). If L is a set of formulas, the L-theory of a structure M is the set of sentences (i.e., formulas without free variables) in L that hold in M.

If M is a relational structure in a signature τ which contains only relational symbols and equality, and ψ is a formula in FO(∃ ∞ , ∃ ℵ 0 , ∃ 2 ℵ 0 ), then the semantics of the above quantifiers are defined as follows.

• M |= ∃ ∞ x ψ if and only if there are infinitely many a ∈ M such that M |= ψ(a).

• M |= ∃ κ x ψ if and only if the set ψ M = {a ∈ M | M |= ψ(a)} has cardinality κ. Theorem 2.2 ([BKR09]) Let n ≥ 1 be an integer. If M is an injec- tively ω-tree-automatic structure then its FO(∃ ∞ , ∃ ℵ 0 , ∃ 2 ℵ 0 )-theory is decid- able. Moreover assume that M = (M, (R M i ) 1≤i≤k
) is a relational structure and an injective ω-tree-automatic presentation of the structure M is formed by a tuple of Muller tree automata (A, (A i ) 1≤i≤k ) with an isomorphism h from L(A) onto M. If ψ(x 1 , . . . , x p ) is a formula in FO(∃ ∞ , ∃ ℵ 0 , ∃ 2 ℵ 0 ), with free variables x 1 , . . . , x p , then the set h -1 [{(a 1 , . . . , a p ) ∈ M p | M |= ψ(a 1 , . . . , a p )}] is a regular subset of L(A) p and one can effectively construct a tree automaton accepting it.

Injectively ω-tree-automatic ordinals

We now give a complete characterization of the class of injectively ω-treeautomatic ordinals.

Theorem 3.1 An ordinal α is an injectively ω-tree-automatic ordinal if and only if it is smaller than the ordinal ω ω ω .

Proof. As usual, an ordinal α is considered as a linear order (α, <). We first asssume that the ordinal (α, <) has an injective ω-tree-automatic presentation, which is given by two Muller tree automata (A, A < ) and an isomorphism h from L(A) onto α. We now distinguish two cases.

First Case. The ordinal α is countable. Then L(A) is countable and then it is proved by Niwinski in [START_REF] Niwinski | On the cardinality of sets of infinite trees recognizable by finite automata[END_REF] that the infinite trees of the tree language L(A) can be represented by finite binary trees and that the corresponding language of finite trees is regular (accepted by a tree automaton reading labelled finite trees). Thus the structure (α, <) is not only injectively ωtree-automatic but also injectively tree-automatic, the relation < being also accepted by a tree automaton reading finite binary trees. But Delhommé proved in [START_REF] Delhommé | Automaticité des ordinaux et des graphes homogènes[END_REF] that the tree-automatic ordinals are the ordinals smaller than the ordinal ω ω ω . Therefore in that case the ordinal α is smaller that the ordinal ω ω ω .

On the other hand it is proved in [START_REF] Colcombet | Transforming structures by set interpretations[END_REF] that every tree-automatic structure is actually injectively tree-automatic and hence also injectively ω-treeautomatic. Thus every ordinal smaller that ω ω ω is also injectively ω-treeautomatic.

Thus the countable injectively ω-tree-automatic ordinals are the ordinals smaller than ω ω ω .

Second Case. The ordinal α is uncountable. Recall that Niwinski proved in [START_REF] Niwinski | On the cardinality of sets of infinite trees recognizable by finite automata[END_REF] that a regular language of infinite binary trees is either countable or has the cardinal 2 ℵ 0 of the continuum and that one can decide, from a tree automaton A, whether the tree language L(A) is countable. Since there is an isomorphism h from L(A) onto the uncountable ordinal α, the tree language L(A) and the ordinal α have cardinality 2 ℵ 0 . In particular, the ordinal α is greater than or equal to the first uncountable ordinal ω 1 . But then the initial segment of length ω 1 of α is definable by the formula φ(x) = ∃ ℵ 0 y(y < x)∨¬∃ ∞ y(y < x). We can now infer from Theorem 2.2 that one can effectively construct some Muller tree automata B and B < such that L(B) ⊆ L(A) and (B, B < ) form an injective ω-tree-automatic presentation of the ordinal (ω 1 , <). Moreover the formula ∀x φ(x) is satisfied in the structure (ω 1 , <) and this can be deduced effectively from the two tree automata B and B < of the injective ω-tree-automatic presentation of the structure. The ordinal (ω 1 , <) being uncountable and injectively ω-tree-automatic it follows as above that it has the cardinality of the continuum and thus that 2 ℵ 0 = ℵ 1 , i.e. that the continuum hypothesis CH is satisfied.

We have seen that we can construct (B, B < ), with an isomorphism h ′ from L(B) onto ω 1 , which form an injective ω-tree-automatic presentation of the ordinal (ω 1 , <). Recall that the usual axiomatic system ZFC is Zermelo-Fraenkel system ZF plus the axiom of choice AC. A model (V, ∈) of the axiomatic system ZFC is a collection V of sets, equipped with the membership relation ∈, where "x ∈ y" means that the set x is an element of the set y, which satisfies the axioms of ZFC. We often say " the model V" in-stead of "the model (V, ∈)". The axioms of ZFC express some natural facts that we consider to hold in the universe of sets, and the axiomatic system ZFC is a commonly accepted framework in which all usual mathematics can be developed. We are here supposed to live in a universe V of sets which is a model of ZFC in which the continuum hypothesis is satisfied. But we know that, using the method of forcing developed by Cohen in 1963 to prove the consistency of the negation of the continuum hypothesis, one can show that there exists a forcing extension V[G] of V which is a model of ZFC in which 2 ℵ 0 > ℵ 1 , see [START_REF] Jech | Set theory[END_REF]. Consider now the pair of Muller tree automata (B, B < ) in this new model V[G] of ZFC. By [START_REF] Niwinski | On the cardinality of sets of infinite trees recognizable by finite automata[END_REF] the tree language L(B) is still uncountable (because this can be decided by an algorithm defined in ZFC and hence does not depend on the ambient model of ZFC) and so it has the cardinality 2 ℵ 0 of the continuum. Moreover the automaton B < defines a linear order < 0 on L(B) such that every initial segment is countable (because this is expressed by the sentence ∀x φ(x) in FO(∃ ∞ , ∃ ℵ 0 , ∃ 2 ℵ 0 ) which was satisfied in the structure (ω 1 , <) in the universe V) . But one can show that this is in contradiction with the negation 2 ℵ 0 > ℵ 1 of the continuum hypothesis. Indeed we can construct by transfinite induction a strictly increasing infinite sequence (x α ) α<ω 1 in L(B) such that each initial segment I α = {x ∈ L(B) | x < 0 x α } is countable. Consider now the union I = α<ω 1 I α . The set I forms an initial segment of L(B) for the linear order < 0 and it has cardinality ℵ 1 < 2 ℵ 0 . Then there exists an element y ∈ L(B) such that x α < y for each α < ℵ 1 . But then the initial segment {x ∈ L(B) | x < 0 y} is uncountable and this leads to a contradiction. Finally we have proved that there are no uncountable injectively ω-tree-automatic ordinals.

The hierarchy of ω n -automatic structures

Let Σ be a finite alphabet, and α be an ordinal; a word of length α (or α-word) over the alphabet Σ is an α-sequence (x(β)) β<α (or sequence of length α) of letters in Σ. The set of α-words over the alphabet Σ is denoted by Σ α . The concatenation of an α-word x = (x(β)) β<α and of a γ-word y

= (y(β)) β<γ is the (α + γ)-word z = (z(β)) β<α+γ such that z(β) = x(β) for β < α and z(β) = y(β ′ ) for α ≤ β = α + β ′ < α + γ; it is denoted z = x • y or simply z = xy.
We recall now the definition and behaviour of automata reading words of ordinal length Definition 4.1 ([Woj84, Woj85, Bed96]) An ordinal Büchi automaton is a sextuple (Σ, Q, q 0 , ∆, γ, F ), where Σ is a finite alphabet, Q is a finite set of states, q 0 ∈ Q is the initial state, ∆ ⊆ Q × Σ × Q is the transition relation for successor steps, γ ⊆ P (Q) × Q is the transition relation for limit steps, and F ⊆ Q is the set of accepting states.

A run of the ordinal Büchi automaton A = (Σ, Q, q 0 , ∆, γ, F ) reading a word σ of length α, is an (α + 1)-sequence of states x defined by: x(0) = q 0 and, for i < α, (x(i), σ(i), x(i + 1)) ∈ ∆ and, for i a limit ordinal, (Inf (x, i), x(i)) ∈ γ, where Inf (x, i) is the set of states which cofinally appear during the reading of the i first letters of σ, i.e.

Inf (x, i) = {q ∈ Q | ∀µ < i, ∃ν < i such that µ < ν and x(ν) = q}
A run x of the automaton A over the word σ of length α is called successful if x(α) ∈ F . A word σ of length α is accepted by A if there exists a successful run of A over σ. We denote L α (A) the set of words of length α which are accepted by A. An α-language L is a regular α-language if there exists an ordinal Büchi automaton A such that L = L α (A).

An ω n -automaton is an ordinal Büchi automaton reading only words of length ω n for some integer n ≥ 1.

Recall that we can obtain regular ω n -languages from regular ω-languages and regular ω n-1 -languages by the use of the notion of substitution. Proposition 4.2 (see [START_REF] Finkel | Locally finite languages[END_REF]) Let n ≥ 2 be an integer. An ω n -language L ⊆ Σ ω n is regular iff it is obtained from a regular ω-language R ⊆ Γ ω by substituting in every ω-word σ ∈ R each letter a ∈ Γ by a regular ω n-1language L a ⊆ Σ ω .

We can obtain the following stronger result which will be useful in the sequel.

Proposition 4.3 Let n ≥ 2 be an integer. An ω n -language L ⊆ Σ ω n is regular iff it is obtained from a regular ω-language R ⊆ Γ ω by substituting in every ω-word σ ∈ R each letter a ∈ Γ by a regular ω n-1 -language L a ⊆ Σ ω , where for all letters a, b ∈ Γ, a = b, the languages L a and L b are disjoint.

Proof. The result follows easily from Proposition 4.2 and the two following facts: (1) The class of regular ω n -languages is effectively closed under finite union, finite intersection, and complementation, i.e. we can effectively construct, from two ω n -automata A and B, some ω n -automata C 1 , C 2 , and C 3 , such that L(C 1 ) = L(A) ∪ L(B), L(C 2 ) = L(A) ∩ L(B), and L(C 3 ) is the complement of L(A). (2) The class of regular ω-languages R ⊆ Γ ω is effectively closed under the substitutions Φ : Γ → P (X), where X ⊆ Γ is a finite alphabet and P (X) is the powerset of X.

We have defined in [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF] the notion of ω n -automatic presentation of a structure and of ω n -automatic structure which is simply obtained by replacing in the above Definition 2.1 the Muller tree automata by some ω n -automata.

Recently Schlicht and Stephan have independently considered finite word α-automatic structures which are relational structures whose domain and relations are accepted by automata reading finite α-words: they define finite α-words as words x of length α over an alphabet Σ containing a special symbol ♦ such that all but finitely many letters of x are equal to ♦. In particular they have proved the following result.

Theorem 4.4 An ordinal is finite word ω n -automatic, where n ≥ 1 is an integer, iff it is smaller than the ordinal ω ω n .

We now extend this result to the class of injectively ω n -automatic structures considered in [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF].

Theorem 4.5 An ordinal is injectively ω n -automatic, where n ≥ 1 is an integer, iff it is smaller than the ordinal ω ω n .

Proof. Notice first that Schlicht and Stephan considered only injective presentations of structures. Thus Theorem 4.4 implies that every ordinal α < ω ω n is injectively ω n -automatic.

On the other hand we proved in [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF] that every injectively ω n -automatic structure is also an injectively ω-tree-automatic structure. Thus it follows from the above Theorem 3.1 that every injectively ω n -automatic ordinal is countable.

But we now prove the following result.

Proposition 4.6 All injectively ω n -automatic countable ordinals are finite word ω n -automatic.

Proof. It suffices to show that if L ⊆ Σ ω n is a countable regular ω n -language then the ω n -words in L can be represented by finite ω n -words. This is clear for n = 1 because a countable regular ω-language L ⊆ Σ ω is of the form L = ∪ 1≤j≤n U i • v ω i , where for each integer i ∈ [1, n] U i ⊆ Σ ⋆ is a finitary language and v i is a non-empty finite word over Σ. Then the general result can be proved by an easy induction on n, using Proposition 4.3. Finally Theorem 4.5 follows from Theorem 4.4 and Proposition 4.6.

Corollary 4.7 The hierarchy of injectively ω n -automatic structures, n ≥ 1, which was considered in [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF], is strict.

Remark 4.8 Schlicht and Stephan asked whether a countable α-automatic structure is finite word α-automatic. The proof of Proposition 4.6 about ordinals extends easily to the general case of (injectively) ω n -automatic structures. Thus any countable (injectively) ω n -automatic structure is finite word ω n -automatic.

Concluding remarks

We have determined the classes of injectively ω-tree-automatic ordinals and of injectively ω n -automatic ordinals, where n ≥ 2 is an ordinal. This way we have also proved that the hierarchy of injectively ω n -automatic structures, n ≥ 1, which was considered in [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF], is strict.

The next step in this research would be to determine the classes of ωtree-automatic and ω n -automatic ordinals, in the non-injective case. Recall that every ω n -automatic structure is a Borel structure, i.e. has a Borel presentation, see [START_REF] Finkel | A hierarchy of tree-automatic structures[END_REF], and thus it follows from a result of Harrington and Shelah in [START_REF] Harrington | Counting equivalence classes for coκ-Souslin equivalence relations[END_REF] that there are no uncountable ω n -automatic ordinals. We conjecture that a countable ω n -automatic structure is actually injectively ω nautomatic. From this it would follow that the ω n -automatic ordinals are also the ordinals smaller than the ordinal ω ω n .

Acknowledgements. We thank the anonymous referees for useful comments on a preliminary version of this paper.


J Rynkiewicz

Estimating the Number of Regimes of Non-linear Autoregressive Models

Keywords: time series, switching regimes, mean square error, asymptotic statistic, models selection, multilayer perceptron

. However, the problem of selecting the number of regimes is far less obvious and hasn't been completely answered yet. When the number of regimes is unknown, identifiability problems arise and, for example, the likelihood ratio test statistic is no longer convergent to a χ 2 -distribution. In this paper, we consider models which allow the series to switch between regimes and we propose to study such models without knowing the form of the density of the noise. The problem we address here is how to select the number of components or number of regimes. One possible method to answer this problem is to consider penalized criteria. The consistency of a modified BIC criterion was recently proven in the framework of likelihood criterion for linear switching models (see Oltéanu and Rynkiewicz 2012). We extend these results to mixtures of nonlinear autoregressive models with mean square error criterion and prove the consistency of a penalized estimate for the number of components under some regularity conditions.

The model -definition and regularity conditions

Throughout the paper, we shall consider that the number of lags is known and, for ease of writing, we shall set the number of lags equal to one, the extension to l time-lags being immediate.

Let us consider the real-valued time series Y t which verifies the following model ( 1)

Y t = F θ 0 X t (Y t-1 ) + ε t ,
where • X t is an iid sequence of random variables valued in a finite space {1, ..., p 0 } and with probability distribution π 0 ;

• for every i ∈ {1, ..., p 0 }, F θ 0 i (y) ∈ F and

F = F θ , θ ∈ Θ, Θ ⊂ R l compact set
is the family of possible regression functions. We suppose throughout the rest of the paper that F θ 0 i are sublinear, that is they are continuous and there exist

a 0 i , b 0 i ∈ R 2 + such that F θ 0 i (y) ≤ a 0 i |y| + b 0 i , (∀) y ∈ R ;
• for every i ∈ {1, ..., p 0 }, (ε t ) t∈Z is an independent centered noise, independent of (Y t ) t∈Z .

We need the following hypothesis which implies, according to [START_REF] Yao | On stability of nonlinear AR processes with Markov switching[END_REF], strict stationarity and geometric ergodicity for Y t :

(HS) p0 i=1 π 0 i a 0 i s < 1
Let us remark that hypothesis (HS) does not request every component to be stationary and that it allows non-stationary "regimes" as long as they do not appear too often.

Estimation of the number of regimes

Let us consider an observed sample {y 1 , ..., y n } of the time series Y t . Then, for every observation y t , the conditional expectation with respect to the previous y t-1 and marginally in X t is

E (Y t | y t-1 ) = p0 i=1 π 0 i F 0 θi (y t-1 ) := g 0 (y t-1 )
As the goal is to estimate p 0 , the number of regimes of the model, let us consider all possible conditional expectation up to a maximal number of regimes P , a fixed positive integer. We shall consider the class of functions

G P = P p=1 G p , G p = g | g (y 1 ) = p i=1 π i F θi (y 1 ) ,
where π i ≥ η > 0, p i=1 π i = 1. For every g ∈ G P we define the number of regimes as p (g) = min {p ∈ {1, ..., P } , g ∈ G p } and let p 0 = p g 0 be the true number of regimes. We can now define the estimate p as the argument p ∈ {1, ..., P } maximizing the penalized criterion

(2) T n (p) = sup g∈Gp E n (g) -a n (p)
where

E n (g) = 1 2 n t=2 (y t -g(y t-1 )) 2
and a n (p) is a penalty term.

Convergence of the penalized mean square estimate For λ > 0, let us define the generalized derivative function :

d λ θ (Y t , Y t-1 ) = e -λ(Y t -g(Y t-1 )) 2 -e -λ(Y t -g 0 (Y t-1 )) 2 e -λ(Y t -g 0 (Y t-1 )) 2 e -λ(Y t -g(Y t-1 )) 2 -e -λ(Y t -g(Y t-1 )) 2 e -λ(Y t -g 0 (Y t-1 )) 2 2 = e -λ ( (Y t -g(Y t-1 )) 2 -(Y t -g 0 (Y t-1 )) 2 ) -1 e -λ((Y t -g(Y t-1 )) 2 -(Y t -g 0 (Y t-1 )) 2 ) -1 2 (1)
and let us define d λ θ -(x, y) = min 0, d λ θ (x, y) . Several statistical and probabilistic notions such as mixing processes, bracketing entropy or Donsker classes will be used hereafter. For parcimony purposes we shall not remind them, but the reader may refer to [START_REF] Doukhan | Mixing: properties and examples[END_REF] and [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF] for complete monographs on the subject.

The consistency of p is given by the next result, which in an extension of [START_REF] Gassiat | Likelihood ratio inequalities with applications to various mixtures[END_REF]:

Theorem 1 : Consider the model (Y t , Y t-1 ) defined by [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF]

and the penalizedlikelihood criterion introduced in (2). Let us introduce the next assumptions :

(A1) a n (•) is an increasing function of p, a n (p 1 )a n (p 2 ) → ∞ when n → ∞ for every p 1 > p 2 and an(p) n → 0 when n → ∞ for every p (A2) the model (Y k , X k ) verifies the weak identifiability assumption

p i=1 π i F i (y 1 ) = p0 i=1 π 0 i F 0 i (y 1 ) ⇔ p i=1 π i δ θi = p0 i=1 π 0 i δ θ 0 i (A3) It exists λ > 0 so that d λ θ , θ ∈ Θ is a Donsker class.

Then, under hypothesis (A1)-(A3), (HS), p → p 0 in probability.

The assumption (A1) is fairly standard for model selection, in the Gaussian case (A1) will be fulfilled by the BIC criterion. Note that the weak identification assumption (A2) does not allowed to use linear regression because the regression functions have to be linearly independents.The assumption (A3) is more difficult to check. First we note: e -λ((Yt-g(Yt-1)) 2 -(Yt-g 0 (Yt-1)) 2 ) -1 2 = e -2λ((Yt-g(Yt-1)) 2 -(Yt-g 0 (Yt-1)) 2 ) -2e -λ((Yt-g(Yt-1)) 2 -(Yt-g 0 (Yt-1)) 2 ) + 1 So, d λ θ is well defined if E e -2λ((Yt-g(Yt-1)) 2 -(Yt-g 0 (Yt-1)) 2 ) < ∞, but

(Y t -g(Y t-1 )) 2 -(Y t -g 0 (Y t-1 )) 2 = (Y t -g 0 (Y t-1 ) + g 0 (Y t-1 ) -g(Y t-1 )) 2 -(Y -g 0 (Y t-1 )) 2 = 2ε(g 0 (Y t-1 ) -g(Y t-1 ) + (g 0 (Y t-1 ) -g(Y t-1 )) 2
where ε = Y tg 0 (Y t-1 ) is the noise of the model. So, if the regression functions are bounded, d λ θ is well defined if λ > 0 exists such that e λ|ε| < ∞ i.e. ε admits exponential moments. Finally, using the same techniques of reparameterization as in Liu and Shao, 2003 or Oltéanu and Rynkiewicz, 2012, assumption (A2) can be shown to be true for mixture of MLP regression models.

Conclusion

We have proven the consistency of BIC-like criteria for estimating the number of components in a mixture of non-linear regression. This result can be shown without knowing the form of the density function of the noise, altought the weak identifiabilty assumption excludes linear regression functions. Finally, a more challenging task may be to get a more precise tuning of penalization term which, according to our result, can be chosen among a wide range of functions.


Cédric Gueguen
email: gueguen@univ-mlv.fr

Abderrezak Rachedi
email: rachedi@univ-mlv.fr

Coverage Extension Based on Incentive Scheduler for Mobile Relaying Nodes in Wireless Networks

Keywords: Incentative Scheduling, Coverage extension, Cooperation, Quality of Service, opportunistic scheduling

In this paper, we propose a new scheduler able to extend the wireless coverage by using an incentive approach for potential mobile relaying nodes. Indeed, the cost of cooperation can be expensive in terms of QoS and energy consumption which do not motivate the nodes to cooperate. Our incentive approach rewards the cooperative nodes. The percentage of cooperation is considered in the QoS management in order to incite the border nodes to cooperate and then to extend the wireless area. Moreover, the monitoring mechanism is proposed to correctly evaluate the cooperation rate of each node. The results show that not only the proposed solution allows the border nodes to cooperate without the negative impact but also enhance the QoS parameters.

I. INTRODUCTION

The basic idea of the coverage extension area in wireless networks is to increase the network connectivity without increasing the infrastructure. This is one of the main applications of the cooperative communications in wireless networks. The coverage extension issue requires the cooperation of border mobile nodes to relay the packets of neighbouring nodes which are located outside of the base-station area. Many researchers dealt with the strategies to find the optimal placement for the relayed nodes in order to guarantee the Quality of Services (QoS) [START_REF] Sadek | Distributed relay-assignment protocols for coverage expansion in cooperative wireless networks[END_REF]. Other works dealt with the optimal number of hops of relayed nodes in wireless networks [START_REF] Florea | The optimal number of hops in infrastructure-based fixed relay networks[END_REF] [START_REF] Sikora | Bandwidth and power-efficient routing in linear wireless networks[END_REF]. However, they assume that the relayed nodes by definition are fixed and cooperative which is not interesting in the case of a dynamic wireless network where the nodes freely move and potentialy selfish. No incentive approach is considered in these works. The relayed nodes must share their throughput with other neighbouring nodes which can impact their own packets transmission. In addition, the energy consumption of the relayed nodes is more important than the one of other classical nodes. They do not only transmit their own packets but also the packets of other neighbouring nodes. Therefore, the user of the potential relayed node can disable the cooperative functionality in order to keep the performance in terms of QoS only for its own transmission.

In this paper, we consider that the mobile relayed nodes are not part of the fixed wireless infrastructure. That is why the incentive strategy for potential mobile relay nodes is very important to consider in the design of the cooperation protocol. Our proposal is based on an incentive approach with QoS consideration for the mobile relayed nodes in order to extend the coverage area. This approach consists in increasing the priorities of the relayed nodes according to their cooperation rate. The idea is to reward the relayed nodes for their cooperation instead of penalized them by the cost of the cooperation. Consequently, the nodes do not have the benefits to select and act as selfish behaviour using their throughput only to transmit their own packets. Our proposed solution called "Coverage Extension based on Incentive scheduling" (CEI) is also build in a cross layer approach. Moreover incentive nodes to cooperate, the physical layer informations are used in order to take advantage of the time, frequency and multiuser diversity and maximize the system capacity close to the Shannon limit. Unlike the existing models, our solution can be widely implemented, evaluated and compared to the existing ressource allocation strategies like classical Round-Robin (RR) [START_REF] Nagle | On packet switches with infinite storage[END_REF] and acknowledged MaxSNR [5][6].

This paper is organized as follows: in Section 2, we present the related works. Section 3 describes the proposed coverage extension protocol based on the incentive approach. The fourth Section presents the obtained simulation results and their analysis. Finally, Section 5 concludes the paper and presents our future works.

II. RELATED WORK

The conventional access methods like Round Robin (RR) [START_REF] Nagle | On packet switches with infinite storage[END_REF] and Random Access (RA) are not adapted to the wireless environment and provide poor throughput. More recently intensive research efforts have been made in order to propose more efficient schedulers: the opportunistic schedulers. They preferably allocate the resources to the active mobile(s) with the most favourable channel conditions at a given time. Taking benefit of multiuser and frequency diversity in order to maximize the system throughput, all these schemes strongly rely on diversity for offering their good performances. One major scheduling algorithm has emerged and appeared in literature as the reference: Maximum Signal-to-Noise Ratio (MaxSNR).

In MaxSNR, priority is given at every scheduling event to the mobiles which have the greatest signal-to-noise ratio (SNR). It allocates the resource at a given time to the active mobile with the greatest SNR [START_REF] Wong | Multiuser OFDM with adaptive subcarrier, bit, and power allocation[END_REF] [START_REF] Wang | An OFDM-TDMA/SA MAC protocol with QoS constraints for broadband wireless LANs[END_REF]. Taking profit of multiuser and frequency diversity, MaxSNR scheduling continuously allocates the radio resource to the mobile with the best spectral efficiency. Consequently, MaxSNR strongly increases the system throughput. Dynamically adapting the modulation and coding allows one to always make the most efficient use of the radio resource and come closer to the Shannon limit.

However MaxSNR does not take into consideration other aspect than throughput and particularly MaxSNR scheduling does not manage priority in order to favour cooperative mobiles. Consequently, no reward is guaranteed to cooperative mobiles. Their supplementary energy consumption and the personal throughput loss are not compensated. These results in a severe penalty for them which do not encouraged cooperative networks and coverage extension.

III. THE INCENTIVE SCHEDULER ALGORITHM

The scheduler grants radio resources to each mobile as a function of its: channel state, currently cooperation ratio, network confidence percentage, traffic backlog. The knowledge of the channel state is supposed to be available at the receiver [START_REF] Li | Channel estimation for OFDM systems with transmitter diversity in mobile wireless channels[END_REF]. The channel attenuations are estimated by the access node based on the SNR. The CEI scheduling algorithm relies on weights that set the dynamic priorities for allocating the resource. These weights are built in order to satisfy two major objectives: to maximise system throughput and to encourage nodes to cooperate.

A. System Throughput Maximization Parameter

The CEI scheduler maximizes the system throughput in a MAC/PHY opportunistic approach. Data integrity requirements of the mobiles are enforced adapting the modulation scheme and the transmission power to the mobile specific channel state. At each scheduling epoch, the scheduler computes the maximum number of bits m k,n that can be transmitted in a time slot of subcarrier n if assigned to a mobile k, for all k and all n. This number of bits is limited by two main factors: the data integrity requirement and the supported modulation orders.

The bit error probability is upper bounded by the symbol error probability and the time slot duration is assumed to be equal to the duration T s of an OFDM symbol [START_REF] Wong | Multiuser OFDM with adaptive subcarrier, bit, and power allocation[END_REF]. The required received power P r (q, k) for transmitting q bits in a Resource Unit (RU) while keeping below the data integrity requirement BER target,k of the service flow of mobile k is a function of the modulation type, its order and the single-sided power spectral density of noise N 0 . For QAM and a modulation order M on a flat fading channel [START_REF] Proakis | Digital Communications[END_REF]:

P r (q, k) = 2N 0 3T s erf c -1 BER target,k 2 2 (M -1), (1)
where M = 2 q and erf c is the complementary error function. P r (q, k) may also be determined in practice based on BER history and updated according to information collected on experienced BER. The transmission power P k,n of mobile k on subcarrier n is upper bounded to a value P max which complies with the transmission Power Spectral Density regulation:

P k,n ≤ P max . (2)
Given the channel gain a k,n experienced by mobile k on subcarrier n (including path loss and Rayleigh fading):

P r (q, k) ≤ a k,n P max . (3)
Hence, the maximum number of bits q k,n of mobile k which can be transmitted on a time slot of subcarrier n while keeping below its BER target is:

q k,n ≤     log2   1 + 3P max × T s × a k,n 2N 0 erf c -1 BER target,k 2 2         . (4)
We further assume that the supported QAM modulation orders are limited so that q belongs to the set S = {0, 2, 4, . . . , q max }. Hence, the maximum number of bits m k,n that will be transmitted on a time slot of subcarrier n if this RU is allocated to the mobile k is:

m k,n = max {q ∈ S, q ≤ q k,n } . ( 5
)
MaxSNR based schemes allocate the resources to the mobiles which have the greatest m k,n values. This bandwidth allocation strategy maximizes the bandwidth usage efficiency but do not encourage nodes cooperation. In order to extend coverage area while preserving the system throughput maximization, a new parameter is added on m k,n which modulates these pure opportunistic resource allocation.

B. Incentative Parameter

The second major objective of the CEI is to incite nodes to participate to frame relay in order to extend the network coverage zone. This is achieved by extending the above crosslayer design to other layers. A new "Incentive Parameter" (IP k ) is introduced based on the current estimation of the cooperation ratio:

IP k = R k D k = D k + i i=0...i=K D ki D k , (6)
where R k is the global amount of data transmitted by the mobile k. It is the sum between D k , the amount of data transmitted to the mobile k for its own requirement and D ki , the amount of data transmit to the mobile k for a mobile i (then these data will be relayed to mobile i by mobile k in the relaying subframe). This information could be directly monitored by the access point, or signalled by each mobile to access point. We also define the cooperation ratio C k as the number of packets that the mobile k is ready to relay for other mobiles when it receives 100 packets for its own consumption, for exemple:

Supposing that there are always packets to relay out of the cell, the IP k will be respectively for these three cases equal to 1, 1.5 and 2. Consequently, based on the ressource allocation on IP k allows to give more priority on mobiles which cooperate to extend coverage zone with frame relaying.

C. Confidence Parameter

We assume that each mobile signals its R k and D k to the access point. Thanks to this information, the CEI scheduler will make the adequate resource allocation rewarding the mobile for its cooperation degree. However in order to block malicious mobiles which could lie on this information, we introduced the last paramater called confidence parameter. The confidence parameter T k depends on the correspondence between the announced cooperative ratio and the observed forwarding ratio. This control is done by a monitor node (in our case the AP or cluser-head). Each T k is varies between 0 and 1 included. When the access point watching on R k and D k correspond to the annonced cooperative ratio, T k is set to 1. On the contrary, when the mobile does not relay the announced amount of data for which it had previously received more priority, its T k is set to 0 for one round of scheduling in order to punish it. This assures a deterrent threat for mobile which would try to mislead the system.

D. Global CEI Algorithm Description

In the allocation process of a given time slot, the priority of a mobile with respect to another is determined by the magnitude of its CEI parameter :

CEI k,n = m k,n × R k D k × T k . ( 7
)
Based on the m k,n and IP k factor, the CEI k,n directly takes into account the channel states and the mobile behavior. Physical layer information is used with m k,n in order to take advantage of the time, frequency and multiuser diversity and maximize the system capacity. Cooperation information is exploited in a weighted system with IP k parameter that introduces dynamic priorities between mobiles for ensuring good awards to mobiles which help to extend the coverage zone. This results in an efficient scheme which guarantees a better network connectivity while avoiding tradeoff with the system capacity.

The T k paramater is an additionnal factor which allows to temperate CEI k,n value function of network confidence. Include T k parameter allows to be resistant to malicious nodes which would lie on its i i=0...i=K D ki . Thanks to this control parameter, no malicious mobiles will retransmit the packets which they have previously received supplementary priorities will be penalised and the good behavior are stimulate.

IV. PERFORMANCE ANALYSIS

In this section we evaluate by simulation the performance of the proposed CEI scheduling and we compare it with the classical Round Robin allocation and the well acknowledged MaxSNR scheduler. We consider four kinds of nodes: the first kind of nodes are selfish and do not relay packets (C k = 0%), the second kind relay a few packets with C k = 10%, the third kind is more cooperative with C k = 50% and the last kind of nodes is really network friendly with a high cooperative ratio of 100%.

A. Simulation setup

We assume that each frame is formed by 128 subcarriers and 5 time slots. The channel gain model on each subcarrier considers free space Path Loss a k and multipath Rayleigh fading α 2 k,n [START_REF] Parsons | The Mobile Radio Propagation Channel[END_REF]:

a k,n = a k × α 2 k,n . (8)
where a k is dependent on the distance between the access point and mobile k and α 2 k,n represents the flat fading experienced by mobile k on subcarrier n. α k,n is Rayleigh distributed with an expectancy equal to unity. Additionally, the maximum transmission power satisfies:

10 log 10 P max T s N 0 × a ref = 31 dB (9)
and the BER target is equal to 10 -3 . With this setting, the value of m k,n is 3 bits when α 2 k,n equals one. In order to be close to the reality, we consider that all mobiles run the same videoconference application. This demanding type of application generates a high volume of data with a high sporadicity and requires tight delay constraints which substantially complicate the task of the scheduler. The traffic load variation is done through increasing the mobile bit rate requirement of each mobile all together.

B. Delay impact

First we focus on the mean mobile packet delay provided by each scheduler function of different traffic loads keeping a special attention on their ability to encourage mobile cooperation with a low guaranteed delay. The obtained results are plotted in figure Fig. 1 with the mean throughput required by each mobile of the cell represented on the abscissa.

Figure 1(a) shows the case of RR with different cooperation ratios of nodes. We remark that the classical RR fails to promote cooperation activities. The RR fairly allocates the RUs to the mobiles without taking into account the effort of the cooperative mobile nodes which share their allocated resources with other nodes located out of the primary access point cell. Consequently, the more cooperative the nodes are, the less resources for their own transmission they have. Moreover, the RR does not take benefit of multiuser diversity which results in a bad utilization of the bandwidth and in turn, poor system throughput. Consequently, unacceptable packet delay is experienced even with relatively low traffic loads.

Figure 1(b) illustrates the obtained results in the case of MaxSNR with different cooperation ratios of nodes. We point out that even if a higher traffic load is supported with an acceptable packet delay, the cooperative nodes are not rewarded and their performance in terms of QoS are less compared to the non-cooperative nodes. Figure 1(c) shows the obtained results in the case of the proposed CEI scheduler. We remark that CEI not only encourages the nodes to cooperate but also to enhance the performance in terms of delay. When the nodes increase their cooperation ratio, the enhancement of their delay is more important. For exemple, the nodes with 100 % as cooperative ratio, have a delay inferior to 100ms when the mean required throughput is less than 3 × 10 5 bps which is not possible with other schedulers. The CEI dynamically and gradually adjusts the relative priorities of the mobiles in order to fairly and adequately reward them according to their relative cooperation ratio. With this approach, sparingly delaying the selfish mobiles, the CEI helps the others and whatever the traffic load, the mobile which provides the best cooperative ratio experiences the lowest packet delay.

C. Relay efficiency impact

Figure 2 illustrates the relay efficiency in terms of the total mean throughput that each scheduling algorithm has allowed to provide out of the cell. We remark that RR provides the worst performances compared to MaxSNR and CEI. MaxSNR allows to relay more packets but it is the CEI which gives the best number of provided throughputs out of the cell. The RR and MaxSNR curves decrease after the peaks due to the best cooperative mobile penalizing when the system capacity is reached. The CEI, according more priority to friendly mobiles, continues to increase the total amount of forwarding throughputs until a high traffic load which corresponds to a high network extension capacity. With this new resource allocation strategy, when the mean required throughput of each mobile is equal to 500 Kbps, the total amount of data transmitted out of the cell in order to extend the coverage area can be increased around 59% compared to the well acknowledged MaxSNR and around 129% compared to the classical scheduling algorithm RR.

V. CONCLUSION

In this paper we have proposed a new incentive approach which encourages nodes to relay neighboors frame. With our proposition, the mobile stays free to cooperate or not but the proposed scheduler CEI sparingly awards participative nodes so that it is more interesting for them to actively contribute to a good network connectivity. This result is a well-balanced resource allocation which allows to increase the network coverage area while never reducing the global system throughput thanks to a combined opportunistic approach. A minimum throughput is guaranteed to all mobiles of the cell and, thanks to its high spectral efficiency, the mean packet delay provided to the selfish mobiles by CEI is close to the best RR performance. These CEI performance results are all the more interesting that a significant priority is given to mobiles which help the network providing a low packet delay and a high personal throughput. In the future work, we plan to introduce services differentiation in our proposed solution.

Fig. 1 .

 1 Fig. 1. Measured mobile mean delay with respect to their cooperation ratio.

Fig. 2 .

 2 Fig. 2. Relay efficiency.

• if mobile k relay no traffic out of the cell, C k equal 0%;• if it is ready to relay 50 packet out of the cell since it receive 100 packets for its own consumption, C k equal 50%; • if the mobile relay as much packet out of the cell that its own received for its own consumption, C k equal 100%.


Didier Henrion
email: henrion@laas.fr

Martine Ganet-Schoeller

Samir Bennani

Measures and LMI for space launcher robust control validation *

Keywords: space launcher control, safety, verification, robustness analysis, measures, moments, polynomials, LMI, convex optimization

We describe a new temporal verification framework for safety and robustness analysis of nonlinear control laws, our target application being a space launcher vehicle. Robustness analysis, formulated as a nonconvex nonlinear optimization problem on admissible trajectories corresponding to piecewise polynomial dynamics, is relaxed into a convex linear programming problem on measures. This infinitedimensional problem is then formulated as a generalized moment problem, which allows for a numerical solution via a hierarchy of linear matrix inequality relaxations solved by semidefinite programming. The approach is illustrated on space launcher vehicle benchmark problems, in the presence of closed-loop nonlinearities (saturations and dead-zones) and axis coupling.

Introduction

This work is carried out within the scope of the project SAFE-V (Space Application Flight control Enhancement of Validation Framework). The objective of this project is to analyse, develop and demonstrate effective design, verification and validation strategies and metrics for advanced guidance navigation and control systems. These new strategies should be implemented in an incremental manner in the traditional validation framework in order to be applicable both for current launchers validation and for future launcher and re-entry vehicle validation.

Control design performance and compliance with system specifications should be demonstrated by a suitable combination of stability analysis and simulation. This validation is made on a wide domain of variation of influent parameters; it includes worst case validation and Monte Carlo simulation for statistic requirements, see e.g. [START_REF] Rongier | Robustness of the Ariane 5 guidance navigation and control[END_REF] for Ariane 5 validation.

Generally speaking, the Monte Carlo method is a numerical integration method using sampling, which can be used, for example, to determine the quantities of interest of a variable of interest for a random input variable such as mean or standard deviation or probability density function. One of the most commonly used methods inside ASTRIUM ST for the confidence interval on a quantile estimation is Wilks formula, see [START_REF] Openturns | Treatment of Uncertainties RisksN Statistics)[END_REF]. This formula of the confidence interval depending on the probability level and the number of samples can be used, for Monte Carlo only, in two ways: to determine the probability and confidence level for the values of samples chosen by the user, or in reverse to determine the number of simulations to be carried out for the values probability and confidence level chosen by the user. Monte Carlo simulations are widely used in traditional validation process and they are very robust because they do not require smoothing assumptions on the model neither importance factor assumptions on the parameters. The random events are simulated as they would occur naturally. The generic Monte Carlo method is however very time consuming, because a great number of simulations to demonstrate extreme level of probability (in general much more than 1/p simulations to estimate the probability of an event of probability p). The number of simulations required may be critical for long duration mission simulation (like launcher flight). In the context of robust control, see [START_REF] Tempo | Randomized algorithms for analysis and control of uncertain systems[END_REF] for related probabilistic approaches. Note that these approaches may miss worst case behavior, especially when the number of parameters is large.

In contrast, worst case validation techniques are aimed at identifying worst case instances and behavior. These techniques include Lyapunov and mu-analysis approaches. Among them, the most relevant for space vehicle applications are:

• results by the GARTEUR group, see Fielding et al. (2002), and the COFCLUO group, see COFCLUO (2007), for recent aeronautic applications and new developments;

• results by ASTRIUM ST in the frame of ESA ITT for Robust LPV for launcher application, see ASTRIUM ST (2009a,b); [START_REF] Ganet-Schoeller | Hybrid Optimization Tools for Worst Case Analysis of Flexible Launcher[END_REF], and for robust stability analysis for ATV industrial validation, see [START_REF] Ganet-Schoeller | Nonlinear and robust stability analysis for ATV rendezvous control[END_REF];

• results by NGC on rendezvous validation process in the frame of the ESA VVAF project, see Di Sotto et al. (2010a,b); [START_REF] Kron | Mu-analysis based verification and validation of autonomous satellite rendezvous systems[END_REF]; [START_REF] Paulino | Worst case and safety analysis tools for autonomous rendezvous system[END_REF], and on re-entry validation, see [START_REF] Kron | The generalized (µ,MU)-iteration illustrated by flexible aircraft robust controller reduction[END_REF].

Many other non industrial applications are available in the literature. Two of them could be considered as major for our applications: a PhD thesis on missile application by [START_REF] Adounkpé | Robustesse dans un cadre non linéaire[END_REF] and the IQC toolbox developed by [START_REF] Scherer | Robustness with dynamic IQCs: an exact state-space characterization of nominal stability with applications to robust estimation[END_REF] with Delft University that was applied for HL-20 re-entry vehicle analysis in [START_REF] Veenman | Analysis of the controlled NASA HL20 atmospheric re-entry vehicle based on dynamic IQCs[END_REF], under ESA contract.

In this paper we propose an original approach to validation of control laws based on measures and convex optimization over linear matrix inequalities (LMIs). The approach can be seen as a blend between worst case validation techniques and statistical simulation techniques. On the one hand, we optimize over worst case admissible trajectories to validate a system property. On the other hand, we propagate along system trajectories statistical information (probability measures) instead of deterministic initial conditions. Moreover, our approach is primal in the sense that we optimize directly over systems trajectories, we do not seek a dual Lyapunov certificate.

First we rephrase our validation problem as a robustness analysis problem, and then as a nonconvex nonlinear optimization problems over admissible trajectories. This is the approach followed e.g. in [START_REF] Prajna | Convex programs for temporal verification of nonlinear dynamical systems[END_REF] where the authors verify or prove temporal properties such as safety (all trajectories starting from a set of initial conditions never reach a set of bad states), avoidance (at least one trajectory starting from a set of initial conditions will never reach a set of bad states), eventuality (all trajectories starting from a set of initial conditions will reach a set of good states in finite time) and reachability (at least one trajectory starting from a set of initial conditions will reach a set of good states in finite time).

Following [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF], we formulate our nonconvex nonlinear trajectory optimization problem as a linear programming (LP) problem on measures, with at most three unknowns: the initial measure (modeling initial conditions), the occupation measure (encoding system trajectories) and the terminal measure (modeling terminal conditions); in some cases a subset of these decision variables may be given. The final time can be given or free, finite or infinite. If all the data are polynomials (performance measure, dynamics, constraints), then we formulate the infinite-dimensional LP on measures as a generalized moment problem (GMP), and we approach it via an asymptotically converging hierarchy of finite-dimensional convex LMI relaxations. These LMI relaxations are then modeled with our specialized software GloptiPoly 3, solved with a general-purpose semidefinite programming (SDP) solver, and system properties are validated from the numerical solutions.

The models we can deal with should be piecewise polynomial (or rational) in time and space. We consider a compact region of the state-space over which the trajectories are optimized. The system dynamics are defined locally on explicitly given basic semialgebraic sets (intersections of polynomial sublevel sets) by polynomial vector fields. The objective function (stability or performance measure) consists of a polynomial terminal term and of the time integral of a polynomial integrand.

Our contribution is to use the approach described in [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] jointly with the corresponding software tool GloptiPoly 3 by [START_REF] Henrion | GloptiPoly 3: moments, optimization and semidefinite programming[END_REF] to address temporal validation/verification of control systems properties, in the applied context of space launcher applications. The use of LMI and measures was already investigated in [START_REF] Prajna | Convex programs for temporal verification of nonlinear dynamical systems[END_REF] for building Lyapunov barrier certificates, and based on a dual to Lya-punov's theorem described in Ranzter (2001). Our approach is similar, in the sense that optimization over systems trajectories is formulated as an LP in the infinite-dimensional space of measures. This LP problem is then approached as a generalized moment problem via a hierarchy of LMI relaxations.

Historically, the idea of reformulating nonconvex nonlinear ordinary differential equations (ODE) into convex LP, and especially linear partial differential equations (PDE) in the space of probability measures, is not new. It was Joseph Liouville in 1838 who first introduced the linear PDE involving the Jacobian of the transformation exerted by the solution of an ODE on its initial condition Liouville (1838), see [START_REF] Ehrendorfer | The Liouville equation and its potential uefulness for the prediction of forecast skill. Part I: theory[END_REF][START_REF] Ehrendorfer | The Liouville equation in atmospheric predictability[END_REF] for a survey on the Liouville PDE with applications in meteorology. The idea was then largely expanded in Henri Poincaré's work on dynamical systems at the end of the 19th century, see in particular (Poincaré, 1899, Chapitre XII (Invariants intégraux)). This work was pursued in the 20th century in [START_REF] Kryloff | La théorie générale de la mesure dans son application à l'étude des systèmes dynamiques de la mécanique non linéaire[END_REF], (Nemystkii and Stepanov, 1947, Chapter VI (Systems with an integral invariant)) and more recently in the context of optimal transport by e.g. [START_REF] Rachev | Mass transportation problems. Volume I: theory[END_REF], [START_REF] Villani | Topics in optimal transportation[END_REF] or [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]. Poincaré himself in [START_REF] Poincaré | L'avenir des mathématiques[END_REF] Section IV) mentions the potential of formulating nonlinear ODEs as linear PDEs, and this programme has been carried out to some extent by [START_REF] Carleman | Application de la théorie des équations intégrales linéaires aux systèmes d'équations différentielles non linéaires[END_REF], see also [START_REF] Lasota | Probabilistic properties of deterministic systems[END_REF], [START_REF] Kowalski | Dynamical systems and Carleman linearization[END_REF], Hernández-Lerma and [START_REF] Hernández-Lerma | Markov chains and invariant probabilities[END_REF] and more recently [START_REF] Barkley | The moment map: nonlinear dynamics of density evolution via a few moments[END_REF], [START_REF] Vaidya | Lyapunov measure for almost everywhere stability[END_REF], [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. For recent studies of the Liouville equation in optimal control see e.g. [START_REF] Kwee | Optimal control using the transport equation: the Liouville machine[END_REF] and [START_REF] Brockett | Optimal control of the Liouville equation[END_REF], and for applications in uncertainty propagation and control law validation see e.g. [START_REF] Mellodge | Uncertainty propagation in abstracted systems. Model abstraction in dynamical systems: applications to mobile robot control[END_REF] and [START_REF] Halder | Dispersion analysis in hypersonic flight during planetary entry using stochastic Liouville equation[END_REF].

Piecewise polynomial dynamic optimization

Consider the following dynamic optimization problem with piecewise polynomial differential constraints

J = inf x(t) h T (T, x(T )) + T 0 h(t, x(t))dt s.t. ẋ(t) = f k (t, x(t)), x(t) ∈ X k , k = 1, 2, . . . , N x(0) ∈ X 0 , x(T ) ∈ X T , t ∈ [0, T ] (1)
with given polynomial dynamics f k ∈ R[t, x] and costs h, h T ∈ R[t, x], and state trajectory x(t) constrained in closed basic semialgebraic sets

X k = {x ∈ R n : g kj (t, x) ≥ 0, j = 1, 2, . . . , N k } for given polynomials g kj ∈ R[t, x].
We assume that the state-space partitioning sets, or cells X k , are disjoint, i.e. all their respective intersections have zero Lebesgue measure in R n , and they all belong to a given compact semialgebraic set, e.g.

X = {x ∈ R n :
x 2 2 ≤ M} for a sufficiently large constant M > 0. Finally, initial and terminal states are constrained in semialgebraic ets

X 0 = {x ∈ R n : g 0j (t, x) ≥ 0, j = 1, 2, . . . , N 0 } ⊂ X and X T = {x ∈ R n : g T j (t, x) ≥ 0, j = 1, 2, . . . , N T } ⊂ X for given polynomials g 0j , g T j ∈ R[t, x].
It is assumed that optimization problem (1) arises from validation of a controlled (closedloop) system behavior. In this problem, the infimum is sought over absolutely continuous state trajectories x(t), but the same methodology can be extended to (possibly discontinuous) trajectories of bounded variation (e.g. in the presence of impulsive controls or state jumps), see [START_REF] Claeys | Measures and LMI for impulsive optimal control with applications to space rendezvous problems[END_REF].

The final time T is either given, or free, in which case it becomes a decision variable, jointly with x(t). Similarly, the initial and terminal constraint sets X 0 and X T are either given, or free, in which case they also become optimized decision variables.

Linear programming on measures

In this section, we formulate optimization problem (1), which is nonlinear and nonconvex in state trajectory x, into a convex optimization problem linear in µ, the occupation measure of trajectory x. This reformulation is classical in the modern theory of calculus of variations and optimal control and it can be traced back to [START_REF] Young | Generalized curves and the existence of an attained absolute minimum in the calculus of variations[END_REF], see also [START_REF] Ghouila-Houri | Sur la généralisation de la notion de commande d'un système guidable[END_REF], [START_REF] Young | Lectures on the calculus of variations and optimal control theory[END_REF] and Gamkrelidze (1975), amongst many others.

Occupation measure

To understand the basic idea behind the transformation and the elementary concept of occupation measure, it is better to deal with the single nonlinear ODE

ẋt = f (x t ) (2)
where x t is a shortcut for x(t), a time-dependent state vector of R n and f : R n → R n is a uniformly Lipschitz map. It follows that the Cauchy problem for ODE (2) has a unique solution x t for any given initial condition x 0 ∈ R n . Now think of initial condition x 0 as a random variable of R n , or more abstractly as a nonnegative probability measure µ 0 with support X 0 ⊂ R n , that is a map from the sigmaalgebra of subsets of X 0 to the interval [0, 1] ⊂ R such that µ 0 (X 0 ) = 1. For example, the expected value of x 0 is the vector

E[x 0 ] = X 0 xµ 0 (dx).
Now solve ODE (2) for a trajectory, or flow x t , given this random initial condition. At each time t, state x t can also be interpreted as a random variable, i.e. a probability measure that we denote by µ t (dx). We say that the measure is transported by the flow of the ODE. The one-dimensional family, or path of measures µ t satisfies a PDE

∂µ t ∂t + div(f µ t ) = 0 (3)
which turns out to be linear in the space of probability measures. This PDE is usually called Liouville's equation. As explained e.g. in (Villani, 2003, Theorem 5.34), nonlinear ODE (1) follows by applying Cauchy's method of characteristics to linear transport PDE (3), see also (Evans, 2010, Section 3.2) for a tutorial exposition. In equation ( 3), div models the divergence operator, i.e.

div(v) = n i=1 ∂v i ∂x i
for every smooth function v. Its action on measures should be understood in the weak, or distributional sense, i.e.

v div(ν) = -Dv • dν
where v is a smooth test function, ν is a vector-valued measure and D stands for the gradient operator. Given a subset T × X in the sigma-algebra of subsets of [0, T ] × X, we define the occupation measure

µ(T × X ) = T µ t (X )dt
which encodes the time-space trajectories x t , in the sense that µ([0, T ] × X ) is the total time spent by trajectory x t in a subset X ⊂ X of the state space. In its integral form,

transport PDE (3) becomes div(f µ) = µ 0 -µ T (4)
where µ T is the terminal probability measure with support X T ⊂ R n . PDE (4) can equivalently be formulated as

X Dv • f dµ = X T vdµ T - X 0 vdµ 0 (5)
for all smooth functions v compactly supported on X. Problem ( 5) is an infinite-dimensional linear system of equations linking occupation measure µ, initial measure µ 0 and terminal measure µ T , consistently with ODE (2).

Finite terminal time

If we apply these ideas to problem (1), encoding the state trajectory x(t) in an occupation measure µ supported on X, we come up with an infinite-dimensional LP problem

J ∞ = inf µ h T (T, x)dµ T (x) + k h(t, x)dµ k (t, x) s.t. k ∂v(t,x) ∂t dµ k (t, x) + k Dv(t, x) • f k (t, x) dµ k (t, x) = v(T, x) dµ T (x) -v(0, x) dµ 0 (x) (6)
for all smooth test functions v, where each occupation measure µ k is supported on set X k and the global occupation measure is

µ = k µ k with normalization constraint µ([0, T ] × X) = T (7)
such that T is a finite terminal time.

In problem (6), final time T , initial measure µ 0 and terminal measure µ T may be given, or unknown, depending on the original optimization problem.

Generalized moment problem

More concisely, LP problem (6) can be formulated as follows:

inf µ k X k c k dµ k s.t. k X k a ki dµ k = b i , ∀i (8)
where the unknowns are a finite set of nonnegative measures µ k , with respective compact semialgebraic supports

X k = {x : g kj (x) ≥ 0, ∀j}. (9)
Note that here we have incorporated time variable t into vector x, for notational conciseness. If all the data in problem ( 6) are polynomials, and if we generate test functions v(x) using a polynomial basis (e.g. monomials, which are dense in the set of continuous functions with compact support), all the coefficients a(x), b(x), c(x) are polynomials, and there is an infinite but countable number of linear constraints indexed by i.

We will then manipulate each measure µ k via its moments

y kα = X k x α dµ k (x), ∀α (10)
gathered into an infinite-dimensional sequence y k indexed by a vector of integers α, where we use the multi-index notation x α = x α 1 1 x α 2 2 . . . LP measure problem (8) becomes an LP moment problem, or GMP, see [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]:

inf y k α c kα y kα s.t. k α a kiα y kα = b i , ∀i
provided we can handle the representation condition (10) which links a measure with its moments. It turns out [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF], Chapter 3) that if sets X k are compact semialgebraic as in (9), we can use results of functional analysis and real algebraic geometry to design a hierarchy of LMI relaxations which is asymptotically equivalent to the generalized moment problem. In each LMI relaxation

J d = inf y c T y s.t. Ay = b M d (y) 0 M d (g kj , y) 0, ∀j, k (11)
we truncate the moment sequence to a finite number of its moments. If the highest moment power is 2d, we call LMI relaxation of order d the resulting finite-dimensional truncation. Matrix M d (y) is symmetric and linear in y, it is called the moment matrix, and it must be positive semidefinite for (10) to hold. Symmetric matrices M d (g kj , y) are also linear in y and they are called localizing matrices. They ensure that the moments correspond to measures with appropriate supports. For d finite, problem ( 11) is a finitedimensional convex LP problem in the cone of positive semidefinite matrices, or SDP problem, for which off-the-shelf solvers are available, in particular implementations of the primal-dual interior-point methods described in [START_REF] Yu | Interior-point polynomial algorithms in convex programming[END_REF].

Convergence

The hierarchy of LMI relaxations (11) generates an asymptotically converging monotonically increasing sequence of lower bounds on LP (6), i.e. J d ≤ J d+1 for all d = 1, 2 . . . and lim d→∞ J d = J ∞ . Obviously, the relaxed LP cost in problem ( 6) is a lower bound on the original cost in problem (1), i.e. J ∞ ≤ J. Under mild assumptions, we can show that actually J ∞ = J, but the theoretical background required to prove this lies beyond the scope of this application paper.

Long range behavior

In the case terminal time T tends to infinity, normalization constraint (7) becomes irrelevant, and the overall mass of occupation measure µ tends to infinity. A more appropriate formulation consists then of normalizing measure µ to a probability measure

π = µ T = k µ k T = k π k so that PDE (4) becomes div(f π) = lim T →∞ µ 0 -µ T T = 0 and LP problem (6) becomes inf π k h(t, x)dπ k (x) s.t. k ∂v(t,x) ∂t dπ k (x) + k Dv(t, x) • f k (t, x) dπ k (t, x) = 0
where test functions v(t, x) satisfy appropriate integrability and/or periodicity properties. Measure π is called an invariant probability measure and it encodes equilibrium points, periodic orbits, ergodic behavior, possibly chaotic attractors, see e.g. [START_REF] Lasota | Probabilistic properties of deterministic systems[END_REF], [START_REF] Diaconis | Iterated random functions[END_REF], Hernández-Lerma and [START_REF] Hernández-Lerma | Markov chains and invariant probabilities[END_REF] and [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF].

Application to orbital launcher control law validation

In this section we report the application of our method to simplified models of an attitude control system (ACS) for a launcher in exo-atmospheric phase. We consider first a one degree-of-freedom (1DOF) model, and then a 3DOF model. The full benchmark is described in (SAFEV, 2011, Section 9) but it is not publicly available, and our computer codes cannot be distributed either.

ACS 1DOF

First we consider the simplified 1DOF model of launcher ACS in orbital phase. The closedloop system must follow a given piecewise linear angular velocity profile. The original benchmark includes a time-delay and a pulsation width modulator in the control loop, but in our incremental approach we just discard these elements in a first approximation. See below for a description of possible extensions of our approach to time-delay systems.

As to the pulsation width modulator, it can be handled if appropriately modeled as a parametric uncertainty.

Model

The system is modeled as a double integrator

I θ(t) = u(t)
where I is a given constant inertia and u(t) is the torque control. We denote

x(t) = θ(t) θ(t)
and we assume that both angle x 1 (t) and angular velocity x 2 (t) are measured, and that the torque control is given by

u(x(t)) = sat(K T dz(x r (t) -x(t)))
where x r (t) is the reference signal, K ∈ R 2 is a given state feedback, sat is a saturation function such that sat(y) = y if |y| ≤ L and sat(y) = L sign(y) otherwise, dz is a deadzone function such that dz(x) = 0 if |x i | ≤ D i for some i = 1, 2 and dz(x) = 1 otherwise. Thresholds L > 0, D 1 > 0 and D 2 > 0 are given.

We would like to verify whether the system state x(t) reaches a given subset X T = {x ∈ R 2 : x T x ≤ ε} of the deadzone region after a fixed time T , and for all possible initial conditions x(0) chosen in a given subset X 0 of the state-space, and for zero reference signals. We formulate an optimization problem (1) with systems dynamics defined as locally affine functions in three cells X k , k = 1, 2, 3 corresponding respectively to the linear regime of the torque saturation

X 1 = {x ∈ R 2 : |K T x| ≤ L}, f 1 (x) = x 1 -K T x the upper saturation regime X 2 = {x ∈ R 2 : K T x ≥ L}, f 2 (x) = x 1 -L
and the lower saturation regime

X 3 = {x ∈ R 2 : K T x ≤ -L}, f 3 (x) = x 1 L .
The objective function has no integral term and a concave quadratic terminal term h T (x) = -x(T ) T x(T ) which we would like to minimize, so as to find trajectories with terminal states of largest norm. If we can certify that for every initial state x(0) chosen in X 0 the final state x(T ) belongs to set included in the deadzone region, we have validated our controlled system.

Validation script

The resulting GloptiPoly 3 script, implementing some elementary scaling strategies to improve numerical behavior of the SDP solver, is as follows: If we want to use this approach to simulate a particular trajectory, in the code we must modify the definition of the initial measure. For example for initial conditions x 1 (0) = 50, x 2 (0) = -1, we must insert the following sequence:

% given moments of initial measure = Dirac at x0 p = genpow(3,d); p = p(:,2:end); % powers theta0 = 50; omega0 = -1; % in degrees y0 = ones(size(p,1),1)*[theta0 omega0]*pi/180; y0 = prod(y0.^p,2);

As previously, the sequence of bounds on the maximum squared Euclidean norm of the final state is constantly equal to 1.0 32 5.3 0.74 0.30 0.21 0.15 0.17 dµ 3 32 5.1 7.1 6.9 6.8 6.9 7.0

This indicates that most of the time (approx. 93%) is spent in the linear regime, with approx. 7% of the time spent in the lower saturation regime, and a negligible amount of time is spent in the upper saturation regime. This is confirmed by simulation, see Figure 1.

Dealing with uncertainty

Finally, note that we can incorporate real parametric uncertainty in the dynamics, if required. Each uncertain parameter must be introduced as an additional state of the system, and this generally makes the dynamics polynomial in the extended state. The parameters must be constrained to an explicitly given compact semialgebraic set, and the overall trajectory optimization problem consists of finding the worst-case uncertain parameter instance.

For illustration, we modify the previous script to cope with uncertainty entering affinely the dynamics. We assume that the (reciprocal of the) inertia is subject to multiplicative relative uncertainty, i.e. we replace occurences of I with 1 1+u I for u a real parameter such that u 2 ≤ U with U > 0 a given threshold, say U = 1 2 . The resulting GloptiPoly 3 script is as follows: Introducing an uncertain parameter amounts to adding a state variable to the model, and this has a significant impact on the overall computational time, as shown in the table below: Note however that for this example there is no effect on the bounds, and the control law is also validated in the presence of uncertainty on inertia.

relaxation order d 1 2 3 4 upper bound J d 1.0 • 10 -5 1.0 • 10 -5 1.0 • 10 -5 1.

ACS 3DOF

The main difficulty with the 3 degree-of-freedom version of the ACS benchmark is the large number of states (4 quaternions for the launcher attitude, 3 angular velocities, and 2 additional states for the quaternion reference signal) and the non-linear (quadratic) dynamics.

Model

The following example illustrates the validation of a given control law following a given constant roll velocity equal to w R 1 = 20 o /s, in the absence of actuator saturation, but in the presence of non-linear coupling between the three axes of the launcher. System dynamics are given by ẋ = f (x) with

x =   q w q R   , f (x) =   f q (x) f w (x) f q R (x)   , and
f q (x) = 1 2 Q(x) 0 w , f w (x) = I -1 (u(x) -Ω(x)Iw), f q R (x) = -1 2 q R 1 w R 1 1 2 q R 0 w R 1 where Q(x) =     q 0 -q 1 -q 2 -q 3 q 1 q 0 -q 3 q 2 q 2 q 3 q 0 -q 1 q 3 -q 2 q 1 q 0     , Ω(x) =   0 -w 3 w 2 w 3 0 -w 1 -w 2 w 1 0   , u(x) = -K D w E -K P q E , K D =   K D 1 0 0 0 K D 2 0 0 0 K D 3   , K P =   0 0 0 0 K P 2 0 0 0 K P 3   , w E =   w 1 -w R 1 w 2 w 3   , q E = 2     q R 0 q R 1 0 0 -q R 1 q R 0 0 0 0 0 q R 0 q R 1 0 0 -q R 1 q R 0     q,
and q ∈ R 4 is the quaternion modeling the attitude, w ∈ R 3 models the angular velocities, q R ∈ R 2 is the reference quaternion to follow, which depends on the constant reference velocity w R 1 ∈ R, q E ∈ R 4 is the quaternion error, and w E ∈ R 3 is the angular velocity error. Both vectors q(t) and q R (t) have constant Euclidean norm since d dt q 2 2 = 2q T q = q T f q (x) = 0 for all x and similarly d dt q R 2 2 = 0, so we enforce the algebraic constraints

q 2 2 = 1, q R 2 2 = 1.
The proportional-derivative control law u(x) is designed so that velocity w 1 (t) follows reference velocity w R 1 , and our validation tasks consists of optimizing over the worst-case trajectory starting at a given initial condition x(0) and maximizing the objective function

J = T 0 (w 1 (t) -w R 1 ) 2 dt
for a given time horizon T . If we can guarantee a sufficiently small upper bound J d on this objective function, for a given relaxation order d = 1, 2, . . ., the control law is validated.

Validation script

We use the following GloptiPoly 3 script (note however that the function G2I to generate the inertia data is not available for the reader): -50, -1100; -50, 45300, -220; -1100, -220, 44100]; [Ge2In,Ip] = G2I(Ig); % controller gains

% Inertia Ig = [ 27500,

Coping with saturations and dead-zones

We can modify the above code to cope with actuator saturation and dead-zone, proceeding exactly as described above for the ACS 1DOF benchmark problem, namely by splitting the trajectory occupation measures into local occupation measures defined on semialgebraic (here polyhedral) cells. In the presence of saturations on the three actuators, this generates 3 3 = 27 local occupation measures. This should not be an issue from the computational point of view, the limiting factor being essentially the size of the largest SDP block, which grows polynomially with the relaxation order, but with an exponent which is the number of variables (here equal to 10).

Following a time-varying reference signal

Finally, if we want to validate a control law to follow a time-varying velocity profile (instead of a constant velocity as above), we must introduce time as an additional variable in the trajectory occupation measure, we must introduce the velocity as an additional state, and we must split the occupation measure into local occupation measures, each of which corresponding to given time-invariant dynamics. Alternatively, we can also model time-varying dynamics, but the dependence on time must be polynomial (which is not the case e.g. if the reference signal is piecewise linear).

Discussion

In this section we discuss the limitations of our approach and we also describe possible extensions.

Limitations

The main limiting factor for our method is the total number of variables (number of states plus number of uncertain parameters), since the computational burden of solving the LMI relaxations grows polynomially as a function of the LMI relaxation order, but the order of the polynomial dependence depends linearly on the number of variables in the measures (here the number of states). So the critical dimension is not too much the number of measures (that is the number of cells partioning the state-space, which corresponds to our models of nonlinearities), or the degree of the polynomials in the cost and/or dynamics, but rather the number of states. Systems with a high number of states can be handled with these techniques, but only when exploiting problem structure and sparsity.

More explicitly, a rough complexity analysis can be carried out as follows. If an LMI relaxation has the simplified form inf y c T y s.t. M d (y) 0 where M d (y) is the moment matrix of a measure of n variables at relaxation order d, then the number of variables in vector y is

N = n + 2d n = (n + 2d)! n! (2d)!
and the size of matrix M d (y) is

M = n + d n = (n + d)! n! d! .
A standard primal-dual interior-point algorithm to solve this LMI at given relative accuracy ǫ > 0 (duality gap threshold) requires a number of iterations (Newton steps) growing as O(M 1 2 log ǫ), whose dependence on M is sublinear, hence almost negligible. In practice, on well-conditioned problems, we observe that the number of Newton iterations is between 5 and 50, see [START_REF] Yu | Interior-point polynomial algorithms in convex programming[END_REF], [START_REF] Vandenberghe | Semidefinite programming[END_REF] and [START_REF] Ben-Tal | Lectures on modern convex optimization[END_REF].

In the real model of computation (for which each addition, subtraction, multiplication, division of real numbers has unit cost), each Newton iteration requires O(N 2 M)+O(N 3 M)+ O(N 2 M 2 ) operations to form the linear system of equations, and O(M 3 ) operations to solve the system and find the search direction. When solving a hierarchy of simple LMI relaxations as described above, the number of variables n is fixed, and the relaxation order d varies, so the dominating term in the complexity estimate grows in O(d 4n ), which clearly shows a strong dependence on the number of variables. Even though the growth of the computational burden is polynomial in the relaxation order, the exponent is 4 times the number of variables. One should however keep in mind that these estimates are (usually very loose) asymptotic upper bounds, and that the observed computational complexity grows much more moderately in practice (at least in the absence of conditioning and numerical stability issues).

Finally, if the LMI relaxation has the form

inf y c T y s.t. M d (y k ) 0, k = 1, 2, . . . , K
where M d (y k ) is the moment matrix of a measure µ k of n variables at relaxation d, and we have K measures, then the above complexity estimate grows in O(Kd 4n ). The impact on the computation burden of the number n of variables in each measure µ k is thus much more critical than the number K of measures. In our target application, n is the number of states and uncertain parameters, K is the number of cells used to model nonlinearities, and a lower bound on the minimum relaxation order d is given by the degree of the polynomial data (dynamics, constraints, objective function).

Accuracy

The original validation problem can be formulated as infinite-dimensional linear programming problems on the dual cones of nonnegative measures and continuous functions, so that there is no hope of finite convergence of finite-dimensional LMI optimization techniques. Convergence is guaranteed only asymptotically. Obviously, the accuracy depends on the speed of convergence of the hierarchy of LMI problems, but this speed is impossible to evaluate a priori. The only guarantee is that we have a monotically increasing sequence of lower bounds on the objective function to be minimized.

As far as solving finite-dimensional LMI problems is concerned, we should emphasize the fact that there is currently no proof of backward stability of implementations of LMI solvers. Moreover, evaluation of the conditioning of a given LMI problem is difficult, and estimates can be obtained only at the price of solving several instances of (slightly modified versions) the original LMI problem. To certify the output of a numerical algorithm, we must simultaneously ensure backward stability of the algorithm and well-conditioning of the problem. None of these two properties can be ensured when solving LMI problems, given the current state-of-the-art in LMI solvers. But this is not specific to our method, and any validation-verification technique which is based on LMI optimization is subject to the same limitations.

Time-delay systems

Our techniques can in principle cope with time-delays, but this requires a significant modification of the approach described in this document. Let us only sketch the main ideas, in the case of an ordinary differential equation with one time-delay

ẋ(t) = f (x(t)) + g(x(t -τ )), ∀t ∈ [0, T ] with boundary conditions x(t) = ξ(t), ∀t ∈ [-τ, 0]
where ξ(t) is a given function recording the state history due to the given delay τ ∈ R.

Instead of transporting a probability measure µ t (dx) supported on X ⊂ R n from initial time t = 0 to terminal time t = T , we must transport the state history in an occupation measure µ t (ds, dx) supported on [-τ, 0] × X for t ∈ [0, T ].

Discrete-time systems

In this paper we deal only with continuous-time systems. Discrete-time systems of the form

x k+1 = f (x k ) (12)
must be handled differently. Denoting by µ k (x) the probability measure transported along dynamical system (12), the discrete-time analogue of Liouville's transport equation

(3) reads µ k+1 (X ) = f -1 (X ) µ k (dx) = I X (f (x))µ k (dx) ( 13
)
where X is any subset of the sigma-algebra of state set X, and I X is the indicator function equal to one in X and zero outside. It follows from ( 13) that the moments of measure µ k+1 can be expressed linearly as functions of moments of measure µ k . Besides this analogy, the resulting discrete-time generalized moment problem differs significantly from its continuous-time counterpart. For this reason, handling discrete-time systems would require an important modification of the approach described in this document.

Conclusion

In this paper we do not follow the mainstream Lyapunov approach to dynamical systems stability and/or performance validation. Lyapunov (1892) originally introduced his approach to conclude about systems behaviour while avoiding explicit computations of system trajectories. In the 1990s it was combined with LMI techniques, see e.g. [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF][START_REF] Yu | Interior-point polynomial algorithms in convex programming[END_REF] to provide numerically quadratic certificates of stability and/or performance for linear and nonlinear systems. Later on, it was extended to more general, but still numerical polynomial sum-of-squares (SOS) certificates, see e.g. [START_REF] Henrion | Positive polynomials in control[END_REF]. By contrast, our approach is genuinely primal, in the sense that we directly optimize numerically over systems trajectories, using measures, moments and LMI techniques. Indeed, since the computation of SOS certificates is numerical anyway, we believe that in the context of validation, LMI techniques should rather be used to optimize systems trajectories, instead of optimizing indirect certificates for these trajectories. Since we are using primal-dual SDP solvers to solve the hierarchy of LMI problems, solving the dual to the problem of optimizing trajectories provides anyway certificates of feasibility or infeasibility, i.e. of stability, instability and/or performance.

As briefly described in §5, the measure/LMI approach can be extended easily to systems with real uncertainties, as soon as the uncertainties enter polynomially in the dynamics, and they are bounded in explicitly given semialgebraic sets (e.g. balls or boxes). See e.g. (Lasserre, 2009, Section 13.1) for more information on robust optimization. As usual, the number of real uncertain parameters should be kept small enough to ensure computational tractability.

Note that in this paper the set of initial conditions is assumed to be given, and when the initial measure is unknown, its support is constrained to be included in the set of initial conditions. In the context of validation/verification it could be relevant to maximize the size of the set of initial conditions, and we are currently investigating this problem from the point of view of occupation measures.

The measure/LMI approach can readily deal with systems with piecewise polynomial (or rational) models, e.g. systems with input/output saturations and dead-zones. It can be seen as a (primal) extension of early attempts of use of LMI/Lyapunov techniques in the context of systems with saturations, see e.g. [START_REF] Henrion | LMI relaxations for robust stability of linear systems with saturating controls[END_REF] and more recently [START_REF] Garulli | Piecewise polynomial Lyapunov functions for global asymptotic stability of saturated uncertain systems[END_REF][START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF]. As briefly discussed in §5, our approach can be extended without major difficulty to time-delay, stochastic, discretetime systems or hybrid system dynamics where transitions between models are ruled e.g. by probability measures, see e.g. [START_REF] Diaconis | Iterated random functions[END_REF], Hernández-Lerma and [START_REF] Hernández-Lerma | Markov chains and invariant probabilities[END_REF] or [START_REF] Barkley | The moment map: nonlinear dynamics of density evolution via a few moments[END_REF].

  pi/180; dz2 = 0.05*pi/180; % deadzone levels thetamax = 50; omegamax = 5; % bounds on initial conditions epsilon = sqrt(1e-5); % bound on norm of terminal condition T = 50; % final time d = input('order of relaxation ='); d = 2*d; % measures mpol('x1',2); m1 = meas(x1); % linear regime mpol('x2',2); m2 = meas(x2); % upper sat mpol('x3',2); m3 = meas(x3); % lower sat mpol('x0',2); m0 = meas(x0); % initial mpol('xT',2); mT = meas(xT); % terminal % dynamics on normalized time range [0,1] % saturation input y normalized in [-1,1] K = -[kp kd]/L; y1 = K*x1; f1 = T*[x1(2); L*y1/I]; % linear regime y2 = K*x2; f2 = T*[x2(2); L/I]; % upper sat y3 = K*x3; f3 = T*[x3(2); -L/I]; % lower set % test functions for each measure = monomials g1 = mmon(x1,d); g2 = mmon(x2,d); g3 = mmon(x3,d); g0 = mmon(x0,d); gT = mmon(xT,d);

IFigure 1 :

 1 Figure 1: Torque input with lower saturation during approx. 7% of the time range.

  • 10 -5 , and in the following table we represent as functions of the relaxation order d the masses of measures µ k , k = 1, 2, 3 which are indicators of the time spent by the trajectory in the respective linear, upper saturation and lower saturation regimes:

	relaxation order d 1	2	3	4	5	6	7
	dµ 1	37 89	92	92	93	93	93
	dµ 2

Acknowledgment

We are grateful to Dimitri Peaucelle and Prathyush Menon for their comments on this work.

* This work was performed in the frame of the SAFE-V (Space Application Flight control Enhancement of Validation Framework) ESA (European Science Agency) TRP (Technology Development Project) under contract number 4000102288. The research of the first author was also supported by research project 103/10/0628 of the Grant Agency of the Czech Republic.

Numerical results

With this Matlab code and the SDP solver SeDuMi 1.3 (see sedumi.ie.lehigh.edu) we obtain the following sequence of upper bounds (since we maximize) on the maximum squared Euclidean norm of the final state:

0.2 0.5 0.7 0.9 number of moments 30 75 140 225

In the table we also indicate the CPU time (in seconds, on a standard desktop computer) and the total number of moments (size of vector y in the LMI relaxation (11). We see that the bound obtained at the first relaxation (d = 1) is not modified for higher relaxations. This clearly indicates that all initial conditions are captured in the deadzone region at time T , which is the box [-2, 2] 10 -1 π 180 × [-5, 5] 10 -2 π 180 ⊃ {x ∈ R 2 : x T x ≤ 10 -5 }. % occupation measure mpol('q',4); mpol('w',3); mpol('qr',2); x = [q;w;qr]; m = meas(x); % initial measure mpol('q0',4); mpol('w0',3); mpol('qr0',2); x0 = [q0;w0;qr0]; m0 = meas(x0); % terminal measure mpol('qt',4); mpol('wt',3); mpol('qrt',2); xt = [qt;wt;qrt]; mt = meas(xt); % dynamics on normalized time range [0,1] Omega = [0 -w(3) w(2); w(3) 0 -w(1); -w(2) w(1) 0]; Qr = [qr(1) qr(2) 0 0; -qr(2) qr(1) 0 0; 0 0 qr(1) qr(2); 0 0 -qr(2) qr( 1)]; dq = 2*Qr*q; Q = [q(1) -q(2) -q(3) -q(4); q(2) q(1) -q(4) q(3); q(3) q(4) q(1) -q(2); q(4) -q(3) q(2) q( 1

Numerical results

We obtain the monotonically decreasing upper bounds J d , d = 1, . . The first LMI relaxation (110 moments of degree up to 2) seems to be unbounded above (the dual LMI problem is infeasible) and hence it conveys no useful information. The second LMI relaxation (770 moments of degree up to 3) and the third LMI relaxation (1430 moments of degree up to 4) are solved with SeDuMi 1.3 in a few seconds. The fourth LMI relaxation is more challenging, with 5720 moments of degree up to 5, and it requires less than one hour of CPU time on a standard PC. We observe however that a useful upper bound on J is already obtained at a low relaxation order (say d = 2 or d = 3) at a low computational cost, and that the computational burden increases significantly for d = 4 without dramatic improvements in the quality of the upper bound.