Keywords: Absolute localization, Intervals analysis, Digital Elevation Map, autonomous vehicle

Absolute position estimation is a key function of autonomous systems and is often required at the mission level. One way the get such information is to use an initial model of the environment provided for instance, by a geographic database or built by an other robot. This study aims at providing an algorithm which matches the output of a sensor with an initial model to estimate the pose of the robot in a guaranteed way, using the intervals analysis framework. By construction, initial models do not perfectly fit the reality and the acquire data set can contains an unknown, and potentially large, percentage of outliers. When the environment is described by a surface, the set membership estimator GOMNE can be used to concurrently estimate the number of outliers and the localisation parameters. However, in the general case, with full 3D representation and partially mapped objects, it can not be used. To cope this issue, a new algorithm called Outer-GOMNE is proposed. By combining intervals methods with local estimation algorithms, it have been applied to the localization problem. After a simulated test case with a 2D map, an experimental validation using real laser data and different 3D models is reported to illustrate the performance of the method. Results are compared with ground truth provided by a differential GPS. Outer-GOMNE is able to robustly enclose the ground truth in a sub-paving (union of non overlapping boxes)

Résumé

La localisation d'un robot mobile dans son environnement est une fonction essentielle des systèmes autonomes. C'est l'une des briques de base permettant d'effectuer des actions de plus haut niveau, tel que la cartographie ou la planification de trajectoire. L'objectif de ce stage était d'élaborer une méthode de localisation absolue, c'est-à-dire sans connaissance a priori de la position du robot, à partir d'un modèle numérique de terrain initial et des données extéroceptives fournies par les capteurs du robot. Cependant, par construction, ces modèles contiennent de nombreuses erreurs qui amènent à estimer, conjointement à la position du robot, le nombre de mesures aberrantes présentes dans les données. Dès lors que le modèle de terrain n'est pas surfacique, les techniques d'estimations du nombre de mesures aberrantes ne fonctionnent plus. Outer-GOMNE a donc été développé pour palier à ce manque. En combinant les méthodes ensemblistes avec des algorithmes d'estimation locaux (optimisation de fonction, Monte-Carlo, ...) il s'est avéré efficace et robuste pour ce type de problème. Une expérience avec des données réelles, et différents modèles de terrain, a été conduite pour évaluer et valider l'algorithme.

Introduction

Metric maps of buildings, urban and natural environments are becoming widely available, and make map-based localization a key function that perfectly complements dead-reckoning and SLAM approaches to solve the overall robot localisation problem. By providing an "absolute" position estimate (actually relative to the map frame), map-based localization can indeed be exploited either as a solution to the kidnapped robot problem, to palliate the inherent drift of dead-reckoning, in SLAM approaches, or in multi-robot contexts where one robot localizes itself within a map built by an other.

Related work

The literature provides various approaches to this problem, which can be characterized by the way the acquired data are matched to the a priori map (data association), and the algorithms that compute the position from the matches (estimation). Skyline localisation approaches match the skyline detected on images with a predicted skyline derived from a 3D model in urban environments [START_REF] Ramalingam | SKYLINE2GPS: Localization in Urban Canyons using Omni-Skylines[END_REF], or from a digital elevation map (DEM) in planetary environments [START_REF] Stein | Map-based localization using the panoramic horizon[END_REF][START_REF] Cozman | Outdoor visual position estimation for planetary rovers[END_REF]. In this latter context, the absence of absolute localization means and the availability of sub-meter accuracy orthoimages and DEMs has motivated the development of various approaches, that register orthoimages built from camera and LIDAR images [START_REF] Sheshadri | Position estimation by registration to planetary terrain[END_REF], correlate locally built DEMs with the global one [START_REF] Hwangbo | Integration of orbital and ground image networks for the automation of rover localization[END_REF], or matches features (peaks) extracted from LIDAR scans and the global DEM [START_REF] Patrick | Long-Range Rover Localization by Matching LIDAR Scans to Orbital Elevation Maps[END_REF]. An other relevant approach is [START_REF] Vandapel | Unmanned ground vehicle navigation using aerial ladar data[END_REF], where spin-images are associated between aerial and ground data. Note that the "pinpoint landing" problem, which consists in localizing a lander within existing orbiter DEMs and orthoimages, also call for similar map-based techniques [START_REF] Trawny | Vision-aided inertial navigation for pin-point landing using observations for mapped landmarks[END_REF].

In urban environments, the building of maps dedicated to vehicles localization can yield excellent precision: for instance [START_REF] Levinson | Map-Based Precision Vehicle Localization in Urban Environments[END_REF] propose algorithms to build an orthoimage of ground reflectivity, which is later used to augment the precision of GPS/IMU/odometry solution. Particle filters are often used as the estimation engine for map-based localisation [START_REF] Fox | Particle filters for mobile robot localization[END_REF]. One of their advantage is that they do not rely on explicit data association: only a likelihood measure is required between the acquired data and the initial map for a given position hypothesis, and the estimation scheme progressively discards wrong hypotheses. Yes this approach need to be complemented with additional motion estimation techniques for the prediction step, and depending on the environment, can require lengthy displacements to provide a precise localisation estimate when starting without any prior information on this estimate.

Motivations and approach

In most of the outdoor and indoor environments, map-based localisation approaches are challenged by the presence of dynamic elements, which generate numerous outliers, and hence considerably hinder the data association process or the likelihood measure of a position hypothesis. The environment dynamics have a variety of frequencies, from fast moving people or vehicles, to slowly seasonal changes or infrastructure evolution, via changes in the illumination conditions. Our aim is to define an approach to estimate the robot position that is robust with respect to these changes. For that purpose, we rely on the following two choices:

• Exploit geometric information: the geometry of the environment is an intrinsic property that does not depend on illumination conditions, neither on the sensors that capture it and their viewpoints. Furthermore, most of initial maps are geometric, and in particular Digital Elevation Maps are easily built with either aerial or ground means;

• Exploit interval analysis, that cast the localisation problem into a set inversion problem.

The literature has shown that interval analysis defines robust estimation solutions, while ensuring integrity: when a estimation is provided, the actual position is guaranteed to lie within the associated bounds.

Yet, the potential large number of outliers, the facts that initial maps are always partial and that the initial position can very poorly known challenge the existing set membership approaches. The main contribution of this work is the proposal of an interval analysis algorithm that can cope with these problem characteristics.

Outline The next section briefly recall the basics of interval artihmetics and set inversion, so as to familiarize the reader with this formalism. Section 3 depicts the GOMNE algorithm, an essential contribution of the interval analysis literature to achieve robust set inversion in the presence of an unknown proportion of outliers. Our extension to this algorithm is presented in section 4, with an extensive analysis of its characteristics in simulation. An extension to the estimate of the whole trajectory is presented in section 5, and finally section 6 presents its application to the robot localization problem using range data and initial Digital Elevation Maps. A discussion concludes the report.

Background On Interval Analysis For Estimation

This section summarizes the basics of interval analysis, for a more detailed introduction see [START_REF] Jaulin | Applied Interval Analysis[END_REF].

Interval Arithmetics

An interval, denoted by [x] is a closed and connected subset of R defined by a lower and upper bound such as:

[x] = x -, x + = x ∈ R, x -≤ x ≤ x + (1)
The width w of an interval is defined by w

([x]) = x + -x -. A box [x] of R n is a Cartesian product of n intervals.
The set of all boxes of R n is denoted IR n . Basic operations on real numbers or vectors can be extended to intervals in a natural way. For two intervals [x] and [y] of R n and an operator ⋄ ∈ {∩, ∪, +, -, * , /}, we define [x] ⋄ [y] as the smallest interval containing all feasible values for x ⋄ y when x ∈ [x] and y ∈ [y]:

[x] ⋄ [y] = [{x ⋄ y|x ∈ [x], y ∈ [y]}] (2)
where [.] denotes the convex hull. For instance:

[-1, 5] ∩ [-3, 2] = [-1, 2] [-3, 3] ∪ [4, 10] = [-3, 10] [-3, 3] -[4, 10] = [-13, 1] (3)
Let f : R n → R m . An inclusion function denoted [f ] of f is defined by:

∀ [x] ∈ IR n , [f ] ([x]) [{f (x) , x ∈ [x]}] (4) An inclusion function [f ] is convergent if, for any sequence of boxes [x](k) lim k→∞ w ([x](k)) = 0 ⇒ lim k→∞ w ([f ]([x](k))) (5)
For any function obtained by the composition of elementary operators such as +, -, *, /, cos, sin, exp, . . . , a natural inclusion function can be built by replacing these operators by their interval counterpart. For instance, let

f (x 1 , x 2 ) = x 1 * cos (x 2 ). The natural inclusion function [f ] is given by: [f ] ([x 1 ], [x 2 ]) = [x 1 ] * cos ([x 2 ]) (6)
Its evaluation is done by applying interval arithmetic:

[f ] [-2, 1], [ π 6 , π 3 ] = [-2, 1] * cos [ π 6 , π 3 ] = [-2, 1] * [ 1 2 , √ 3
2 ] = [- √ 3, √ 3 2 ]
(7)

Constraint Propagation

A contractor C is an operator from IR n to IR n such as:

C ([x]) ⊂ [x] (contractance) [x] ⊆ [y] ⇒ C ([x]) ⊆ C ([y]) (monotonicity) (8)
A set X is consistent with the contractor C (denoted as X ∼ C) if for all [x], we have

C ([x]) ∩ X = [x] ∩ X (9)
Let X be a set defined with a constraint C such as an equality or an inequality. The contractor C out associated with C makes it possible to remove part of a arbitrary initial box [x] which are not consistent with C without removing any feasible values. This contractor is called the outer contractor. If X has a non empty volume ( there exist points in X which do not belong to the border of X), based on the De Morgan rules, it is also possible to define the inner contraction which removes parts of [x] which belong to the X set. For instance, a curve such as a line has an empty volume because it is only described by its border (X = δX) and no inner contractor can be defined. On the contrary, it can be defined for a disk or a ring. These two example are depicted on figure 1. When available, these two contractors have complementary actions and from the pair {C in , C out }, a separator [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF] can be built. In this case the separator S associated with the set X is an application such as:

S : IR n -→ IR n × IR n [x] -→ {[x in ], [x out ]} = {C in ([x]), C out ([x])}
We define the remainder δS of a separator S as

δS = C in ([x]) ∩ C out ([x]
). δS is a contractor on the border of the set (see figure 2).

   X = {(x, y) ∈ R 2 , x 2 + y 2 ≤ 1} X = {(x, y) ∈ R 2 , x 2 + y 2 > 1} δX = {(x, y) ∈ R 2 , x 2 + y 2 = 1} X X δX X δX = X X = {(x, y) ∈ R 2 , f (x, y) = 0} X = {(x, y) ∈ R 2 , f (x, y) = 0}
Figure 1: The disk, described by an equality admits an inner contractor. On the contrary, only an outer contractor can be build for the curve.

X

[x]

C out ([x]) [x] C in ([x]) δS ([x])
Figure 2: Illustration of a separator on three different initial intervals. The outer contractor removed the blue dashed area and the red dashed area was removed by the inner contractor. Their combined action leads to the yellow box.

Set Inversion

Let f be a function from R n to R m (possibly non-linear) and let Y be a subset of R m such as for instance a sub-paving (finite union of non-overlapping boxes). The Set Inversion is the characterization of :

X = {x ∈ R n |f (x) ∈ Y} = f -1 (Y) (10)
For any Y ⊂ R m and for any function f admitting a convergent inclusion function [f ] (.), two regular sub-pavings X -and X + can be obtained with an arbitrary precision using the algorithm SIVIA (Set Inversion Via Interval Analysis) [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF] such that :

X -⊂ X ⊂ X + (11)
SIVIA is a branch-and-bound algorithm with several implementations. It takes an initial box [x 0 ] large enough to contain the solution, and applies on it contractors (or separators) to remove parts which can be classified into X or X. If the width of the resulting box is larger than a threshold ǫ the box is bisected (cut into 2 sub-boxes) and the algorithm is run on each new box. The complexity of the algorithm is exponential with respect to the dimension n of the initial box, and terminate after less than ( w([x 0 ]) ǫ + 1) n iterations. A recursive version of the algorithm is given in algorithm 1 using separators and, another using a list and only an outer contractors is given in algorithm 2.

Algorithm 1: SIVIAS(in: [x], S, ǫ, InOut: X -, X + ) 1 {[x in ], [x out ]} = S ([x]); 2 Store [x in ]\[x] into X -and also into X + ; 3 [x] = [x in ] ∩ [x out ]; 4 max ← a 1 ; 5 if w ([x])) ≤ ǫ) then 6 store [x] into X + 7 else 8 bisect [x] into [x 1 ] and [x 2 ]; 9 SIVIAS([x 1 ], S, ǫ, X -, X + ); 10 SIVIAS([x 2 ], S, ǫ, X -, X + ); Algorithm 2: SIVIA(in: [x 0 ], C out , ǫ; Out:X + ) 1 L ← [x 0 ]; 2 while L = ∅ do 3 Pull [x] from L; 4 [x] = C out ([x]); 5 if w ([x])) ≤ ǫ) then 6 store [x] into X +

Relaxed Intersection

Parameters estimation using interval analysis is more efficient if the number of sets involved is greater than the number of unknowns. However the presence of outliers will yield an empty solution set. To solve this issue, the relaxed-intersection [START_REF] Jaulin | Robust set membership state estimation ; application to underwater robotics[END_REF] has been introduced.

Consider N sets X 1 , . . . , X n of IR n . The q-relaxed intersection denoted by {q} X i is the set of all x ∈ R n which belong to all X i 's, except q at most. Figure 3 illustrates this notion for N = 6 and q = 2, 3, 4.

Since the relaxed-intersection can be written as a combination of unions and intersections, it is monotonic with respect to X i and q:

q 1 , q 2 ∈ N * , q 1 ≤ q 2 ⇒ {q 1 } X i ⊆ {q 2 } X i (X 1 ⊂ Y 1 , . . . , X N ⊂ Y N ) ⇒ {q} X i ⊂ {q} Y i (12)
Computing the exact q-intersection of p n-dimensional boxes has a complexity in O (p n ). However, using the Marzullo's algorithm ( [START_REF] Chew | Masking failures of multidimensional sensors[END_REF]), an approximation can be computed in O (np log(p)) by calculating the intersection along each dimension independently. A study of the relaxedintersection can be found in [START_REF] Carbonnel | Q-intersection Algorithms for Constraint-Based Robust Parameter Estimation[END_REF] including a better approximation in O np 2 . From a set of N contractors C 1 , . . . , C N associated with sets X 1 , . . . , X N , the contractor C q is defined by algorithm 3.

Algorithm 3: C q (in: q, C 1 , . . . , C N , inout: [x]) 1 for k from 1 to N do 2 [x k ] = C k ([x]); 3 [x] = [x] ∩ {q} i∈{0,...,N } [x i ];
Knowing q, the relaxed-intersection is mainly used inside the set inversion algorithm instead of the classical intersection to take into account at most q outliers.

Guaranteed Outlier Minimal Number Estimator

In actual cases, the percentage of outliers is not known and needs to be estimated concurrently with the parameters. The Guaranteed Outlier Minimal Number Estimator (GOMNE ) algorithm, introduced in [START_REF] Jaulin | Guarateed Robust nonlinear parameter Bounding[END_REF] and [START_REF] Jaulin | Guaranteed robust nonlinear estimation with application to robot localization[END_REF], has been designed for this purpose.

GOMNE Algorithm

Let's take N sets (X 1 , . . . , X N ) defined by constraints such as, for instance, f i (x) ∈ [y i ]. Let j be a cost function which returns the number of measurements inconsistent with x.

j : R n -→ [[1, N ]] x -→ card {i|f i (x) / ∈ [y i ]} (13)
The GOMNE aims at characterizing the set S * of feasible parameters which minimizes the number of outliers :

S * = arg min x∈R n j (x) (14)
Using the relaxed intersection, the set of all parameters that are consistent with at least Nq measurements can then be seen as a set inversion :

S q = {x ∈ R n |j(x) ≤ q} = j -1 ([0, q]) (15)
GOMNE runs a finite sequence of set inversions with different values of q in order to enclose S q * between two sub-pavings S - q * and S + q * in a guaranteed way. With an initial box [x 0 ], it starts with q = 0 and calculates S - q and S + q by running SIVIA([x 0 ],C q ,ǫ). If the outer sub-paving S + q is empty, there does not exist any solution consistent with at least Nq data. The number of outliers q is incremented and SIVIA is run again. If S + q is not empty and S - q is, a smaller value of ǫ is used to find a non empty inner set or an empty outer set. In the first case, the algorithm ends and q * = q, otherwise SIVIA is run again with a greater value of q. The full GOMNE algorithm is given in listing 4. A minimal value ǫ 0 of ǫ has been introduced in order to avoid infinite loop when the interior of S * is empty which could happen when, for instance, S * is a singleton.

Algorithm 4: GOMNE(in: [x 0 ] , ǫ 0 , ǫ min , out: q * , S - q * , S + q * ) 1 q = -1 ; 2 repeat 3 ǫ = ǫ 0 ; q = q + 1; 4 repeat 5 S - q , S + q = SIV IA([x 0 ],q,ǫ); 6 ǫ = ǫ/2; 7 until S - q = ∅ or S + q = ∅ or ǫ < ǫ min ; 8 until S - q = ∅;

Illustration : Simple Map-Based Localization

A simple localization problem within an initial 2D map illustrates how GOMNE can be used and its limitations. In order to build an inner and outer contractor, the map must divide the space into two complementary sets, with non-zero volume: M (gray in figure 4a) and M. Measurements will belong to the border of M. In our case, the map is composed of a non-convex polygon with a hole, which for instance represents a room. The robot R is equipped with a laser sensor which returns the distance

d i ∈ [d - i , d + i ] between R and the nearest obstacle in a given direction (α i ∈ [α i -ǫ α , α i + ǫ α ]
). We also denote by a i ( resp. b i ) the nearest (resp. furthest) impact point associated with the distance d - i (resp. d + i ). The parametrization is shown in figure 4b. The goal of this problem is to estimate the position x = (p x , p y ) of the robot on the map. Its heading θ is assumed to be known and is not estimated in this example.

Let f (x, α, d) be the function which translates the point x by a distance d in the direction given by α + θ.

Let f i -(x) = f x, α i , d - i (resp. f i + (x) = f x, α i , d + i
) be the function which transforms the x into a i (resp. b i ). The set X i of feasible configurations consistent with the i th measurement is the set of all x such that the point a i belongs to M and b i belongs to M given by:

X i = {x ∈ R n |f i -(x) ∈ M and f i + (x) / ∈ M} ie. X i = f -1 i -(M) ∩ f -1 i + M (16)
The complementary set X i is defined by :

X i = f -1 i -M ∪ f -1 i + (M) (17)
With the contractor arithmetics it is possible to build the separator

S = {C in , C out } such as C in ∼ X i and C out ∼ X i .
The GOMNE algorithm is then applied to find the solution set X * which minimizes the number of outliers:

X * = {q * } i X i (18)
Results of the set inversion are shown in figure 5 with two values of bounds for α i . In the top right sub-figure, the value of α i is assumed to be precisely known (ǫ α = 1e -4 rad.) and the resulting solution set has an inner part. On the contrary, in the bottom right sub-figure, the interval for α is larger (ǫ α = 1e -2 rad.) and the inner part disappears. This situation could occur frequently if a sonar, which has a larger emission cone than a laser, is for instance used. In this case, the GOMNE algorithm fails to find a non empty X -even if a smaller value of ǫ is used. Maybe, the condition for x to belong to X i is too restrictive because the equation ( 16) must be satisfied for all α i in [α i ]. For that purpose the projection algorithm [START_REF] Chabert | Contractor Programming[END_REF] should be used to find the set X

∩[α i ] i
defined by:

X ∩[α i ] i = x ∈ R n |∃α i ∈ [α i ], f -1 i -(M) ∩ f -1 i + M (19)
But this requires additional bisections on α i and hence increases the complexity of the algorithm. However, depending on the problem, it may be useless to precisely characterize the inner set. For instance, with the localization problem, in order to ensure a guarantee on the solution, only X + is taken into account. Finally, the inner contractor may not exist and the GOMNE algorithm can not be used. This is main reason why a new method is introduced. With only an outer contractor, the classical GOMNE approach fails to find a guaranteed result because the stopping condition (S -= ∅) can not be satisfied. To cope this lack of constraint, we propose to use the local information given by j (equation 13) to find the smallest majorant of q * as possible. The outer contractor will remove inconsistent parts of the initial space while a local method looks for q * . Let's take a set of N measurements, each of them defining a set X i . Let C i out be the outer contractor associated with the set X i . The method is based on the following propositions: proposition 1: Based on the definition of the outer contractor, we have :

∀[x] ∈ IR n , C i out ([x]) = ∅ ⇔ [x] / ∈ X i ( 20
)
which means that if C i out ([x]) = ∅, [x]
is not consistent with the i th measurement. In the special case where [x] is a singleton denoted {x}, we define the following function :

R n → [0, 1] x → µ i (x) = 1 if C i out ({x}) = ∅ 0 otherwise (21)
proposition 2: An explicit expression for j : R n → [[0, N -1]] is given by :

j (x) = N i=1 µ i (x) (22)
remark 1: Punctual values are used in j instead of intervals in order to limit the pessimism during the evaluation process. remark 2: Instead of having additional information, properties of j such as continuity, variations, monotonicity are unknown: j can only be evaluated at a given value.

Properties: Based on the construction of j and on the monotonicity of the relaxed intersection we have :

S * = ∅ ⇔ ∃x ∈ R n , j (x) = q * (23) ∀x ∈ R n , S q * ⊆ S j(x) (24)
∀x / ∈ S * , j (x) > q * (25)

With equation [START_REF] Ramalingam | SKYLINE2GPS: Localization in Urban Canyons using Omni-Skylines[END_REF], the minimum of j can always be found if the search algorithm has no time limit. Otherwise, with the property [START_REF] Radu | 3D is here: Point Cloud Library (PCL)[END_REF], it can be stopped at any time and the guarantee is preserved. With the last property, algorithms can be run safely with any x ∈ R n . It is now possible to define a contractor for [q]. Given a subset A of R n :

C j A ([q]) = [q] ∩ [0, min x∈A j(x)] (26)
By construction C j A is monotonic, but not always minimal, depending on the result of the minimization algorithm used. If the minimum is very hard to find, q * will be overestimated.

To minimize the function j several strategies can be used. A trivial way is to randomly test N points in A and keep the minimal value. The probability of finding the minimum increases when N tend towards infinity or when the volume of A approaches towards 0. A smarter way is to consider j as a likelihood function and use the particle filter framework to estimate S * . Alternatively, the issue of finding the best value of j could be seen as a black-box optimization problem because no information are available on j. We choose to use the NOMAD library [START_REF] Digabel | Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm[END_REF][START_REF] Abramson | The NOMAD project[END_REF] which is a C++ implementation of the Mesh Adaptive Direct Search (MADS) algorithm [START_REF] Audet | A progressive barrier for derivative-free nonlinear programming[END_REF][START_REF] Audet | Mesh adaptive direct search algorithms for constrained optimization[END_REF] designed for constrained optimization of black-box functions.

Algorithm

The initial GOMNE algorithm is modified in order to take into account the new way to estimate q * . At the beginning, [q] = [q -, q + ] is unknown and ranges between 0 and N -1. Similarly to the GOMNE algorithm a set inversion is performed. If the resulting sub-paving S + is empty, there are at least q outliers and q -= q + 1 else the contractor C j S + is called to find a better value for q + . The algorithm stops when q -≥ q + , and then q * = q + . Remarks:

• The returned value is not always optimal and depends on the result of the minimization.

However, if S * is small, the initial search space is also small and an optimal value is found quickly. So, after the optimization, if no improvement of the value of q is found (line 8), a smaller accuracy (ǫ) is used to reduce the size of boxes which composed the S + . Then SIVIA is called until S * becomes empty or ǫ < ǫ lim where ǫ lim is chosen by the user.

• When the number of outliers is greater than 50% testing each value of q could be very slow. To solve this issue, when S * is empty, q can be incremented by ∆q and the best value which satisfies q ->= q + will be found quickly. However, this value will not always be optimal, but satisfies q +q * ≤ ∆q.

• At any time, the algorithm could be stopped and the solution set S + q + will contain the solution in a guaranteed way.

Algorithm 5: Outer-GOMNE(in: X 0 , ǫ 0 , ǫ lim , out:

q + ) 1 q -= 0 , q + = N -1, ǫ = ǫ 0 ; 2 while (q + > q -) do 3 S + = SIV IA([x 0 ], C q -, ǫ); 4 if (S + = ∅ or ǫ < ǫ lim ) then 5 q -= q -+ 1, ǫ = ǫ 0 6 else 7 q = min x∈S + j(x); 8 if (q < q + ) then 9 q + = q 10 else 11 ǫ = ǫ/2
Results of this algorithm are shown in the next section. A compromise needs to be found between the time allocated to set inversion and to optimization.

Test Case and Results

Introduction

Let's take a robot R in an environment composed of lines and circles, which could for instance represent trees and buildings. An initial map containing only the border of obstacles is provided but in this case, no inner and outer parts can be defined (the volume of M is null). To our knowledge, no inner contractor has been found for this kind of situation and only an outer contractor is available. The robot is equipped with a laser sensor which returns the distance d i between R and the nearest obstacle in a given direction (α i ). Let y i = (d i , α i ) a measurement. The pose x = (p x , p y , θ) of R must satisfy :

g (x, y i ) ∈ M ( 27
)
where g is the function which translates the point (p x , p y ) by a distance d in the direction given by α + θ defined by:

g : R 3 × R 2 -→ R 2 (p x , p y , θ, d i , α i ) -→ q x q y = d i . cos (α i + θ) sin (α i + θ) + p x p y
To ease the reading, g(x, y i ) is denoted by g i (x). The map is represented by a binary occupancy grid and the constraint g i (x) ∈ M is implemented with the image contractor (see 4.3.2 for a short introduction), which is very efficient for this kind of application. In addition, the localization algorithm has to deal with erroneous measurements and detection of unmapped objects such as people, trees or cars. The situation is depicted in figure 6. The solution set X is then defined by:

X = {q * } g -1 i (M) (28)
where q * = min x∈R 3 j(x).

Map contractor

Unstructured environments can be described by an occupancy grid represented as a binary image i in which the value of each pixel is 1 if it is occupied, 0 otherwise. The image contractor introduced in [START_REF] Sliwka | Using set membership methods for robust underwater robot localization[END_REF][START_REF] Guyonneau | Méthodes ensemblistes pour la localisation en robotique mobile[END_REF] can then be used. This contractor is based on the summed area table (also known as integral image), widely used in computer vision, defined by:

I (x, y) = x ′ ≤x y ′ ≤y i x ′ , y ′ (29)
Using the integral image, the number of 1-valued pixel contained in any rectangular region can be computed in four operations. Let φ be the function which returns the number of occupied cells in a given box

[x] = [x - 1 , x + 1 ] × [x - 2 , x + 2 ]
aligned on the grid. Let A(x -, y -), B(x + , y -), C(x + , y + ), D(x -, y + ) be the coordinates of the rectangle corners, we have :

IN 2 -→ N φ ([x]) -→ I(A) + I(C) -I(B) -I(D) (30)
On figure 7a

             y - 1 = max(x ∈ [x 1 ], φ([x - 1 ; x] × [x 2 ]) = 0) y + 1 = min(x ∈ [x 1 ], φ([x; x + 1 ] × [x 2 ]) = 0) y - 2 = max(x ∈ [x 2 ], φ([x 1 ] × [x - 2 ; x]) = 0) y + 2 = min(x ∈ [x 2 ], φ([x 1 ] × [x; x + 2 ]) = 0) (31)
The min and max can be computed using a dichotomy, which has a logarithmic complexity.

Results

A trajectory of 20 configurations is generated, for each of which the number and the position of unmapped objects changes (case of a very dynamic environment). All different configurations are shown in figure 11b. The heading θ th is assumed to be known with an accuracy of ± 30 degrees while their position is completely unknown. The optimization procedure use MADS algorithm [START_REF] Audet | Nonsmooth optimization through mesh adaptive direct search and variable neighborhood search[END_REF] to escape local minima. For each box of the resulting sub-paving S + an optimization is run with a maximum of 100 tries. Finally the ǫ used in SIVIA is ǫ = (1m, 1m, 3 • ). The algorithm is run on a Intel(R) Core(TM) i5-2450M CPU at 2.50GHz and the computing time for each situation is shown in figure 8a. The width of the box which enclosed S q * and the error between its center and x th are given in table 1.

Outer-GOMNE retrieves, in a guaranteed way, the set S q * . Nevertheless, due to undetected and spurious outliers, the true solution x th sometimes do not belong to the resulting solution set S q * . This is the case for the 4 t h configuration in figure 9a. To avoid this, like the original GOMNE algorithm, an arbitrary number of undetected outliers r can be added to the returned value q * to characterize the set S q * +r , which is more likely to contain x th . However, when the proposition of outliers is higher than 50%, another configuration consistent with a set of outliers may exist in the initial box. This is the case of the 15 th situation in figure 9a, localized as shown in figure 9b. An additional constraint such as "a beam can not cross the map" can be added to avoid this situation but, in the general case, without additional information, such situation 1: Results of each configuration. The size of the union of all boxes of the sub-paving is given and the error between the theoretical position and the center of this union of the sub-paving. Bold lines are configuration inconsistent with the theoretical position.

can not be avoided. Results of the set inversion with a smaller initial box are shown in figure 8b with r = 1. Here the computing time is lower and the algorithm never fails.

During our tests, Outer-GOMNE was able to cope with more than 70% outliers. Its main drawback, like its original version, is that it looks for the set which minimizes the number of outliers and, if only static information are used, this set may not contain the true solution. To address this issue, the fact that successive configuration constitute a trajectory can be considered. 5 Extension to Trajectories

Principle and Algorithm

The previous method can be applied as contractor for tubes ( which is an interval of functions [START_REF] Bethencourt | Solving non-linear constraint satisfation problems involving time-dependant functions[END_REF]) to retrieve the trajectory of the robot by using its dynamics or any independent relative motion estimation measure such as provided by odometers or inertial navigation systems. Trajectories contain much more information which can deal with symmetric and ambiguous situations, especially when the true solution does not belong to S q * +r . Let's take a system which provides N k outputs at each time step k described by :

x k = f (x k-1 , u k ) + w k y k = g k (x k ) ( 32
)
where

y k = (y 1 k , . . . , y N k k ) and g k (x k ) = g 1 1 (x k ), . . . , g N k k (x k ) . A state noise w k is assumed to belong to the box [w(k)].
Each output contains q k outliers. Let p k = (x k , q k ), the vector of parameters to be estimated for time k.

Constraint propagation can be used to estimate the trajectory z = (p 0 , . . . , p N k ) of the system. Let C ik be the outer contractor associated to the constraint

g i k (p k ) ∈ [y i k ]. We define the contractor C k which contract a box [p k ] = [x] × [q k ] : C k : IR n × IN -→ IR n × IN ([x], [q]) -→      [q] = [q] ∩ [0, min x∈[x] j(x)] [x] = [x] ∩ q + C ik ([x]) (33)
The algorithm 5 inspired from [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF] contracts an unknown trajectory [z] = ([x 0 , . . . , x N ]). For each time step k, [x k ] is predicted from the previous [x k-1 ] using equation (32). Then, C k is called to estimate the number of outliers and contract [x k ] with this value. Using interval arithmetic, a backward propagation is then done. This contractor procedure can be called several times until it has no effect on [z] (the fixed point is then reached).

Algorithm 5: C traj (in: C k , f inout :[z] = (p 0 = [x 0 ] × [q 0 ], . . . , p k = [x k ] × [q k ])) 1 for k from 1 to N do 2 [x k ] = f ([x k-1 ]); 3 C k ([p k ]) see equation 33; 4 for k from N-1 to 0 do 5 [x k ] = [x k ] ∩ f -1 (x k+1 )
In the end, the trajectory can be retrieved only with contractions and optimizations. However, if one step becomes inconsistent with the other the resulting trajectory will be empty. To avoid this, the relaxed intersection can be used on the trajectory space. Let z i be the trajectory corrected only with the output of the system at time k. The final trajectory is now defined by

z = qz z i (34) w 0 ([x])(m) w 0 ([y])(m) w 0 ([θ])(

Test Case

Consider the whole trajectory of the previous test case. The observation model remains the same and the motion model is defined by:

x k = x k-1 + u k + w k (35)
where u k is the movement of the robot between time t k-1 and t k . An additive noise of more or less 10% of the input vector value is used for the prediction :

[w k ] = [-0.1 * u k , 0.1 * u k ].
We allow also 10% of spurious states ( q z = 2). Table 2 shows the computing time needed to retrieve a 65 meter trajectory long with different sizes of the initial box, and different numbers of run used by the optimization algorithm. Figures 11a and11b show the estimated trajectory. Last steps are less accurate because less information constrains the position. Instead of processing the whole trajectory in one run, it is possible to apply the algorithm on a sliding temporal window and compute the trajectory while the robot is moving. With a 10 time steps window, and a small box containing the true position, the computing time is less than 12 seconds. Outer-GOMNE can be used instead of C k .

Localization in a Digital Elevation Map

DEM's represent a surface z = f (x, y) on a regular Cartesian grid where each cell contains the value of the elevation z. Depending on the data and algorithms used to build such a structure, additional information such as a standard deviation and extrema values for z can be encoded in each cell. DEM's are widely used, either in geographic data systems because they present a good compromise between expressiveness, simplicity and compactness1 , or in robotics, where they are easily built from range data and exploited to assess traversability, to plan itineraries or trajectories, or to assess visibilities between positions.

If a DEM is a nearly perfect structure to represent surfaces, it fails to represent verticals and overhangs, and this has numerous consequences, especially when it comes to associate range data acquired by the robot to the map. Furthermore, the presence of unmapped elements in the environment, and the fact that some objects may only be partially mapped make the localization challenging. This section presents how Outer-GOMNE can cope with these issues, using range data acquired by the robot.

Initial Models

Map Structures

With a DEM, the space can be split into two sets: the one over the surface (zf (x, y) > 0) and the one below (zf (x, y) < 0). In such cases, an inner and an outer contractor can be defined, and the same formalism of section 3.2 can be applied to localize the robot using range data, using the classical GNOME approach to cope with the presence of outliers. DEMs are particularly adapted to open terrains for which a surfacic representation is well suited, or in more complex (non-surfacic) terrains, in the case where the range data used to localize the robot are gathered in conditions similar to the data from which the DTM has been built.

For instance, if a DEM is built from an Unmanned Aerial Vehicle (UAV) using a downwards oriented sensor (see figure 12b), range data acquired by other UAVs using a downwards oriented sensor can be well matched with the DEM structure. But on the contrary, range data acquired by a ground robot will hardly match such a map: e.g. the tree canopy is mapped by the UAV, whereas the ground robot only perceives tree trunks and lower branches and leaves. From the point of view of a UAV, the representation as a surface (convex world) is a good approximation. However from the point of view of the ground robot, the elevation between two adjacent cells is assumed to be continuous, which leads to a poor representation of verticals and overhangs.

As a compromise, the DEM can be voxelized, by filling the gaps between two different elevations of neighboring cells. The figure 13 shows the resulting environment model after applying this process to the DEM of figure 12b. Note that this representation also allows the use of an image contractor extended to the 3D case. But if this copes for verticals such as walls, it also clearly introduces erroneous information: the extension to the ground of the tree canopy does not fit the reality and hides the tree trunks for instance.

To better solve this issue, a volumetric representation using a three-dimensional voxel grid can be used: non surfacic elements are then accurately represented. However, the building of such models has a high computational cost and requires a lot of memory. And still, some objects may not be entirely mapped in the model: e.g. sometimes only parts of buildings, cars, trees, are modeled and the resulting model does not define closed connected components, as can be seen in the model of figure 14, where only walls perceived by the robot are modeled -and not the entire buildings. The hypothesis of a convex world then does not hold, and no inner contractor are available: in such a situation, Outer-GOMNE needs to be used.

Test Maps

.

For our localisation tests, three different models of the same area (training camp of Caylus, figure 12a) are used. The first model, denoted "atLaas", is built from a set of 3D scans acquired by a Velodyne Lidar mounted on a ground robot (see figure 14). The second and third models are built from data acquired by an UAV equipped with a camera (15a) and a Lidar (figure 15b), and are respectively denoted "UAV-Vision" and "UAV-Lidar".

The voxelization is done by marking as occupied all cases such as z = f (x, y) and filling gaps. The formulation is :

i(x i , y i , z i ) =        1 if z i ∈   min x∈x i-1 ,x i y∈y i-1 ,y i f (x, y), max x∈x i-1 ,x i y∈y i-1 ,y i f (x, y)   0 otherwise (36)
6.2 Data Pre-processing Before using the Outer-GOMNE algorithm, the range data acquired by the ground robot Lidar need to be preprocessed. Indeed, the complexity of the algorithm is linear with respect to the number of measurements, and a 360 • panoramic scan indeed contains up to 288.000 points. Also, in the context of our localization problem where only p x , p y , ϕ need to be retrieved, measurements which are likely to correspond to the ground are removed. First, the density of the point cloud dramatically decreases with the distance to the sensor: the initial point cloud (figure 16a) is down-sampled to homogenize its spatial density (see figure 16b). To select out the points corresponding to the ground, the normal vector of each point is estimated using a circular neighborhood support. We assume that a point belongs to the ground if the angle between the vertical z axis and the estimated normal vector is lower than a given value (based on the analysis of histograms of normal vectors, the threshold is set to 30 • ). The resulting point cloud is shown in figure 16c. Finally, based on the distance between points, an euclidean clustering is performed to split the cloud onto spatially coherent sub-clouds (figure 16d). Table 3 indicate the time required by these processes, which are implemented with the Point Cloud Library [START_REF] Radu | 3D is here: Point Cloud Library (PCL)[END_REF].

Algorithm Extension to 3D Data

The algorithms presented in section 4 can readily be extended to the case of localization in the 3D world, as only the contractor on the DEM and the constraint g (equation 27) need to be adapted.

Let y i Bs = (x i , y i , z i ) be a point acquired by the sensor in its own frame (B s ). The goal is to estimate parameters x = (p x , p y , p z , θ, ψ, φ) T of the transformation matrix T B 0 →Bs (x) between yet a way to select relevant measurements has to be defined. 17.

Table 5 and figure 18 illustrate the localisation results of two other configurations A and B, with the three initial DEMs considered. In figure 18, the resulting sub-paving is displayed in yellow, the red triangle is the estimated position of the robot, and the range data that match the DEM in green, while the estimated outliers are in red. A straightforward observation is that using the aerial DEMs, Outer-GNOME can deal with the numerous model errors, that generate up to 59% outliers.

Figure 18a and 18b show the final situation when the "atLaas" DEM is used. As for the former trajectory, the robot is very accurately localized: the range data are indeed acquired in conditions close to the acquisition conditions to built the model, and contain less than 10% of outliers.

With the "UAV-Lidar" DEM, Outer-GOMNE succeeds to find the pose with a much higher proportion of outliers (figure 18c and18d). The "UAV-Vision" DEM contains a scale error and it is nearly impossible for the scan to match two parallel walls in the same time. As a result, the estimated sub-paving is made of three disconnected components, which correspond to different data pairing. Figure 18f show the data points reprojected considering one of this component: they are well aligned with walls of the building on the top, but other points are inside the middle building. Consequently, the number of outliers increases. An other component of the sub-paving, for which the data points are shown in figure 18e, corresponds to an unfeasible pose because the constraint used allows scan to pass through walls.

Outer-GOMNE seems to be robust to gross model errors and able to face complex situation. Of course, the better the initial model is and the shorter the initial box is, the faster and more accurate the localization is. 7 Discussion and Hints for Future Work

One of the main issue of map-based localization is the presence of outliers, that mainly come from the dynamic elements of the scene and limitations of the map structure. Being robust with respect to these outliers is of course essential: set membership approaches to this problem have shown to have such a robustness, the GOMNE algorithm being able to estimate the position parameters even in the case of an unknown bound on the outliers. However GOMNE requires an inner contractor, which can not be provided when the initial model is partial and exhibits non-closed components -which is the case in most operational applications. To cope with this, we proposed Outer-GOMNE, that combines a set inversion algorithm with an optimization method. With experimental results, we have shown that Outer-GOMNE can be applied to the absolute localization problem with Digital Elevation Maps of a 3D environment as initial models, even though the DEM data structure is not well suited to represent verticals and overhangs, and thus causes gross map errors.

Various improvements can be considered at the algorithmic level. In particular, monitoring the bisection process can save considerable time: a too small value of ǫ indeed yields numerous useless bisections in the SIVIA algorithm.

Yet Outer-GOMNE has similar drawbacks as GOMNE, as it and only looks for the solution set which minimizes the number of outliers. Most of the time, this set is the good one, but in scarce cases the true solution may not belong to this set. This is one of the main difference with the Bayesian estimation framework, which, by using infinite density of probability, is less precise but never fails. An arbitrary number of undetected outliers can be used to cope with this issue. But more interestingly, the incremental estimation of the whole trajectory will bring higher robustness.

7 else 8 bisect

 78 [x] and push into L the two resulting boxes

Figure 3 :

 3 Figure 3: q-relaxed intersection of 6 sets for q=2 (red), q=3 (green), q= 4 (blue) , q= 5 (yellow)

Figure 4 :Figure 5 :

 45 Figure 4: Initial situation, the laser sensor returns 56 measurements with 15 outliers (dashed lines)

Figure 6 :

 6 Figure 6: The map M , in black, is composed of lines and circles. The opened polygon does not define a closed area and consequently the map has an empty volume. Dashed red circles are unmapped objects which generated false detection. Consistent measurements are in green, while outliers are in red.

  the box contains two black pixel, so φ ([x]) = 2. From an initial box [x 0 ], the image contractor aims at finding the smallest box [x] included in [x 0 ] which contains exactly the same number of 1-valued pixels, ie. φ([x 0 ]) = φ([x]). Denote by C the image contractor. Consider [x] ∈ IN 2 and C([x]) = [y]. We have:

Figure 7 :

 7 Figure 7: Image contractor

Figure 8 :

 8 Figure8: Computing time needed to find S q * with different initial estimates, sorted by the percentage of outliers. The configuration's number is in white. Resulting sub-pavings, which does not contain the true solution, are marked with a red cross.

  Actual position (green robot) and estimated position (red robot) of the 15 th configuration Close-up of the estimated position

Figure 9 :

 9 Figure9: Result of the 15 th ambiguous situation. In figure9a, the initial configuration (green robot) has 73 outliers. Outer-GOMNE finds another configuration (in yellow) with only 67 outliers (figure9b). The constraint used allows beams to cross the map and leads an erroneous estimation.

  Recovered trajectory with the dynamic algorithm (b) Projection on the (x,y) plan

  (a) Aerial orthoimage of a test site (b) DEM

Figure 12 :

 12 Figure12: Aerial view of a 80m x 80m test site, and DEM built from a UAV with a downwards looking Lidar. Note that buildings roofs and tree canopy are disconnected from the ground surface.

Figure 13 :

 13 Figure 13: View of the DEM of figure 12b after its "voxelisation".

Figure 14 :

 14 Figure 14: The ground robot with a Velodyne Lidar in the site of figure 12a, and the corresponding "atLaas" DEM (0.1m resolution).

Figure 15 :

 15 Figure 15: DEMs of the test area of figure 12a generated by UAV data. Note that the DEM built from images contains gross errors, due to the absence of visual features matches: the elevation has been interpolated by the photogrametric algorithms to the ground elevation for most of the trees.

Figure 16 :

 16 Figure 16: Point cloud filtering: The raw point cloud (16a) is first down-sampled (16b), and based on the normal direction, points which are likely to belong to the ground are removed (16c). Then, the remaining point cloud is split into smaller but spatially coherent subsets (16d).

  Positions estimated for 8 configurations.

  ) = 50x50x60 • w(x) = 10x10x60 • w(x) = 10x10x2 • (b) Computation time with different sizes of the initial boxes.

Figure 17 :

 17 Figure 17: Localisation results using the "atLaas" DEM as the initial map (here shown on the "UAV-Vision" DEM), and execution times as a function of the size of the initial box. Light colored bars correspond to the time needed by the set inversion algorithm.

  (a) Situation A on atLaas DEM (b) Situation B on atLaas DEM (c) Situation A on UAV-Lidar DEM (d) Situation A on UAV-Lidar DEM (e) Situation B on UAV-Vision DEM -first component of the sub-paving (f) Situation B on UAV-Vision DEM -second component of the sub-paving

Figure 18 :

 18 Figure 18: Localisation results on the 3 different DEMs considered.

Table 2 :

 2 Computing time needed to retrieve the trajectory with different parameters

	rad) number of run time(s)

Table 4 :

 4 Characteristics of the resulting sub-paving (with an initial box of 50m x 50m x 60 • ) for the 8 position estimates shown figure

	1	1174	105 (9%)	0.78	0.81	82.3	0.20	9.8
	2	1154	34 (2%)	0.41	0.55	50.7	0.20	3.6
	3	786	35 (4%)	0.56	0.55	53.9	0.16	18.2
	4	738	85 (11%)	0.67	0.79	56.3	0.03	6.3
	5	805	73 (9%)	0.80	0.75	67.0	0.09	1.0
	6	912	89 (9%)	0.81	0.67	56.3	0.20	6.3
	7	1038	84 (8%)	0.76	0.75	80.0	0.10	18.2
	8	851	45 (5%)	0.63	0.68	67.0	0.28	24.8

Table 5 :

 5 Results of Outer-GOMNE for two positions estimated on three different DEMs.

			situation A			situation B
			414 measurements		293 measurements
		atLaas UAV-Vision UAV-Lidar	atLaas UAV-Vision UAV-Lidar
	N outliers	26 (6%)	245 (59%)	222 (53%)	10 (3%)	158 (54%)	102(34%)
	time(s)	40.33	122	63.1	7.68	105	8.38
	w([p x ]) m	0.72	0.80	1.11	0.36	3.67	0.81
	w([p y ]) m	0.89	0.97	1.39	0.37	4.54	0.76
	w([θ]) rad	0.071	0.064	0.200	0.060	0.04	0.035

For instance, most commercial aerial mapping systems generate DEMs along with orthoimages

 16a.

the global frame B 0 and the sensor frame B s (ϕ, θ, ψ respectively denote the yaw, pitch and roll). Like in the 2D test cases, the point y i B 0 must belong to the map, so the constraint g is now:

The localization problem is then formalized as a robust set inversion :

where Outer-GOMNE can be used to concurrently determine x and q * . Based on the generalized notion of integral image in higher dimensional space (see [START_REF] Tapia | A note on the computation of high-dimensional integral images[END_REF]), the image contractor can be used with a voxel grid which represent the map.

Results

The approach has been evaluated in different configurations using two range data sets acquired in two different days, with ground truth positions provided by a centimeter accuracy RTK-GPS. One data-set has been used to produce the "atLaas" DEM (figure 14), while the second data-set is used to test the localisation algorithm in this DEM and the two "UAV-Vision" and "UAV-Lidar" DEMs derived from UAV data of the site.

Localisation results on the "atLaas" DEM for 8 robot configurations are shown in figure 17a, execution time with different sizes of the initial box are given in figure 17b (measured on a Intel(R) Core(TM) i5-2450M CPU at 2.50GHz), and the characteristics of the resulting sub-paving are given in table [START_REF] Audet | Mesh adaptive direct search algorithms for constrained optimization[END_REF].

In each case, in accordance with the theory, the integrity of the result is satisfied and the precision of the retrieved position is smaller than one meter. Table 4 also shows the errors err xy and err θ , which have been computed by comparing the ground truth with the center of the box that bounds the resulting sub-paving.

Computing times are higher than for the 2D test case: on the one hand, the function involved in the set inversion is more complex and introduces more pessimism: hence the contraction is less efficient and more bisections are made. To contraction process can be improved by using more constraints on the pose of the robot. On the other hand, the complexity of the algorithm is linear with respect to the number of measurements used. We use the point cloud given by the segmentation algorithm which contains on average 1000 data points (see table 4 column 2). Using less measurements would speed up the algorithm and especially the optimization part -


A Seuret
email: aseuret@laas.fr

F Gouaisbaut

Hierarchy of LMI conditions for the stability analysis of time delay systems

Keywords: Bessel inequality, Jensen's inequality, stability analysis, time-delay systems, Lyapunov-Krasovskii functionals

Assessing stability of time-delay systems based on the Lyapunov-Krasovskii functionals has been the subject of many contributions. Most of the results are based, first, on the design of more and more involved class of functionals and, finally, on the use of the famous Jensen's inequality. In contrast with this design process, the present paper aims at providing a generic set of integral inequalities which are asymptotically non conservative and then to design functionals driven by these inequalities. The resulting set of stability conditions forms a hierarchy of LMI which is competitive with the most efficient existing methods (delay-partitioning, discretization and sum of squares), in terms of conservatism and of complexity. Finally, some examples show the efficiency of the method.

Introduction

In the past two decades, a large number of papers appeared on stability of time delay systems [START_REF] Gu | Stability of Time-Delay Systems[END_REF][START_REF] Niculescu | Delay Effects on Stability. A Robust Control Approach[END_REF][START_REF] Sipahi | Stability and stabilization of systems with time delay[END_REF]. Several research lines have been explored with success, regarding the case of a linear system with a single delay. A first idea consists in the study of the characteristic equation associated with the delay system (see [START_REF] Gu | Stability of Time-Delay Systems[END_REF][START_REF] Sipahi | Stability and stabilization of systems with time delay[END_REF] and references therein). Even though these methods provide necessary and sufficient condition for stability, extensions to time-varying delays or robustness issue lead to unavoidable difficulties. Another interesting method relies on the modeling of the original delay system as a closed loop between a nominal system and the delay operator, which is embedded into a suitable uncertainty. Then robust analysis tools like Small gain Theorem [START_REF] Gu | Stability of Time-Delay Systems[END_REF][START_REF] Zhang | Stability of time-delay systems: Equivalence between Lyapunov and scaled small-gain conditions[END_REF], IQCs [START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF] or Quadratic separation [START_REF] Gouaisbaut | Bessel inequality for robust stability analysis of time-delay system[END_REF] allow to conclude. On the other side, one of the most fruitful fields of research in this area relies on the exhibition of Lyapunov-Krasovskii Functionals (LKF) for functional differential equations. Unlike the case of ordinary differential equations, the introduction of delay hugely complicates the task and was therefore the source of many theoretical developments [START_REF] Fridman | New Lyapunov-Krasvoskii functionals for stability of linear retarded and neutral type systems[END_REF][START_REF] Kolmanovskii | Stability of some linear systems with delay[END_REF]. For a single and constant delay, a Lyapunov functional candidate, called complete LKF [START_REF] Kharitonov | Lyapunov-Krasovskii approach to the robust stability analysis of time-delay systems[END_REF] leads to a necessary and sufficient stability condition. Nevertheless, the parameters, which compose the complete LKF, i.e. the matrix functions depend on the solution of a differential equation with boundary conditions and is recognized to be numerically difficult to handle especially for high dimensional systems [START_REF] Mondie | Complete Lyapunov-Krasovskii functional with a given cross term in the time derivative[END_REF]. A number of investigations then turns to approximating the parameters of the complete LKF. Hence, the problem has been then partially solved by considering a discretization scheme proposed originally by [START_REF] Gu | Stability of Time-Delay Systems[END_REF]. The interval of integration for the quadratic forms are partitioned into smaller delay domain and the structure of the matrix functions are chosen to be piecewise continuous within these intervals. It results a numerically tractable Linear Matrix Inequality (LMI) problem, which is efficient on some examples. More recently, this approximation method has been improved by considering polynomial-like parameter of arbitrary degree. By considering polynomial relaxations, the original problem of finding an LKF may be recasted into an SOS optimization scheme [START_REF] Peet | Positive forms and stability of linear time-delay systems[END_REF]. All these two approximation methods provide with only sufficient conditions, but potentially tends to be necessary if the discretization step or the polynomial degree are refined [START_REF] Peet | On the conservatism of the sum-of-squares method for analysis of timedelayed systems[END_REF]. On the other side, many papers have been concentrated on the research of a simpler LKF whose parameters are optimized via an LMI setup. The goal is often to study the trade-off between computational complexity of the resulting LMI and the pessimism induced by the simplification of the LKF. In the literature, two mainstream methods, often coupled in practice, have been followed. The first one proposed by [START_REF] Ariba | An augmented model for robust stability analysis of time-varying delay systems[END_REF][START_REF] Gouaisbaut | Delay-dependent stability analysis of linear time delay systems[END_REF][START_REF] Han | A discrete delay decomposition approach to stability of linear retarded and neutral systems[END_REF][START_REF] Shao | New delay-dependent stability criteria for systems with interval delay[END_REF][START_REF] Sun | Improved delay-range-dependent stability criteria for linear systems with time-varying delays[END_REF][START_REF] Sun | Delay-dependent stability and stabilization of neutral time-delay systems[END_REF] is to consider a simpler structure for the LKF than the previously mentioned complete LKF which is applied to some higher dimensional comparison systems. The additional degrees of freedom come from the consideration of extrastates like higher derivatives of the states [START_REF] Ariba | An augmented model for robust stability analysis of time-varying delay systems[END_REF][START_REF] Kim | Note on stability of linear systems with time-varying delay[END_REF], delayed states [START_REF] Gouaisbaut | Delay-dependent stability analysis of linear time delay systems[END_REF][START_REF] Gu | Stability of Time-Delay Systems[END_REF], integral delayed terms [START_REF] Sun | Delay-dependent stability and stabilization of neutral time-delay systems[END_REF]. Nevertheless, it is not so clear that considering some particular extrastate will be more accurate than choosing another. The second one remains in the use of technical lemmas to provide tight inequalities [START_REF] He | Further improvement of free-weighting matrices technique for systems with time-varying delay[END_REF][START_REF] Jiang | Delay-dependent robust stability for uncertain linear systems with interval time-varying delay[END_REF][START_REF] Kim | Note on stability of linear systems with time-varying delay[END_REF][START_REF] Park | Reciprocally convex approach to stability of systems with timevarying delays[END_REF]. Among them, Jensen inequality has been the core of many interesting improvements [START_REF] Briat | Convergence and equivalence results for the Jensen's inequality -application to time-delay and sampled-data systems[END_REF]. Combined with some delay partitioning approaches [START_REF] Gouaisbaut | Delay-dependent stability analysis of linear time delay systems[END_REF][START_REF] Han | A discrete delay decomposition approach to stability of linear retarded and neutral systems[END_REF], it leads to efficient criteria. More recently, several papers have employed a new inequality based on Wirtinger Lemma [START_REF] Liu | Wirtinger's inequality and Lyapunov-based sampled-data stabilization[END_REF][START_REF] Seuret | On the use of the Wirtinger's inequalities for time-delay systems[END_REF][START_REF] Seuret | Wirtinger-based integral inequality: Application to time-delay systems[END_REF]. It allows a significant reduction of the conservatism with a slight modification of the usual LKF. However, it seems, at the first sight, complicated to refine these techniques. In this paper, we concentrate on this second method and introduce a set of new integral inequalities, which encloses the Jensen inequality [START_REF] Gu | A further refinement of discretized Lyapunov functional method for the stability of timedelay systems[END_REF] and the Improved Wirtinger-based inequality [START_REF] Seuret | Wirtinger-based integral inequality: Application to time-delay systems[END_REF] as special cases. These inequalities introduce extra-states functions representing the projection of the delay state x t on the set of Legendre polynomials of degree less than a prescribed upper-bound N . Guided by these inequalities and these extra-states, we thus propose a new LKF. The resulting analysis leads to a hierarchy of sufficient LMI stability conditions for the pointwise delays and delay range cases. Finally, some academic examples show the effectiveness of our method even for non trivial examples. In all the cases, the numerical results tend to the analytical bound with a reasonable numerical burden compared to efficient techniques based on the discretization, delay partitioning or SOS optimization.

Notations: Throughout the paper R n denotes the ndimensional Euclidean space with Euclidian norm | • |, R n×m is the set of all n × m real matrices. The notation P ≻ 0, for P ∈ R n×n , means that P is symmetric and positive definite. The sets S n and S + n represent, the set of symmetric and symmetric positive definite matrices of R n×n , respectively. The set of continuous functions from an interval [-h, 0] ⊂ R to R n which are, consequently, square integrable is demoted as space C. . Moreover, for any square matrix A ∈ R n×n , we define He(A) = A + A T . The matrix I represents the identity matrix of appropriate dimension. The notation 0 n,m stands for the matrix in R n×m whose entries are zero and, when no confusion is possible, the subscript will be ommitted. For any function x : [-h, +∞) → R n , the notation x t (θ) stands for x(t + θ), for all t ≥ 0 and all θ ∈ [-h, 0]. The notation k l refers to the binomial coefficients given by k! (k-l)!l! .

Problem formulation

Consider a linear time-delay system described by:

ẋ(t) = Ax(t) + A d x(t -h), ∀t ≥ 0, x(t) = φ(t), ∀t ∈ [-h, 0], (1)
where x(t) ∈ R n is the state vector, φ is the initial conditions and A and A d , are constant matrices. The delay is assumed to be constant. In the context of the stability analysis of such systems using LKF, several types of functionals have been provided in the literature. All of them are composed by the sum of several typical terms [START_REF] Gu | Stability of Time-Delay Systems[END_REF]. Among them, one of the most relevant components which was introduced in [START_REF] Fridman | A descriptor system approach to H ∞ control of linear time-delay systems[END_REF] is an integral quadratic term given by

V (x t ) = 0 -h 0 s ẋT t (θ)R ẋt (θ)dθds,
where ẋt (θ) = ẋt (t + θ) represents the state of the timedelay system and h > 0 the delay. The positivity of such functionals is ensured by R ≻ 0. This class of Lyapunov-Krasovskii terms has been widely used in the literature mainly because the computation of its time derivative leads to conditions which depend explicitly on the value of the delay h. Indeed, differentiating this term with respect to the time variable t leads to

V (x t ) = h ẋT (t)R ẋ(t) - 0 -h ẋT t (s)R ẋt (s)ds. ( 2
)
This term is relevant to ensure the negativity of V (x t ) because of the negative contribution of the second term.

In order to transform (2) into a suitable LMI setup, this integral term should be expressed appropriately in terms of x t (0) and x t (-h). This is made possible thanks to the Jensen inequality (see for instance [START_REF] Gu | Stability of Time-Delay Systems[END_REF]). Naturally, it is likely to entail some inherent conservatism and several works have been devoted to the reduction of such a gap (see [START_REF] Briat | Convergence and equivalence results for the Jensen's inequality -application to time-delay and sampled-data systems[END_REF] and the references therein). Recently, in [START_REF] Seuret | On the use of the Wirtinger's inequalities for time-delay systems[END_REF][START_REF] Seuret | Wirtinger-based integral inequality: Application to time-delay systems[END_REF], an alternative inequalities called Wirtingerbased inequalities is employed in order to obtain a more accurate bound for this integral and therefore to reduce the conservatism of the resulting stability conditions. It has been shown that an inequality was able to encompass the Jensen's inequality and leads to less conservative results without increasing drastically the computation complexity of the derived LMI's. The present paper aims at extending the results provided in [START_REF] Seuret | On the use of the Wirtinger's inequalities for time-delay systems[END_REF][START_REF] Seuret | Wirtinger-based integral inequality: Application to time-delay systems[END_REF] by considering an alternative method based on the Bessel inequality and polynomial approximations. Thanks to this setup, a set of new integral inequalities is provided as showed in following section.

3 New integral inequalities

Legendre polynomials

In the following, a brief recall of the Legendre polynomials and their relevant properties is proposed.

Definition 1 The Legendre polynomials considered over the interval [-h, 0] are defined by

∀k ∈ N, L k (u) = (-1) k k l=0 p k l u + h h l . with p k l = (-1) l k l k+l l
.

The set of Legendre polynomials {L k , k ∈ N} forms an orthogonal sequence with respect to the inner product:

f, g = 0 -h f (t)g(t)dt, ∀f, g ∈ C. (3)
Hence, the Legendre polynomials described in Definition 1 satisfy the following properties:

Property 2 P1 Orthogonality:

∀(k, l) ∈ N 2 , 0 -h L k (u)L l (u)du = 0, k = l, h 2k+1 , k = l. (4) P2 Boundary conditions: ∀k ∈ N, L k (0) = 1, L k (-h) = (-1) k . P3 Differentiation: Lk (u) = 0, k = 0, k-1 i=0 (2i+1) h (1 -(-1) k+i )L i (u), k ≥ 1.

Bessel-Legendre inequalities

Based on the Legendre polynomials, the following lemma is derived.

Lemma 3 Let x ∈ C and R ∈ S + n and h > 0. Define I(x, h) as follows

I(x, h) := 0 -h x T (u)Rx(u)du.
Then, the inequality

I(x, h) ≥ 1 h N k=0 (2k + 1)Ω T k RΩ k , (5)
holds for all N ∈ N, where

Ω k = 0 -h L k (u)x(u)du, k = 0, . . . N . Proof : Consider a function x in C, a matrix R in S + n and h > 0. Define the function z N by z N (u) = x(u) - N k=0 2k + 1 h Ω k L k (u)
Clearly, z N is in C and it represents the approximation error between x and its projection to the polynomial set {L k , k = 0, . . . , N } with respect to the inner product (3). The quantity

I(z N , h) = 0 -h z T N (u)
Rz N (u)du exists and the orthogonal property of the Legendre polynomials yields

I(z N , h) = 0 -h x T (u)Rx(u)du -2 N k=0 2k+1 h 0 -h L k (u)x(u)du T RΩ k + N k=0 2k+1 h 2 0 -h L 2 k (u)du Ω T k RΩ k .
Finally by recalling that

Ω k = 0 -h L k (u)x(u)du and 2k+1 h 2 0 -h L 2 k (u)du = 2k+1 h , it yields I(z N , h) = I(x, h) - N k=0 2k+1 h Ω T k RΩ k .
Finally, the inequality ( 5) is obtained by noting that I(z N , h) > 0 since R ≻ 0. Notice that this proof is related to the Bessel inequality applied to x ∈ C equipped with the inner product (3). Hence, inequality (5) will be recalled latter on as Bessel-Legendre (B-L) Inequality. ♦

A suitable corollary for the stability analysis of time-delay systems

As it was mentioned in the introduction, the problem is often to derive a lower bound of 0 -h ẋT (u)R ẋ(u)du. The next corollary addresses this particular problem.

Corollary 4 Let x be such that ẋ ∈ C, R ∈ S + n and h > 0. Then, the integral inequality

I( ẋ, h) ≥ 1 h ξ T N N k=0 (2k + 1)Γ N (k)RΓ N (k) ξ N , (6)
holds, for all integer N ∈ N, where

ξ N =          [x T (0) x T (-h)] T , if N = 0, x T (0) x T (-h) 1 h Ω T 0 . . . 1 h Ω T N -1 T , if N > 0, Γ N (k) =        I -I , if N = 0, I (-1) k+1 I γ 0 N k I . . . γ N -1 N k I , if N > 0. γ i N k = -(2i + 1)(1 -(-1) k+i ), if i ≤ k, 0, if i > k.
and where Ω k is defined in Lemma 3.

Proof : According to Lemma 3 to the order N , the following inequality holds

I( ẋ, h) ≥ 1 h N k=0 (2k + 1) Ωk R Ωk , (7)
where Ωk = 0 -h L k (u) ẋT (u)du, for all k = 0, 1, . . . , N . An integration by parts ensures that, for all k ≥ 0

Ωk = L k (0)x(0) -L k (-h)x(-h) - 0 -h Lk (u)x(u)du.
Thanks to properties P2 and P3 of the Legendre polynomials, the following expression is derived

Ωk = x(0) -(-1) k x(-h) - k-1 i=0 γ i N k Ω i = Γ N (k)ξ N . (8
) Replacing Ωk by its expression using the matrices Γ N (k), k = 1, . . . , N leads to [START_REF] Fridman | A descriptor system approach to H ∞ control of linear time-delay systems[END_REF], which concludes the proof. ♦ Remark 1 An interpretation of the B-L inequality in the context of robust analysis is provided in [START_REF] Gouaisbaut | Bessel inequality for robust stability analysis of time-delay system[END_REF].

Remark 2 Considering the B-L inequality with N = 0 allows retrieving the celebrated Jensen Inequality, which has been widely use in the literature. The case N = 1 has been already proposed in [START_REF] Seuret | Wirtinger-based integral inequality: Application to time-delay systems[END_REF]. The proof of this inequality was derived by application of the Wirtinger inequality. Unfortunately, this first approach did not suggest such a generic method to generate more inequalities. This corresponds to one of the main achievement of the present paper.

4 Application to the stability analysis of systems with a discrete delay

Systems with constant and known delay

We present in this sub-section a first stability result for time-delay systems, which is based on the use of the B-L inequality developed in the previous section. The following stability theorem is provided by the use of Corollary 4 with an arbitrary N .

Theorem 5 For a given integer N and a constant delay h, assume that there exist a matrix P N ∈ S (N +1)n and two matrices S, R ∈ S + n such that the LMI

Θ N (h) = P N ≻ 0, if N = 0, P N + 1 h diag(0, S N -1 ) ≻ 0, if N > 0, Φ N (h) = Φ N 0 (h) -     Γ N (0) . . . Γ N (N )     T R N     Γ N (0) . . . Γ N (N )     ≺ 0, (9)
hold, where Γ N (k), for all k = 0, . . . , N , are defined in Corollary 4 and

Φ N 0 (h) = He G T N (h)P H N + SN + h 2 F T N RF N , SN = diag{S, -S, 0 N n }, R N = diag{R,
F N = A A d 0 n,nN , , G N (h) = I 0 n 0 n,nN 0 nN,n 0 nN,n hI nN , H N = F T N Γ T N (0) Γ T N (1) . . . Γ T N (N -1) T .
Then the time delay system (1) is asymptotically stable for the constant delay h.

Proof : Guided by the B-L inequality (6) and the signals involved, we consider the following extra-states xN (t) defined by:

xN (t) =        x t (0) 0 -h L 0 (s)x t (s)ds . . . 0 -h L N -1 (s)x t (s)ds        , if N ≥ 1 and x0 (t) = x t (0), if N = 0.
(t) = H N ξ N (t), (10)
where

ξ N (t) =           x t (0)
x T t (-h)

1 h 0 -h L 0 (s)x t (s)ds . . . 1 h 0 -h L N -1 (s)x t (s)ds           , N ≥ 1,
and, if N = 0, ξ T 0 (t) = x T t (0) x T t (-h) . It appears that this augmented system is the interconnection of the original delay system and a LTI system defined by the states 0 -h L k (s)x t (s)ds, for k = 0, . . . , N -1. It is also worth mentioning that the only delayed term in ( 10) is x t (-h). Then, a natural choice for the LKF is

V N (x t , ẋt ) = xT N (t)P N xN (t) + t t-h x T (s)Sx(s)ds +h t t-h t θ ẋT (s)R ẋ(s)dsdθ, (11)
On a first hand, following the procedure provided in [START_REF] Gu | Stability of Time-Delay Systems[END_REF], the condition S ≻ 0 allows applying Lemma 3 to the second term of V N to give a more accurate lower bound of the functional. In order to be consistent with the definition of xN , Lemma 3 is considered with the order N -1. It thus yields

V N (x t , ẋt ) ≥ xT N (t)Θ N (h)x N (t) +h t t-h t θ ẋT (s)R ẋ(s)dsdθ.
Then the positive definiteness of V N results from the conditions S ≻ 0, R ≻ 0 and Θ N (h) ≻ 0. This also implies that there exists a sufficiently small ǫ 1 > 0, such that Θ N (h) ≻ ǫ1I 0 0 0 . It follows that V N (x t , ẋt ) ≥ ǫ 1 |x t (0)| 2 . Furthermore, there exists a sufficiently large scalar λ > 0 such that P N ≺ λdiag(I, I, 3I, 5I, . . . , (2N -1)I). It thus holds

V N (x t , ẋt ) ≤ λ|x t (0)| 2 + λ N -1 i=0 (2i + 1)Ω T i Ω i + t t-h x T (s)Sx(s)ds +h t t-h t θ ẋT (s)R ẋ(s)dsdθ.
Thanks to Lemma 3, we obtain

V N (x t , ẋt ) ≤ λ|x t (0)| 2 + t t-h x T (s)(λhI + S)x(s)ds +h t t-h t θ
ẋT (s)R ẋ(s)dsdθ, which guarantees that there exists a scalar ǫ 2 > 0, such that

V N (x t , ẋt ) ≤ ǫ 2 |x t | 2 h , for all t > h, where xt = xt ẋt . Then it holds ǫ 1 |x t (0)| 2 ≤ V N (x t , ẋt ) ≤ ǫ 2 |x t | 2 h . (12)
Consider now the derivative of V N , for all t ≥ h. We obtain

VN (x t , ẋt ) = 2x T N (t)P N ẋN (t) + x T t (0)Sx t (0) -x T t (-h)Sx t (-h) + h 2 ẋT t (0)R ẋt (0) -h 0 -h ẋT t (s)R ẋt (s)ds. ( 13
) By noting that xN (t) = G N (h)ξ N (t), ẋN (t) = H N ξ N (t) and ẋt (0) = F N ξ N (t), it yields VN (x t , ẋt ) = ξ T N (t)Φ N 0 (h)ξ N (t) -hI( ẋ, h) (14)
Finally, applying Corollary 4 to the order N and injecting the resulting inequality into ( 14) leads to VN (x t , ẋt ) ≤ ξ T N (t)Φ N (h)ξ N (t). Hence, if the LMI (9) are satisfied, there exists a scalar ǫ 3 > 0 such that Φ N (h) ≺ -ǫ3I 0 0 0 . We finally obtain

VN (x t , ẋt ) ≤ -ǫ 3 |x t (0)| 2 , ∀t ≥ h. ( 15
)
The end of the proof is taken from the proof of Theorem 1 from [START_REF] Fridman | Robust sampled-data stabilization of linear systems: An input delay approach[END_REF]. Integrating [START_REF] Jiang | Delay-dependent robust stability for uncertain linear systems with interval time-varying delay[END_REF] we have

V N (x t , ẋt ) -V N (x h , ẋh ) ≤ -ǫ 3 t h |x s (0)| 2 ds (16)
and, hence, (12) yields [START_REF] Hale | Introduction to Functional Differential Equations[END_REF] p. 168) and ẋ, defined by the right-hand side of (1), satisfies

ǫ 1 |x t (0)| 2 ≤ V N (x t , ẋt ) ≤ V N (x h , ẋh ) ≤ ǫ 2 |x h | 2 h . Since |x h | h ≤ c 1 |φ| h , c 1 > 0 (cf.
| ẋh | h ≤ c 2 |φ| h , c 2 > 0, we obtain that |x t (0)| 2 ≤ V N (x h , ẋh )/ǫ 1 ≤ c 3 |φ| 2 h , c 3 > 0.
Hence, ( 1) is stable. To prove asymptotic stability we note that, for any initial condition φ, x is uniformly continuous on [0, ∞) (since ẋ defined by the right-hand side of ( 1) is uniformly bounded). Moreover, ( 16) yields that |x t (0)| 2 is integrable on [h, ∞). Then, by Barbalat's lemma, x t (0) → 0 as t → ∞. Consequently, if the LMI of Theorem are satisfied, the delay system (1) is asymptotically stable for the constant delay h. ♦ Remark 3 Taking N = 0 in Theorem 5 allows retrieving one of the most classical delay-dependent stability conditions based on Jensen's inequality and LMI [START_REF] Gouaisbaut | A note on stability of time delay systems[END_REF]. Additionally, choosing N = 1 leads to the stability conditions from [START_REF] Seuret | Wirtinger-based integral inequality: Application to time-delay systems[END_REF] 4

.2 Remark on the choice of the LKF

A comment on the LKF, V N , and its relation with the class of functionals studied in [START_REF] Gu | A further refinement of discretized Lyapunov functional method for the stability of timedelay systems[END_REF][START_REF] Peet | Positive forms and stability of linear time-delay systems[END_REF] is highlighted here. Indeed by considering the functional [START_REF] Gu | Stability of Time-Delay Systems[END_REF] and by defining the polynomial matrix D(s) = diag(0 n , L 0 (s)I, L 1 (s)I, . . . , L N (s)I) and the matrices

P = I 0 nN,n T P N I 0 nN,n , S(s) = S Q(s) = I 0 nN,n T P N D(s), R(s, ξ) = D T (s)P N D(ξ).
Therefore, the functional V N can be rewritten as

V N (x t , ẋt ) = x T (t) P x(t) + 2x T (t) 0 -h Q(s)x t (s)ds + 0 -h 0 -h x T t (s)R(s, ξ)x t (ξ)dsdξ + 0 -h x T t (s)S(s)x t (s)ds +h 0 -h 0 θ ẋT t (s)R ẋt (s)dsdθ. ( 17
)
The three first terms of V N are similar to the ones employed in [START_REF] Peet | Positive forms and stability of linear time-delay systems[END_REF] and in [START_REF] Gu | A further refinement of discretized Lyapunov functional method for the stability of timedelay systems[END_REF]. In [START_REF] Peet | Positive forms and stability of linear time-delay systems[END_REF], the degree of freedom comes from the degree of the polynomial matrices Q(s), S(s) and R(s, ξ), denoted as D p . In [START_REF] Gu | A further refinement of discretized Lyapunov functional method for the stability of timedelay systems[END_REF], the degree of the polynomial is always 1 but the degree of freedom comes from the degree of discretization, denoted latter on as D d . A first difference with respect to these two approaches is that in our setup, the polynomial matrix S(s) is constant. The consequence is that our method requires less parameters to define the functional when increasing D p or D d . Another difference relies on the last integral quadratic term of V N which depends on ẋt . Finally, the previous theorem does not need to enter into the sum of squares framework which generally requires the use of additional decision relaxation variables when testing the stability conditions.

Delay range stability

In Theorem 5, the delay is supposed to be perfectly known. Therefore, Theorem 5 only ensures the stability of the delay system for the pointwise delay h. This subsection extends this first result by considering that the delay h is unknown but belongs to a prescribed interval [h 1 , h 2 ]. We aim therefore at providing a criterion which ensures stability for all constant delays in this pocket.

Theorem 6 For a given integer N and an uncertain constant delay h ∈ [h 1 h 2 ], assume that there exist a matrix P N ∈ S (N +1)n and S, R ∈ S + n such that

Θ N (h 2 ) ≻ 0, Φ N (h 1 ) ≺ 0, Φ N (h 2 ) ≺ 0, (18)
hold, where Θ N (h) and Φ N (h) are defined in Theorem 5.

Then the time delay system (1) is asymptotically stable for all constant delays h in the interval [h 1 , h 2 ].

Proof : Assume that the LMI (18) hold and consider any delay h in the interval [h 1 , h 2 ]. We aim at proving that the two LMI conditions of Theorem 5 are satisfied of the delay h. On the first hand, the positive definiteness of S directly implies that Θ N (h) ≻ Θ N (h 2 ), which thus guarantees the positive definiteness of Θ N (h).

On the other hand, the application of the Schur Complement to the term h 2 F T N RF N with the following decomposition (hRF N ) T R -1 (hRF N ) leads to an affine and, consequently convex LMI with respect to h. Then, this convexity property thus guarantees that

Φ N (h) h -h 1 h 2 -h 1 Φ N (h 2 ) + h 2 -h h 2 -h 1 Φ N (h 1 ).
Finally the two last LMI of ( 18) ensure that the right hand side of the previous inequality is definite negative, which, by transition, ensures that Φ N (h) ≺ 0. ♦

Hierarchy of LMI stability conditions

This section is devoted to proving that the previous stability conditions form a hierarchy of LMI conditions. This is formulated in the following theorem based on the stability conditions of Theorem 5.

Theorem 7 For any time delay system (1), define the set H N by

H N := h ∈ R + s.t. Θ N (h) ≻ 0, Φ N (h) ≺ 0, P N , S(N ) ≻ 0, R(N ) ≻ 0
Then, it holds

H N ⊂ H N +1 , ∀N ≥ 0.
Proof : Let N ∈ N. If H N is empty, the inclusion is trivial. If H N is not empty, then consider an element h ∈ H N . From the definition of H N , there exist symmetric matrices P N , S(N ) ≻ 0 and R(N ) ≻ 0 such that Θ N (h) ≻ 0 and Φ N (h) ≺ 0. Taking advantages of the construction of the LKF [START_REF] Gu | Stability of Time-Delay Systems[END_REF], we suggest the matrices

P N +1 = P N 0 0 ǫI , S(N + 1) = S(N ) = S, R(N + 1) = R(N ) = R,
where ǫ is a scalar to be chosen. Then the matrix Θ N +1 (h) can be rewritten as

Θ N +1 (h) = Θ N (h) 0 0 2N +1 h S + ǫI .
Since Θ N (h) ≻ 0 and S ≻ 0, there exists a sufficiently small ǫ such that Θ N +1 (h) ≻ 0. Moreover it holds

H N +1 = H N 0 N n,n Γ N +1 (N ) , G N +1 (h) = G N (h) 0 N n,n 0 n,N n hI , SN+1 = SN 0 N n,n 0 n,N n 0 n , F N +1 = F N 0 n .
From these expressions, the matrix Φ N +1 (h) can be expressed using the matrix Φ N (h) as follows

Φ N +1 (h) = I 0 Γ N +1 (N + 1) T Φ N (h) 0 0 -(2N + 3)R × I 0 Γ N +1 (N + 1) + ǫHe hΓ T N +1 (N ) 0 I .
Since Φ N (h) ≺ 0, R ≻ 0 and by noting that the matrix

I 0 Γ N +1 (N + 1)
is non singular, the first term of the previous expression is negative definite. It implies that there exists a sufficiently small ǫ for which Φ N +1 (h) ≺ 0 and Θ N +1 (h) ≻ 0. Therefore, h belongs to H N +1 . Finally, since h is any element of H N , it implies that H N ⊂ H N +1 . ♦ Theorem 7 proves that, the stability conditions provided in Theorem 5 at the order N + 1 delivers, at least, the same result the same condition taken at the order N . Moreover, since Theorem 5 only provides sufficient stability condition, the set H N , for a given N ∈ N represents an inner approximation of the stability pockets. However, the previous theorem does not prove that the conditions of Theorem 5 will converge to the analytical bounds of the delay. A similar development shows that Theorem 6 also forms a hierarchy of stability conditions.

Examples

The purpose of the following section is to illustrate on academic and non trivial examples how the inequalities given in Section 3 lead to a relevant reduction of conservatism in the stability condition.

Example 1

Consider the linear time-delay system (1) with the matrices

A = -2 0 0 -0.9 , A d = -1 0 -1 -1 . ( 19
)
This system is a well-known delay dependent stable system, that is the delay free system is stable and the maximum allowable delay h max = 6.1725 can be easily computed by delay sweeping techniques. The results are reported in Table 1. Many recent papers give the same result since they are intrinsically based on the same LKF and use the same bounding cross terms technique i.e. Jensen inequality. Some papers [START_REF] Sun | Improved delay-range-dependent stability criteria for linear systems with time-varying delays[END_REF], [START_REF] Sun | Delay-dependent stability and stabilization of neutral time-delay systems[END_REF] which use an augmented Lyapunov, based on the addition on a triple Integral term on the LKF can go further but with a numerically increasing burden, compared to our proposal. The partitioning approach proposed by [START_REF] Han | A discrete delay decomposition approach to stability of linear retarded and neutral systems[END_REF] based on the discrete delay decomposition gives an upperbound which tends to the analytical value even if the numerical complexity remains important. The robust approach [START_REF] Kao | Stability analysis of systems with uncertain time-varying delays[END_REF] gives a very good upper-bound with a similar computational complexity than our present result. The discretized LKF proposed by [START_REF] Gu | A further refinement of discretized Lyapunov functional method for the stability of timedelay systems[END_REF] as well as the sum of square optimization scheme developed by Peet et al [START_REF] Peet | Positive forms and stability of linear time-delay systems[END_REF] give a delay upperbound very closed to the maximum allowable delay with an increasing numerical complexity.

Example 2

This example is provided to illustrate Theorem 7. Note that this system has not been studied in the literature of time delay system using the Lyapunov-Krasovskii Theorem and LMI conditions. It is extracted from the dynamics of machining chatter [START_REF] Zhang | Stability of time-delay systems: Equivalence between Lyapunov and scaled small-gain conditions[END_REF][START_REF] Sipahi | Stability and stabilization of systems with time delay[END_REF] and is given by

ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t), with A =      0 0 1 0 0 0 0 1 -10 10 0 0 5 -15 0 -0.25      , B =      0 0 1 0      , C =      1 0 0 0      T .
A delayed static output feedback controller is proposed:

u(t) = -Ky(t) + Ky(t -h),
where K is the gain of the controller and h is an unknown constant delay. The resulting dynamics is thus modeled We have also proved that the set of stability conditions forms a hierarchy of LMI indexed by the polynomial degree N , in the sense that increasing N reduces the conservatism of the proposed method. Future works will include the study of the asymptotic necessity of this approach and also the case of time-varying delays.

.

  For any function f ∈ C, the norm |f | h refers to sup θ∈[-h, 0] |f (θ)|. diag(A, B) stands for the diagonal matrix A 0 0 B

  3R, . . . , (2N + 1)R}, S N = diag{S, 3S, . . . , (2N + 1)S},

Fig. 1 . 7 Conclusions

 17 Fig. 1. Stability region in the plan (K, h), obtained using Theorem 5 for N = 0, . . . , 7.by a time-delay system:ẋ(t) = A 0 x(t) + A 1 x(t -h),with A 0 = A -BKC and A 1 = BKC. The results are summarized in Figure1which shows the stability regions in the (K, h) plane. The blue region represents the instability region which have been calculated using a griding over K along with the allmargin function of the Control Toolbox of Matlab c . Then, Theorem 5 provides inner approximations of the stability region delimited by colored curves. The curve N = 0 perfectly detects the independent of the delay stability region (for K ≤ 0.3) as well as a first delay dependent stability pocket. It corresponds to the maximal allowable delay when Jensen's Lemma is used when establishing the stability criterion. Taking N = 1 allows retrieving the same results as in[START_REF] Seuret | Wirtinger-based integral inequality: Application to time-delay systems[END_REF] which uses Wirtinger's Lemma. Clearly, increasing N (N = 0, 1, . . . , 7) allows reducing the pessimism

  The augmented vector xN is composed by the instantaneous state x t (0) and the projections of the state function x t to the N first Legendre polynomials. Following the proof of Corollary 4 and equation (8), an integration by parts allows expressing the time derivative of xN as follows ẋN

Table 1

 1 Results for Example[START_REF] Kolmanovskii | Stability of some linear systems with delay[END_REF] for constant delay h. The degree of discretization D d and the degree of the polynomial Dp are defined in subsection 4.2.

	Theorems	hmax	number of variables	Theorems	hmax	number of variables
	[9, 14, 28, 31]	4.472	1.5n 2 + 1.5n	[10], D d = 1	6.053	7.5n 2 + 3.5n
	[16]	6.1107	1.5n 2 + 9n + 9	[10], D d = 2	6.165	10.5n 2 + 4.5n
	[2]	5.120	7n 2 + 4n	[10], D d = 3	6.171	14.5n 2 + 4.5n
	[31]	5.02	18n 2 + 18n	[10], D d = 4	6.171	20.5n 2 + 5.5n
	[18]	4.97	69n 2 + 5n	[25], Dp = 1	5.19	7n 2 + 3n
	[1]	5.120	6.5n 2 + 3.5n	[25], Dp = 2	5.90	12.5n 2 + 4.5n
	[30]	5.30	8.5n 2 + 3.5n	[25], Dp = 3	6.10	21n 2 + 6n
	[13],D d = 2	5.71	4n 2 + 2n	Th.5, N = 0	4.472	1.5n 2 + 1.5n
	[13],D d = 3	5.96	6.5n 2 + 2.5n	Th.5, N = 1	6.059	3n 2 + 2n
	[13],D d = 4	6.05	10n 2 + 3n	Th.5, N = 2	6.168	5.5n 2 + 2.5n
	[26]	5.901	3n 2 + 2n	Th.5, N = 3 6.1725	9n 2 + 3n