Benoît

Cette thèse a été réalisée en collaboration avec le département Drones Navals et Guerre des Mines de la Direction Générale de l'Armement -Techniques Navales sur le site de Brest et l'Ensta Bretagne. Je souhaite remercier l'ensemble de ces acteurs pour m'avoir fourni la possibilité de réaliser ces travaux de recherche. Mes premières pensées vont bien évidemment aux personnes qui m'ont encadré, conseillé et soutenu.

En premier lieu, je tiens à remercier Luc Jaulin pour avoir été un excellent directeur de thèse. Ces travaux n'auraient pu être aussi aboutis sans son aide précieuse, sa grande disponibilité et nos nombreuses discussions ! Je remercie également les membres de mon jury pour le temps qu'ils m'ont consacré ainsi que l'attention et l'intérêt qu'ils ont montrés lors de la relecture de mon manuscrit et ma présentation : leurs corrections, questions et remarques pertinentes et constructives m'ont permis d'améliorer ce manuscrit et de valoriser mon travail.

En particulier, je tiens à remercier mes rapporteurs : Sylvie Galichet, pour avoir apporté un point de vue issu de la communauté de la logique flou et Warwik Tucker pour son approche théorique, Simon Lacroix d'avoir présidé ce jury et Vincent Creuze, Fabrice Le Bars, Nicolas Drogi, Alain Bertholom d'avoir endossé les rôles d'examinateurs.

Je remercie également mes collègues de DGA Tn, qui m'ont apporté une vision pragmatique sur les problématiques de positionnement sous-marin et en particulier toute l'équipe du groupe Drone navale et Guerre des mines de Brest. Pour la partie théorique, je remercie vivement les différents personnels, doctorants, de l'Ensta Bretagne qui m'ont aidé durant ces 3 années. En particulier, je tiens à remercier Simon R., Thomas L. M., qui m'ont accueilli dans leur bureau lors de mes passages à l'Ensta Bretagne.

Je tiens à remercier tout particulièrement mes relecteurs Quentin B, Simon R. pour leur travail minutieux.

A titre personnel, je tiens à remercier l'ensemble de ma famille et en particulier mes parents et beaux-parents pour m'avoir encouragé durant mes nombreuses heures de travail. Et un grand merci pour le pot de thèse bien sûr ! Enfin, et par-dessus tout, je souhaite remercier ma femme Charlotte pour sa tendresse et son soutien indéfectible durant ces 3 années.

List of Acronyms

Introduction

This thesis focuses on new methods, based on Interval Analysis, able to deal with constraints that involve uncertain sets. Many engineering problems, such as robot localization, global motion path-planning or control, can be modeled by constraint networks. A Constraint Network (CN) is composed of a set of variables that belong to some domain and which must satisfy a given number of elementary facts and rules: so-called constraints. In classical formulations [Mackworth, 1977], variables are vectors of real numbers and domains are subsets of R n .

In this manuscript, we consider a larger class of problems that involve subsets of R n as variables too. As the main application, we deal with the Simultaneous Localisation And Mapping (SLAM) problem in unstructured environments, i.e that cannot be described by geometrical features. Sensor readings will not be depicted by a model that can be described by parameters of R but directly by shapes or subsets of R n .

The uncertainties will then be handled, in a set-membership approach, by an interval of sets called Thick Set.

Motivation and background

Autonomous vehicles are used in support of humans to perform dull, dirty and dangerous tasks. The more dangerous and challenging the environment, the more autonomous they must be.

To be fully autonomous, vehicles must be able to perform long-term and long-range missions with little or no human intervention. They must position themselves with meaningful and reliable bounds on their positioning uncertainty that can be used for decision-making purposes. For instance, in the context of seabed mapping using an Autonomous Underwater Vehicle (AUV), the positioning uncertainty must be monitored, in order to assess that the mission requirements are successfully achieved. This can concern the assessment of the covered area by the vehicle before returning to the mother-ship, or that the positioning errors on payload data are acceptable with respect to the mission requirements. If not, the initial plan must be modified and specific actions need to be taken.

To localize itself, a modern autonomous vehicle relies on an Inertial Navigation System (INS) that merges proprioceptive measurements (acceleration, velocity, attitude, depth, . . . ) in order to continuously determine, via dead reckoning methods, the position, orientation and velocity of the vehicle. INS provides a good short-term accuracy but accumulates errors that grow unbounded with time and distance traveled. Consequently, this drift needs to be fixed with, for instance, data from Global Navigation Satellites Systems (GNSS). However, even if these systems can provide a centimeter position accuracy, they are not always available and remain sensible to jamming and spoofing. Another solution is to use physical features of the environment observed by embedded sensors, such as camera, LIDAR or sonar, to provide an onboard, real-time estimated, positioning solution.

Mapping

Environments, structured or not, are composed of features which are gathered in a map. In this thesis, structured environments are taken to mean terrains that can be described by geometric structures such as lines, segments, landmarks, . . . . For instance, man-made environments are generally structured. On the contrary, unstructured environments are composed of shapes that cannot be parameterized. This is the case of underwater environments since there are few natural occurring geometrical landmarks.

It exists several methods to depict the environment, but according to [Thrun, 1998] two main representations are commonly considered: topological and metric maps.

A topological map refers to a map that has been simplified so that only vital information remains and unnecessary details have been removed. Only the neighborhood relation between places is maintained. The tube map of the Paris underground is an example of such representation where only the connections between stations are depicted without considering scales and distances.

On the contrary, metric maps assign coordinates to objects of the environment that allows computing distances or angles between them. This representation can be classified into two categories [Thrun, 2002], which are referred to as feature-based map and location-based map.

Feature-based maps

Feature-based maps are composed of a set of features together with their Cartesian location. The map has a parametric structure where the features are points, segments, corners, or any other parametric shapes [START_REF] Castellanos | Multisensor fusion for simultaneous localization and map building[END_REF]. In the case of the SLAM problem, feature parameters are included among state variables [START_REF] Dissanayake | A solution to the simultaneous localization and map building (SLAM) problem[END_REF], Montemerlo et al., 2003] and are estimated concurrently with the pose of the vehicle. Probabilistic methods such as Kalman filtering, Bayesian estimation or particle filtering [Thrun, 2002] can be used. Interval-based methods have also been proved to be efficient to solve the feature-based SLAM [Le [START_REF] Bars | Interval SLAM for underwater robots; A new experiment[END_REF], Jaulin, 2016].

However, geometrical features are spatially localized and do not include information on their environment. A particular attention must be paid when solving data association ambiguities between measurements and features [START_REF] Neira | Data association in stochastic mapping using the joint compatibility test[END_REF]. Loop closure methods [START_REF] Ho | [END_REF]Newman, 2006, Aubry et al., 2013], that aim at finding places that have been previously visited, can provide a prior data association. This critical point made the localization of featured map a challenging problem in complex and unstructured environments like these faced in outdoor or underwater domains.

Location-based maps

On the contrary, location-based maps offer a label to any location in the world. A classical location-based map representation is known as occupancy map [Elfes, 1987] (also called pose-based map). They assign to each point of the world an occupancy value (a Boolean number or a probability of occupancy) that specifies whether or not a pose is occupied by an obstacle. For computational reasons, occupancy maps are most of the time based on a regular grid. Each cell of the grid stores the probability, between 0 and 1, to be occupied by an obstacle. The main difficulty with a grid-based representation is to correctly handle both the sensor and the pose uncertainties when updating the grid. This is especially the case when the map is only perceived through sensors, such as sonars [Thrun, 2003] or cameras, that partially observe the environment. To handle this difficulty properly, we introduce the notion of shape.

Shape.

A shape is a subset of R q , where q is smaller than the dimension of the workspace. In this thesis, a shape will always be a subset of dimension 2. A shape corresponds to a part of the environment that is captured by some sensors of the robot. As shapes may be uncertain, we will also introduce shape intervals (a particular case of thick sets) to represent this uncertainty in a set-membership framework. A shape separates the environment into two complementary parts, associated with physical features that can be measured by vehicle's sensors.

For examples, these features can be:

• the part of the space that can be proved to be free of obstacle. For instance, from range measurements acquired by laser or sonar sensors;

• contour map from Digital Elevation Map (DEM) that can be used to partition the space based on a given elevation value;

• segmented images based on texturing information. For instance, from data collected by a side scan sonar, which provides an acoustic image of the seabed, it is possible to partition the space based on homogeneity properties [START_REF] Picard | Potential of the intrinsic dimensionality for characterizing the seabed in the ATR context[END_REF] or a seabed classification [START_REF] Leblond | Use of classification and segmentation of sidescan sonar images for long term registration[END_REF].

Given a particular feature, labeled as a class A, the notion of uncertain shape, called shape interval, allows us to classify the space into three categories: [START_REF]Density probability function of random variable x and[END_REF] parts that can be classified for sure as A, (2) those that can be classified as not A and ( 3) parts where nothing can be said, called the penumbra. From the raw image, Figure 1.1a, a segmentation algorithm separates the mud (darker area) and non-mud area. The boundaries, between these two classes, cannot be detected accurately and are assumed to be enclosed in the penumbra. The area outside the image, which has not been observed, also belongs to the penumbra. The second part of this manuscript deals with the interval shape (or thick set) representation. The raw image 1.1a is segmented and interpreted as an uncertain shape 1.1b.

Shape based localization and mapping

From the previous set-based map representation, the localization problem needs to be defined. From a theoretical point of view, a localization problem can be described by the following set of equations.

    
ẋ(t) = f (x(t), u(t)) (evolution equation) Z(t) = g(x(t), M) (observation equation) x(0) ∈ X 0 (initial state) (1.1) where t ∈ R is the time, x ∈ R n is the state vector of the vehicle, u ∈ R m is the input vector (in general associated with proprioceptive sensors), f : R n × R m → R n is the evolution function. The set M is the map which is a subset of R q , where q is the dimension of the workspace (two or three in practice). The map M corresponds to a simplified representation of the environment. The extension to maps composed of several classes is considered as trivial and is not treated in our formalism, for clarity. The set Z(t) ⊂ R r , where r ≤ q, is a shape that represents a view of the map M, is collected by some exteroceptive sensors (for instance a camera, a sonar, a telemeter) and expressed in the robot frame. In practice, due to the uncertainties, the measurement process does not provide us Z (t) directly, but an interval shape [Z] (t) containing Z (t). The fact that the measurement corresponds to a shape instead of a vector is not classical in the SLAM community and is related to the location-based SLAM formalism. The treatment of the shape-SLAM is the main goal of this thesis. The observation function g : R n × R q → R r links the measured shape Z(t) to the map M. The initial state at t = 0 is known to belong to an initial domain X 0 ⊂ R n . The functions f and g may be non-linear and uncertain.

With respect to the unknowns of the problem, three cases can be considered: the initial localization, the position tracking, and the Simultaneous Localization And Mapping problem.

The initial localization problem corresponds to the situation where only the map M is available and the initial state must be recovered. In the probabilistic context, in the specific case where the shape Z(t) is a vector of R q or geometric features, this problem has been successfully addressed with Particle Filter (PF) methods such as Monte Carlo Localization (MCL) [START_REF] Dellaert | Monte Carlo localization for mobile robots[END_REF] or Rao-blackwellised Particle Filter (RBPF). However, proving the convergence and robustness of these algorithms remains a challenging task. On the contrary, interval methods [Jaulin et al., 2001a, Clérentin et al., 2008, Gning, 2006, Langerwisch, 2014] have been proved to be able to deal efficiently with this situation. The thesis [START_REF] Guyonneau | A visibility information for multi-robot localization[END_REF] provides an overview of set-membership methods applied to global localization.

The pose tracking problem arises when both the map M and the initial condition X 0 are known. It consists of keeping track of the vehicle's current position by matching sensor readings with the map. It can be regarded as a special case of the initial localization problem where the search space is restricted. Probabilistic methods, such as the Kalman filter [Kalman, 1960] in the linear case, or the particle filter [START_REF] Jensfelt | Pose tracking using laser scanning and minimalistic environmental models[END_REF] in the non-linear case, have been used for a long time. However, due to changes in the environment or sensor reading uncertainties, the tracked solution can be lost. These algorithms may diverge and need to be initialized with a new starting position, that needs to be found.

Finally, the situation where M needs to be estimated corresponds to the SLAM problem [Leonard et al., 1992]. For an autonomous robot, moving in an unknown environment, the SLAM problem is to build a map while simultaneously using it to compute its location. Note that the map is built relatively to the initial state which, if undetermined, is fixed arbitrarily. A significant amount of work has been completed around this topic, still subject to further research [START_REF] Leonard | Mapping Partially Observable Features from Multiple Uncertain Vantage Points[END_REF], Thrun, 2002].

Underwater SLAM [START_REF] Hidalgo | Review of underwater SLAM techniques[END_REF] and terrain based navigation [START_REF] Melo | Survey on advances on terrain based navigation for autonomous underwater vehicles[END_REF] remain challenging topics for autonomous underwater vehicles in long-term operations due to the limitations of subsea localization sensors and real-time payload exploitation. For these reasons, set-membership methods will be used since shape measurements are difficult to be handled by probabilistic methods.

A vector of R can be easily modeled by a random variable, but it is not the case for a subset of R n . Some formulations of the SLAM problem have been proposed using Random Finite Sets [START_REF] Molchanov ; Molchanov | Theory of Random Sets[END_REF], Mullane et al., 2011] but remains limited to punctual features.

Probabilistic and set-membership approaches

In this thesis, the set-membership approach has been chosen to deal with localization issues. The following paragraph highlights the main difference between probabilistic methods and set-membership ones.

Probabilistic methods are largely used in mathematics and engineering and has been proved to be efficient. A random vector x of R n is described (or encoded) using a density probability function: π : R n -→ R + , such as R n π(x)dx = 1. The variable is contained in a set X = {x ∈ R n | π(x) = 0}, called support, which also provides a description of how this variable is spread over X.

With a set membership approach, the variable is only represented by a set such as x ∈ X. The domain X can enclose inconsistent values, but all consistent values for x should be inside of X. Consequently, the set-membership representation is poorer than the probabilistic approach but needs little knowledge about the statistical properties of variables and allows to handle a larger class of problems.

The following example explains the difference between the two approaches.

Example 1.1. Let us take three random variables x, y, z of R linked with the constraint z = x + y. From statistical properties on x and y, we aim at deducing those of z with probabilistic and set membership methods.

Probabilistic approach. Assume that x and y are uniform such as π(x) = 1/5 over [0,5] and π(y) = 1/2 over [START_REF]Density probability function of random variable x and[END_REF]3]. Computing z requires the knowledge of the joint probability density for x and y. Assume for instance that x and y are independent, so the density of the couple is uniform with π(x, y) = 1/10 over the box [0,5] × [START_REF]Density probability function of random variable x and[END_REF]3]. The Probability Density Function (PDF) for z is given by: which gives us: [START_REF]Density probability function of random variable x and[END_REF]3], (1.4) 3,6], (1.5) π(z) = 1 10 (-z + 8) if z ∈ [6,8], (1.6) π(z) = 0 if z ≥ 8.

π(z) = R π x (x) • π y (z -x)dx (1.2)
π(z) = 0 if z ≤ 1, (1.3) π(z) = 1 10 (z -1) if z ∈ [
π(z) = 1 5 if z ∈ [
(1.7)

The function π(z) is depicted on Figure 1.2. Note that the integral calculus is generally much harder even in the particular cases where x and y are Gaussian or uniform.

Set-membership approach. The uncertainty associated to x et y is given by intervals: [x] = [0, 5] and [y] = [START_REF]Density probability function of random variable x and[END_REF]3]. The random variable z = x + y is contained in the interval [z] = [0 + 1, 5 + 3] = [START_REF]Density probability function of random variable x and[END_REF]8]. So the set membership approach leads to an easier computation with less hypotheses.

Note that representing a random vector x by a set X does not mean that its density function π(x) is uniform over X. It only supposes that x is in X or, in other words, that the support of π(x) is included in X. In the Example 1.1, the random variable z is described by the interval [START_REF]Density probability function of random variable x and[END_REF]8] but is not uniform on [START_REF]Density probability function of random variable x and[END_REF]8].

As illustrated by the Example 1.1, the computation of π(z) is possible because the PDF has been assumed to be uniform and with a linear dependency between x, y and z. When the relation between variables is non-linear, Equation 1.2 becomes hard to solve in the general case. Moreover, the independence between variables must always be assumed mainly because the join law between variables is rarely available or because it does not correspond to the reality. It frequently happens when dealing with data from real sensors. However, when probabilistic properties are well-known, for instance when the problem is linear Gaussian, probabilistic methods are efficient and outperform largely set-membership methods.

In other cases, strong hypothesis must be done over the probabilistic properties of the problem and no proof can be done on the results. Filters, such as extended Kalman filters or particle filters, can converge to false or inconsistent values without raising any error. For this last point, interval methods allow to correctly, with realistic hypothesis, handle uncertainties.

In this manuscript, probabilistic methods will not be presented or compared.

With respect to the bounded error hypothesis, set-membership methods provide reliable and guaranteed results. They can be used concurrently with probabilistic approaches to find the best estimate that lays into a guaranteed solution set. In other words, the hard part of strongly non-linear problems, such as data association or initial localization, can be made using interval analysis. Once done, probabilistic approaches can be used to efficiently find the best solution that lies in a guaranteed domain of validity.

Contributions and Outlines

The main contribution of this thesis is to consider that both the map and observations can be described by shapes, associated with physical features of the environment. We aim at tackling both initial localization and SLAM problems through a formalism that is not commonly used in the literature. For that purpose, the notion of interval shape is introduced, with its associated arithmetics and the definition of the uncertain set inversion problem. A specific contractor is developed to solve the so-called shape registration and carving problem.

Another contribution concerns a way to build minimal contractors and separators consistent with equations of the form y = f (x). It allows us to introduce minimal operators associated with the change from polar to cartesian coordinates when range and bearing measurements are collected.

In this manuscript, we chose to only focus on the shape-based SLAM. For this reason, two main contributions, not directly related to the SLAM problem, have been placed in Appendix A and B. The first one is dedicated to a specific algorithm that solves the uncertain transformation of an interval shape when the uncertain function can be enclosed by an interval of functions. The contribution is dedicated to the characterization of the area explored by a robot when the dimension of the visible space at time t is smaller than that of the workspace. Both of them have required a huge amount of work and helped us to mature the concepts needed to properly formalise and solve the Shape SLAM problem.

This thesis is decomposed into two parts. The first one, from Chapters 2 to 4, is dedicated to the classical interval analysis used to solve problems involving variables that are vectors of R n . The second part, composed with Chapters 5 and 6, focuses on constraints that involve shapes as variables.

In Chapter 2, basic notions of interval analysis and contractors are introduced. A special focus on the separator operator is done in this chapter. This operator, introduced in [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF], allows simple and efficient set manipulation and characterization.

Contributions and Outlines

In Chapter 3, we introduce a new theorem that can be used to build minimal contractors consistent with equations, and another new theorem to derive an optimal separator from a minimal contractor. As an application, we focus on the polar constraint associated with the change between Cartesian coordinates and Polar coordinates. We illustrate our method on the localization problem of an actual underwater robot where both range and goniometric measurements of landmarks are collected. The material of this chapter is published in [Desrochers and Jaulin, 2016a].

In Chapter 4, a new type of constraint is introduced, that links shapes via a transformation depending on some parameters. Based on the Polar separator, introduced in Chapter 3 that allows handling sectoral sensor, a first initial localization problem on an unstructured map is presented. This chapter makes the transition from classical constraint networks to shape-valued constraint networks where variables become shapes. This work is published in [Desrochers and Jaulin, 2017c].

Chapter 5 deals with intervals of shapes, which is a particular type of Thick Sets. Shape intervals provide an interval based representation of uncertain shapes. An uncertain shape Y is bracketed between an inner shape Y -and an outer shape Y + . After defining operations that can be applied to shape intervals, tools used to deal with constraint networks involving shapes as variables are introduced. A special focus is done on the set inversion problem involving an uncertain function that depends on parameters. The shape registration constraint, presented in Chapter 4 is then extended to the uncertain case.

Chapter 6 deals with the SLAM problem where the observations are not based on scalar values or equations, but are subsets of R n . Two applications are considered.

In the first one, the initial localization problem presented in Chapter 4 is extended to the context of the SLAM, where the map needs also to be estimated. In the second application, measurements taken by an imaging sensor are depicted by an uncertain shape. A SLAM method, based on this set-membership representation is then illustrated. This chapter summarizes the targeting application of this thesis.

Chapter 7 concludes this thesis.

Appendix A deals with the set inversion problem X = f -1 (Y) in the case where the function f : R n → R m and the set Y are both uncertain. The uncertainty is treated under the form of intervals. More precisely, for all x, f (x) is inside an interval of functions [f ](x) and the uncertain set Y is depicted by a thick set. The introduction of new tools such as thick intervals and thick boxes will allow us to propose an efficient algorithm that handles the uncertainty of sets in an elegant and efficient manner. Some elementary test-cases that cannot be handled easily and properly by existing methods show the efficiency of the approach. This work is published in [START_REF] Desrochers | Thick set inversion[END_REF].

Appendix B deals with the guaranteed characterization of the part of the space that has been explored by a robot. The main difficulty of the problem is to take into account the uncertainty associated with the trajectory and the fact that the dimension of the visible space at time t may be smaller than that of the workspace. An example involving an experiment made with an actual underwater robot is presented in order to illustrate the efficiency of the approach (published in [Desrochers and Jaulin, 2017b]).

PyIbex: an Interval Library

To support and capitalize the developments done during this thesis, a Python library called pPyIbex has been developed. The source code is available on GitHub at https://github.com/benEnsta/pyIbex.

Based on the contractor programming library ibex-lib, it aims at providing a python interface to basic functionality, such as interval, contractor, separator, and adds extra-modules dedicated to tools introduced in this thesis.

PyIbex is presented as a multi-platform Python module. It was used by IAMOOC, a MOOC on Interval Analysis with application to parameter estimation and robot localization.

The documentation and examples are given in http://benensta.github.io/pyibex

The figures of the manuscript are generated using VIBes [START_REF] Drevelle | VIBes: A Visualizer for Intervals and Boxes[END_REF].

Part I Interval Analysis

Chapter 2

Interval Analysis Contents

Introduction

Our choice is to address robotic problems with set-membership methods. Basically, the solution of the problem, such as the estimating of the position of the robot, is characterized by a set defined by non-linear constraints. This set needs to be represented by an over approximation which can be manipulated by the computers. As illustrated in Figure 2.1, in a set membership framework, several approximations can be found in the literature such as zonotopes or polyhedral enclosures [Combastel, , Walter and Piet-Lahanier, 1989], ellipsoids [START_REF] Rokityanski | Application of ellipsoidal estimation to satellite control[END_REF], intervals or subpavings [Jaulin and Walter, 1993c]. In this thesis, we will focus on the last representation which is related to interval analysis. Interval analysis has become over the past few years a strong alternative to traditional probabilistic approaches [START_REF] Thrun | The Graph-SLAM Algorithm With Applications to Large-Scale Mapping of Urban Structures[END_REF] to solve complex non-linear systems of equations.

Interval analysis makes it possible to characterize a set of R n defined by constraints. For instance, the set may correspond to:

• the set of all parameters that are consistent with some interval measurements [START_REF] Lévêque | Vehicle Localization from Inaccurate Telemetric Data: A Set-Inversion Approach[END_REF], Gning and Bonnifait, 2006, Kreinovich et al., 1997, Gning et al., 2013, Chabert and Jaulin, 2009a],

• to the set of all configuration vectors such that a robot does not meet any obstacles [START_REF] Porta | A space decomposition method for path planning of loop linkages[END_REF], Jaulin, 2001], to the set of all parameter vectors of a controller such that the closed loop system is stable [START_REF] Wan | A numerical approach to design control invariant sets for constrained nonlinear discrete-time systems with guaranteed optimality[END_REF], Didrit et al., 1995],

• to set all calibration parameters [Daney et al., 2006, Ramdani and[START_REF] Ramdani | [END_REF] • to attractors of dynamical systems [Tucker, 1999].

More formally, a mathematical problem can be formulated by means of a Constraint Network (CN) where a set of variables, that belong to some domains, must satisfy simultaneously a given number of elementary facts and rules: so-called constraints [Mackworth, 1977].

A CN is defined as a triple X, D, L where

X = {x 1 , . . . , x n } is a set of variables, D = {X 1 , . . . , X n }
is a set of the respective domains of values, and

L = {L 1 , . . . , L m } is a set of constraints.
The problem to be considered here is to find the smallest sub-domains of the X i for variables x i consistent with all constraints. Note that, in this part, variables are real numbers or vector of R n and domains are intervals, boxes, or sets. Later, the class of variables to be handled will be extended: they could correspond to shapes and their domains will be shape intervals.

Remark 2.1. In the literature, there exists a confusion between Constraints Networks (CN) and Constraint Satisfaction Problem, or CSP for short. Stress that a CSP is composed of a CN (which is the model) and a question. Solving a CSP is done by providing a solution of the corresponding CN, which is not exactly what we want to do here. Instead, we want to find an enclosure of all solutions of the CN.

The purpose of this chapter is to familiarize the reader with interval analysis concepts. In the first section, elementary notions on set theory are introduced. Then, in Section 2.3, basics of interval arithmetic are provided. In order to get an inner and an outer approximation of sets defined with constraints, contractors and separators are presented in Section 2.4. With respect to the applications of this manuscript, common separators associated to images, subpavings or polygons are introduced in Section 2.5. We will also give two applications illustrating how these operators can be combined to solve robotics problems.

Set Theory

The fundamental notion that is used at the heart of this thesis is the notion of set. A set is a well-defined collection of distinct objects. The objects that make up a set (also known as the set's elements or members) can be anything: numbers, people, letters of the alphabet, graphs, other sets, and so on. The set theory, that aims at studying sets, was developed at the end of the 19th century by George Cantor and Richard Dedekind [Cantor, 1847]. In this manuscript, sets are denoted in blackboard bold font (mathbb). The set of real numbers is denoted R.

A B

(a) A ∩ B. A B (b) A ∪ B. A B (c) A \ B.

Operations

Set theory begins with a fundamental binary relation that defines whether or not an element a belongs to a set A, denoted by a ∈ A. From that relation, we are mainly considering sets composed of real numbers defined from predicates or constraints. For instance, the set of positive real numbers, commonly denoted by R + is defined by:

R + = {x ∈ R | x ≥ 0} . (2.1)
Two special sets must be considered: the empty set denoted by ∅ that contains no element and the universe, denoted Ω that contains all elements. Basic operations on sets used in this manuscript are now introduced.

The complement of a set A, denoted A is:

A = {x ∈ Ω | x / ∈ A} . (2.2)
Given two sets A and B, the inclusion between sets is defined by: .3) The equality between two sets is then defined by: .4) where ∧ stands for the logical and operator. The logical or is denoted ∨. Complex sets are defined from binary operations such as: In this thesis, a special intersection called the relaxed-intersection [Jaulin, 2009] will be used. It allows solving inconsistent CN by relaxing a small number of constraints. In the context of parameter estimation, it makes it possible to be robust to a given number of outliers (erroneous measurements).

A ⊂ B ⇔ ∀x ∈ A, x ∈ B. ( 2
A = B ⇔ (A ⊂ B ∧ B ⊂ A) . ( 2
Intersection A ∩ B = {x | x ∈ A ∧ x ∈ B}. Union A ∪ B = {x | x ∈ A ∨ x ∈ B}. Difference A \ B = {x | x ∈ A ∧ x / ∈ B} = A ∩ B. Cartesian product A × B = {(x, y) | x ∈ A ∧ y ∈ B}.
Consider N sets X 1 , . . . , X N of R 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 2.3 illustrates this notion for N = 6 and q = 2, 3, 4. The relaxed intersection can be formulated in term of unions of intersections. For three sets X 1 , X 2 , X 3 , the 1-relaxed intersection is defined by:

{1} X i = (X 1 ∩ X 2 ) ∪ (X 1 ∩ X 3 ) ∪ (X 2 ∩ X 3 ).
(2.5) Figure 2.3: q-relaxed intersection of 6 sets for q=2 (red), q=3 (green), q=4 (blue) , q=5 (yellow)

Set image

Given a function f : A → B and A 1 ⊂ A, the direct image by f of A 1 is defined by:

f (A 1 ) = {f (x) ∈ B | x ∈ A 1 } (2.6)
The reciprocal image by f of B 1 ⊂ B is then:

f -1 (B 1 ) = {x ∈ A | f (x) ∈ B 1 } (2.7)
A canonical problem for which interval analysis has been proved to be efficient is the Set Inversion Problem. Let f be a function from R n to R m (possibly non-linear)

and Y be a subset of R m . The Set Inversion problem aims at characterizing the reciprocal image of Y by f :

X = {x ∈ R n |f (x) ∈ Y} = f -1 (Y) (2.8)

Link with interval analysis

Elementary sets can be defined by constraints and combined in order to build more complex ones. The interval analysis framework provides algorithms and operators, such as contractors and separators, that makes it possible to compute, combine and represent sets of R n that are defined by constraints. The following sections introduce these tools.

Interval Analysis

Interval Analysis was first developed to quantify the error on numerical computations [Moore, 1966]. In computers, real numbers are usually represented by floats with limited significant digits. This limitation leads to a small error that can drastically propagate and increase along successive operations. For instance, the decimal number "0.1" cannot be represented exactly in 32-bits floating-point pricision (8 bits exponent and 23 bits mantissa) in a computer. But it can be enclosed between two numbers such as:

0.0999999940395355224609375 < 0.1 < 0.100000001490116119384765625
Interval methods make it possible to enclose numerical and physical uncertainty, providing a rigorous theoretical framework. In this thesis, we will not focus on this low level uncertainty, which is absorbed by interval libraries such as filib++ [START_REF] Lerch | FILIB++, a Fast Interval Library Supporting Containment Computations[END_REF], gaol [Goualard, 2008],. . . , but we will use interval analysis to manipulate the uncertainty on input parameters. For instance, a range sensor returns a distance, with an error, that will be represented by an interval. This section summarizes the basics of interval analysis. For a more detailed introduction, see [START_REF] Jaulin | Applied Interval Analysis[END_REF].

Interval arithmetics

Intervals

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

[x] = x -, x + = x ∈ R, x -≤ x ≤ x + .
(2.9)

[x] [y] = [{x y|x ∈ [x], y ∈ [y]}] (2.10)
where [.] denotes the convex hull. The convex hull of a set A ⊂ R is the smallest interval that contains A.

Example 2.1. We have

[-1, 5] ∩ [-3, 2] = [-1, 2] (2.11) [-3, 3] ∪ [4, 10] = [-3, 10] (2.12) [-3, 3] -[4, 10] = [-13, -1] (2.13) w([-2, 9]) = 11
(2.14)

Boxes

An interval vector [x] of R n , or a box for short, is a Cartesian product of n intervals.

The set of all boxes of R n is denoted by IR n .

[ .15) In this thesis, the width of a box [x] ∈ IR n is defined by: 2.16) This definition differs from the classical one which considers the length of the largest interval along all dimensions:

x] = [x 1 ] × [x 2 ] × • • • × [x n ]. ( 2
w([x]) = n i=1 w([x i ]) 2 (
w ∞ ([x]) = max 1≤i≤n w([x i ]).
(2.17) The first definition is preferred because, for two boxes [x] and [y], it satisfies the following relation :

1 2 3 4 5 1 2 3 4 [x] = [1, 4] × [1, 3] x 1 x 2 (a) Box 1: w ∞ ([x]) = 3, w([x]) = √ 13 1 2 3 4 5 1 2 3 4 [x] = [1, 4] × [1, 2] x 1 x 2 (b) Box 2: w ∞ ([x]) = 3, w([x]) = √ 10
[x] ⊂ [y] [x] = [y] ⇒ w([x]) < w([y]). (2.18)
which is not true from w ∞ . Consequently, it will be used to compare the accuracy of the approximation given by the proposed algorithms. This is illustrated in Figure 2.4 where w and w ∞ are compared.

Basic operations on real numbers or vectors can be extended to intervals by considering interval computations on each component of the box and to vector operations such as dot or cross products.

Example 2.2. For instance, with 5,5] we have:

[x] = [1, 2] × [2, 9] and [y] = [-1, 2] × [-
[x] + [y] = [0, 4] × [-3, 14] (2.19) [x] • [y] = [-47, 49] (2.20) w([x]) = √ 50 (2.21) w ∞ ([x]) = 7 (2.22)
where • denotes the product.

Inclusion function

One of the purposes of interval analysis is to provide, for a large class of functions f , a fast and accurate way to compute a box which encloses the direct image of a box by f . This operator is called an inclusion function of f and is now introduced.

Definition and properties

Consider f : R n → R m . For any subset X of R n , the image set of f is defined by:

f (X) = {y ∈ R m | ∃x ∈ X, y = f (x)} . (2.23) An inclusion function [f ]:IR n → IR m of f : R n → R m is a function which satisfies: ∀ [x] ∈ IR n , f ([x]) ⊂ [f ] ([x]) . (2.24)
Thanks to interval arithmetic, and its implementation, the evaluation of an inclusion function provides a robust and reliable result taking into account round-off errors, imprecision on input data, numerical truncation, etc. To illustrate the notion of inclusion function, let us take

f from R 2 to R 2 and a box [x] ∈ IR 2 . The image set f ([x]
), in gray on Figure 2.5, can have any shape, be non convex and even disconnected. The inclusion function

[f ] of f makes it possible to compute a box [f ]([x]) guaranteed to contain f ([x]
).

An inclusion function

[f ] of f is said to be inclusion monotonic if: [x] ⊂ [y] ⇒ f ([x]) ⊂ f ([y]). (2.25) Consequently, [f ] is minimal if for any [x], [f ]([x]
) is the smallest box that contains f ([x]). The minimal inclusion function for f is unique and is denoted [f] * .

Enclosing the image set by a box axis-aligned introduces pessimism that can be seen as the set of extra points that belong to [f ]([x]) but not to the image set f ([x]).

The two main causes of pessimism are the multiple occurrences of variables in the expression of f , that can be balanced by formal simplification [START_REF] Araya | Exploiting common subexpressions in numerical CSPs[END_REF] and the wrapping effect. This over-estimation can be reduced [START_REF] Chabert | Computing the Pessimism of Inclusion Functions[END_REF] by using specific inclusion functions, multi-precision arithmetic [START_REF] Revol | Motivations for an Arbitrary Precision Interval Arithmetic and the MPFI Library[END_REF], Taylor extensions [START_REF] Berz | Verified Integration of ODEs and Flows Using Differential Algebraic Methods on High-Order Taylor Models[END_REF]], Bernstein expansions for polynomials [Garloff, 1985], etc. .26) [f] is said to be thin if, for any degenerate interval vector

An inclusion function

[f] is convergent if, for any sequence of nested boxes [x](k) lim k→∞ w ([x](k)) = 0 ⇒ lim k→∞ w ([f ]([x](k))) = 0. ( 2
[x] = {x} we have [f ]([x]) = f ({x}).
In other words, the image of a singleton should be a singleton. This is not always the case. In Chapter 5 and Appendix A, we will consider specific methods to deal with thick functions. These thick (i.e., non thin) functions associate to an element of R n , not a point, but a box of IR m .

x 1

x 2

[x 1 ] [x 2 ] [x] = [x 1 ] × [x 2 ] y 1 y 2 f([x]) [f] * ([x]) [f]([x]) f (x)

Natural inclusion function

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.

Example 2.3. The natural inclusion function

[f ] of f (x 1 , x 2 ) = x 1 • cos (x 2 ) is: [f ] ([x 1 ], [x 2 ]) = [x 1 ] • [cos] ([x 2 ]) (2.27)
Its evaluation is done by applying interval arithmetic. For instance:

[f ] [-2, 1], [ π 6 , π 3 ] = [-2, 1] • [cos] [ π 6 , π 3 ] = [-2, 1] • [ 1 2 , √ 3
2 ] = [- √ 3, √ 3 2 ].
(2.28)

If f is continuous, its natural inclusion function is monotonic, convergent and thin, but generally not minimal.

Contractors and Separators

In order to characterize efficiently the domain of the variables involved in a CN [Ceberio andGranvilliers, 2000, Araya et al., 2008], interval analysis algorithms rely on operators called contractors and separators. These operators are able to remove parts of a box that do not satisfy the constraints. This section introduces these two operators.

Contractors

This contractor operator, associated with a constraint, is able to reduce, or contract, an initial box by removing parts that does not satisfy the constraint.

Definition

A contractor C is an operator IR n → IR n (see e.g., [Di Loreto et al., 2007, Chabert andJaulin, 2009b]) such that

C([x]) ⊂ [x] (contractance) [x] ⊂ [y] ⇒ C([x]) ⊂ C([y]).
(monotonicity) (2.29)

A set X is consistent with the contractor C (we will write

X ∼ C) if for all [x], we have C([x]) ∩ X = [x] ∩ X. (2.30)
As a consequence, a contractor C X associated with the set X is an operator able to contract a box of R n without removing a single point of the subset X of R n to which it is associated. This is illustrated by Figure 2.6. We define the inclusion between two contractors C 1 and C 2 as follows: .31) Two contractors C and C 1 are equivalent (we will write C ∼ C 1 ) if we have: We define the negation ¬C of a contractor C as follows

C([x]) C * ([x]) X [x] [x] ∩ X ¬C([x])
C 1 ⊂ C 2 ⇔ ∀ [x] ∈ IR n , C 1 ([x]) ⊂ C 2 ([x]). ( 2
X ∼ C ⇔ X ∼ C 1 . ( 2
¬C([x]) = {x ∈ [x] | x / ∈ C([x])} (2.34)
Note that, as illustrated by the blue dashed area on Figure 2.6, ¬C([x]) is not a box in general, but a union of boxes. As a consequence, the negation of a contractor is not a contractor. This is a problem which will motivate later the introduction of separators.

Building contractors

Contractors can be built for a large class of constraints. When constraints are defined by equations and inequalities, such that f (x) ∈ [y], a contractor based on an inclusion function [f] of f can be easily built. Given a box [x], it is defined by:

C([x]) = ∅ if [f]([x]) ∩ [y] = ∅ [x]
otherwise (2.35) This binary contractor contracts to the empty set every box such as

[f]([x])∩[y] = ∅.
Of course, this contractor is not minimal.

An efficient improvement is the forward/backward contractor also called HC4-revise [Messine, 2008, Benhamou et al., 1999]. The principle is to evaluate f(x) using interval arithmetic (forward step) but instead of only considering the result of the intersection with [y] it retro-propagates the information to the initial domains. This is illustrated by the following simple example.

Example 2.4. Given three variables x ∈ [x], y ∈ [y] and z ∈ [z], the contractor for the constraint z = x -y is given by:

C -:    [z] [x] [y]    =    [z] ∩ ([x] -[y]) forward step [x] ∩ ([z] + [y]) backward step [y] ∩ ([x] -[z])    (2.36)
For instance C -( [4,5], [0,3], [-2, 2]) = ( [4,5], [2,3], [-2, -1])

Contractors algebra

Contractors can be combined to solve systems of constraints. If C 1 and C 2 are two contractors, we define the following operations s [Chabert and Jaulin, 2009b]:

(C 1 ∩ C 2 )([x]) = C 1 ([x]) ∩ C 2 ([x]), (2.37) (C 1 C 2 )([x]) = C 1 ([x]) C 2 ([x]), (2.38) (C 1 • C 2 )([x]) = C 1 (C 2 ([x])) , (2.39) C ∞ 1 ([x]) = C 1 • C 1 • C 1 • • • • • C 1 ([x]). (2.40)
where is the union hull defined by .41) Based on this formalism, primitive contractors can be combined in order to build more complex ones. Generally, these elementary contractors are built from dedicated algorithms.

[x] [y] = [x] ∪ [y] . ( 2
Example 2.5. In order to illustrate the use of contractors, we consider the example of the goniometric localization of a robot. This problem consists in using angles measured between beacons, whose position is known, and the robot for localization.

Let us consider the situation, depicted on Figure 2.7 with two beacons m 1 and m 2 .

We call bearing the angle θ i between the vehicle and the beacon m i , measured with respect to the x-axis direction. The measurements, taken with an uncertainty ∆θ of 3 Recall that two vectors u, v of R 2 are colinear if their determinant is equal to zero, in other words if det(u, v) = 0. Thus, for each beacon m i = (m i 1 , m i 2 ) , we have the relation:

det m i 1 -x 1 m i 2 -x 2 , cos(θ i ) sin(θ i ) = 0. (2.42)
The set of feasible positions for the vehicle is defined by:

X i = x ∈ R 2 | ∃θ i ∈ [θ i ], (m i 1 -x 1 ) • sin(θ i ) -(m i 2 -x 2 ) • cos(θ i ) = 0 . (2.43)
These sets are depicted on Figure 2.7. Using a forward/backward propagation algorithm, an elementary contractor C i consistent with X i can be built. With two beacons, the position belongs to the intersection of the X i and the contractor C consistent with the solution set X is then:

X = X 1 ∩ X 2 (2.44) ∼ C = C 2 • C 1 (2.45)
When we use contractors, we are looking for boxes as small as possible. For that

m 1 m 2 X 2 X 1 θ 2 θ 1 ∆θ -6 -4 -2 0 2 4 6 8 -6 -4 -2 0 2 4 6 8 Figure 2.7: Goniometric localization purpose, C is iteratively composed, denoted C ∞ = C • C • • • • • C,
in order to reach the fixed point. Since all contractors used are monotonic, it can be proved that the iterated composition C ∞ will always converge to the greatest box [z], such that C ∞ ([z]) = [z], regardless the order of the contractors in composition [START_REF] Montanari | Constraint relaxation may be perfect[END_REF]. 2.8 illustrates, step by step, the propagation process. On Figure 2.8a, the contraction C([x 0 ]), depicted by the orange box, is decomposed. The blue dashed area has been removed by C 1 , the cyan one by C 2 . The yellow box, on Figure 2.8b, shows the contraction done by applying C again on the orange box. This process is repeated until no more contraction happens anymore. The final result is shown on Figure 2.8c. Now, in our case, with respect to the configuration of angles θ 1 and θ 2 , the propagation process provides an enclosure of the solution set that is pessimistic. If a thinner approximation is required, a classical solution to reduce this pessimism is to bisect (split into two parts) this box, and appling C ∞ again on each sub-boxes. This process is illustrated in Figure 2.8c. This procedure can be recursively used until the width of resulting subboxes are smaller than a given threshold. It produces the subpaving (union of non-intersecting boxes) that encloses X as depicted in Figure 2.8d This algorithm, which classifies the space, will be detailed in Section 2.4.2. The main issue with this method is to choose the right number of bisections and to characterize the pessimism introduced by the approximation. Since boxes, which are completely included in X, are never contracted and always bisected, their accumulations (see the yellow boxes in Figure 2.8d) lead to a huge computational cost. This justifies the introduction of separators, that aim at providing an inner approximation of the solution set, and the need for developing minimal contractors in order to avoid as much as possible bisections. Note also that the formulation of Equation (2.43) cannot be used to get an inner approximation of X and will be reformulated in Chapter 3 in order to build a dedicated contractor.

Given an initial box [x

0 ] = [-4, 6] 2 , Figure
m 1 m 2 [x 0 ] C2 • C1([x0]) -6 -4 -2 0 2 4 6 8 -6 -4 -2 0 2 4 6 8 (a) C 1 ([x 0 ]) m 1 m 2 C([x 0 ]) C • C([x 0 ]) -6 -4 -2 0 2 4 6 8 -6 -4 -2 0 2 4 6 8 (b) C 1 • C 2 • C 1 ([x 0 ]) m 1 m 2 C ∞ ([x 0 ]) C ∞ ([x 1 ]) C ∞ ([x 2 ]) -6 -4 -2 0

Separators

Given a set X, the contractor consistent with X only provides an outer approximation. If X has a non-empty volume (there exists points in X which do not belong to its border), it can be useful to prove that a part of the initial space belongs effectively to the solution set. For that, the complementary contractor consistent with X needs to be considered. In order to characterize an inner and outer approximation of the solution set, we introduce the notion of separator.

Definition

A separator S associated with the set X is an application such as

S : IR n -→ IR n × IR n [x] -→ ([x in ], [x out ]) (2.46)
with the following properties

(i) [x] = [x out ] ∪ [x in ] (ii) [x out ] ∩ X = [x] ∩ X (iii) [x in ] ∩ X = [x] ∩ X (2.47)
A separator can also be seen as a pair of contractors {S in , S out } such that, for all [x] ∈ IR n , we have

S in ([x]) ∪ S out ([x]) = [x] (complementarity). ( 2

.48)

A set X is consistent with the separator S (we will write X ∼ S), if X ∼ S out and X ∼ S in . (2.49) where X = {x | x / ∈ X}. This notion of separator is illustrated by Figure 2.9.

X [x] S out ([x]) [x] S in ([x]) ¬S in ¬S in X Figure 2
.9: Illustration of a separator applied on different initial boxes. The outer contractor removes the blue dashed area and the red dashed area is removed by the inner contractor

We define the inclusion between two separators S 1 and S 2 as follows

S 1 ⊂ S 2 ⇔ S in 1 ⊂ S in 2 and S out 1 ⊂ S out 2 .
(2.50)

A separator S is minimal if S 1 ⊂ S ⇒ S 1 = S. (2.51)
It is trivial to check that if S is minimal then the two contractors S in and S out are both minimal.

Separator algebra

The algebra for separators is a direct extension of the contractor algebra [Chabert and Jaulin, 2009b]. The main difference is that the contractor algebra does not allow any non monotonic (or decreasing) operation. It means that, if a contractor C is defined by an expression E of other contractors C i , then we always have center (1,3) (4,23) radius [3, 6] [1, 5] [3, 4] Let us take S i = {S in i , S out i }, we could define the following operations:

∀i, C i ⊂ C i ⇒ E (C 1 , C 2 , . . . ) ⊂ E C 1 , C 2 , . . . . ( 2
S 1 ∩ S 2 = {S in 1 ∪ S in 2 , S out 1 ∩ S out 2 } (intersection) S 1 ∪ S 2 = {S in 1 ∩ S in 2 , S out 1 ∪ S out 2 } (union) {q} S i = {m-q-1} S in i , {q} S out i (relaxed intersection) S 1 \S 2 = S 1 ∩ S 2 . (difference) (2.54)
These operations are used to combine constraints. More details about separators are given in [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF].

The following example shows how separators can be used to manipulate sets defined with equations.

Example 2.6. Consider the set X defined by:

X = (X 1 ∩ X 2 ) ∪ (X 3 \ X 1 ) . (2.55)
where X 1 , X 2 , X 3 are three rings with center and radius defined in table 2.1. For each set, a separator S i = C i , C i can be built using a forward/backward propagation method. Using the separator algebra, the separator S = C, C consistent with X is:

S = (S 1 ∩ S 2 ) ∪ S 3 ∩ S 1 .
(2.56)

It will automatically apply the De Morgan rules:

C = (C 1 ∩ C 2 ) ∪ C 3 ∩ C 1 C = C 1 ∪ C 2 ∩ C 3 ∪ C 1 (2.57)
Combined with a Paver (see next section), the separator S can be used to characterize the set defined by Equation (2.55).

Separators allow us to deal efficiently with set operations by manipulating the inner and outer approximation of sets.

Paver

Separators are used by pavers to characterize sets defined by constraints or set operations. A paver is a branching algorithm which calls the separator S to classify ), the part of [x] that is proved to be inside X, into X -and also into X + .

Step 5 computes ∂S ([x]) by intersecting [x out ] and [x in ]. If this box is too small (i.e. with a width smaller than ε), it is stored inside X + and will not be studied anymore. Otherwise, it is bisected, generally along its largest dimension, at Step 7 and stored into L waiting to be processed. After completion of the algorithm, we have the enclosure X -⊂ X ⊂ X + .

(2.58) Remark 2.2. . For the implementation, the resulting paving can be represented by a binary tree (i.e., each node has two sons or is a leaf). The binary tree is said to be minimal if for any node i 1 (not the root) with brother i 2 and father j, we have

(i) [x in ](i 1 ) = ∅, [x out ](i 1 ) = ∅ (ii) [x in ](j) ∩ [x out ](j) = ([x in ](i 1 ) ∩ [x out ](i 1 )) ([x in ](i 2 ) ∩ [x out ](i 2 )) , ( 2
.59) From a given tree generated by the paver, it is possible to simplify it into a minimal tree without changing the approximation for X. This simplification can be done by scanning all nodes of the tree upward from the leaves to the root. If for a given node, properties (i) or (ii) of Equation (2.59) are not satisfied, then the two brother nodes i 1 , i 2 with the father j are reunited into a single node. The procedure has a complexity which is linear with respect to the number of nodes of the tree. It Pull [x] from L;

4: {[x in ], [x out ]} = S([x]); 5: Store [x] \[x in
] into X -and also into X + ;

6:

[x] = [x in ] ∩ [x out ]; 7: if w ∞ ([x]) < ε then 8: store [x] in X + 9:
else 10:

bisect [x] and push into L the two resulting boxes 11:

end if 12: end while will allow us to drastically reduce the number of elements of the subpaving. An illustration is provided on Figure 2.11. Figure 2.11: Simplification of the paving generated by the paver; left: before simplification, right: after. By convention, red boxes are proved to be inside the solution set, blue boxes are outside, and yellow boxes are undetermined.

Remark 2.3. In the worst case, the complexity of the algorithm is exponential with respect to the dimension n of the initial box, and terminates after less than ( w([x 0 ]) +1) n iterations. However, in practice, the number of iterations is O(A 1 n-1 ) where A is the area of the accumulation zone (in yellow) [Jaulin, 1994].

Remark 2.4. The area painted in red by the paver is not

[x in ] (or [x out ]) but [x]\[x in ]
which corresponds to the part of the initial box which has been removed by the inner contractor.

Example 2.7. Consider the set

X = x ∈ R 2 , (x 1 -2) 2 + (x 2 -2.5) 2 ∈ [1, 1.5] 2 (2.60)
which corresponds to a ring centered in (2, 2.5). The minimal contractor C X consistent with X can be built using a forward-backward constraint propagation. Figure 2.12: Pavings associated to Example 2.7, Left: paving obtained using the contractor, Right: paving obtained using the separator. Red boxes belong to X (the ring); blue boxes are outside of X. No conclusion can be given on the yellow boxes.

Using a paver, this contractor is then used to obtain an outer approximation for X. This is illustrated by Figure 2.12 (left) where C X removes parts of the space outside X (painted blue). But due to the consistency property, see Equation (2.30), C X has no effect on boxes included in X. A box partially included in X cannot be eliminated and is bisected, except if its length is larger than a given value ε.

A separator consistent with X makes it possible to get an inner approximation.

We are now able to quantify the pessimism introduced by the set inversion and to prove the existence of solutions. In a computational point of view, the paver becomes faster since the accumulation area is of dimension 1, instead of 2.

Inversion of separators

As it has been introduced previously, a set can be represented by its associated separator. The inversion of a separator will make it possible to get a separator for the inverse image of a set by a function. The inverse of a set Y ⊂ R n by a function f : R n → R m is defined as

X = f -1 (Y) = {x | f (x) ∈ Y} .
(2.61)

The function f can be a translation, rotation, dilation, projection, or any other function. If C Y is a contractor for Y, a contractor C X for X can be defined using a generalization of the forward-backward contractor as shown in [Chabert, 2013].

The contractor C X is called the inverse of C Y by f and we write

C X = f -1 (C Y ). If
S Y is a separator associated with a set Y, we define the inverse of the separator S Y as follows

f -1 (S Y ) = f -1 (S in Y ), f -1 (S out Y ) . (2.62) Theorem 2.1. The separator f -1 (S Y ) is a separator associated with the set X = f -1 (Y), i.e., f -1 (Y) ∼ f -1 (S Y ) .
(2.63)

Proof. We have : [Chabert, 2013])

Y ∼ {S in Y , S out Y } ⇔ Y ∼ S out Y , Y ∼ S in Y (see 2.49) ⇒ f -1 (Y) ∼ f -1 (S out Y ), f -1 Y ∼ f -1 (S in Y ) (see
⇔ f -1 (Y) ∼ f -1 (S out Y ), f -1 (Y) ∼ f -1 (S in Y ) ⇔ f -1 (Y) ∼ {f -1 (S in Y ), f -1 (S out Y )} (see ((2.49))) ⇔ f -1 (Y) ∼ f -1 (S Y ) (see (2.63)) (2.64)
which terminates the proof.

The following example shows how the inversion of a separator can be used.

Example 2.8. Consider the polygon M defined on Figure 2.13a, for which a separator can be built (see Section 2.5.1 for more details) and a rotation of an angle of π 4 denoted R π 4 . Figure 2.13b depicts the set R π 4 (M) computed using Theorem 2.1. Figure 2.13c shows how separators can be combined to compute the set R π 4 (M) ∩ M.

Dedicated Contractors and Separators

A contractor can be used to represent an information we have on a vector of R n . This information can be an equation or an inequality, but it can also come from

S [∆x] = C ∆S ([x]) ? ? [x 2 ] [x 1 ] (a) Boundary contraction using C ∆S S [∆x] = C ∆S ([x]) T S (x 1 ) = true T S (x 2 ) = f alse [x in ] [x out ] (b) Subboxes classification Figure 2
.14: Illustration of boundary-based separator. An initial box is contracted using C ∆S (left). On each removed part of the initial box, a point is tested and, by continuity, the subbox is classified (right). Note that, for the illustration, both [x in ] and [x out ] have been enlarged. a database. For instance, it can be the map of the surrounding environment or the shape of an object. This is an important point in our SLAM context where we have to combine equations (coming from observations or the evolution of the robot) with some other knowledge coming from a database. We will illustrate the principle by building some dedicated separators associated to data sets.

Boundary-based separator

Sometimes, isolating explicitly the corresponding inner and outer contractor associated to a set S is not so easy. However, if a contractor on the boundary of S and a test which is able to prove that a given point is inside or outside S are available, it is still possible to build a separator using the following boundary-based approach. Given a contractor C ∆S on the border of set S, and a test T S such that, for .66) This method is now applied in order to build a separator associated to a polygon.

x ∈ R n T S (x) = true if x ∈ S f
S ∆S ([x]) = {[x in ], [x out ]} = {[∆x] ∪ [x 2 ], [∆x] ∪ [x 1 ]} . ( 2

Polygon separator

Two-dimensional structured environments, such as rooms or buildings, can be described by a set of segments that forms a closed polygon, convex or not, with eventually holes inside. When a polygon is convex, this set is described by a system of linear inequalities that can be solved with linear solvers. When not, this method cannot be applied and the boundary-based approach can be used instead. We will now introduce how to build a separator for general polygons of R 2 .

Contractor on the border

Consider an oriented segment [a, b] where (a, b) ∈ R 2 . The point m is in the segment if it satisfies the two following equations det (b -a, a -m) = 0 min (a, b) ≤ m ≤ max (a, b) .
(2.67)

A general polygon P is composed of N oriented segments. The border ∆P of the polygon satisfies the following constraint:

∆P = m ∈ R 2 , ∃i ∈ [[1, N ]], m ∈ [a i , b i ] . (2.68) Let us take C a i ,b i as a contractor for the segment [a i , b i ].
The contractor for ∆P is:

C ∆P = N i=1 C a i ,b i . (2.69)
This contractor is illustrated on Figure 2.15a. As expected, it only focuses on the border of the polygon.

Remark 2.5. Because the union of minimal contractors is minimal (see Proposition 3.1), C ∆P is a minimal contractor for the border of the polygon P.

Test: point in polygon

To identify if a point is inside a polygon, we will use the Winding Number [Krantz, 2012] which represents the total number of times that a curve travels counterclockwise around the point. The winding number depends on the orientation of the curve, and is negative if the curve travels around the point clockwise. Let us take a polygon

P with vertices V 1 , V 2 , . . . , V n = V 1
and m a point not on the border of P. The winding number is defined by: So, as shown in Figure 2.16, if m is outside P we will have wn(m, P) = 0, otherwise if m is inside, wn(m, P ) = 1.

wn(m, P) = 1 2π n i=1 θ i = 1 2π n i=1 arccos (V i -m).(V i+1 -m) (V i -m) (V i+1 -m) (2.
Remark 2.6. Faster algorithms that avoid the use of trigonometric functions can be found in [O'Rourke, 1998].

Finally, we build a separator S P for the polygon P. Figure 2.15b shows the set X (in red) of all points inside the polygon.

Since this approach is generic, this type of separator can be extended to any set for which its border can be described by a union of elementary curves. The main difficulty is to find an efficient test to check whether a point belongs or not to the set.

Separator on subpavings

Data sets can also be stored in the form of a subpaving, i.e. a union of nonoverlapping boxes generally organized as a binary tree, which represents an inner and an outer approximation of a set. From a given subpaving, a separator (or contractor) can be defined. Considering that each node N j has two sons N i 1 and

N i 2 , Algorithm 2.2 allows to separate a box [x] into [x in ] and [x out ] in a recursive manner.
The complexity of such an algorithm is O(log(n) + m) where n is the number of nodes and m is the number of nodes intersecting the initial box [de Berg et al., 2008].

Algorithm 2.2 separate(in: N , [x], out: [x in ], [x out ]) 1: if [x] = ∅ or N is a leaf then Return {[x in (j)] ∩ [x], [x out (j)] ∩ [x]} 2: end if 3: [x] = [x] ∩ [x in (j)] ∩ [x out (j)] 4: {[x in 1 ], [x out 1 ]} = separate(N i 1 , [x]) 5: {[x in 2 ], [x out 2 ]} = separate(N i 2 , [x]) 6: Return {[x in (j)] ∪ [x in 1 ] ∪ [x in 2 ], [x out (j)] ∪ [x out 1 ] ∪ [x out 2 ]}
Note that any subpaving built using a paver can be simplified into a minimal subpaving (see Equation (2.59)) and used as a separator. Consequently, the results of long-running algorithms can be saved efficiently for future use.

Separator from images

An unstructured dataset X can be described by an occupancy grid, represented as a binary image, denoted i, in which the value of each pixel is 1 if it belongs to X, 0 otherwise. The image contractor and separator are now introduced.

Image contractor

The image contractor was first introduced in [Sliwka, 2011, Guyonneau et al., 2013]. This contractor is based on the summed area table (also known as integral image), denoted I, widely used in computer vision [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF] and defined by:

I (n 1 , n 2 ) = n 1 ≤n 1 n 2 ≤n 2 i (n 1 , n 2 ) .
(2.71)

Using the integral image, the number of 1-valued pixels 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

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

IN 2 -→ N φ ([n]) -→ I(A) + I(C) -I(B) -I(D) (2.
]) = φ([n]).
Denote by C i the image contractor associated to the image i.

Consider [n] ∈ IN 2 and C i ([n]) = [m] that satisfies: L :              m - 1 = max(x ∈ [m 1 ], φ([m - 1 ; m] × [m 2 ]) = 0) m + 1 = min(x ∈ [m 1 ], φ([m; m + 1 ] × [m 2 ]) = 0) m - 2 = max(x ∈ [m 2 ], φ([m 1 ] × [m - 2 ; x]) = 0) m + 2 = min(x ∈ [m 2 ], φ([m 1 ] × [m; m + 2 ]) = 0) (2.73)
The min and max can be computed using a dichotomy, which has a logarithmic complexity. To apply a continuous contractors on a data set described by a geo-referenced image, world continuous coordinates must be converted into discrete image ones. Let x 0 be the coordinate of the top-left corner and = ( x , y ) the horizontal and vertical pixel size. The change between image coordinate [n] to workspace coordinates

n 1 n 2 A B C D [n] (a) φ function evaluation n 1 n 2 A B C D A' B' C' D' C([n]) (b) Contraction of the box
[x] is [x] = x 0 + [n] • ({0} ∪ { }) (2.74)
Remark 2.7. Because of the image convention, with vertical axis oriented downward, the value of y is negative and

({0} ∪ { }) is equal to [0, x ] × [ y , 0].
Remark 2.8. When a set X is digitized as an image, one valued pixel must describe an over approximation of X in order to keep the consistency property of the contractor. Consequently, the contractor is minimal with respect to the image but not with respect to X itself.

Remark 2.9. The image integral can be generalized to the n-dimensional case [Tapia, 2011]. In particular, the image contractor has been extended to the threedimensional case in [START_REF] Desrochers | Setmembership approach to the kidnapped robot problem[END_REF] to handle spacial occupancy grids, built from a 3D dense point clouds.

Image separator

A separator can be built using two complementary image contractors. Figure˜2.18 shows how, from a digitized set X, images representing an upper and lower approximations of X can be used to build the separator.

X ∆X

X - n 2 n 1 (a) Digitization of set X n 2 n 1 (b) Image used by S out n 2 n 1 (c) Image used by S in Figure 2
.18: Image approximation of the set X, orange dashed . Black pixels, which are completely included in X, belong to a lower approximation X -. The upper approximation X + = X -∪ ∆X contains black pixels and those which intersect the boundary of X. The outer contractor will be built using the left image while the inner contractor will be built using the complementary of X -(right one).

Remark 2.10. The image separator, as well as the image contractor, has a very high speed of execution regardless the size of the input box. It only involves sums and differences of integer numbers. However, data are discretized on a regular grid and a trade off needs to be found between the storage space, the size of area, and the resolution of the grid. On the contrary, a subpaving can have better accuracy on the border of a set with a smaller memory footprint but with an high access cost.

Remark 2.11. The separators on images and subpaving are minimal with respect to the dataset. But, they are not minimal when approximated by an image or a subpaving. As illustrated by Figure 2.18a, the true border of the set, in orange, is enclosed by the gray pixels that introduce pessimism.

Differential contractor on tube

In the formalism presented in the introduction, the evolution of the state of the robot x(t), t ∈ R, is assumed to be constrained by a differential equation such as

ẋ(t) = f (x(t)), u(t)) . (2.75)
where u(t) is the input of the system. The temporal evolution of a variable x is described by a trajectory, denoted by x(.) : R → R n . When this trajectory is uncertain, interval analysis makes it possible to enclose it in a Tube, which is an interval of trajectories.

The contractors and separators, introduced previously, are mainly dedicated to constraints that can be applied at a given time t. Now, we present briefly the tools needed to deal with an uncertain trajectory, constrained by a differential equation.

A tube [x](.) as defined in [Le Bars et al., 2012, Bethencourt andJaulin, 2014], is an interval of trajectories [x -(.), x + (.)] such that ∀t, x -(t) ≤ x + (t). A trajectory x(.) belongs to the tube [x](.) if ∀t, x(t) ∈ [x](t). Figure 2.19 provides an illustration of a tube. It is important not to make the confusion between x(.) which is a trajectory and x(t) which is a vector of R n .

[x](.)

x + (.)

x -(.)

x(.) ), in grey, enclosing a trajectory x(.)

[x] t x(t 1 ) t 1 Figure 2.19: A tube [x](.
In this thesis, tubes are mainly used to propagate a contraction, done at a given time step, to the whole trajectory. For that purpose, a dedicated contractor, denoted C d dt , has been developed in [START_REF] Rohou | Guaranteed computation of robot trajectories[END_REF]. It is associated to the following constraint:

     ẋ(.) = u(.) x(.) ∈ [x](.) u(.) ∈ [u](.) (2.76)
The contractor C d dt allows us to contract the tube [x](.) with respect to its derivative u(.) enclosed by a tube [u](.).

Remark 2.12. In the initial definition proposed in [START_REF] Rohou | Guaranteed computation of robot trajectories[END_REF], the contractor C d dt requires the derivative of the tube as an argument. It is written as

[x](.) := C d dt ([x](.), [u](.)).
(2.77)

Now, since the tube [u](.) cannot be contracted with this contractor (the proof can be found in the corresponding article), we choose to simplify the notation by omitting the last argument. The tube [u](.) enclosing the derivativeis handled as an implicit parameter of the contractor.

Example 2.9. In order to illustrate how tubes are used, we consider the following example:

         ẋ(t) = u(t) z 1 = x(t 1 ) x(0) ∈ [-3, 3] z 1 ∈ [2.2, 2.8] (2.78)
where u(.) is assumed to belong to a bounded tube [u](.), depicted in Figure 2.20b, and defined by:

∀t > 0, [u](t) = 4 • sin(t -5) + (t -3.3) • [-0.1, 0.1].
(2.79)

In our robotics applications, [u](.) is built from proprioceptive measurements, such as the speed or the angular velocity of the robot. 3,3]. Only trajectories which are not consistent with the differential equation and the initial conditions are removed.

](.) defined by ∀t > 0, [x 0 ](t) = [-∞, ∞] and [x 0 ](0) = [-
In the same way, a local information, such as x(t 1 ) = z 1 , can be propagated using C d dt to the whole trajectory. This is illustrated by the dark grey tube on Figure 2.20a. The light grey tube is contracted in order to be consistent with the observation [z 1 ]. Note that contracting a tube at a given time step is not as trivial as doing [x](t 1 ) := [x](t 1 ) ∩ [z](t 1 ). It requires the use of a specific contractor which is not presented here. More details on tubes can be found in [Rohou, 2017]. In this thesis, when we manipulate tubes, we use the open-source library called tubex-lib available at:

t 1 [z](t 1 ) -5 5 10 t [x](.) (a) Tube [x](.) -4 -2 2 4 t [u](.) (b) Tube [u](.)
• http://www.simon-rohou.fr/research/tubex-lib

Applications

We now consider two classical applications that will illustrate how interval analysis, and mainly separators, can be used to solve robotics problems. The first application, taken from [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF], deals with path-planning problems. The second one deals with the initial localization problem using a scanning laser rangefinder on a map defined by a polygon.

Application to path planning

The goal of path planning is to find a collision-free path for a robot in a given space with obstacles. The issue of path planning in a known environment has been addressed since many years (see, e.g., [Laumond, 1986, Lozano-Perez, 1981[START_REF] Yap ; Ó'dúnlaing | A "retraction" method for planning the motion of a disc[END_REF], Koditschek, 1987]) and can easily be combined with set-membership techniques to take into account some uncertainties [START_REF] Ceccarelli | Set Membership Localization and Map Building for Mobile Robots[END_REF]. Most approaches are based on the concept of configuration space (C-space). Each coordinate of the C-space represents a degree of freedom of the object. An example of such robots are industrial robots which are kinematic chains of adjacent links are connected by n prismatic or rotary joints, each with one degree of freedom. The positions and orientations of each link can be characterized by n real numbers, which are the coordinates of a single n-dimensional point in the C-space [Lozano-Perez, 1983]. The feasible configuration space M is the subset of the C-space corresponding to feasible configurations of the robot. Partitioning

1 x 1 x 2 f 2 (x) f 3 (x) Y Figure 2
.21: Two-dimensional wire loop game. For a feasible configuration (x 1 , x 2 ), as represented here, the gray segment should cross the wire. Is it possible to follow the circular path ? the C-space with subpavings makes it possible to solve the problem using graph algorithms [Brooks andLozano-Perez, 1985, Jaulin, 2001].

The objective of this application (extracted from [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF]) is to show how separators could solve the path planning problem on a simple example which is a 2D version of the wire loop game. This game involves a metal loop on a handle and a length of curved wire (see Figure 2.21). The player holds the loop in one hand and attempts to guide it along the curved wire without touching it. In our 2D version of this game, the player is an articulated robot with two rotary joints and the loop is a segment. The curved wire corresponds to the boundary of a set Y with an inside (grey in the figure) and an outside part. The length of the first and second arms are 4 and 2, respectively. The length of the loop is 1. The feasible configuration space is

M = {(x 1 , x 2 ) | f 2 (x) ∈ Y and f 3 (x) / ∈ Y} = f -1 2 (Y) ∩ f -1 3 Y , ( 2.80)
where

f (x) = 4 cos x 1 sin x 1 + cos (x 1 + x 2 ) sin (x 1 + x 2 ) , ∈ {2, 3} . (2.81)
If S Y is a separator associated with Y (built as explained in Section 2.5.1), then a separator for M is

S M = f -1 2 (S Y ) ∩ f -1 3 S Y . (2.82)
The paver is able to approximate the feasible configuration space as illustrated by Figure 2.22a. A graph-based method is able to compute a feasible path (see the black path of the figure) which corresponds to a solution of our wire loop game.

Remark.

For our example, the set to be characterized is given by

M = f -1 2 (Y) ∩ f -1 3
Y and the separator is obtained by the same expression

S M = f -1 2 (S Y ) ∩ f -1 3 S Y .
If we separate this expression, we get the contractor counterpart

S in M = f -1 2 (S in Y ) ∪ f -1 3 (S out Y ) S out M = f -1 2 (S out Y ) ∩ f -1 3 (S in Y ).
(2.83)

The separator formalism thus makes it possible to get directly the separator from the expression of the set that we want to characterize. This was not possible using contractor algebra which does not allow any decreasing operation. With contractors, the user has to rewrite the complementary expression. Moreover, the separator algebra requires less computations since a part of the common computation made for S in M and S out M is factorized.

Simple static map-based localization

A simple localization problem within an initial 2D map illustrates how separator operations can be used to compute the position of a robot. 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˜2.23a) 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 . At a given time step, N measurements are acquired with outliers. An outlier is defined as a measurement that is not consistent with the bounded error model used. The full situation is depicted on Figure 2.23a.

We also denote by a i (sresp. b i ) the nearest (resp. furthest) impact point associated with the distance

d - i (resp. d + i ).
The parametrization is shown in Figure˜2.23b. 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 and the uncertainty on α i is assumed to be negligible.

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 (2.84)
Given S M associated with M, the separator S i consistent with X i is:

S i = f - i (S M ) ∩ f + i (S M ) (2.85)
Without outliers, the feasible solution set must be consistent with the intersection of all X i . However, when dealing with outliers, this intersection is often empty. To solve this issue, the q-relaxed intersection is used to compute the set of feasible positions compatible with at least N -q measurements.

X = {q} i X i ∼ {q} i S i (2.86)
In this test case, 56 ranges are used with an uncertainty of +/-0.1 meter. Results of the set inversion are shown in Figure 2.24 with q = 15. Because measurements are redundant and accurate, the solution set is small.

Finding automatically the number of allowed outliers is a hard task. Algorithms such as GOMNE (Guaranteed Outlier Minimal Number Estimator) can be used [Jaulin et al., 2002]. Roughly, a series of set inversion are performed with an increasing value of q until the solution set becomes not empty. An application with 3D data can be found in [START_REF] Desrochers | Setmembership approach to the kidnapped robot problem[END_REF].

Conclusion

In this chapter, the interval analysis theory was introduced in the context of robotics localization. Many solution sets that may correspond to the position of a robot can be decomposed into a combination of unions, intersections, projections or inversions of elementary sets. Using operators such as contractors and separators, these sets can be characterized using interval analysis. For that purpose, constraint propagation and bisections mechanisms are combined inside an algorithm, called a paver, in order to get an inner and an outer approximation of a solution set. Now the performance of the paver relies on the efficiency of the elementary operators that are used. We introduced some dedicated contractors and separators associated to datasets which has been illustrated on two robotics examples. In the next chapter, another contractor and separator will be introduced to deal with the polar constraint, often met in localization when angle and bearing are measured.

Introduction

In the previous chapter, interval analysis tools have been introduced. These tools, such as contractors and separators, are used to solve rigorously complex problems involving uncertainties and nonlinear equations [Ceberio andGranvilliers, 2000, Kreinovich et al., 1997]. Coupled with a Paver, they allow computing an inner and an outer approximation of a set defined by constraints.

Contractor programming relies on a catalog of elementary contractors. Most of the time, these elementary contractors are built using interval arithmetic [Moore, 1966]. Then, by combining all these elementary contractors, we can construct a more sophisticated contractor consistent with the solution set of the problem we want to solve. The principle can be extended to separator programming [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF] in order to compute an inner and an outer approximation of the solution set. Now, combining contractors introduces a pessimism which has to be balanced by additional bisections performed by the paver. For more efficiency, it is important to extend the catalog by adding some new specific contractors.

The following example motivates briefly the advantage of using a dedicated and minimal contractors or separators.

Example 3.1. The set X = f -1 ([y]) with [y] = [3,5] and:

f (x) = x 2 1 + |2x 2 | -x 1 (3.1)
can be characterized using two different separators. The first one S is built from the natural inclusion function of f . Because the variable x 1 appears twice, neither this inclusion function is minimal nor S. The second separator S * is built after changing the expression of f into:

f * (x) = (x 1 - 1 2 ) 2 - 1 4 + |2x 2 | (3.2)
in which multiple occurrences of x 1 have been removed. This leads to a minimal inclusion function and to the minimal separator. Figure 3.1 shows the resulting subpaving using S (left) and S * (right). The minimal behavior of S * can be observed by the fact that each box touches the border of X. The number of bisections is also smaller when using a minimal separator. In this case, 839 bisections have been required to get the left subpaving while the right one requires 593 bisections only.

The first outer approximation of the minimal separator, depicted by dashed white boxes, is more accurate. This point is critical when dealing with high dimensional problems, for which bisections need to be avoided. The more efficiently a separator is, the better the final approximation will be.

To be as efficient as possible, minimal separators, and therefore contractors, must be used when available. As illustrated by the previous example, they lead to 56 All sides touch the border of X. This is not the case with S where the pessimism has to be balanced with additional bisections.

Building Minimal Contractors

an accurate outer approximation and reduce the number of bisections needed to compute an inner and an outer approximation of the solution set. This justifies the need for developing specific operators dedicated to our problems.

In this chapter, we propose some new theorems in order to build more easily optimal contractors/separators consistent with equations often used, for instance, in the field of robotics [START_REF] Kieffer | Guaranteed Mobile Robot Tracking Using Interval Analysis[END_REF], Daney et al., 2006, Langerwisch and Wagner, 2012]. As an application, we will consider the polar constraint associated to the change of coordinates between Cartesian and polar form [START_REF] Candau | Complex interval arithmetic using polar form[END_REF]]. This constraint is essential for solving the localization of robots when both goniometric and/or distance measurements are available [Colle andGalerne, 2013, Di Marco et al., 2001]. Robotics applications, given in Section 3.4, will illustrate the use of the polar separator.

Building Minimal Contractors

Building minimal contractors for sets, defined by inequalities, can sometimes be done using interval based methods [START_REF] Araya | Exploiting common subexpressions in numerical CSPs[END_REF]. In the special case where each variable occurs only once in the expression and when all involved operators are continuous, for instance with the constraint a + sin (b + c • d) = 0, a simple interval evaluation followed by a backward propagation in the syntactic tree of the constraint provides a minimal contraction [START_REF] Benhamou | Revising hull and box consistency[END_REF]. When the constraint is monotonic with respect to all variables, then again, it is possible to reach the minimality [START_REF] Chabert | Hull Consistency Under Monotonicity[END_REF]. In this section, we propose some new results that will allow us to extend the class of constraints for which we can provide a minimal contractor.

Transformations of contractors and separators

We now present some results obtained in [Sliwka, 2011] and [Herréro, 2006] about the symmetries and the minimality of contractors. These results will be used in the next part to get the minimal contractor for specific constraints.

Proposition 3.1. If C 1 and C 2 are the two minimal contractors consistent with X 1 and

X 2 then C 1 C 2 is the minimal contractor consistent with X 1 ∪ X 2 .
Proof. The minimal contractor consistent with

X 1 ∪ X 2 is (X 1 ∪ X 2 ) ∩ [x] = (X 1 ∩ [x]) ∪ (X 2 ∩ [x]) ((A ∪ B) ∩ C = (A ∩ C) ∪ (B ∩ C)) = X 1 ∩ [x] ∪ X 2 ∩ [x] ( A ∪ B = A ∪ B ) = X 1 ∩ [x] X 2 ∩ [x] ( [x] ∪ [y] = [x] [y] ) = C 1 ([x]) C 2 ([x]) (minimality of C 1 and C 2 ) (3.
3) which terminates the proof.

Definition 3.1. A bijective function

f : R n → R n is box-conservative if for all A ⊂ R n , f ( A ) = f (A) . (3.4) Proposition 3.2. If f is box conservative so is f -1 .
Proof.

f -1 ( A ) = f -1 ( f • f -1 (A) ) (f is bijective) = f -1 • f ( f -1 (A) ) (f is box conservative) = f -1 (A) . (3.5) Example 3.2. A rotation from R 2 to R 2 of angle α is box-conservative iff α = k • π 2 , k ∈ Z. Definition 3.2. If f : R n → R n
is a bijective function, we define the image by f of a contractor as follows:

f (C X ) = f • C X • f -1 . (3.6)
This new definition will make it possible to extend the contractor algebra to more complex operations. For instance (f

0 f 1 f 2 ) (C X ) defines the following contractor [x] → C X ([x]) f • C X • f -1 ([x]) f 2 • C X • f -2 ([x]) . (3.7)
This contractor is consistent with the set X ∪ f (X) ∪ f • f (X), as shown at least partly by the new following proposition.

Proposition 3.3. Define a set X for which we have a minimal contractor

C X . If f is box-conservative, then f (C X ) is the minimal contractor for f (X). Proof. The minimal contractor for f (X) is [x] → f (X) ∩ [x] . Now, f (X) ∩ [x] = f • f -1 ( f (X) ∩ [x] ) (f is bijective) = f ( f -1 • f (X) ∩ f -1 ([x]) ) (f -1 is box conservative) = f ( X ∩ f -1 ([x]) ) (f -1 • f (X) = X) = f (C X (f -1 ([x]))) (minimality of C X ) = (f (C X )) ([x]) (definition of f (C X ) ) (3.8) Thus f (C X ) = f • C X • f -1 is the minimal contractor for f (X). Corrolary 1. If f is box conservative and if I is the identity function, from Proposition 3.3, the minimal contractor for the set X∪f (X) is (I f ) (C X ).
Proof. The minimal contractor consistent with X∪f (X) is

(X ∪ f (X)) ∩ [x] = (X ∩ [x]) ∪ (f (X) ∩ [x]) ((A ∪ B) ∩C = (A ∩ C) ∪ (B ∩ C)) = X ∩ [x] f (X) ∩ [x] ( A ∪ B = A B ) = C X ([x]) f (C X ) ([x]) = (C X f (C X )) ([x]) = ((I f ) (C X )) ([x]) . (3.9) Corrolary 2. If f is box-conservative and if S X = {S in X , S out X } is the minimal separator for X, then the minimal separator for f (X) is f (S X ) = f • S in X • f -1 , f • S out X • f -1 . (3.10)
Proof. This is a direct consequence of Proposition 3.3.

Example 3.3. Let us consider the set X defined in Example 2.7, and the following box-conservative function: From a minimal separator for the set X as defined in Example 2.7, Corollary 2 allows us to obtain minimal separators consistent with the sets X ∪ f (X) and

f (x) = √ 3 -x 2 x 1 . (3.11) -2 -1 0 1 2 3 4 -2 -1 0 1 2 3 4 (a)X ∪ f (X) -2 -1 0 1 2 3 4 -2 -1 0 1 2 3 4 (b)X ∪ f (X) ∪ f 2 (X)
X ∪ f (X) ∪ f 2 (X).
Corresponding subpavings are depicted on Figure 3.2.

Minimal contractors

We consider here an equation of the form y = f (x). The following new theorem defines a way to build minimal contractors for this equation. Theorem 3.1. The minimal contractor consistent with

S = {(x, y) | y = f (x)}, where f : R n → R p , is C S [x] [y] = [x] ∩ f -1 ([y]) [y] ∩ f ([x]) . (3.12) Proof. Define x =i = (x 1 , . . . , x i-1 , x i+1 , . . . , x n ) and [x =i ] = [x 1 ] × • • • × [x i-1 ] × [x i+1 ] × • • • × [x n ]
where the [x i ]'s are the interval components of the box [x]. The optimal contraction for x i is

{x i ∈ [x i ] | ∃x =i ∈ [x =i ] , ∃y ∈ [y] , y = f (x)} (3.13) = x i ∈ [x i ] , ∃x =i ∈ [x =i ] , x ∈ f -1 ([y]) (3.14) = proj x i [x] ∩ f -1 ([y]) . (3.15)
where the projection operator proj x i will be defined in Chapter 4. Now, since for any subset of R n :

A = proj x 1 (A) × • • • × proj xn (A) , (3.16) we get [x] ∩ f -1 ([y]) = proj x 1 [x] ∩ f -1 ([y]) × • • • × proj xn [x] ∩ f -1 ([y]
) .

(3.17) Let us apply the same reasoning with y i . The optimal contraction for y i is

{y i ∈ [y i ] | ∃y =i ∈ [y =i ] , ∃x ∈ [x] , y = f (x)} (3.18) = {y i ∈ [y i ] | ∃y =i ∈ [y =i ] , y ∈ f ([x])} (3.19) = proj y i ([y] ∩ f ([x])) . (3.20) Thus, [y] ∩ f ([x]) = proj y 1 ([y] ∩ f ([x])) × • • • × proj yp ([y] ∩ f ([x])) . (3.21)
As a consequence, C S corresponds to the minimal contractor consistent with S.

Polar contractor

The results given in the previous section are applied here to build a minimal contractor for the polar set defined by:

Π = p = (x, y, ρ, θ) ∈ R 4 | (x, y) = π(ρ, θ) (3.22)
where

π ρ θ = ρ cos θ ρ sin θ (3.23) is the polar function. Define Π 0 = [p 0 ] ∩ Π. (3.24) with [p 0 ] = R + × R + × R + × 0, π 4 . (3.25)
On [p 0 ], we have

x y = π ρ θ ⇔ ρ = √ x 2 + y 2 θ = atan( y x ) (3.26) i.e., ρ θ =π -1 x y = √ x 2 + y 2 atan( y x )
.

(3.27) From Theorem 3.1, the minimal contractor for Π 0 is

C Π 0 [x] × [y] [ρ] × [θ] = [x] × [y] ∩ π ([ρ] × [θ]) [ρ] × [θ] ∩ π -1 ([x] × [y]) (3.28) Figure 3.3 illustrates the contraction of five different boxes [x] × [y] × [ρ] × [θ].
The light gray pies is the initial domain for [ρ] and [θ] while the dark gray pies are the resulting domains obtained after contraction. For instance, the box

[x 1 ] × [y 1 ] × [ρ] × [θ] on the left is contracted into the box [x 1 ] × [y 1 ] × [ρ 1 ] × [θ 1 ] on the right.
The minimal contractor C Π for Π can be deduced from C Π 0 using the following proposition.

Proposition 3.4. Define the following symmetries

σ 1 : (x, y, ρ, θ) → y, x, ρ, π 2 -θ (3.29) σ 2 : (x, y, ρ, θ) → (x, -y, ρ, -θ) (3.30) σ 3 : (x, y, ρ, θ) → (-x, y, ρ, π -θ) (3.31) σ 4 : (x, y, ρ, θ) → (x, y, -ρ, π + θ) (3.32) γ : (x, y, ρ, θ) → (x, y, ρ, θ + 2π) .
(3.33)

A minimal contractor for Π is

C Π = i γ i • (I σ 4 ) • (I σ 3 ) • (I σ 2 ) • (I σ 1 ) (C Π 0 ) . (3.34)
Proof. By composing transformation functions, the initial domain, restricted to

R + ×R + ×R + ×[0, π 4 ] is extended to R 4 .
We have

C Π 0 ∼R + × R + × R + × 0, π 4 ∩ Π (3.35) (I σ 1 ) C Π 0 ∼R + × R + × R + × 0, π 2 ∩ Π (3.36) (I σ 2 ) • (I σ 1 ) C Π 0 ∼R + × R × R + × - π 2 , π 2 ∩ Π (3.37) i∈{3,2,1} (I σ i ) C Π 0 ∼R × R × R + × [-π, π] ∩ Π (3.38) i∈{4,3,2,1} (I σ i ) C Π 0 ∼R × R × R × [-π, π] ∩ Π (3.39) i γ • i∈{4,3,2,1} (I σ i ) C Π 0 ∼Π. (3.40)
The minimality is a consequence of the fact that all transformations are boxconservative.

Building Minimal Separators

The previous section introduced a new way to build a minimal contractor related to constraints which are built by composition of elementary constraints and boxconservative transformations. We now extend these results to build minimal separators.

Minimal separators

Consider the set

X = {x | f (x) ∈ Y} = f -1 (Y) , (3.41)
where f is a function mapping R n into R m . We assume here that Y is a subpaving (a finite union of boxes) which may overlap, i.e.,

Y = i [y i ] . (3.42)
Since Y is a subpaving, its complementary set Y is also a subpaving. For instance, in Figure 3.4, we have [START_REF]Density probability function of random variable x and[END_REF]3]) .

Y = [1, 2] × [1, 3] ∪ [3, 4] × [1, 3] Y = (R × [3, ∞]) ∪ (R × [-∞, 1]) ∪ ([-∞, 1] × [1, 3]) ∪ ([2, 3] × [1, 3]) ∪ ([4, ∞] × [
= f -1 (Y), is S ([x]) = S in , S out ([x]) =      [y]∈Y ∂ x C ([x] , [y]) , [y]∈Y ∂ x C ([x] , [y])      . ( 3

.46)

Proof. To prove that S is minimal, it suffices to prove that the two contractors S in and S out are minimal. Let us first prove that S out is minimal. Define

X = {x | f (x) ∈ Y} = f -1 (Y) . (3.47) For a given box [x], the minimal contractor yields X∩ [x] . Now X∩ [x] = f -1 (Y) ∩ [x] (definition of X) = f -1 ( [y]∈Y [y])∩ [x] (definition of Y) = [y]∈Y f -1 ([y]) ∩ [x] (f -1 (A ∪ B) = f -1 (A) ∪ f -1 (B) ) = [y]∈Y (f -1 ([y])∩ [x]) ((A ∪ B) ∩ [x] = (A ∩ [x]) ∪ (B ∩ [x])) = [y]∈Y f -1 ([y])∩ [x] ( A ∪ B = A ∪ B ) = [y]∈Y ∂ x C ([x] , [y]) (minimality of C [y] ) = [y]∈Y ∂ x C ([x] , [y]) (definition of ) = S out ([x]) .
Let us now prove that S in is minimal. Define

X = {x | f (x) / ∈ Y} = f -1 Y . (3.48)
A reasoning, similar to the first part of the proof, gives us

X∩ [x] = [y]∈Y ∂ x C ([x] , [y]) = S in ([x]
) , (3.49) which terminates the proof.

Remark 3.1. This theorem shows that getting the optimal set inversion mostly depends on the function f and not on Y. As a consequence, it is worthwhile to spend time to compute optimal contractors for the constraint f (x) = y in order to derive minimal separator for f (x) ∈ Y.

Polar separator

We consider the problem of approximating the projection of Π (see Equation (3.22)) along the (ρ, θ) plan defined by:

[ρ] × [θ] ∩ π -1 ([x] × [y]) . (3.50)
Let C Π be the minimal polar contractor consistent with Π. From Theorem 3.2, the minimal separator consistent with this set is

S [x],[y] π ([ρ], [θ]) =      [x 1 ]×[y 1 ]∈[x]×[y] ∂ ρθ C Π ([x 1 ], [y 1 ], [ρ], [θ]) , ∂ ρθ C Π ([x], [y], [ρ], [θ])      , (3.51) where [x] × [y] is a subpaving corresponding to the complement of the box [x] × [y].
Here, [x], [y] are the parameters of the separators and the contractions operate on the box [ρ] × [θ]. The same reasoning applied on ρ and θ, concludes that the minimal separator consistent with the set:

[x] × [y] ∩ π ([ρ] × [θ]) (3.52)
is given by

S [ρ],[θ] π -1 ([x], [y]) =      [ρ 1 ]×[θ 1 ]∈[ρ]×[θ] ∂ xy C Π ([x], [y], [ρ 1 ], [θ 1 ]) , ∂ xy C Π ([x], [y], [ρ], [θ])      . (3.53) -3 -2 -1 1

Application to Localization

Sensors returning range and goniometric measurements of a given landmark, such as sonar or LIDAR, are commonly used in robotics for localization. When measurements are related to some landmarks, the problem can be modeled using the polar constraint [START_REF] Guyonneau | A visibility information for multi-robot localization[END_REF]. Now, we propose examples of localization to illustrate the efficiency of the polar separator.

Static localization

Localization with one landmark

Consider one landmark at a known position m = (m 1 , m 2 ). One robot at position x = (x 1 , x 2 ) is able to measure the distance y 1 to m and the direction y 2 of m in the local reference frame with an interval accuracy. In this example, we take Figure 3.6: Approximation of X using the minimal separator (Left) and by using a classical forward-backward separator (Right). The initial domain is X 0 = [-6, 4] 2 .

[y] = [4,6] × -π 12 , π 4 . The set X of feasible x is:

X = {x ∈ R 2 | ∃y ∈ [y] ,m -x = π(y)} = {x ∈ R 2 | m -x ∈ π([y])} = g m • π ([y]) . (3.54)
where π is defined in Equation (3.23) and

g m (x) = m -x. (3.55)
which is box-conservative. From Corollary 2, the minimal separator S m , associated with landmark m, consistent with X is

S m = g m • S [y] π -1 . (3.56)
where S

[y]

π -1 is the minimal separator defined in Equation (3.53). Figure 3.6 shows the results of the set inversion using our minimal separator (left) and a forward-backward separator (right). Note that the new Theorem 3.2 has to be used in order to be able to get an inner approximation, even with the forward/backward contractor.

Example with several landmarks

Assume now that we have 3 landmarks m i . The measurements [y i ] and the position of the landmarks are given in Table 3.1. The feasible set X is now

X = i g m i • π ([y i ]) .
(3.57) The corresponding separator S is

S = i g m i • S [y i ] π -1 . (3.58)
The paving obtained using S is shown in Figure 3.7. The white dashed box correspond to the result of the first call of the forward/backward separator and is bigger than the one obtained with S.

Underwater localization

We now illustrate the efficiency of the polar separator on a data-set extracted from a mission performed by the actual AUV (Autonomous Underwater Vehicle) VAMA (see Figure 3.8) in the Road-stead of Brest (France, Brittany), July 18, 2015. We focus on the transit phase of the mission where the AUV follows a set of waypoints in order to reach its mission area. During this phase, ranges and bearings

y k = (y k 1 , y k 2 )
between the ship at position m = (m 1 , m 2 ) and the AUV at time k are periodically measured using an ultra short baseline (USBL) and then sent to the robot through an acoustic communication. We assume that no communication delays exists. Moreover, thanks to the pressure sensor, the robot can measure its depth and the localization problem can be easily projected on the two-dimensional horizontal plane.

An estimated reconstitution of the part of the mission made by VAMA that will be used here is depicted on Figure 3.9. The motion of the robot is assumed to be described by the discrete-time state equation:

x(k + 1) = ϕ k (x(k)) = x(k)+ cos ψ(k) -sin ψ(k) sin ψ(k) cos ψ(k) • dt • v(k), (3.59)
where x(k) corresponds to the 2D coordinates of the center of the robot at time k expressed in an absolute inertial frame, ψ(k) is the heading and v(k) is the horizontal speed vector of the robot expressed in its own coordinate system. The speed and heading are measured using a MEMS IMU (Xsens Mti) and are assumed to be known with an uncertainty of ±0.05 m.s -1 and ±1 • . The range is measured with an accuracy of ±1% and the bearing with an accuracy of ±2 • . As a state estimator, we propose to compute the set X k of feasible states x consistent with the last ī measurements. It is given by

X k = x ∈ R 2 | ∀i ∈ {0, ī}, ∃y ∈ [y i ] ,m k-i -ϕ k:k-i (x) = π(y k-i ) (3.60)
where ϕ k 1 :k 2 is the flow defined as follows

ϕ k 1 :k 2 = ϕ k 2 -1:k 2 • • • • • ϕ k 1 :k 1 +1 if k 2 ≥ k 1 ϕ k 1 :k 2 = ϕ -1 k 2 :k 1 otherwise.
(3.61)

We have

X k = i∈{0, ī} x ∈ R 2 | ∃y ∈ [y k-i ] ,m k-i -ϕ k:k-i (x) = π(y) = i∈{0, ī} x ∈ R 2 | m k-i -ϕ k:k-i (x) ∈ π([y k-i ]) = i∈{0, ī} ϕ k-i:k • g m k-i • π ([y k-i ]) . (3.62)
where g m k-i is the function associated with the position of the k -i th landmark, defined in Equation 3.55. The corresponding separator is

S X k = i∈{0, ī} ϕ k-i:k • g m k-i • S [y k-i ] π -1 . (3.63)
For an implementation of the observer, Equation (3.63),

S [y k-i ] π -1
can either be the minimal polar separator or a separator based on a classical forward/backward propagation. For a comparison, let us first apply the separator S X k on the initial box [x 0 ] = [-∞, ∞] 2 for all k without bisection and with ī = 0. As expected, no inner contraction occurs. The diameters of boxes [x k ], which enclose X k , are depicted in Figure 3.10 in blue. It also shows the diameter of the boxes [x k ] that are obtained with a forward/backward separator (green). The peak, between t = 40s and t = 60s, arises when the measured angles are close to π 4 . We conclude that the forward/backward separator is indeed more pessimistic, which is consistent with the fact that the polar contractor is minimal (see also the video in [Desrochers and Jaulin, 2016a]). Figure 3.11 shows an approximation of X k for ī = 5 and t = 55 sec. Since there is an acoustic measurement every 5 sec and since the sampling time is δ = 0.1 sec, a synchronization of all the data (not discussed here) had to be done. The correspondence between k and t is k = 5 • t.

To be robust with respect to outliers, we may allow the observer to relax on the time window of length ī at most q outliers [START_REF] Desrochers | Setmembership approach to the kidnapped robot problem[END_REF]. The solution set X k becomes:

X k = {q} i∈{0, ī} ϕ k-i:k • g m k-i • π ([y k-i ]) (3.64)
where the q-relaxed intersection [Jaulin andWalter, 2002, Carbonnel et al., 2014] is used. The associated separator becomes

S X k = {q} i∈{0, ī} ϕ k-i:k • g m k-i • S [y k-i ] π -1 . (3.65)
The result obtained using this observer for q ∈ {0, 1, 2} are illustrated on a video given in [Desrochers and Jaulin, ]. where x is the unknown state vector, y an exteroceptive measurement vector, u a proprioceptive measurement vector. It is known that for all t, u(t) ∈ [u] (t) and y(t) ∈ [y] (t). The set M is a map, assumed to be known. Equivalently, the observation constraint can then be written as

∃m(t) ∈ M | m(t) = g(x(t), y(t)).
If we take into account the uncertainty on y, we can say that the state vector x(t) is consistent with the interval measurement [y] (t) if

∃m(t) ∈ M, ∃y(t) ∈ [y] (t) | m(t) = g(x(t), y(t)).
The existential quantification ∃m(t) ∈ M highlights the requirement of solving the so-called data association problem which aims at finding which point of the map is associated with the measurement vector y. For us, M is composed with finite number of isolated points and we deal with the initial localization problem on a field of point landmarks that are indistinguishable. All measurements have the same aspect and cannot be associated directly with a particular point of the map.

This problem frequently arises when acoustic sensors are used to detect underwater environmental features [START_REF] Fallon | Relocating underwater features autonomously using sonar-based SLAM[END_REF].

Remark. The observation constraint differs from the shape-based observation constraint given by Z(t) = g(x(t), M) as given in Chapter 2. For the shape observation, the measurement is a shape Z whereas here it is a vector y. This is consistent with the classical feature-based localization. We have chosen here to consider the more classical feature-based localization to illustrate the use of the polar contractor described in this chapter. Later, it will be used in problems involving shape measurements.

Problem

We consider a robot moving on a plane, the motion of which is described by the state equation

ẋ(t) = f (x(t)), u(t)) =    cos(ψ(t)) -sin(ψ(t)) sin(ψ(t)) cos(ψ(t)) • v(t) ω(t)   
The state vector is x = (p x , p y , ψ), where p = (p x , p y ) is its position and ψ is its heading. The input vector is u = (v x , v y , ω), where v = (v x , v y ) is the horizontal speed in the vehicle frame, measured for instance with a Doppler Velocity Log (DVL) (in case of an underwater robot), and ω(.) is the angular velocity measured by gyroscopes.

At some times t i ∈ T, the robot collects the range-bearing vector y

(t i ) = (ρ(t i ), ϕ(t i )) T to a landmark m(t i ) = (m x (t i ), m y (t i ))
T which belongs to the map M, composed of a collection of georeferenced points. This leads to the following

constraint ∃m ∈ M, m = g (x(t i ), y(t i )) (3.67) with g(x(t i ), y(t i )) = p x (t i ) p y (t i ) + ρ(t i ) • cos(ψ(t i ) + ϕ(t i )) sin(ψ(t i ) + ϕ(t i ) . (3.68)
The localization problem is thus described by the following set of equations:

         ẋ(t) = f (x(t), u(t)) m(t i ) = g(x(t i ), y(t i )) m(t i ) ∈ M x(0) ∈ X 0 (3.69)
where the trajectory x(.), and the landmark associated to the measurements taken at time t i both need to be estimated. It can be decomposed into the following constraint network:

                   (i) ẋ(.) = f (x(.), u(.)) (ii) m(t i ) ∈ M (iii) a i = m(t i ) -p(t i ) (iv) α i = ψ(t i ) + ϕ(t i ) (v) a i = ρ • cos(α(t i )) sin(α(t i )) (3.70)
where a contractor can be defined for each constraint. In particular, the contractor C d dt (see Section 2.5.4) will be used to contract the tube [x](.) with respect to (i). The map can be depicted by a subpaving, or by an image, for which a contractor introduced in Section 2.5.2 can be used. For any box [x] ∈ IR 2 , this contractor returns the smallest box which contains all landmarks included in [x]. And the constraint (v) corresponds to the polar one where the minimal contractor, introduced in Section 3.2.3 can be used. The other contractors can be built as explained in Example 2.4.

Test-case

Consider an AUV starting its mission with a huge position uncertainty. This can happen during a dive in deep water [McPhail, 2009] or when, for discretion purposes, a long-range underwater transit phase is required to reach the working area.

For operational reasons, no external positioning system, such as acoustic beacons or USBL, are deployed. We assume that a part of the mission area has been previously mapped during a previous survey and this area is large enough to be reached by the AUV. The corresponding map M describing this area is modeled by a set of 280 point landmarks.

Our robot performs a small mission pattern as depicted in Figure 3.12. It senses its environment using a imaging forward-looking sonar oriented toward the seabed, the scope of which is represented by the blue pie. Every three seconds, it is able to measure the distance and bearing between its pose to some landmarks that range between 10 and 70 meters. The positions of the detected landmarks are depicted by green dots. The 90 red segments represent the measurements. Note that only a small number of mapped landmarks have been detected.

Only the 2D position of the robot and its heading need to be estimated since other state variables (roll, pitch, altitude, depth) are directly measured.

Assumptions. For simplicity, the fact that only landmarks inside the pie are seen by the AUV is not taken into account by the method. Moreover, in an underwater context, the detected landmarks cannot be distinguished from the others, since for instance, two different rocks can have the same aspect and dimension. Moreover, the landmark detection process is sensitive to change in the point of view of the sensor, some landmarks of the map cannot be seen during the survey. Thus, no reliable data associations based on the shape of the landmarks can be assumed. Moreover successive measurements corresponding to the same landmark could be associated by the sonar tracking system. Again, we consider this matching as non reliable and it will not be used for the localization. This example aims at providing a practical illustration of how the constraint propagation methods can be used to find the set of feasible initial states of the vehicle and to solve the data association problem without using a combinatorial algorithm. Once these two issues are solved, any classical localization method, such as an EKF, can be used to get a more accurate estimate of the trajectory. The filter is initialized with a reliable starting point and only updated with observations which are correctly associated.

Results

Given [x(0)], . . . , [x(k max )] a set of boxes enclosing the state trajectory of the robot. All boxes are initialized to R 2 × [-π, π]. The initial heading is assumed to be known with an accuracy of ±10 degrees. Figure 3.15 shows the final trajectory obtained after 12 iterations in less than 1 minute on i7-5600U CPU@2.60GHz.

In Table 3.3, for each iteration, the times needed to contract the whole trajectory are given. This time is constant at each iteration. During the constraint propagation process, the thinner the trajectory, the smaller the number of landmarks contained in [m](t i ), and vice versa. As an indicator, Columns 3 (resp. 4) of Table 3.3 shows the minimal (resp. maximal) number of landmarks contained in [m](t i ) among all measurements. The last column corresponds to the number of correct associations, i.e, when [m](t i ) contains a single landmark. Figure 3.13 shows diameters of boxes along the trajectory for different iterations which illustrates the constraint propagation process. The constraint propagation methods is shown to be powerful in situations involving a huge number of possible data associations. In comparison, existing method [START_REF] Fallon | Relocating underwater features autonomously using sonar-based SLAM[END_REF] often meet difficulty when both the initial position and the data associations are unknown.

Conclusion

Contractor-based techniques are particularly attractive when solving engineering applications, due to the fact that they can handle and propagate uncertainties in a context where the equations of the problem are non-linear and non-convex. Now, the performances of paving methods are extremely sensitive to the accuracy of the contractors that are associated with the equations. One of the equations often met in practice is the polar equation which links Cartesian to polar coordinates.

In this chapter, we proposed two new theorems that could help to build more easily minimal contractors and separators consistent with some specific constraints. These theorems allowed us to build the minimal contractor and a minimal separator for the important polar constraint, which was not done before, to our knowledge.

The efficiency of these new operators and their ability to get an inner and outer approximation of the solution set was illustrated on the problem of the localization of a robot when both goniometric and range measurements are collected. In the next chapter, the polar contractor will be intensively used in the context of the shape-based localization that fits the formalism presented in the introduction of this thesis.

Introduction

In the previous chapter, a dedicated separator, associated with the polar constraint, has been introduced. It allows characterizing a set defined by angles and distances. The efficiency of this separator has been illustrated with localization problems in structured environments composed of landmarks. The variables of these problems are the landmark positions, the pose of the vehicle and the measurements ρ and φ. These variables are all real numbers. The domain which contains feasible values was computed using the polar separator. The objective of this previous chapter was to illustrate how minimal contractors / separators could be built for a constraint that is important in our localization context. Now, the existing catalog of primitive contractors built for dedicated geometrical constraints should be extended to incorporate all constraints that could be useful in a robotics context.

A B

f (p, A) Implementing these specific primitive contractors will allow us to be more efficient by limiting the decompositions (which yields pessimism).

In this chapter, we will assume that both the map M and measurements Y i are subsets of R q , or shapes, where q is the dimension of the map (2 or 3 in practice). We also assume that shapes are perfectly known; i.e. they are thin sets. In the next chapters, this will not be the case since the shape will become uncertain, and we will introduce the notion of interval shape.

This chapter deals with a new type of constraint, called shape-registration. Given two shapes A ⊂ R n , B ⊂ R m and p ∈ R p linked by a constraint L that depends on a parameter p ∈ R p , the shape-registration problem aims at finding the set P of parameter vectors consistent with L. For instance, it may correspond to finding the parameter vectors p ∈ P of a non linear transformation f such that f (p, A) ⊂ B. This is illustrated by Figure 4.1.

These problems are related to the registration of points cloud, generally done by Iterative Closest Points methods [START_REF] Pomerleau | A Review of Point Cloud Registration Algorithms for Mobile Robotics[END_REF], images [Zitova and Flusser, 2003] and surfaces [START_REF] Audette | An algorithmic overview of surface registration techniques for medical imaging[END_REF]]. But we are in a set-membership context.

As the main illustration, the localization problem of a robot with sonar range-finder in an unstructured environment is handled. This problem is considered as difficult due to the fact that the sonar returns a measurement under the form of an impact point inside an emission cone. This specific type of measurement makes the problem partially observable. Moreover, our environment is not represented by geometric features such as segments or disks, but by a shape which cannot be translated into analytical constraints. Now, as shown in Figure 2.5, an unstructured map can be cast into a contractor form which allows us to use contractor/separator algebra.

Here, we propose first to use a separator-based method to perform a reliable simulation necessary to generate realistic data (see [START_REF] Taha | Acumen : An Open-source Testbed for Cyber-Physical Systems Research[END_REF] for a survey on reliable simulation). Then, once these data have been generated, we consider the inverse problem, i.e., the robot localization with large-cone sonar measurements in an unstructured map. This problem has never been considered yet, to our knowledge, at least in an unstructured environment (see e.g., [START_REF] Jaulin | Set inversion for χ-algorithms, with application to guaranteed robot localization[END_REF], Kieffer et al., 1999, Jaulin et al., 2002, Lévêque et al., 1997, Colle and Galerne, 2013] in the case where the map is made with geometrical features). We will also show the link with Minkowski operations and propose separator counterparts for these operations.

Section 4.2 presents the concept of shape registration and shows that it can be interpreted as a set projection problem. The principle is illustrated on a localization problem. Section 4.3 proposes to formulate the Minkowski operations as a specific shape registration problem where the transformations correspond to translations. Section 4.4 illustrates the application of the Minkowski operation to the problem of localization of a robot in an unstructured environment. Section 4.5 concludes the chapter.

Shape Registration as a Projection

Many sets that are defined with quantifiers can be described in terms of projection, inversion, complement, and composition. The projection operator is now defined with its associated separator.

Projection of separators

Definition 4.1. Given X ⊂ R n and Y ⊂ R p . If Z = X × Y, then the projection of a subset Z 1 of Z onto X (with respect to Y) is defined as:

proj X (Z 1 ) = {x ∈ X | ∃y ∈ Y, (x, y) ∈ Z 1 } . (4.1)
Figure 4.2 illustrates the projection of a set along the X dimension.

When a separator is available for Z 1 , it is possible to define a separator consistent with proj X (Z 1 ) based on the following theorem.

proj x (S)([x]) =    y∈[y] ∂ x S in ([x], [y]), y∈[y] ∂ x S out ([x], [y])    (4.2)
Theorem 4.1. Given a separator S associated to Z, we have: Proof. To prove 4.3, we must prove that the outer contractor is consistent with proj x (Z) and the inner with its complement. We have:

proj x (Z) ∼ proj x (S). (4.3) Y X Z 1 proj X (Z 1 )
proj x (Z) = {x | ∃y ∈ [y], (x, y) ∈ Z} (4.4) = y∈[y] {x, (x, y) ∈ Z} (4.5) ∼ y∈[y] ∂ x S out (4.6)
Similarly, since Z ∼ S in , we have:

proj x (Z) = {x | ∃y ∈ [y], (x, y) ∈ Z} (4.7) = {x | ∀y ∈ [y], (x, y) ∈ Z} (4.8) = y∈[y]
x, (x, y) ∈ Z (4.9)

∼ y∈[y] ∂ x S in (4.10)
that ends the proof.

The implementation of the projection of separators is based on Proj-Union and Proj-intersection algorithms described in [Chabert and Jaulin, 2009b]. Roughly, the initial domain [y] is split into a list of small intervals [y i ]. For the outer contraction, each sub-domain [x] × [y i ] is contracted using the outer contractor

∂ x S out ([x], [y i ])
and the hull of the results is returned. For the inner contraction, a sequence of points y i are sampled from [y i ] and ∂ x S in is enforced on [x]×{y i }. The intersection of the results is returned.

In practice, many tricks are used to limit the number of separations. A specific implementation that tried to take advantage of the complementarity of the inner and outer contractions of S has been done. The implementation provided in pyIbex fits the generic separator formalism which can be used with border-based separators.

It does not rely on explicit inner and outer contractors. Note that the proposed projection algorithm assumes that set Z has a non-empty volume in order to define an inner contractor.

Shape registration

In this section, we formulate the shape registration as a set projection. Consider a function

f : R n × R p -→ R m (a, p) -→ f (a, p) (4.11)
For a given p ∈ R p , A ⊂ R n , B ⊂ R m we shall use the following notations:

f (A, p) = {b|∃a ∈ A, b = f (a, p)} f -1 (B) = {z = (a, p) |∃b ∈ B, b = f (a, p)} (4.12)
Consider the set defined by :

P = {p ∈ R p | f (A, p) ⊂ B}. (4.13)
The vector p corresponds to a parameter vector associated to a transformation. A transformation vector p is consistent if after transformation of A it is included inside B. We have

f (A, p) ⊂ B ⇔ ∀a ∈ A, f (a, p) ∈ B ⇔ ¬∃a ∈ A, f (a, p) ∈ B ⇔ ¬∃a ∈ A, (a, p) ∈ f -1 B ⇔ ¬∃a, (a, p) ∈ A × R p ∧ (a, p) ∈ f -1 B . (4.14) As a consequence, P = proj p {(A × R p ) ∩ f -1 B }. (4.15)
Therefore, if we have separators S A , S B for A, B then a separator S P for P can be obtained using the separator algebra [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF]. It is given by

S P = proj p {(S A × S R p ) ∩ f -1 S B }. (4.16)
Combining this separator with a paver, we are able to obtain an inner and outer approximation of P.

Remark 4.1. The same method applies to build the separator consistent with the sets: 4.18) or

P = = {p ∈ R p , f (A, p) ∩ B = ∅} (4.17) = p ∈ R p , f (A, p) ⊂ B (
P = = {p ∈ R p , f (A, p) ∩ B = ∅} (4.19) = P = (4.20)

Example

A robot at position (0, 0) is inside an environment defined by the map

M = {x ∈ R 2 | x 1 < 5 or x 2 < 3}. (4.21)
It emits an ultrasonic sound in the cone with angles π 4 ± π 24 . For a simulation purpose, we want to compute the distance d returned by the sonar. This distance corresponds to the shortest distance inside the emission cone to the complement of the map:

d = inf{d | f (S 1 , d) ∩ M = ∅}, (4.22)
where f (x, d) = d • x is the scaling function and S 1 is the unit cone defined by Recall that for this type of set, a minimal separator has been introduced in Chapter 3. Equivalently, we have

S 1 = {(x, y) | x 2 + y 2 <
d = sup{d | f (S 1 , d) ⊂ M}. (4.24)
To solve our problem, we first characterize the set:

D = {d | f (S 1 , d) ⊂ M}. ( 4

Minkowsky Sum and Difference

The elementary inclusion shape registration defined by Equation (4.13) has been illustrated in the previous section. Our aim is to derive efficient shape registration algorithms, based on Minkowski operations, to deal with the special case where the transformations are translations.

Minkowski operations are used in morphological mathematics to perform dilation or inflation of sets. As it will be shown in Section 4.4, it can also be used for a localization purpose. Efficient algorithms (see [Shoji, 1991]) have been proposed to perform Minkowski operations with sets represented by subpavings. In this section, we show that Minkowski sum and difference can be seen as a shape registration problem. This will allow us to build separators for these Minkowski operations. In the next paragraphs, Minkowski difference and addition are defined with their associated separators. This will be always true for our applications.

Minkowski difference

S B S A = proj p {(S A × S R n ) ∩ f -1 S B }. (4.28)
where f (p, a) = a + p, is a separator for B A.

Proof. Computing the Minkowski difference can be seen as a specific shape registration problem where f (p, a) = a + p , i.e., the transformation corresponds to a translation of vector p. As a consequence, Proof. We have:

B A = proj p {(A × R p ) ∩ f -1 B }. ( 4
A ⊕ B = {p | ∃a ∈ A, ∃b ∈ B, p = a + b} = {p | ∃a ∈ A, ∃b ∈ B, p -a = b} = {p | (p -A) ∩ B = ∅} = {p | (p -A) ∩ B = ∅} = {p | (p + (-A)) ⊂ B}(see Eq 4.27) = B -A.

Localization in an Unstructured Environment

Consider a robot R at position p = (p 1 , p 2 ) in an unstructured environment described by the set M. We assume that the heading θ of R is known with a good accuracy (for instance, by using a compass) and does not need to be estimated.

The robot is equipped with several sonars which return the distance between the robot and the map with respect to their emission cone. This section deals with the localization of the robots using a shape registration method. Several authors have already studied this problem using interval analysis [START_REF] Lévêque | Vehicle Localization from Inaccurate Telemetric Data: A Set-Inversion Approach[END_REF], Meizel et al., 2002, Colle and Galerne, 2013] but in an environment made with segments.

Each sensor emits an acoustic wave in its direction α i which propagates inside a cone of half angle γ corresponding to the aperture of the beam. By measuring the time lag between the emission and the reception of the wave, reflected by the map, an interval

[d i ] = [d - i , d + i
] contains the true distance d i to the nearest obstacle which lies in the scope of the sensor can be obtained. The situation is depicted in Figure 4.6. The area swept by the wave between 0 and d i is free of obstacles whereas the map is hit by the wave at distance d i . Define

S i = {(x, y) | x 2 + y 2 < d - i and atan2(y, x) ∈ [α i -γ, α i + γ]} ∆S i = {(x, y) | x 2 + y 2 ∈ [d i ] and atan2(y, x) ∈ [α i -γ, α i + γ]} (4.33)
The set S i is called the free sector and ∆S i is called the impact pie. These sets are depicted on Figure 4.6. The set of all feasible positions P consistent with [d i ] is

di γ αi (p1, p2) M M S i ∆S i (p1, p2) M M
P(i) = {p ∈ R 2 | (p + S i ) ⊂ M and (p + ∆S i ) ∩ M = ∅} (4.34) = (M S i ) ∩ (M ⊕ -∆S i ).
With several measurements [d i ] the set of all positions consistent with all data is s

P = i (M S i ) ∩ (M ⊕ -∆S i ). (4.35)
Denote by S M , S S i , S ∆S i separators for M, S i , ∆S i . Then a separator for P is

S P = i (S M S S i ) ∩ (S M ⊕ -S ∆S i ). ( 4

.36)

As an illustration, consider the situation described by Figure 4.7 (left), where a robot collects 6 sonar data. The width of the intervals corresponding to the range measurement is ±1m. The first measurement corresponding to i = 1 is painted in green. Figure 4.9 (left) corresponds to an approximation of the set P, obtained using a paver on with the separator S P . It corresponds to the set of positions for the robot that all six impact pies ∆S i intersect the outside of the map M and all six free sectors are inside M. A zoom of the solution set is given in Figure 4.9 (right). The computing time is 127 sec. and 205 boxes have been generated. Note that obtaining an inner approximation of the solution set was not possible using existing approaches that are not based on separators.

Conclusion

Separator-based techniques are particularly attractive when solving engineering applications, due to the fact that they can handle and propagate uncertainties in a context where the equations of the problem are non-linear and non-convex. Now, this chapter extended the class of constraints by considering the shape registration problem. For instance, the classical set inversion problem can be described by the following constraint network where f ∈ F(R n , R m ) and Y is a known set that comes from input data. Now, to deal with unstructured environments, a new constraint which links sets has been introduced

     Variable: x ∈ R n Constraint: f (x) ∈ Y Domain: X ⊂ R n
     Variable : p ∈ R p Constraint : f p (X) ⊂ Y Domain: P ⊂ R p (4.38) where f ∈ F(R n × R p , R m
) is a parametric function and X, Y are two known sets.

As illustrated by the application, this formulation allows us to deal with unstructured data sets (the shapes) that represent a part of the environment (the map) or a measurement (area obstacle-free). Sensor readings are not described by scalars, with an interval as a domain, but by a shape. The goal of the next part is to address the case where the shapes X and Y are variables of our Constraint N etwork. This happens with the SLAM problem inwhich the map is unknown. This issue will be discussed more deeply in Chapter 6.

Part II Shape Valued Constraint Network

domains. They have been used to compute the set of feasible values that satisfy a given number of constraints and have been applied to parameter estimation and robot localization. Now, we will consider the case where variables are subsets of R n . This typically appears in engineering applications when arbitrary shapes (i.e. that cannot be parameterized) are involved or when the model used is uncertain. It corresponds, for instance, to the reconstruction of a three-dimensional object from photos [START_REF] Bottino | Reconstructing 3D Objects from Silhouettes with Unknown Viewpoints: The Case of Planar Orthographic Views[END_REF], mapping an environment from sonar measurements [Leonard and Durrant-Whyte, 1992] or SLAM [Jaulin, 2011]. For us, the uncertainty will be treated under the form of an interval of sets, called Thick Set, or Interval Shape in a robotics context. More precisly, in the context of the shape-based SLAM formalism, presented in the introduction of this manuscript, the observation equation takes the following form Y = f (p, X).

(5.1)

where g : R p × R n → R m , p is an element of R p , and X (resp. Y) is a subset of R n (resp. R m ). All variables, and in particular X and Y, of the equation are uncertain, and need to be estimated. This is the goal of this section which aims at providing tools and methods needed to solve the canonical problem, that we called the shape registration and carving problem, defined in Equation 5.1.

After recalling some notions on lattices, the first section introduces the concept of intervals of sets called Thick Set, or shape interval. This new notion will be used as domains for set-valued variables (or shapes). Section 5.3 extends the notion of separators to characterize thick sets. The canonical problem of set inversion is generalized in Section 5.4 for dealing with uncertain sets, and uncertain functions. Finally, Section 5.5 deals with contractors on thick sets and will be used to refine the approximation of an uncertain set defined by constraints.

Thick Sets

This section recalls more advanced notions on lattices and intervals that will be used to define the thick sets. Behind the provided mathematical background, this section aims to underline the fact that the thick sets theory is mainly an extension of interval analysis. Now intervals are in an infinite dimensional space corresponding to the space of shapes.

Lattices

Lattices. Interval methods can be applied as soon as the set of domains for the variables has a lattice structure [START_REF] Apt | The essence of constraint propagation[END_REF]. A lattice (E, ≤) is a partially ordered set, closed under least upper and greatest lower bounds [START_REF] Davey | Summary for Policymakers[END_REF].

The least upper bound of x and y is called the join and is denoted by x ∨ y. The greatest lower bound is called the meet and is written as x ∧ y. Let us now give three examples.

• The set (R n , ≤) is a lattice with respect to the partial order relation given by x ≤ y ⇔ ∀i ∈ {1, . . . , n} , x i ≤ y i . We have x ∧ y = (x 1 ∧ y 1 , . . . , x n ∧ y n ) and x ∨ y = (x 1 ∨ y 1 , . . . , x n ∨ y n ) where x i ∧ y i = min (x i , y i ) and x i ∨ y i = max (x i , y i ) .

• The set (F, ≤) of the functions that map R to R is a lattice with respect to the partial order relation given by

f ≤ g ⇔ ∀t ∈ R, f (t) ≤ g (t) . We have f ∧ g : t → min {f (t) , g (t)} and f ∨ g : t → max {f (t) , g (t)} .
• The set IR of closed intervals, as introduced by Moore [Moore, 1966], is a complete lattice with respect to the inclusion ⊂. The meet corresponds to the intersection (generally denoted by ∩) and the join corresponds to the interval hull (generally denoted by ). For instance [START_REF]Density probability function of random variable x and[END_REF]4] ∩ [2, ∞] = [2,4] and [1, 4] [8, 9] = [START_REF]Density probability function of random variable x and[END_REF]9] .

(5.2)

A lattice E is complete if for all (finite or infinite) subsets A of E, the least upper bound ∧A and the greatest lower bound ∨A belong to E. When a lattice E is not complete, we can add two elements corresponding to ∧A and ∨A to make it complete. For instance, the set R is not a complete lattice whereas R = R ∪ {-∞, ∞} is. By convention, for the empty set, we set ∧∅ = ∨E and ∨∅ = ∧E.

Intervals. A closed interval (or interval for short) [x] of a complete lattice E is a subset of E which satisfies [x] = {x ∈ E | ∧ [x] ≤ x ≤ ∨[x]} . Both ∅
[x] ∧ [y] = [x -∨ y -, x + ∧ y + ] [x] ∨ [y] = [x -∧ y -, x + ∨ y + ] . (5.
3)

The Remark 5.1. The bracket notation is here used to denote an interval. The brackets can be interpreted as an operator which associates to an unknown variable x, an interval domain [x] which contains it. This operator is used when solving Constraint Satisfaction Problems with intervals [Davis, 1987]. When we define spontaneously an interval [x], then at the same time we define implicitly a variable x enclosed by

[x]. X X + X -

Thick set

Denote by (P(R n ), ⊂) the powerset of R n equipped with the inclusion ⊂ as an order relation. The set P(R n ) is a complete lattice with respect to ⊂. The meet operator corresponds to the intersection and the join to the union. A thick set

[X] of R n is an interval of (P(R n ), ⊂). If [X]
is a thick set of R n , there exist [START_REF] Yao | Towards More Adequate Representation of Uncertainty: From Intervals to Set Intervals, with the Possible Addition of Probabilities and Certainty Degrees[END_REF] two subsets of R n called the lower bound X -and the upper bound

X + such that [X] = [X -, X + ] = {X ∈ P(R n ) | X -⊂ X ⊂ X + } . (5.4)
As illustrated on Figure 5.1, a thick set partitions R n into three zones: the clear zone X -, the penumbra X + \X -and the dark zone R n \X + .

Remark 5.2. For simplicity, the clear zone is also denoted by X in = X -, the dark zone by X out = X + and the penumbra by X ? = X + \ X -.

Remark 5.3. The convention used to draw sub-paving associated with a thick set is the following: red boxes belong to X -(clear zone), orange boxes belong to the penumbra and blue boxes are outside X + (dark zone).

A thick set [X] is a sub-lattice of (P(R n ), ⊂), i.e., if A ∈ [X], B ∈ [X], then A ∩ B ∈ [X] and A ∪ B ∈ [X]
. The set of thick sets of R n will be denoted by

IP(R n ). If for the thick set [X] = [X -, X + ],
we have X -= X + then [X] is said to be thin. It corresponds to a singleton in P(R n ) or equivalently a classical subset of R n .

Remark 5.4. Other possibilities to represent vagueness exist in the literature such as Fuzzy Sets [START_REF] Dubois | Possibility Theory[END_REF] or Fuzzy Interval used in fuzzy control [Galichet andFoulloy, 1995, Boukezzoula et al., 2004]. The closest definition, provided by Equation (5.4), is the one used to define the Rough Set [Pawlak, 1992].

In the Rough-set model, a given set is represented by a pair of ordinary sets called the lower and upper approximations. The approximation space is constructed based on an equivalence relation defined by a set of attributes. On the contrary, the thick set representation, also called i-set in [Jaulin, 2012] or interval-set [START_REF] Yao | Towards More Adequate Representation of Uncertainty: From Intervals to Set Intervals, with the Possible Addition of Probabilities and Certainty Degrees[END_REF], Gervet, 1997], is simply defined as an interval of sets. According to [START_REF] Yao | Comparison of Rough-set and Interval-set Models for Uncertain Reasoning[END_REF], rough-set and thick set models are different extensions of set theory. More precisely with the rough-set theory, the domain of a random variable is described by a membership function which is also uncertain. We can talk about a modelization of the uncertainty of the uncertainty. On the contrary, with the thickset model, the variable is a subset of R n and its domain is handled with the classical interval formalism.

The rough-set model introduces two additional set-theoretic operators based on an equivalence relation on the universe while the thick set model extends the standard set-theoretic operators like intersection, union or complement.

Operations

Some operations which can be applied on thick sets are now defined. Since thick sets are intervals of sets, two types of operators can be considered: extension of classical set operations to elements of thick sets and operations on intervals of sets.

The first class extends binary operators on subsets of R n which correspond to classical operations on vectors of R n (in the Minkowsky sense). If is a binary operator in R n (such as +, -, the multiplication when n = 1, the vector product when n = 3), this extension is:

[A] [B] = {A B, A ∈ [A] , B ∈ [B]} . (5.5)
The second class is composed of binary operators such as ∈ {∪, ∩, ×, \, . . . } where × is the Cartesian product, \ is the restriction (or trim) operator, for subsets of R n that do not correspond to any extension of operators in R n .

Given a collection {X i } i∈I of subsets of R n , we denote by {X i, i ∈ I} the smallest thick set which encloses all X i , i ∈ I. We have:

{X i, i ∈ I} =   i∈I X i , i∈I X i   . (5.6)
It is possible to extend the operators to thick sets as follows:

[A] [B] = {C, ∃A ∈ [A] , ∃B ∈ [B] , C = A B} . (5.7)
From the monotony of the following operators, we have:

[A] ∩ [B] = A -∩ B -, A + ∩ B + (5.8) [A] ∪ [B] = A -∪ B -, A + ∪ B + (5.9) [A] × [B] = A -× B -, A + × B + (5.10) [A] \ [B] = A -\ B + , A + \ B - (5.11) [A] ⊕ [B] = A -⊕ B -, A + ⊕ B + (5.12) If f is a function from R n to R m the image of the thick set [A] ⊂ P(R n ) by f is defined as: f ([A]) = f (A), A ∈ A -, A + = f (A -), f (A + ) .
(5.13)

For example, the following illustrates the two types of intersections that can be defined:

[X] ∩ [Y] = {Z ∈ P(R n ) | Z = X ∩ Y, X ∈ [X] , Y ∈ [Y]} [X] [Y] = {Z ∈ P(R n ) | Z ∈ [X] , Z ∈ [Y]}. (5.14)
The first ∩ corresponds to the extension to IP(R n ) of the intersection in P(R n ) whereas the second corresponds to the intersection in IP(R n ). Therefore .15) As shown in [Jaulin, 2012], we have:

X ∈ [X] , Y ∈ [Y] ⇒ X ∩ Y ∈ [X] ∩ [Y] Z ∈ [X] , Z ∈ [Y] ⇒ Z ∈ [X] [Y] . ( 5
[X] ∩ [Y] = [X -∩ Y -, X + ∩ Y + ] [X] [Y] = [X -∪ Y -, X + ∩ Y + ] .
(5.16)

The same type of relations applies to the union

[X] [Y] = [X -∩ Y -, X + ∪ Y + ] [X] ∪ [Y] = [X -∪ Y -, X + ∪ Y + ] (5.17)
Comparisons between ∩, ,∪ and are depicted in Figure 5.2.

Thick Separators

To characterize a thin set using a paver, we may use a separator inside a paver.

As for operations on sets, separators can be immediately generalized to thick sets.

x 2

x 1 -3 -2 -1 1 2 3 -3 -2 -1 1 2 3 X x 2 x -3 -2 -1 1 2 3 -3 -2 -1 1 2 3 Y X Y x 2 x 1 -3 -2 -1 1 2 3 -3 -2 -1 1 2 3 X ∩ Y x 2 x -3 -2 -1 1 2 3 -3 -2 -1 1 2 3 X Y X ∩ Y X Y x 2 x 1 -3 -2 -1 1 2 3 -3 -2 -1 1 2 3 X ∪ Y x 2 x -3 -2 -1 1 2 3 -3 -2 -1 1 2 3 X Y X ∪ Y X Y Figure 5
.2: Illustration of the two different types of operators on thick sets. Now, the penumbra as a non-zero volume and, for efficiency reasons, it is important to avoid any accumulation of the paving deep inside the penumbra.

Definition

A thick separator [S] for the thick set [X] is an extension of the concept of separator to thick sets. As illustrated by Figure 5.3, a thick separator is a 3-uple of contractors S in , S ? , S out such that, for all

[x] ∈ IR n S in ([x]) ∩ X in = [x] ∩ X in S ? ([x]) ∩ X ? = [x] ∩ X ? S out ([x]) ∩ X out = [x] ∩ X out (5.18)
As for separators, which are built in practice from a pair of two complementary contractors, thick separators can be built from a pair of classical separators. Given separators S -= S - in , S - out consistent with X -and S + = S + in , S + out consistent with X + , the thick separator [S] consistent with [X] is:

[S] = S - out , S - in ∪ S + out , S + in (5.19)
Proof. This is a direct consequence of definitions of S -, and S + . We have

S - out ∼ X -= X in (5.20) S + in ∼ X + = X out (5.21) S - in ∪ S + out ∼ X -∩ X + = X ?
(5.22)

The following example illustrates how a thick set can be defined from a grayscale image and characterized using a thick separator.

Example 5.1. We consider the problem of characterizing a thick set [X] defined from the ternary images depicted in Figure 5.4a. We assume that white pixels belong to X -, white and gray pixels belong to X + and black pixels do not belong to [X]. Using the image separator introduced in Section 2.5.3 on page 44, it is possible to build a thick separator [S] for [X]. The result of the characterization of [X] is depicted in Figure 5.4b.

Algebra

Algebra for thick separators is defined in the same manner than what has been done for contractors [Chabert and Jaulin, 2009b] or for separators [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF]. The main motivation is to provide methods to compute with

X + X - S in ([x]) S out ([x]) S ? ([x])
Figure 5.3: Thick separator. S out is consistent with X + , S in with X -and S ? with X ? = X + \ X -. thick sets. We only focus on the intersection of thick separators but a similar reasoning can be done to define other operations, introduced in Section 5.1, such as the union, restriction, squared intersection, exclusive-or, or relaxed intersection [START_REF] Desrochers | Relaxed intersection of thick sets[END_REF].

Intersection Consider two thick sets [X] = [X -, X + ] and [Y] = [Y -, Y + ] with their associ- ated thick separators [S X ] = S in X , S ? X , S out X and [S Y ] = S in Y , S ? Y , S out Y .
A thick separator for the thick set

[Z] = Z -, Z + = [X] ∩ [Y] = X -∩ Y -, X + ∩ Y + . (5.23) is [S Z ] = S in Z , S ? Z , S out Z = S in X ∩ S in Y , S ? X ∩ S in Y S ? X ∩ S ? Y S in X ∩ S ? Y , S out X S out Y .
(5.24)

Proof. We have

Z in = Z - = X -∩ Y - = X in ∩ Y in Z ? = Z + \Z -= X + ∩ Y + \(X -∩ Y -) = X + ∩ Y + ∩ X -∩ Y - = X + ∩ Y + ∩ X -∪ Y - = X in ∪ X ? ∩ Y in ∪ Y ? ∩ X out ∪ X ? ∪ Y out ∪ Y ? = X ? ∩ Y in ∪ X ? ∩ Y ? ∪ X in ∩ Y ? Z out = Z + = X + ∩ Y + = X + ∪ Y + = X out ∪ Y out .
From the separator algebra, we get that a contractor for

Z in is S in Z = S in X ∩ S in Y , a contractor for Z out is S out Z = S out X S out
Y and a contractor for Z ? is

S ? Z = S ? X ∩ S in Y S ? X ∩ S ? Y S in X ∩ S ? Y .
(5.25)

Remark 5.5. In the case where [S X ] = S - X , S + X and [S Y ] = S - Y , S + Y , the separator consistent with Equation (5.2) is

[S Z ] = S - X ∩ S - Y , S + X ∩ S + Y .
(5.26)

The proof is a direct consequence of Equation (5.9). We have

S - X ∩ S - Y ∼ X -∩ Y -and S + X ∩ S + Y ∼ X + ∩ Y + .

Uncertain Set Inversion

The canonical problem of the set inversion X = f -1 (Y) is well formalized and many algorithms, such as the ones presented in Chapters 2 and 3, are commonly used to solve it. In this section, we extend this problem to the thick set case where both Y and f are uncertain. A first theorem is given in order to provide a generic formulation of the uncertain set inversion problem.

However, in the general case, this problem is hard to solve. In order to provide explicit methods for solving this problem, assumptions on the nature of f have to be made. When f depends on an uncertain parameter vector p ∈ R p , a method for dealing with the set inversion problem is proposed in this section. It allows handling uncertain nonlinear transformations such as rotations of an angle which belongs to a given interval. The introduced method mainly relies on the projection of separator, presented in Section 4.2.1, which can be slow when dealing with a high-dimensional parameter space. On the contrary, when f is assumed to belong to an interval of functions, a dedicated and more efficient method is introduced in Appendix A. Note also that this appendix, published in [START_REF] Desrochers | Thick set inversion[END_REF], contains additional examples which have motivated the need for developing specific algorithms.

Problem statement

Given a function f : R n → R m and a set Y ⊂ R m , set inversion aims at bracketing from inside and outside the set .27) This formalism has been used for more than 20 years with interval methods to solve problems in bounded-error parameter estimation [START_REF] Vehi | Using Interval Methods for Control Systems Design in the Parameter Space[END_REF], robot localization [START_REF] Meizel | Initial localization by set inversion[END_REF], Colle and Galerne, 2013, Drevelle and Bonnifait, 2013] and robust control [START_REF] Malti | Stability and resonance conditions of elementary fractional transfer functions[END_REF], Herrero et al., 2005, Jaulin and Walter, 1993a]. Most interval algorithms for set-inversion alternate some interval tests or contractions [Chabert and Jaulin, 2009b] to certify that a box (i.e., a Cartesian product of intervals) is inside or outside the solution set X and bisect the boxes for which no conclusion can be reached.

X = f -1 (Y). ( 5
In the case where both f and Y are uncertain, a relaxation of the resulting uncertain constraints can be performed by adding quantifiers as made in [Goldsztejn,

Theorem 5.1. Given F ⊂ F(R n , R m ) and the thick set [Y] = [Y -, Y + ],
the smallest thick set which encloses all sets X such that

∃f ∈ F, ∃Y ∈ [Y] | X = f -1 (Y). (5.31) is the thick set [X] = [X -, X + ]
where

X -= f ∈F f -1 (Y -), X + = f ∈F f -1 (Y + ).
(5.32)

Proof. Denote by {X i } i∈I the set of all solutions of Equation ( 5.31). The smallest thick set [X] containing {X i } i∈I is the thick set

[X] =   i∈I {X i }, i∈I {X i }   . (5.33) Now, for each X i ∈ {X i } i∈I , ∃f i ∈ F, ∃Y i ∈ [Y] | X i = f -1 i (Y i ).
(5.34)

Thus, [X] is given by

  f ∈F Y∈[Y] f -1 (Y) , f ∈F Y∈[Y] f -1 (Y)   .
(5.35)

Now, since f -1 is monotonic with respect to the inclusion ⊂, we get

[X] =   f ∈F f -1 (Y -), f ∈F f -1 (Y + )   .
(5.36)

Remark 5.6. Theorem 5.1 provides the exact formulation of the thick set inversion problem and defines the corresponding two sets we want to compute, i.e., the two sets X -and X + as defined by Equation (5.32). The main difficulty is to get an inner approximation of the penumbra X + \X -.

Existing interval methods can still be used to deal with this type of uncertainty, but they accumulate inside the penumbra. This accumulation makes these methods inefficient, since they spend most of the computation time to test tiny boxes that are deep inside the penumbra, without proving that they are effectively in the penumbra. Consequently, specific algorithms have to be developed.

As it is illustrated by the next example, when F contains a single function f from R n to R m , from Equation (5.13), the smallest thick set solution is [X] = f -1 (Y -), f -1 (Y + ) , (5.37) where classical set inversion algorithms can be used. .38) The result of the set inversion

f (x) = x 1 sin(x 1 + x 2 ) x 2 cos(x 2 -x 1 ) . ( 5
[X] = f -1 ([Y]
) is given by Figure 5.6. The thick separator, used by the paver, is built from a pair of separators consistent with the two classical set inversion problems:

f -1 (Y -) and f -1 (Y + ).
When F is a set of functions that depend on a parameter, elimination of quantified parameters can be done using projection algorithms (see Section 4.2.1 on page 81). This case is discussed in the next section. However, it requires additional bisections along the parameter space that become costly when dealing with high dimensional parameters.

On the contrary, when F can be described by an interval of functions, called a Thick Function, a new algorithm has been developed to compute the solution set without bisecting inside the parameter space. This method is discussed in Appendix A.

Set inversion with parametric function.

We consider here the case where f depends on a parameter vector p ∈

[p] ⊂ R p .
The function should thus be written as f (x, p), but for simplicity of notation, we will often write f p (x). For simplicity, the set inversion problem will be written as:

[X] = f [p] ([Y]) (5.39)
From Theorem 5.1, the bounds of the smallest thick set [X -, X + ] are defined by:

X -= p∈[p] f -1 p (Y -) = x|∀p ∈ [p], f p (x) ∈ Y - (5.40) X + = p∈[p] f -1 p (Y + ) = x|∃p ∈ [p], f p (x) ∈ Y + .
(5.41)

Using the notion of projection introduced in Section 4.2, it is possible to build separators S -, consistent with X -, and S + consistent with X + . Then the thick separator consistent with [X] can be defined and used inside a paver to characterize the solution set.

The following example illustrates the set inversion with the uncertain parametric rotation of a box of R 2 . It compares a sampling method used to guess the shape of the solution set, the classical set inversion algorithm which accumulates in the penumbra and the proposed set inversion using separators.

Example 5.4. Let R θ be a rotation of an angle θ ∈ [ π 6 , π 3 ] and [y] = [0, 3] × [4,5]. The solution of the set inversion problem:

[X] = R [θ] ([y]) (5.42)
is depicted in Figure 5.7 with and without the use of projections. The classical interval methods, on Figure 5.7b, accumulate in the penumbra and add a lot of pessimism. Both the red and blue areas are smaller than the true one, drawn in Figure 5.7c.

Shape Registration and Carving

Previous sections have defined the notions of thick sets and algorithms for solving the uncertain set inversion problem. As defined in the introduction, when a set is used as a variable of our constraint network we call it a shape. It belongs to a domain that we called an interval shape. In this section, the problem of the shape registration, presented in Chapter 4, is extended to the case where the shapes are not thin anymore and are uncertain. This problem is defined by the following .7: The result of the set inversion using a sampling method along the parameter space is depicted in Figure 5.7a. For each θ ∈ [θ], the image of the red box by the rotation of an angle -θ is drawn using a transparency effect. It makes it possible to guess sets X -, X + . Figure 5.7b: the result with classical methods where the penumbra is not characterized. Figure 5.7c: the result of the set inversion using projection algorithms. shape-valued constraint network :

         f (A) = B A ∈ [A] B ∈ [B] f ∈ F (5.43)
Two sub-problems have to be considered in order to contract all domains of the involved variables. The registration problem which aims at finding the smallest set of functions, included in F, compatible with the constraint

f (A) = B, A ∈ [A] and B ∈ [B]
. And, the carving problem which intends at refining the thick set approximation of [A] and [B] for a given set of functions F.

For that purpose, the notion of contractor on thick sets is introduced in this section with some elementary contractors dedicated to Equation 5.43. Then the shape registration problem, presented in Section 4.2.2, will be extended to the thick set case. Finally, this section ends with an example which illustrates the method.

Contractor on thick sets

Since Thick Sets are intervals of sets, it is possible to define contractors for them in the same way as introduced in Section 2.4.1.1.

A contractor on a thick set is an operator from IP(R n ) to IP(R n ) such as :

C([X]) ⊂ [X] (contractance) [X] ⊂ [Y] ⇒ C([X]) ⊂ C([Y]) (monotonicity) (5.44)
The effect of a contraction on a thick set is illustrated on Figure 5.8. A contraction of the lower bound corresponds to an inflation of the under approximation (Fig. 5.8c).

Similarly, decreasing the upper bound of [X] contracts the upper approximation (Fig. 5.8b). Since thick sets are interval of P(R n ), it is only possible to move its bounds and have an outer contraction.

X X + IP(R n ) X - (a) initial domain X X + IP(R n ) X - (b) lower contraction X X + IP(R n ) X - (c) upper contraction Figure 5
.8: Illustration of contractor applied on a thick set. Recall that a thick set is an interval of P(R n ). The thin set X, delimited by the dashed line, is a singleton in this space.

Shape carving

The specific primitive contractors on thick sets are now presented. They will allow to contract shape domains of Equation (5.43).

Constraint A ⊂ B

Since the equality between two sets can be decomposed into two inclusion constraints, we consider the following proposition:

Proposition 5.1. Given A ∈ [A], B ∈ [B]
, the minimal contractor associated with the constraint A ⊂ B is:

C ⊂ [A] [B] = [A] ([B] \ ([∅, R n ])) [B] ([A] ∪ ([∅, R n ])) (5.45)
or equivalently:

C ⊂ [A] [B] = [A -, A + ∩ B + ] [B -∪ A -, B + ] (5.46)
Proof. The proof can be found in [Jaulin, 2012].

This contractor is illustrated in Figure 5.9.

B A (a) [A] and [B]

A (b) [A] ([B] \ [∅, R n ]) B (c) [B] ([A] ∪ [∅, R n ])
Figure 5.9: Contraction with respect to the constraint A ⊂ B.

Constraint f (A) = B

We have the following theorem:

Proposition 5.2. Given A ∈ [A], B ∈ [B]
, the minimal contractor associated with the constraint f (A) = B where f is bijective is :

C f [A] [B] = [A] f -1 ([B]) [B] f ([A]) = [A -∪ f -1 (B -), A + ∩ f -1 (B + )] [B -∪ f (A -), B + ∩ f (A + )] (5.47)
Proof. The proof can be found in [Jaulin, 2012].

Finally, using the thick set inversion algorithm, introduced in Section 5.4 , the previous contractor can be extended to the case where

f is uncertain, i.e., f ∈ F ⊂ F(R n , R m ). Proposition 5.3. Given A ∈ [A], B ∈ [B], f ∈ F the minimal contractor associated with the constraint f (A) ⊂ B is : C F [A] [B] = [A] F -1 ([B]) [B] F([A]) (5.48)
where

F -1 ([B]
) is the result of the set inversion as formalized in Section 5.4.

Proof. The proof is a direct consequence of Theorem 5.1 and Proposition 5.1.

In this section, only contractors needed to define the shape carving problem have been considered. Other contractors, such as those associated to the constraint A ∩ B = ∅, can be found in [Jaulin, 2012], but are not needed for solving the Shape-based SLAM presented in the next chapter.

Shape registration and carving

The shape registration problem, introduced in Chapter 4.2, is now extended to the uncertain case by taking into account that the shapes involved in the transformation are not thin anymore. Given two shapes linked by a function f that depends on a parameter vector p ∈ R p , the following propositions make it possible to contract the domain for the variables (parameters ans shapes).

Constraint f p (A) ⊂ B

The following proposition allows to build a contractor, and even a separator, able to characterize an outer approximation of the set of feasible parameters associated to the constraint f p (A) ⊂ B.

Proposition 5.4. Given f : R n × R p → R m , A ∈ [A] and B ∈ [B],
we have: .54) where . denotes the convex hull.

{p ∈ R p , f (A, p) ⊂ B} ⊂ p ∈ R p , f (A -, p) ⊂ B + (5.49) Proof. For A ∈ [A], and B ∈ [B], {p, f (A, p) ⊂ B} = {p, ∀a ∈ A, f (a, p) ∈ B} (5.50) B⊂B + ⊂ p, ∀a ∈ A, f (a, p) ∈ B + (5.51) A -⊂A ⊂ p, ∀a ∈ A -, f (a, p) ∈ B + (5.52) ⊂ p, f (A -, p) ⊂ B + (5.53) As a consequence, given p ∈ [p] the contractor C f ⊂ associated to the constraint f (A, p) ⊂ B, A ∈ [A], B ∈ [B] is : C f ⊂ ([p], [A] , [B]) = [p] ∩ proj p (A -× R p ) ∩ f -1 (B + ) . ( 5
Remark 5.7. In the case where f p is bijective, the constraint f p (A) = B can be decomposed into

f p (A) ⊂ B f -1 p (B) ⊂ A (5.55)
for which Proposition 5.1 can be used twice. This leads to the following contractor:

C f = ([p], [A] , [B]) = C f ⊂ ([p], [A] , [B]) ∩ C f -1 ⊂ ([p], [B] , [A]).
(5.56)

5.5.3.2 Constraint f p (A) ∩ B = ∅
The disjoint constraint is also extended to the uncertain case. It has been used in the application of Chapter 4 to handle the impact area of sonar sensors.

Proposition 5.5.

Given f : R n × R p → R m , A ∈ [A] and B ∈ [B],
we have:

{p ∈ R p , f (A, p) ∩ B = ∅} ⊂ p ∈ R p , f (A + , p) ∩ B + = ∅ (5.57) Proof. Given A ∈ [A]
, and B ∈ [B], we have:

{p ∈ R p , f (A, p) ∩ B = ∅} = {p, ∃a ∈ A, f (a, p) ∈ B} (5.58) B⊂B + ⊂ p, ∃a ∈ A, f (a, p) ∈ B + (5.59) A⊂A + ⊂ p, ∃a ∈ A + , f (a, p) ∈ B + (5.60)
As a consequence, the contractor

C f = associated to the constraint f (A, p) ∩ B = ∅, A ∈ [A], B ∈ [B] is : C f = ([p]) = [p] ∩ proj p (A + × R p ) ∩ f -1 (B + ) .
(5.61)

Application

As a preamble of applications presented in Chapter 6, the problem of Shape Registration And Carving (SRAC) is now considered. This situation corresponds, for instance, to a situation where a shape has partially been seen from two different positions. The transformation to move, from one point of view to another, depends on some parameters to be estimated. The goal of the shape registration and carving is to reduce the set of feasible parameters (registration) and contract observed shapes (carving).

Given two shapes A ∈ [A], B ∈ [B] linked by a bijective function f p that depends on a parameter vector p ∈ [p], the SRAC problem is described by the following SVCN: The contraction procedure is given by Algorithm 5.1. The transformation between these two views is defined by: 2,2]. Figure 5.11 shows the result of the set inversion using a separator.

         f p (A) = B A ∈ [A] B ∈ [B] p ∈ [p]
Algorithm 5.1 SRAC(in: f , inout: [p], [A], [B]) 1: [p] := C f = ([p], [A] , [B]) (see Proposition 5.4) 2: [A] [B] := C f [p] [A] [B] ( see
f p (x) = p 1 0 + cos p 2 -sin p 2 sin p 2 cos p 2 • x (5.63) The initial domain for [p] is [0, 4] × [-
Remark 5.8. From Equation (5.4), only an outer approximation of the set of feasible parameters can be computed. However, due to the huge pessimism introduced by f , bisections are needed. The use of the separator, which provides an inner approximation, allows us to save computation time. It avoids useless bisections inside the solution set.

After the contraction, we obtained [p] = [2.65113, 2.90901] × [-1.14117, -0.95238] which contains the theoretical value p * = (2.8, -1) used to generate the images. This result is used to contract the thick set approximation. Figure 5.10 shows these sets after the contraction.

-6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 (a) [A] -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 (b) [B] -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 (c) [A] f -1 [p] ([B]) -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 (d) [B] f [p] ([A])

Conclusion

This chapter introduced the notion of interval of sets, called thick set. The thick set object is used as domains associated with set-valued variables (or shapes) that appear in our shaped-valued constraint network.

As an elementary constraint, the uncertain set inversion problem has been formalized. When the function can be parametric, a generic method that relies on the projection operator has been provided. This method can handle any non-linear function, but at a high computational cost.

Then, it has been used to extend the shape registration problem, introduced in the previous chapter, to the case where the shapes are also uncertain and need to be estimated. For that purpose, the carving problem, that aims at contracting their domains, has been introduced. Based on the notion of contractors on thick sets, which has been defined, a specific contractor applied to this problem was introduced. This contractor will be intensively used as a basic operator by applications presented in Chapter 6. However, computing the set inversion with an uncertain function can be resource demanding and specific methods need to be developed in order to avoid bisections as much as possible. This is the goal of Appendix A, which introduces algorithms able to deal with the set inversion problem for which the uncertain function is described by an interval of functions. This chapter concludes the presentation of the theoretical tools that have been developed. These tools will be used in a robotics context in order to solve the shape-based SLAM problem. Note that Appendix B provides another application of interval shapes. It focuses on a method able to deal with the guaranteed characterization of the part of the space that has been explored by a robot when the dimension of the visible space at time t may be smaller than that of the workspace.

Introduction

The previous chapters have introduced all tools and algorithms needed to deal with the targeted applications of this manuscript. The main goal of this thesis is to provide an interval-based method to handle the SLAM problem in an unstructured environment (or Shape-based SLAM). In Chapter 4, the localization problem was first solved, but the map was assumed to be known without uncertainty. To tackle this issue, Chapter 5 introduced the notion of thick set, and defined a method to deal with the shape registration and carving problem. For the particular case of the translation, Appendix A introduces a method for solving the thick set inversion problem.

In the first section, a generic formulation of the shape-based SLAM is presented with a focus on how shapes are built from measurements. With respect to how a sensor reading is translated into a shape, two categories of interval shapes can be defined. For the first one, the shape provides an inner and an outer approximation of a part of the map. This allows us to formulate, in Section 6.3, the shape-based SLAM with inter-temporal constraints. As a direct application, we propose an efficient algorithm able to deal with the bathymetric SLAM problem. On the contrary, using sensors that only provide an inner approximation of the map, such as those used in Chapter 4, a different formulation will be needed. This formulation, called the Dig-SLAM problem, is introduced in Section 6.4. Section 6.5 concludes this chapter.

Problem Statement

In this section, the general formalism of the SLAM problem where both the map and measurements are shapes is given. Based on the thick-set contractors introduced in Chapter 5, a generic resolution method is then presented. The last part of this section focuses on the definition and the representation of physical features of the environment by shapes.

Formalism

A generic pose-based simultaneous localization and mapping problem can be described by:

     ẋ(t) = f (x(t), u(t)) (evolution equation) Z(t) = g x(t) (M) (map equation) x(0) = x 0 (initial position) (6.1)
where t ∈ [t] ⊂ T is the time, x ∈ R n is the state vector, u ∈ R m is the input vector (in general associated with proprioceptive sensors), f : R n × R m -→ R n is the evolution function and x 0 is the initial state at t = 0. The set M is a subset of R q , where q is the dimension of the map (two or three in practice). The shape Z(t) ⊂ R r is an exteroceptive measurement collected by the robot (for instance by a camera, sonar, telemeter, . . . ), and expressed in the robot frame. The function g : R n × R q -→ R r is the map function, which links the observed local shape Z(t) to the global map M with respect to the state of the vehicle x(t). For a given state vector x(t), g is assumed to be bijective i.e

∀m ∈ R q , g x(t) • g -1 x(t) (m) = m (6.2)
Note that, most geometrical transformations, such as rotations, translations or dilations, are bijective functions that can be parametric.

In order to fit the set-membership formalism used in this thesis, the SLAM problem is handled by the following hybrid constraint network

                     Variables: x(•), M, Z(t i ) Constraints: 1. ẋ(•) = f (x (•) , u (•)) 2. Z(t i ) = g x(t i ) (M) Domains: [x](•), [M], [Z] (t i ) (6.3)
where measurements are taken at time t i ∈ T, and inputs of the system u(t) are assumed to be known and to belong to a bounded tube [u](.). Recall that the set-valued variables Z(t i ) are shapes. Their domains, defined by thick sets, are called thick shapes.

Generic resolution

Solving the SLAM problem implies generally to estimate both the trajectory and the map. The CN of Equation ( 6.3) is composed of two constraints for which contractors C d dt , C SRAC have been defined.

The first contractor C d dt , defined in Section 2.5.4 on page 47, contracts the tube [x](t) with respect to the evolution equation ẋ(t) = f (x(t), u(t)) and the initial condition x 0 .

Then, the Shape Registration And Carving contractor C SRAC , introduced in Section 5.5.4 on page 114, is used to contract[x](t), the map, and the shape interval [Z(t)]. Algorithm 6.1 provides a solution for solving the SLAM problem.

Note that this formulation is generic and, with respect to the nature of the measurements, needs to be adapted. Example of data that can be described by shapes are now given.

Algorithm 6.1 Explicit shape-based SLAM

1: [M] = [∅, R q ] 2: [x](0) = x 0 3: repeat 4: [x](.) := C d dt ([x](.)) 5:
for t i in T do 6:

C SRAC ([x](t i ), [Z] (t i ), [M]) 7:
end for 8: until no more contraction on [x](.)

P (R q ) Z = [Z -, Z + ]
(a) Double-sided observed shape

P (R q ) Z = [Z - , R q ]
(b) One-sided observed shape

Shape-based measurements

As it has been said in the introduction of this manuscript, the map representation, based on subsets of R n , fits into the location-based map framework. Each position of the world is assumed to belong, or not, to a set (or shape) associated with a physical feature of the environment.

This shape is perceived by the embedded sensors of the robot. From raw data, we suppose that an algorithm is able to characterize an initial domain depicted by a thick shape. This initial domain has to take into account the uncertainty of the data acquisition process and the scope of the sensor. Now, with respect to the map and the geometry of the sensor, two types of thick shapes need to be considered. The double-sided observed shape, for which we measure both an inner and outer approximation, and the one-sided observed shape, where only an inner approximation can be measured. These shapes are described in Figure 6.1

Double-sided observable shape

A double-sided observable shape can be defined when both M and M can be measured. More formally, its interval shape approximation [Z] = [Z -, Z + ] satisfies: (6.4) This is illustrated by Figure 6.2. The map M is defined by the dashed area and the scope of the sensor, denoted by V, is depicted by the dotted black rectangle. We assume that from sensor's data, it is possible to get an inner approximation Z - (in red) of M and Z + (in blue) of M that ranges in the scope of the sensor. Since robustly detecting the transitions between M and M can be a challenging task, these transitions are enclosed inside the penumbra. The area, outside the scope of the sensor that has not been observed also belongs to the penumbra.

Z -⊂ M Z + ⊂ M
R 0 M V R 0 Z + Z + Z - Figure 6
.2: Double-sided observable shape. On the left-hand side, the map M is depicted by the dashed black lines. The black rectangle delimits the scope of the sensor which is the area that can be sensed by the vehicle. On the right-hand side, the thick shape [Z] which encloses the map i.e Z -⊂ M ⊂ Z + . The transition between M and M belongs to the penumbra as well as all spaces that have not been observed, i.e. outside the scope of the sensor.

For instance, an airborne lidar measures distances between its carrier and the ground. These measurements, merged together into a local Digital Elevation Map (DEM), split the space into three areas:

• The shape Z -that represents the space between the plane and the ground which is free of obstacles.

• The shape Z + that corresponds to the parts of the space below the ground which are assumed to be occupied.

• The shape Z + that is composed of unobserved areas and the uncertainty of the bottom detection at the interface air/ground.

In this case, the DEM is represented as a shape of R 3 .

In Section 6.3, a new algorithm based on inter-temporal constraints will be introduced. It will take advantage of the fact that both M and M can be observed.

One-sided observed shapes

With partially observed shapes Z, the interval shape approximation only provides an inner view of the map, i.e., Z ∈ [Z -, R q ] . It mainly happens when dealing with occupancy maps that separate free from occupied area sensed with ranging sensors. From distance readings, it is only possible to measure the free area between the vehicle and the nearest obstacle that lies within the scope of the sensor. Consequently, nothing can be inferred from the nature of the map outside this area. Figure 6.3 illustrates this concept with a pair of ultrasonic sensors. The regions of the working space that have not been sounded belong to the penumbra. Note that, the areas of the impact point, denoted ∆Z, have to be used for localization purpose.

This situation corresponds to the application introduced in Chapter 4. The Section 6.4 will extend this application to the SLAM problem. Two sonar measurements are shown. For each of them, the impact point lies in a sector ∆Z. Since its position cannot be determined precisely, it makes this type of measurement partially observable. The thick shape associated to the measurements is drawn on the right side. The sounded part of the space, in red, is free from obstacles and belongs to Z -. Since nothing can be said about the nature of the environment sensed at a distance greater than the one measured, all the remaining space belongs to the penumbra.

Shape-based SLAM

The previous section introduced a general method to handle the SLAM problem with shapes. The proposed algorithm updates sequentially the trajectory, the map and the measurement shapes. However, building and updating large maps can be difficult and requires a lot of computing time. It involves solving many set inversion problems for which a specific algorithm has been introduced in Appendix A.

Sometimes, the map generated by the SLAM algorithm, is not a required output of the problem. Only the pose estimation is expected. In this situation, when shapes can be fully observed, the formalism of Equation ( 6.1) can be simplified.

An inter-temporal formulation can be used to remove the map as a variable of the problem.

Inter-temporal SLAM

An inter-temporal constraint [Le Bars et al., 2012, Bethencourt andJaulin, 2014] links two observations of the same object, taken at two different times. It can be somehow related to Graph SLAM methods [START_REF] Thrun | The Graph-SLAM Algorithm With Applications to Large-Scale Mapping of Urban Structures[END_REF] where reformulations and substitutions are performed to reduce the complexity of the problem.

The main idea behind this approach is to only consider transformations between two measurements, instead of between a measurement and the map.

Given two times (t 1 , t 2 ) ∈ T 2 , the map equation can be changed into

g x(t i ) (Z(t i )) = M g x(t j ) (Z(t j )) = M ⇔ g x(t i ) (Z(t i )) = g x(t j ) (Z(t j )) (6.5)
Now, in order to fit the shape transformation framework, we assume that the relation between the shapes Z(t i ) and Z(t j ) can be expressed by a function h p ij that depends on a parameter p ij . This parameter links the states of the system between times t i and t j . This link is defined by a function ψ. From Equation (6.5), we have:

Z(t i ) = g -1 x(t i ) • g x(t j ) (Z(t j )) (6.6) ⇔ Z(t i ) = h p ij (Z(t j )) p ij = ψ(x(t i ), x(t j )) (6.7)
For instance, in the translation invariant case where g x (z) = z + x, we have

g -1 x(t i ) • g x(t j ) (z) = x(t j ) -x(t i ) + z (6.8)
which becomes:

h p ij (z) = p ij + z (6.9) p ij = φ(x(t i ), x(t j )) = x(t j ) -x(t i ) (6.10)
The SLAM problem is transformed into the following CN: .11) where (t i , t j ) ∈ T 2 are times of measurements. Algorithm 6.2 describes the propagation process used to solve the CN. As before, the trajectory is contracted using C d dt (line 4). For all pairs (t i , t j ) ∈ T 2 , a two-step contraction process is used. Line 6, C SRAC contractor is used to contract [p ij ] with respect to constraint 2. of CN (6.11). Line 7, the contractor C φ is used to propagate this information to the involved states. It can be built using classical forward-backward propagation methods. For instance, the contractor defined in Example 2.4 can be used for Equation (6.10). for (t i , t j ) in T 2 do 6:

                         Variables: x(•), M, Z(t i ), p ij Constraints: 1. ẋ(•) = f (x (•) , u (•)) 2. Z(t i ) = h p ij (Z(t j )) 3. p ij = φ(x(t i ), x(t j )) Domains: [x](•), [M], [Z] (t i ), [p ij ] ( 6
C SRAC ([p ij ], [Z(t i )] , [Z(t j )]) 7: C ψ ([x(t i )], [x(t j )], [p ij ]) 8:
end for 9: until no contraction on [x](.) This is illustrated by Figure 6.4 where two shapes are registered. Roughly, a given parameter vector is feasible if all red areas do not intersect the blue ones. The contractor based approach will guarantee that no feasible parameter vector will be removed. Consequently, this contractor can be applied to any pair (t i , t j ) ∈ T 2 , without taking particular precaution.

Bathymetric SLAM

In order to illustrate the behavior of Algorithm 6.2, we consider a simulated test case that involves an AUV able to sense its environment using a multi-beam echo sounder. The vehicle made a survey pattern at a speed of 4 knots (around 2 m.s -1 ) and covered an area of 600m 2 in 40 minutes. The AUV is equipped with an INS coupled with a DVL that provides speed and heading data every 0.1 seconds. Every 25 seconds (around 40 meters), the last 30 seconds (75m) of bathymetric data are processed in order to build a local Digital Elevation Map (DEM). This DEM is segmented based on one or several elevation thresholds. The short-time positioning drift during this time window is assumed to be negligible with respect to data uncertainty.

The attitude (roll, pitch and heading), depth and altitude of the AUV are directly measured. Their uncertainties are taken into account when the DEM is built. Consequently, only the horizontal position of the robot needs to be estimated. If this hypothesis is not acceptable, the method presented in Appendix B can be used to build a visibility mask in order to only take into account reliable data.

Figure 6.5 provides a simulation of the AUV trajectory on an unknown map. This map, associated to bathymetric data, splits the space into two complementary sets.

A DEM can be seen as a shape of R 3 , which separates points above the seabed from those which are below. Now, in order to illustrate the proposed method, we choose to discretize it into a series of shapes of R 2 : M 1 , M 2 , . . . , M m based on level sets extracted from the original DEM. The SLAM problem is then described by the following equations:

     ẋ(t) = f (x(t), u(t)) (evolution equation) Z j (t) = g(x(t), M j ) (observation equation) x(0) ∈ X 0 (initial state)
Now, combining several shapes associated to different level sets can be trivially done using the constraint network formalism. We consider only the shape associated with one level set. For instance, the white area corresponds to points with an elevation value that ranges between 10 and 13 meters. The black one corresponds to elevation values that are outside this interval. An example of a measured thick shape is given in Figure 6.5b.

x(t 2 )

x(t 1 )

x(t i ) This SLAM problem is described by the following equations: 6.12) where the state vector x = (p x , p y ) is the horizontal position of the AUV. The inputs are its horizontal speed vector v(t) and its headings ψ(t). The quantities v(t) and ψ(t) are measured every 0.1 seconds with an error of 0.02 m.s -1 and 0.05 • respectively. The initial state is assumed to be known and taken at x 0 = (0, 0) T . Since the heading of the vehicle is known with a good accuracy, only translations are considered between the observation shapes and the map. This is modeled by the following system of equations:

     ẋ(t) = cos ψ(t) -sin ψ(t) sin ψ(t) cos ψ(t) • v(t) t ∈ R Z(t i ) = M -x(t i ) t i ∈ T ⊂ R (
                         Variables: x(•), M, Z(t i ), p ij Constraints: 1. ẋ(•) = f (x (•) , u (•)) 2. Z(t i ) = Z(t j ) + p ij 3. p ij = x(t j ) -x(t i ) Domains: [x](•), [M], [Z] (t i ), [p ij ] (6.13)
where Algorithm 6.2 can be used.

The propagation process is now illustrated. The initial tube, which encloses the true trajectory is depicted in red on Figure 6.7. It was obtained by integrating the differential equation, from the initial vector x 0 , using the contractor C d dt . The diameter of the tube is shown on Figure 6.7b. Since, only the proprioceptive measurements are used, its diameter grows almost linearly with time. In order to reduce the position drift, all pairs (t i , t j ) ∈ T are registered.

The improvement obtained after a pair registration is illustrated with the two measurements, acquired by the two red AUV in Figure 6.5. These measurements, acquired at times t 1 = 4 min 15 s and t 2 = 26 min 30 s are depicted in Figure 6.6. The contractor C SRAR is first used to estimate the translation vector between these two thick shapes. The values of the translation parameter vector before and after the contraction are given in Table 6.1. Then, the contractor C -is enforced on the states at t 1 , and t 2 . Finally, the contractor C d dt is called to propagate the contractions on [x](t i ) and [x](t j ) to the whole trajectory. The result is shown by the green tube on Figure 6.7.

After 10 iterations Algorithm 6.2 stopped in less than 5 minutes. The resulting tube is shown on Figure 6.8 with its diameter. The position uncertainty remains bounded during the mission which is the expected result. 6.5, are registered using the contractor C SRAC . Then the contractor C -is used to contract the states at t 1 and t 2 . Note that only a contraction along the y-axis has been done which is consistent with the shapes used.

Discussion

The proposed method provides a way to deal with the SLAM problem where the observations and the map can be observed from double-sided (inside and outside). Thanks to the thick shape model, an efficient algorithm, based on inter-temporal constraints, has been proposed. It allows localizing the vehicle without explicitly building the map.

In the simulation, only the shape constraint, that links a given map representation with the measurement shape, has been used. Since the formalization of the problem fits the constraint network framework, additional constraints can be combined to improve the efficiency of the constraint propagation process.

Another improvement is to avoid testing all the pairs (t i , t j ). Using the thick set inversion algorithm, it is possible to define a test, based on the scope of the sensor, which checks if two shape-measurements overlap or not. For two visibility sets V(t i ) and V(t j ) and a translation vector [p ij ], the following thick set is considered:

[X] = [V] (t i ) ∩ ([V] (t j ) + [p ij ]) (6.14)
If X + = ∅, the two measurements do not overlap and the pair (t i , t j ) should not be tested anymore. Else if X -= ∅, they overlap for sure. This information can be used, as in [Jaulin, 2016], to merge and associate data. In the last case, where [X] = [∅, X + ], the translation is too uncertain to conclude and the pair must be tested.

If the measurement shapes [Z](t i ) are accurate enough, they do not need to be contracted. Consequently, only the transformation between two thick shapes needs to be estimated, which simplifies the contractor C SRAR . This hypothesis is done for this test case since we considered that bathymetric measurements are reliable and accurate.

Finally, the proposed method needs to be validated with real data. With respect to the targeted application, the specific process, that converts raw sensor data to thick shapes, needs to be developed.

Dig SLAM

The previous section deals with double-sided observable shapes that provide an inner and an outer observation of the map. In this section, the case of sensors that provide one-sided shape observation of the environment is discussed. This corresponds to the extension of the localization problem in an unstructured environment, introduced in Chapter 4.4, to the SLAM context.

Formulation of the Dig-SLAM

We consider the problem taken from [Jaulin, 2011], where a robot is moving in its environment and uses an omnidirectional range sensor to sense the free area around him. At a given time t, this sensor measures two sets: the free zone, denoted by Z(t), which represents the area in the scope of the sensor free of obstacles, and the impact zone, denoted ∆Z(t) which is the border of Z(t) and contains the impact point. Note that measurements from most rangefinders (laser, infrared, ultrasound-based) can be cast into such two sets. Some examples are given in Figure 6.9.

This problem fits the formalism defined in this thesis and is modeled by:

     ẋ(t) = f (x(t), u(t)) (evolution equation) Z(t) = g x(t) (M) (observation equation) x(0) ∈ X 0
(initial state) (6.15) where M depicts the set of points in the working space which are free of obstacles and Z(t) is the space dug by the measurement taken at time t. With respect to the uncertainty of the sensor, the measurement is changed into an interval shape

[Z] = [Z -, Z + ] with Z + = R 2 .
This corresponds to a one-sided observable shape, defined in Section 6.2.3.2.

The function g : R 2 × R 2 → R 2 is a transformation that links the local observation Z(t), expressed in the robot frame, to the map. For instance, in the translation invariant case, g is defined by: g x(t) (m) = mx(t). (6.16) Now, in order to take into account the fact that the impact zone ∆Z intersects the map, the following constraint needs to be considered: (6.17) In order to use the contractors previously introduced, the Dig-SLAM problem is formalized by the following constraint network:

∆Z ∩ g x(t) (M) = ∅.
                         Variables: x(•), M, Z(t i ) Constraints: 1. ẋ(t) = f (x (t) , u (t)) 2. Z(t) = g x(t) (M) 3. ∆Z(t) ∩ g x(t) (M) = ∅ Domains: [x](•), [M], [Z] (t) = [Z(t), R q ], [∆Z] (t) = [∅, ∆Z(t)] (6.18) d(t) Z - ∆Z - d(t) Z - ∆Z - d(t) Z - ∆Z - (a) (b) (c)
Figure 6.9: Examples of measurements which can be converted into interval shapes

[Z] = [Z -, Z + ] with Z + = R 2 .
The dotted area represents the free space Z -and the light gray the impact zone ∆Z(t). (a) an omni-directional rangefinder, where the impact point lies on a circle and the free space is a disk. (b) an ultrasonic rangefinder where Z -is a sector and ∆Z(t) is an arc. (c) represents the area swept by a laser rangefinder pointing perpendicularly to the path of the robot where methods introduced in Appendix B can be used. The set ∆Z(t) is depicted by a curve.

where 1. describes the evolution of the robot with respect to its inputs using the evolution function.

2. corresponds to the measured free zone that must be inside the map.

3. corresponds to the impact zone that must intersect the border of the map.

However, because of the nature of [Z] and the fact that [M] needs to be estimated, Equation ( 2) of ( 6.18) can only be used to contract the inner approximation of the map. The following development justifies the need of considering the impact zone.

Given [M] = [∅, R q ] the initial map, a measurement [Z] (t) = [Z(t), R q ] and a box [x](t). The contractor C F defined in Proposition 5.3 on page 112 applied on [Z] and [M] leads to :

C g [x](t) [Z] [M] = [Z] g [x](t) ([M]) [M] g -1 [x](t) ([Z] (t)) (6.19) =   Z ∪ g [x](t) (∅) , R q ∩ g [x](t) (R q ) ∅ ∪ g -1 [x](t) (Z) , R q ∩ g -1 [x](t) (R q )   (6.20) = [Z, R q ] g -1 [x](t) (Z) , R q (6.21)
Consequently, only M -can be contracted and M + will remain undefined. This approximation is called on [Jaulin, 2011] the dig zone and represents the part of the free space that has been proved free of obstacles by the vehicle.

The registration contractor, defined by Proposition 5.4 needs an outer approximation of the map in order to contract the parameter's domain. This is why, Equation 3of CN (6.18) needs to be used with Proposition 5.5 on page 114. We define the contractor

C = g [x](t i ), ∆Z(t i ), M -= [x](t i ) ∩ proj p (∆Z(t i ) × R q ) ∩ g(M -) (6.22)
The following algorithm can be used for solving the Dig-SLAM problem: Line 5, the map is updated and Line 6, the position of the robot is contracted by taking into account the fact that ∆Z intersects M.

Algorithm 6.3 Dig SLAM(in: T, Z(t i ), ∆Z(t i ), inout: [x](.)) 1: [x](0) = x 0 2: repeat 3: [x](.) := C d dt ([x](.)) 4:
for t i in T do 5:

M

-:= M -∪ g -1 [x](t i ) (Z) 6: [x](t i ) := C = g ([x](t i ), ∆Z, M -) 7:
end for 8: until no contraction on [x](.) Note that the map needs to be explicitly built in order to contract the position of the vehicle. This algorithm is illustrated with the following application.

Application: Range only SLAM in an unstructured environment

From an academic point of view, the omnidirectional range-only posed-based SLAM can be seen as a canonical problem. It is the simplest and most significant of a large class of SLAM problems that are difficult to solve. The development of tools to solve properly and efficiently this problem will be useful to solve many other SLAM problems.

As an illustration, consider a mobile robot described by the following range-only SLAM equations:

     ẋ1 (t) = u 1 (t) cos(u 2 (t)) ẋ2 (t) = u 1 (t) sin(u 2 (t)) y(t) = d(x(t), M) (6.23)
The inputs of the system are the speed u 1 and the heading u 2 of the robot. The measurement y(t) corresponds to the closest distance of the robot to the map which could have been obtained by using an omni-directional sonar with angular aperture After three iterations, no significant contraction happens. The maximal pose uncertainty is obtained when the robot is on the left side of the map. The minimum is obtained when it is near its starting position.

Discussion

This method is efficient to solve the SLAM problem where partially observed shapes are used to build the map and localize the map.

In the translation invariant case, the localization process is based on a dilation of the map by ∆Z. If M -contains small holes, as it is the case in Figure 6.12, these holes grow which limits the efficiency of the contraction process. The computation of the inner approximation of the map must be done carefully.

Arbitrary shapes that are depicting the free space can be considered and, in a theoretical point of view, the heading can be estimated but with a higher computing cost.

Conclusion

In this chapter, the Simultaneous Localization and Mapping in an unstructured environment has been tackled using shapes. To our knowledge, this is new in a SLAM context and this contitutes the major contribution of this thesis. This problem is known as hard to solve due to the fact that no geometrical features can be extracted from the measurements. We proposed a formalism that allows us to handle rigorously nonlinearity and sensor uncertainties while keeping the guarantee and integrity of the results, with respect to the initial hypothesis.

In order to illustrate the proposed methods, two test-cases have been considered as examples. In the first one, the shape-based SLAM has been demonstrated on a simulated test-case based on bathymetric data. This type of datas allows to build a double-sided shape for which an inter-temporal-based approache has been used. The second test-case illustrated the dig SLAM formalism which involves sensors that are able to measure the parts of the working space which are free of obstacles. This problem is an extension of the application presented in Chapter 4 in the SLAM context.

The thick shape representation allows us to efficiently represent sensor uncertainties and the limit of the scope of the observation. This last point allows a paire-wise registration without preliminary associations. Merge different sensor.

The proposed method provides a generic abstraction of sensor reading that remains strongly bounded with physical features of the environment. The main difficulty is to be able to find robust algorithms to transform raw data into thick shapes.

Outliers can eventually be handled using the extension of the relaxed intersection on thick sets that has been presented in [START_REF] Desrochers | Relaxed intersection of thick sets[END_REF]. The constraint network must be reformulated in order to efficiently take advantage of this operator. feasible parameters of the transformation which links two sets that we called shapes. This primitive constraint is at the center of our algorithms. A first localization application in an unstructured environment is then tackled using the polar separator (published in [Desrochers and Jaulin, 2017c]). Both the map and measurements are defined by shapes which, in this chapter, are assumed to be known.

In order to handle the shape uncertainty, the notion of intervals of sets, called a thick set, has been introduced in Chapter 5. Even if the notion of thickset is not new in the literature, a main contribution of this thesis is its exploitation for solving robotics problems. At the heart of the thick set theory, lies the notion of the penumbra. The penumbra characterizes the uncertainty that comes from the initial hypothesis of the problem. It allows us to prove that, for a box [x] inside the penumbra, nothing can be proved for any of its elements x with respect to the solution set. This was illustrated in Appendix B where points that belong to the penumbra are seen for some feasible trajectories but not for all. Without uncertainty of the position of the vehicle, the penumbra would not have existed.

When dealing with sensor readings, the penumbra makes it possible to classify parts of the research space for which not enough information is available to conclude anything. Either because the area has not been sensed or because the process that builds the thick set from measurements cannot classify this part of the space. This can be the case with automatic segmentation algorithms.

The arithmetic and operations that can be performed on thick sets were also introduced. Another work which has not been included in the manuscript concerns the relaxed intersection of thick set presented in [START_REF] Desrochers | Relaxed intersection of thick sets[END_REF].

The shape registration problem, introduced in Chapter 4 was extended to the uncertain case, were variables of the constraint network are shapes and their domains are thick sets, that we called thick shapes.

Based on the notion of contraction on thick sets, the shape carving problem has been defined. Given two thick shapes linked by an uncertain function, it aims at contracting the thick shapes using shape contractors. This contraction was done by solving an uncertain set inversion problem. When the uncertain function depends on parameters, a first method based on projection has been provided.

With respect to the SLAM problem, this constraint is mainly used to build the map or filter the measurements.

Finally, Chapter 6 deals with our targeted application. A SLAM algorithm that involves shapes has been introduced. The formalism, based on thick shapes, allows us to handle efficiently unstructured environments. It has been illustrated through two examples.

Open-source library

A significant contribution of this work is the development by the author of the new open-source library Pyibex, freely available at:

• http://www.ensta-bretagne.fr/desrochers/pyibex This project gathers all the elementary tools presented in this document. The reader will be able to process the simulated examples and build its own solvers for the resolution of shape-based problems.

Overall Prospects

We are convinced that the proposed methods are relevant to address the SLAM problem in complex environments. But, in order to be used, they require a way to be tried and diffused. For that purpose, the development of the pyIbex library will be continued in order to make the use of interval analysis in robotics applications easier.

The shape registration problem relies mainly on the use of the projection algorithm. Its efficiency must be improved in order to be able to handle more complex transformations, such as the rotation.

The key point toward real applications is the processing chain that converts raw sensor data to interval shapes. This chain needs to robustly and reliably classify sensor outputs taking into account the notion of the penumbra. For a given set, the classical representation binaries the world. A point must belong or not to the shape. On the contrary, the notion of interval shape allows us to handle all data, which cannot be reliably classified inside or outside a measurement shape. For instance, with side scan sonar images, shadows are classified in the penumbra.

The interest of the shape carving process has not been fully demonstrated in this manuscript. The shape carving method may be used to perform pattern recognition. For instance, an object partially defined by an interval shape can be matched with another one. The interval shape representation allows to only consider the reliable parts of the initial object which are compared.

Outliers could also be taken into account by using the thick relaxed intersection. This needs to be correctly formalized and validated on actual experiments. The chapter is organized as follows. Section A.3 presents the new notion of thick intervals and thick boxes to be used for solving the thick set-inversion problem. Section A.4 generalizes the classical set-inversion algorithm to the thick case by introducing the new notion of thick inclusion function. Section A.5 illustrates the principle of the method on five test cases with one involving an actual underwater robot. Section B.6 concludes the chapter.

The results in this chapter are extracted from [START_REF] Desrochers | Thick set inversion[END_REF]. Thick objects (functions, intervals, sets) are denoted with a double brackets notation . .

A.2 Problem Statement

Notation. We introduce a specific notation involving the quantifier ∀ when dealing with thick sets. Given two thick sets A and B , we define: As explained in Section 5.4,when [f ] can be parametric, projection algorithm can be used but required additional bisection along the domain of the parametric variable. Although generic, this method induces an extra computational cost and failed to handle efficiently high dimensional problems. To avoid this, additional notions such as thick intervals are introduced in the following section.

( A ⊂ B ) ∀ ⇔ ∀A ∈ A , ∀B ∈ B , A ⊂ B ( A ⊂ B ) ∀ ⇔ ∀A ∈ A , ∀B ∈ B , A ⊂ B ( A ∩ B = ∅) ∀ ⇔ ∀A ∈ A , ∀B ∈ B , A ∩ B = ∅ ( A ∩ B = ∅) ∀ ⇔ ∀A ∈ A , ∀B ∈ B , A ∩ B = ∅.

A.3 Thick Intervals

Denote by IR the set of all intervals of R. A thick interval x (see, e.g., [Chabert and Jaulin, 2009a]) is a subset of IR which can be written under the form

x = [x -] , [x + ] = {[x -, x + ] ∈ IR | x -∈ [x -] and x + ∈ [x + ]} . (A.12)
Here, [x -] , [x + ] are two intervals containing the lower bound x -∈ R and the upper bound .14) called the subset bound and the supset bound of x then

x + ∈ R of an uncertain interval [x -, x + ]. If we define the two intervals of R [x ⊂ ] = [x -, x + ] ∈ IR | x -∈ [x -], x + ∈ [x + ] , (A.13) and [x ⊃ ] = [x -, x + ] ∈ IR | x -∈ [x -], x + ∈ [x + ] , ( A
x ⊂ {[x] ∈ IR | [x ⊂ ] ⊂ [x] ⊂ [x ⊃ ]} , (A.15) with an equality if [x -] ∩ [x + ] = ∅.
As a consequence, a thick interval is not necessarily a thick set: it is more precise or equivalently it is narrower. This can be explained using the endpoints diagram [Kulpa, 1994] (see Figure A.2) where an interval is seen as a point of R 2 . For instance, to the interval [START_REF]Density probability function of random variable x and[END_REF]7], we associate the point with coordinates [START_REF]Density probability function of random variable x and[END_REF]7). The degenerated intervals, such as [2,2], all belong to the diagonal. This representation provides a geometrical representation of the relation between intervals. For instance,

[x] ⊂ [y] if [y] is at the top left of [x]
. The intersection between two intervals (or the interval hull) is obtained by taking the bottom-right corner (or the top left corner) of the smallest box which encloses the two interval points. For instance, if [x] = [START_REF]Density probability function of random variable x and[END_REF]4] and [y] = [2,5], the enclosing interval box is painted red. The top left interval is [x] ∪ [y] = [START_REF]Density probability function of random variable x and[END_REF]5] and the bottom-right interval is [x] ∩ [y] = [2,4]. This red box corresponds to the thick box a = [START_REF]Density probability function of random variable x and[END_REF]2], [4,5] . The subset and supset bounds are [a ⊂ ] = [2,4] and [a ⊃ ] = [START_REF]Density probability function of random variable x and[END_REF]5].

In this figure, the orange polygon corresponds to the thick interval b = [3,7], [6,8] . The subset and supset bounds are [b -] = ∅ and [b + ] = [3,8]. As illustrated by the gray zone of Figure A.2 (left), a subset-supset representation adds pessimism when

[x] x 1 x 2 y 1 y 2 [f ] a [f ](a) [f -]([x]) [f + ]([x]) f -(a)
f + (a)

+ . The boxes [f -] ([x]) and [f + ] ([x]) contain f -(a) and f + (a). Therefore, the box [f ] (a) is inside the thick box f ([x]) = [f -] ([x]) , [f + ] ([x]) .

A.4.3 Algorithm

Algorithm A.2, named ThickSivia (Thick Set Inversion Via Interval Analysis), provides an approximation of the solution of the thick set inversion problem X = [f ] -1 ( Y ). The input of this algorithm are (1) the box [x] which is assumed to be large enough to contains X + , the upper bound of X , (2) an accuracy ε > 0, (3) the thick inclusion function f , and ( 4) thick set Y . The output is an approximation of the thick set X = X -, X + . The algorithm decomposes the initial box [x] into four non-overlapping subpavings: (1) The inner subpaving L clear which contains boxes that have been proved to be inside the clear zone X -, (2) the outer subpaving L dark which contains boxes that have been proved to be outside X + (i.e., inside the dark zone), (3) the subpaving L penumbra which contains boxes that have been proved to be inside the penumbra X + \X -and ( 4) the subpaving made with boxes that have been rejected (for which nothing is known and with a width smaller than the desired level of precision ε).

Remark A.2. In the algorithm, we have several tests on thick boxes that require comparisons on thick boxes as introduced in Section A.3. Consider for instance the test (

f ([x]) ∩ Y + = ∅)
∀ and assume that Y + is made with boxes {[y] (1), . . . , [y] ( k)} (this will be the case in the applications presented in Section A.5). Since we have

f ([x]) ∩ Y + = ∅ ∀ ⇔ ∀k, ( f ([x]) ∩ [y] (k) = ∅) ∀ .

Algorithm A.2 Thick set inversion algorithm: ThickSivia

Input: [x], ε, Y , f Output: L clear , L penumbra , L dark 1: L = {[x]}, L clear = ∅, L penumbra = ∅, L dark = ∅ 2: while L = ∅ do Unstack L into [x] 3: if ( f ([x]) ⊂ Y -) ∀ then push [x] into L clear 4: else if ( f ([x]) ∩ Y + = ∅) ∀ then push [x] into L dark 5: else if ( f ([x]) ∩ (Y + \Y -) = ∅) ∀ then push [x] into L penumbra 6:
else if width([x]) > ε then bisect [x] perpendicularly to its largest side and push the two resulting boxes in L 7:

end if 8: end while Our test amounts to checking that ( a ∩

[b] = ∅) ∀ where a = f ([x]) and [b] = [y] (k). Therefore, ( a ∩ [b] = ∅) ∀ ⇔ ∃i ∈ {1, . . . , n}, ( a i ∩ [b i ] = ∅) ∀ (see Proposition A.2) ⇔ ∃i ∈ {1, . . . , n}, b - i -[a i ] + ⊂ R -∧ [a i ] --b + i ⊂ R - (see Proposition A.1) ⇔ ∃i ∈ {1, . . . , n}, b - i -lb [a i ] + ≤ 0 ∧ (ub ([a i ] -)) -b + i ≤ 0 where lb([a i ] + )
if ( f ([x]) ⊂ Y -) ∀ or ( f ([x]) ∩ (Y + \Y -) = ∅) ∀ . Algorithm A.3 Test if ( f ([x]) ∩ Y + = ∅) ∀ Input: a = f ([x]), Y + = {[y] (1), . . . , [y] ( k)} 1: for k = 1 to k do [b] = [y] (k) 2: for i = 1 to n do If b - i -lb [a i ] + > 0 ∨ (ub ([a i ] -)) -b + i > 0 then Return False. End. 3: end for 4: end for 5: Return True i 1 -1 i 1 0 2 1 2 3 2 3 ∞ 1
Table A.1: results of Formula A.31 for different values of i.

A.4.4 Properties Termination

The algorithm always terminates in less than λ = with([x])

ε dim(x)
iterations, where [x] it the input box. The number λ is huge and corresponds to the worst case situation where all tests fail and the algorithm returns X clear = ∅, X penumbra = ∅, X dark = ∅. In practice, as shown in [START_REF] Walter ; Jaulin | Guaranteed nonlinear parameter estimation from bounded-error data via interval analysis[END_REF] the number of

iterations is O A 1 ε dim(x)-1
, where A is the area of the accumulation zone which is composed here with the boundary of the penumbra. Now, we can find some examples where the number of iterations is larger than The interval is not classified by our algorithm if a i + a -b < 0, i.e., if a i < w ([a, b]). Since in our algorithm, for all interval [a, b], we have ε < w ([a, b]) (otherwise, the interval is eliminated), we conclude that if a i < ε, or equivalently, a < ε 1 i the box cannot be eliminated. Thus, the number of remaining yellow boxes η is of order ε The Table A.1 illustrates this formula for different values of i.

We tried to get an experimental validation of this formula, and we obtained the following figure which is consistent with the convergence formula. Now, even for The box [c] is proved to be inside the dark zone, i.e., outside X + .

Theorem A.3. For a given ε, our algorithm provides three lists L clear (ε),L penumbra (ε) and L dark (ε). Take a point a. For ε sufficiently small we have

(i) [f ](a) ⊂ int (Y -) ⇒ a ∈ L clear (ii) [f ](a) ⊂ int Y + ⇒ a ∈ L dark (iii)
[f ](a) ∩ int (Y + \Y -) = ∅ ⇒ a ∈ L penumbra (A.37)

where int(A) denotes the interior of the set A [Berge, 1963].

Proof. The proof is by contradiction. Assume that for all ε the box containing a is never classified. It means that there exists a sequence of boxes [x] (k) converging to a such that none of the three tests is satisfied for all [x] (k). As a consequence, we get that if either (i), (ii) or (iii) is satisfied then a will be classified inside one of the three lists. Moreover, from Equation (A.32), we get that a will be classified on the right list.

A.5 Test-Cases

This section provides five test-cases to illustrate the efficiency of our method. All these test cases solve a thick inversion problem X = [f ] -1 ( Y ). In the figures, all red boxes are shown to be inside the clear zone X -; all blue boxes are inside the dark zone, i.e., outside X + ; all orange boxes are proved to be inside the penumbra. Nothing is known for the small yellow boxes. As we can see on this figure, the penumbra is better (i.e., without any uncertain boxes inside) characterized with a lower/upper bound representation for the thick intervals.

We compared the computing time (on a processor i5-2520M@2.50GHz) and the number of bisections with the traditional approach (which does not characterize the penumbra) and our method. We get the table below. We observe that when ε is small, classical methods are much less efficient due to the fact many bisections take place inside the penumbra. This observation is also valid for all five test-cases considered in this paper. In the clear zone (red), the robot is able to communicate with all marks. everywhere inside the penumbra. Using thick interval arithmetic, we are able to conclude for an orange box that there is no need to bisect it.

A.6 Conclusion

This chapter deals with the set-inversion problem X = f -1 (Y) in the case where both f and Y are uncertain, i.e., f belongs to the interval of functions [f -, f + ] and Y belongs to a thick set, i.e., an interval of sets Y = Y -, Y + . After introducing the new notions of thick intervals and thick boxes, a new algorithm for set inversion has been proposed. It is able to compute a thick solution set X = X -, X + containing all feasible solution sets.

From the computational point of view, thick intervals allow us to have a better understanding of the uncertainty. For instance, for the set inversion problem, we are able to detect that a box is included in the penumbra X + \X -. In this penumbra, we can conclude that any bisection would be useless. This could not have been detected using classical intervals. As a consequence, the accumulation zone (i.e., the part of the search space where tiny boxes are still bisected) for thick interval based algorithms has a zero volume, since it corresponds to the boundary of the penumbra. Using classical intervals instead, we could obtain similar results, but the accumulation zone would correspond to the whole penumbra, which has a nonzero volume. As a result, a large part of the computational burden made by traditional interval algorithms is done on a part of the search space which has no influence on the final result.

get an inner and an outer characterization of the reachable space or in [Jaulin, 2011] for SLAM involving unstructured maps.

The problem to be considered here is the characterization of the explored zone Z of a robot. The set Z is defined by

(i) ẋ = f (x, u) , x (0) = x 0 , u (t) ∈ [u] (t) (ii) Z = t≥0 V (x (t)) (B.1)
where (i) describes the evolution of the robot and (ii) defines the explored zone Z. In the state equation (i) of the robot, x ∈ R n is the state vector and u is the input vector. In order to take into account some state noise, we assume that, for all t, a box [u] (t) which contains u (t) is available. Moreover, the initial condition x 0 is assumed to be known. For each t, a scanner on the robot is able to observe a part of its environment. More precisely, for each t, there exists a subset of the environment V (x (t)) ⊂ R q , q ∈ {2, 3} that is visible by the robot. This set is called the visible set (see, e.g., [START_REF] Guyonneau | A visibility information for multi-robot localization[END_REF] in the context of localization). Note that our robot does not use any exteroceptive sensors for state estimation, for localization or to control its trajectory. It is only able to estimate its position from the proprioceptive sensors u with a given bounded accuracy. Equivalently, we consider that the scan sensors of the robot are used to collect data for exploration that will only be analyzed after the mission by some human users. The following examples illustrate how the visible sets could be defined in practice.

Example 1. Consider a robot, the pose of which is x = (x, y, θ), moving in a plane. This robot is able to scan the environment up to a distance of 3 meters in front of it inside a cone of ± π 6 rad. Therefore, the visible set V (x) contains all z ∈ R 2 which satisfy the following inequalities (z 1 -x) 2 + (z 2 -y) 2 ≤ 9 cos θ • (z 1 -x) + sin θ • (z 2 -y) ≥ cos π 6 (B.2) Example 2. We consider the robot of Example 1 except that now, this robot is able to scan only points that are exactly on its right at a distance in [2,3]. Then V (x) is the set of all z ∈ R 2 such that of all z that have been seen by the robot at times t 1 and t 2 for the two trajectories.

To characterize Z -, it is not sufficient to compute the union of all certainly visible sets ; we also need to add all cross intersections. For instance, whatever is the true trajectory among x a or x b , we are certain that all points of the dark gray set (at the intersection between the two opposite pies) in the figure have been explored. These cross intersections are not taken into account by existing methods [Drevelle et al., 2013].

Now, for sweep exploration, the dimension of the visible sets is q -1 and the certainly visible sets x∈X (t) V (x) become empty. All the information about Z - can thus only be obtained from these cross intersections, i.e., by taking into account different t and different trajectories all together or equivalently, by considering the right-hand side of (B.7) for the characterization of Z -.

B.4 Characterization of the Explored Zone

We now focus our attention on the sweep exploration where visible sets have the form V (x) = {z ∈ R q | ϕ (z, x) = 0 and ψ (z, x) ≤ 0} . (B.8)

In this formula, ϕ : R p × R n → R is the visibility function which is assumed to be continuous and ψ : R p × R n → R models the scope of the sensor. Examples 1, 2, 3 of Section B.2 provide an illustration of these functions.

This section proposes a new algorithm able to bracket the explored zone between two subpavings (union of non-overlapping boxes). More precisely, the workspace will be partitioned into 4 types of boxes: (a) the clear boxes that have been proved to be inside Z -, (b) the dark boxes that have been proved to be outside Z + , (c) From (B.10,ii), we get that there exists t 0 ∈ [t 1 , t 2 ] such that ϕ (z, x (t 0 )) = 0 ∧ ψ (z, x (t 0 )) ⊂ R - which is inconsistent with the left-hand side of (B.17 Assume now that the left-hand side of (B.18) is also true. From (B.10,ii), we get that for a specific t 0 ∈ [t

B.5 Experiment

This section illustrates our sweep exploration method on a real experiment. November 2015, a 46 minutes mission has been performed in the Road-stead of Brest (France, Brittany) with the underwater robot Daurade which realized a classical survey pattern composed of a set of parallel tracks with an altitude of about 10 meters. This robot (see Figure B.5) has been built by ECA robotics and used by DGA Tn (Direction Général de l'Armement -Techniques Navales) and SHOM (Service Hydrographique et Océanographique de la Marine) for performing REA (Rapid Environment Assessment) studies. REA is intended to survey the environmental conditions of a particular location in order to identify any existing or potential dangers. In the counter mine warfare context, attention is focused on mapping the sea floor with acoustic sensors. Daurade is equipped with a Side Scan Sonar (Klein 5500) used to detect potential mines. With this sonar, data are recorded on a line perpendicular to the path of the sensor and images are formed by putting side by side these lines. Characterizing the zone seen by the sonar enters inside the framework of sweep exploration.

For the navigation, Daurade relies on an inertial central (Phins II IXBlue) coupled with a DVL (Doppler Velocity Log), which returns the speed with respect to the ground. Once under the water, no GPS data are available and the estimated position of the robot drifts with the time. A key point of REA missions is to guarantee that the area of interest has been totally explored without any gap. outliers (which mainly occur on exteroceptive sensors), (ii) to data association errors, or (iii) to any information resulting from an interaction between the robot and the environment. The method can also be used before the mission in order to guarantee that a given area will all be covered by the robot, whatever the noise on the proprioceptive sensors are. We only need to assume that the state noise is bounded with known bounds in order to be able to enclose the trajectory x (•) inside a tube [x] (•).

To solve the problem, we had to introduce new interval tests to classify boxes as clear (certainly explored), dark (certainly unexplored) or in the penumbra. This classification requires the development of new quantifier elimination procedures, based on four new theorems. The efficiency of the resulting algorithm has been validated on an actual experiment made by an autonomous underwater robot.

Résumé

Cette thèse étudie le problème de la localisation et de la cartographie simultanée (SLAM), dans des environnements non structurés, c'est-à-dire, qui ne peuvent pas être décrits par des équations ou des formes géométriques. Ces types d'environements sont souvent rencontrés dans le domaine sous-marin.

Contrairement aux approches classiques, l'environnement n'est pas modélisé par une collection de descripteurs ou d'amers ponctuels, mais directement par des ensembles. Ces ensembles, appelés forme ou shape, sont associés à des caractéristiques physiques de l'environnement, comme par exemple, des textures, du relief ou, de manière plus symbolique, à l'espace libre autour du véhicule.

D'un point de vue théorique, le problème du SLAM, basé sur des formes, est formalisé par un réseau de contraintes hybrides dont les variables sont des vecteurs de R n et des sous-ensembles de R n . De la même façon que l'incertitude sur une variable réelle est représentée par un intervalle de réels, l'incertitude sur les formes sera représenté par un intervalle de forme. La principale contribution de cette thèse est de proposer un formalisme, basé sur le calcul par intervalle, capable de calculer ces domaines. En application, les algorithmes développés ont été appliqués au problème du SLAM à partir de données bathymétriques recueillies par un véhicule sous-marin autonome (AUV).

Mots-clés: analyse par intervalles, localisation, SLAM, robotique

Abstract

This thesis deals with the simultaneous localization and mapping (SLAM) problem in unstructured environments, i.e. which cannot be described by geometrical features. This type of environment frequently occurs in an underwater context.

Unlike classical approaches, the environment is not described by a collection of punctual features or landmarks, but directly by sets. These sets, called shapes, are associated with physical features such as the relief, some textures or, in a more symbolic way, the space free of obstacles that can be sensed around a robot.

In a theoretical point of view, the SLAM problem is formalized as an hybrid constraint network where the variables are vectors and subsets of R n . Whereas an uncertain real number is enclosed in an interval, an uncertain shape is enclosed in an interval of sets.

The main contribution of this thesis is the introduction of a new formalism, based on interval analysis, able to deal with these domains. As an application, we illustrate our method on a SLAM problem based on bathymetric data acquired by an autonomous underwater vehicle (AUV).

Keyworks: SLAM, AUV, interval analysis, localisation

  Figure 1.1: The raw image 1.1a is segmented and interpreted as an uncertain shape 1.1b.

Figure 1 . 2 :

 12 Figure 1.2: Density probability function of random variable x and y.

Figure 2 . 1 :

 21 Figure 2.1: Examples of different set-membership enclosures of a set X.

Figure 2 . 2 :

 22 Figure 2.2: Intersection, union and difference of two sets.

Figure 2 .

 2 Figure 2.2 illustrates the intersection, union and difference of two sets.

Figure 2 . 4 :

 24 Figure 2.4: Boxes representation and comparison. The two boxes have the same width (w ∞ ) even if the area of the right one is smaller.

Figure 2 . 5 :

 25 Figure 2.5: The image of the box [x] by the function f , which is not necessarily a box, is depicted by the gray area. The dotted area represents the pessimism induced by [f ] with respect to [f ] * . The function f is thin because the image f (x) of x is a point.

Figure 2 . 6 :

 26 Figure 2.6: Contractor consistent with the set X. The dashed blue area has been removed by the contractor. The minimal contractor is denoted C * and returns the smallest box enclosing [x] ∩ X, colored in red.

Figure 2 . 8 :

 28 Figure 2.8: Illustration of the constraint propagation process. Results of contraction are depicted by red boxes.

Algorithm 2 . 1

 21 Paver(in: [x], S, out: X -, X + ) 1: L := {[x]} ; 2: while L = ∅ do 3:

Figure 2 .

 2 Figure 2.13: Illustration of the inversion of separators. The polygon in (a) is rotated of an angle of π 4 in (b) and intersected with M in (c).

Figure 2 .

 2 Figure 2.15: The union of contractors on segments is applied to characterize the border of the polygon (a). A test is then applied on a point on each box removed by the contractor to classify the box inside or outside the polygon (b).

Figure 2 .

 2 Figure 2.16: In figure (a) m is inside P so T (m) = true and in figure (b) m is outside P and the test is false

Figure 2 .

 2 Figure 2.17: Image contractor

Figure 2 .

 2 Figure 2.20a shows, in light gray, the result of C d dt applied on a tube [x 0](.) defined by ∀t > 0, [x 0 ](t) = [-∞, ∞] and [x 0 ](0) = [-3, 3]. Only trajectories which are not consistent with the differential equation and the initial conditions are removed.

Figure 2 .

 2 Figure 2.20: Illustration of C d dt . The left figure shows the tube [x](.) after the first call of C d dt (light gray) and after using with [z](t 1 ) (dark gray).

2

 2

  Figure 2.22: . A path, corresponding to one solution of the wire loop game, is represented. The corresponding motion in the world space is depicted on Figure 2.22b.

  (a) Initial situation (b) Modeling and parameters

Figure 2 .

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

Figure 3 . 1 :

 31 Figure 3.1: Resulting subpaving of Example 3.1 with S (left) and S * (right).The dashed white box encloses the solution set and is obtained after the first call of the separator. It represents the outer approximation of the solution set and is more accurate with S * . All sides touch the border of X. This is not the case with S where the pessimism has to be balanced with additional bisections.

Figure 3 . 2 :

 32 Figure 3.2: The transformation of a minimal separator by a box-conservative function is minimal.

Figure 3 . 3 :

 33 Figure 3.3: Illustration of the polar contractor. Left: before contraction, Right: after contraction.

Definition 3 . 3 .Figure 3 . 4 :Theorem 3 . 2 .

 333432 Figure 3.4: Y in gray and Y can be represented by the union of boxes

Figure 3 . 5 :

 35 Figure 3.5: Action of the two separators defined as the projection of the polar contractor. The initial domains for (x, y) and (ρ, θ) are represented by the large rectangle and the large pie. Red boxes and pies belong to solutions sets while those in blue do not

Figure 3 . 7 :

 37 Figure 3.7: Paving obtained using the polar separator. The dashed white box is the result of the first call of a forward/backward separator

Figure 3 .

 3 Figure 3.8: VAMA (Véhicule Anti Mine Autonome), an AUV owned by DGA TN.

Figure 3 . 9 :

 39 Figure 3.9: Part of the trajectory of VAMA. The robot changes its direction when it reached one way-point.

Figure 3 .

 3 Figure 3.10: Diameter of the boxes[x k ] obtained by the minimal polar separator (blue) and the forward/backward separator (green). The green peak can be explained by the fact that the forward/backward contractor is less efficient when the angle between the measured angles are almost π 4 .

Figure 3 .

 3 Figure 3.11: Localization of the robot for t = 55 sec. Each of the blue pies corresponds to all feasible positions of the robot associated with one measurement.

Figure 3 .

 3 Figure 3.12: The simulated environment for initial localization. The trajectory of the AUV is depicted by the blue lines. Its starting point is drawn by the red dot. The map is composed of 230 landmarks represented by black dots.

Figure 3 .

 3 Figure 3.13: w([p]) with respect to the iteration number.

Figure 3 .

 3 Figure 3.14: w([θ]) with respect to the iteration number. The peaks in the graph are due to the wrapping effect when the vehicle turns.

Figure 3 .

 3 Figure 3.15: Final trajectory with a correct estimation of the initial position and landmarks association. The true trajectory, in blue, belongs to the red tube.

1 Figure 4 . 1 :

 141 Figure 4.1: Example of a shape registration problem associated with the constraint f (p, A) ⊂ B. The feasible set P corresponds to the parameters of the function f such as the transformation of the shape A is included in B.

Definition 4 . 2 .

 42 Consider a separator S([x], [y]) = {S in ([x], [y]), S out ([x],[y])}, the projection of S is defined by :

Figure 4 . 2 :

 42 Figure 4.2: Set projection

Figure 4 . 3 :

 43 Figure 4.3: The map M is represented by the white space outside the hatched area while the unit pie S 1 is painted in gray. The red pie represents the lower upper bound of D which almost touches the border of the map.

Definition 4 . 3 .Proposition 4 . 1 .

 4341 Given two sets A ∈ P(R n ), B ∈ P(R n ), the Minkowski difference [Najman and (Eds), 1000], denoted , is defined by B A = {p | A + p ⊂ B}. (4.27) Given two separators S A and S B for A and B. The Minkowski difference of two separators, defined by

Example 4 . 1 .Figure 4 . 4 : 4 .Proposition 4 . 2 .

 4144442 Figure 4.4: Minkowski difference of the disk B and the rectangle A

Example 4 . 2 .Figure 4 . 5 :

 4245 Figure 4.5: Minkowski sum of a square B and a triangle A computed using a set membership algorithm.

Figure 4 . 6 :

 46 Figure 4.6: Sensor model used by the robot. The map M, in white, corresponds to parts of the space without obstacles.

Figure 4 .

 4 7 (right) corresponds to an approximation of the set M S 1 , obtained using the separator S M S S 1 with a paver.

Figure 4 .

 4 Figure 4.8 (left) corresponds to an approximation of the set M ⊕ -∆S 1 , obtained using a paver with the separator S M ⊕ -S ∆S 1 . It corresponds to the set of positions for the robot such that the impact pie ∆S 1 intersects the outside of the map M.

Figure 4

 4 Figure 4.7: Left: a robot which collects 6 sonar range measurements. All free sectors S i are included in the map M (in white) and the impact pies ∆S i in yellow intersect M. Right: set of all positions for the robot consistent with the fact that the free sector S 1 (in green on the left picture) is inside M.

Figure 4 .

 4 Figure 4.8 (right) corresponds to the set (M S 1 ) ∩ (M ⊕ -∆S 1 ), i.e., it contains the position consistent with both ∆S 1 and S 1 .

Figure 4 Figure 4

 44 Figure 4.8: Left: Positions for the robot consistent with the impact pie ∆S 1 . Right: Positions consistent with the free sector S 1 and the impact pie ∆S 1 .

  meet [x] ∧ [y] is called the intersection and will be denoted by [x] ∩ [y] . The join [x] ∨ [y] is the interval hull, denoted by [x] [y]. It should not be confused with the classical union ∪ of two intervals.

Figure 5 . 1 :

 51 Figure 5.1: Example of thick set. The dark zone is outside the X + .

Figure 5 . 4 :

 54 Figure 5.4: Characterization of [X] using a paver defined by a grayscale image.

Figure 5 . 6 :

 56 Figure 5.6: Left: [Y]. Right: [X] = f -1 ([Y]).

Example 5 . 3 .

 53 Consider [Y] as depicted in Figure5.6, and the function f : R 2 → R 2 defined by:

Figure 5

 5 Figure 5.7: The result of the set inversion using a sampling method along the parameter space is depicted in Figure5.7a. For each θ ∈ [θ], the image of the red box by the rotation of an angle -θ is drawn using a transparency effect. It makes it possible to guess sets X -, X + . Figure5.7b: the result with classical methods where the penumbra is not characterized. Figure5.7c: the result of the set inversion using projection algorithms.

  can be built to contract the feasible parameter box [p]. Then Proposition 5.1 is used to contract [A] and [B] with respect to f [p] .

Proposition 5. 1 ) 5 . 5 .

 155 Example Consider the two thick sets [A] and [B], depicted in Figure5.10a and 5.10b, that are approximations of the same ellipse but from different positions.

Figure 5 .

 5 Figure 5.10: Shape registration and carving. The initial set is approximated by 5.10a and 5.10b. Using the results of the carving process, the parameter box [p] (see Figure 5.11) has been contracted using Proposition 5.4. Then Proposition 5.1 is used to contract [A] (5.10c) and [B] (5.10d).

Figure 5 .

 5 Figure 5.11: Contractions in the parameter space. The red box is [p] = [2.65113, 2.90901] × [-1.14117, -0.95238]

  . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Generic resolution . . . . . . . . . . . . . . . . . . . . . 6.2.3 Shape-based measurements . . . . . . . . . . . . . . . . 6.3 Shape-based SLAM . . . . . . . . . . . . . . . . . . . . . 124 6.3.1 Inter-temporal SLAM . . . . . . . . . . . . . . . . . . . 6.3.2 Bathymetric SLAM . . . . . . . . . . . . . . . . . . . . 6.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Dig SLAM . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.4.1 Formulation of the Dig-SLAM . . . . . . . . . . . . . . 6.4.2 Application: Range only SLAM in an unstructured environment . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6 . 1 :

 61 Figure 6.1: Observed thick shapes.

Figure 6 . 3 :

 63 Figure 6.3: Occupancy based representation. On the left-hand side, the map is represented by the white area and obstacles are depicted in dashed black lines. Two sonar measurements are shown. For each of them, the impact point lies in a sector ∆Z. Since its position cannot be determined precisely, it makes this type of measurement partially observable. The thick shape associated to the measurements is drawn on the right side. The sounded part of the space, in red, is free from obstacles and belongs to Z -. Since nothing can be said about the nature of the environment sensed at a distance greater than the one measured, all the remaining space belongs to the penumbra.

Algorithm 6 . 2

 62 Shape-based SLAM 1: [x](0) = x 0 2: [p ij ] := [-∞,

Figure 6 . 4 :

 64 Figure 6.4: Example of the inter-temporal registration between two measurements.

Figure 6 . 5 :

 65 Figure 6.5: The map in black and white classifies the area into two complementary classes. The AUV starts at the center of the area and makes a survey pattern. On each yellow triangle, a measurement is taken. An example is shown on the right sub-figure which corresponds to the yellow AUV at time t i . Note that the area which has not been sensed belongs to the penumbra. The scope of the sensor is depicted by the red rectangles. The measurements of the two red AUV, have been taken as an example to illustrate the inter-temporal contraction.

Figure 6 . 6 :

 66 Figure 6.6: Thick shapes taken at times t 1 = 4 min 15 s and t 2 = 26 min 30 s.

  diameter, the vertical lines correspond to times t 1 and t 2 used to illustrate the algorithm.

Figure

  Figure 6.7: Illustration of the contraction between timest 1 and t 2 . The initial tube, in blue, is contracted into the red one.

  Final tube which encloses the true trajectory.

  Final situation. The initial diameter of the tube is in red.

Figure 6 .

 6 Figure 6.8: Final result after contraction

Figure 6 .

 6 Figure 6.10: Trajectory of the robot. Feasible positions are inside the map, delimited by the polygon. Some measurements are depicted by black circles and correspond to the distance to the nearest obstacle around the robot. The framebox is [-10, 10] × [-10, 10].

Figure 6 .Figure 6 .Figure 6 .

 666 Figure 6.11: Illustration of the contraction of [x](t), for t = 1230 sec. Left: M -with a measurement ∆Z(t) depicted by the dashed circle. Right: The box [x](t) before contraction is depicted in green. The dashed green area has been removed by the contraction using Equation 3 of 6.24. The subpaving corresponds to M -⊕-(∆Z(t)). The black box corresponds to [x](t) after calling C d dt

Figure A. 1 :

 1 Figure A.1: Left: Classical interval methods accumulate on the thick boundary (the penumbra). Right: the method we propose here will allow a fast treatment of the penumbra. The frame box is [-2, 4] × [-2, 4] and the black box corresponds to [a] .

(A. 5 )

 5 Thick function. Based on the notion of lattices, introduced in Section 5.2.1, the definition of the intervals of function, called a thick function is now presented.

Figure A. 5 :

 5 Figure A.5: The thick inclusion function f ([x]) encloses all boxes [f ] (a) where a ∈ [x].

1 .

 1 Consider for instance the problem of solving f (x) < 0 wheref (x) = x i . Take an inclusion function: [f ]([x]) = [x] i + [x] -[x]. Assume that [x] = [a, b] where a > 0. We have:[f ] ([a, b]) = a i , b i + [a -b, b -a] = a i + a -b, b i + b -a (A.29)

Figure

  Figure A.7: Tests used for the thick set inversion. The box [a] is proved to be inside the clear zone X -; The box [b] is proved to be inside the penumbra X + \X -. The box [c] is proved to be inside the dark zone, i.e., outside X + .

  (i) Since for all k, ( f ([x] (k)) ⊂ Y -) ∀ is false, there exists a sequence [y] (k) ∈ f ([x] (k)) such that ([y] (k) ⊂ Y -) is false. Now, since f is a convergent thick inclusion function for [f ], d H ([y](k), [f ](a)) → 0. Since int (Y -) is an open set, we cannot have [f ](a) ⊂ int (Y -) . (ii) Since for all k, ( f ([x])(k) ∩ Y + = ∅)∀ is false, using the same reasoning as for (i), we get that we cannot have[f ](a) ⊂ int Y + . (iii) Since for all k,( f ([x])(k) ∩ (Y + \Y -) = ∅)∀ is false, again, we conclude that we cannot have ([f ](a) ∩ int (Y + \Y -) = ∅).

Remark A. 3 .

 3 Figure A.8: Left: With classical intervals; Center: with subset-supset based thick intervals; Right: with lower-upper bounds based thick intervals

Test-case 1 .

 1 Thick translation. Consider the thick set Y = [Y -, Y + ] where Y -, Y + are two disks with center (0, 0) and with radius r -= 1 and r + = 2, respectively. Consider the thick function [f ] (x) = x -[v] where [v] = [0.7, 1.3] × [-0.02, 0.02]. Figure A.8 represents an approximation of X = [f ] -1 ( Y ) using three types of intervals: the classical intervals (left), the thick intervals with a subset-supset representation (center) and the thick intervals defined by lower-upper interval bounds (right). These results have been obtained with ε = 0.1 and the frame box corresponds to [-2, 4] × [-2, 4]. These figures have been generated in 2.1sec for the left figure; 0.21sec for the centered figure and 0.19 for the right figure.

Figure

  Figure A.11: Thick set X representing the communication region. The dark zone (blue) corresponds to position where the robot cannot communicate with all marks.In the clear zone (red), the robot is able to communicate with all marks.

Example 3 .

 3 Figure B.1: The space is partitioned into three zones : Z + is the dark zone, Z -∩ Z + is the penumbra and Z -is the clear zone

Figure B. 4 :

 4 Figure B.4: This configuration for the tubes [ϕ] (•) and [ψ] (•) allows us to conclude that [z] is dark

). Theorem B. 4 .

 4 Given x(•) ∈ [x](•), z ∈ [z] and t 1 , t 2 ≥ 0. We have: ϕ (z, x(t 1 )) .ϕ (z, x(t 2 )) ≤ 0 ∧ ∀t ∈ [t 1 , t 2 ], ψ (z, x(t)) ≤ 0 ⇒ ¬T dark ([x] (•) , [z]) . (B.18)Proof. The proof is by contradiction. Assume thatT dark ([x] (•) , [z]) is true, i.e., ∀t ≥ 0, 0 / ∈ [ϕ] ([z] , [x] (t)) ∨ [ψ] ([z] , [x] (t)) ⊂ R + Since x(•) ∈ [x](•) and z ∈ [z],we get ∀t ≥ 0, 0 = ϕ (z, x (t)) ∨ ψ (z, x (t)) > 0.

Figure B. 6 :

 6 Figure B.6: Online estimated trajectory (plain) and true trajectory (dotted line)

  An interval [x] containing a single element is called a degenerate interval or a singleton. It is denoted with brackets {x}. The width w of an interval is defined by w ([x]) = x + -x -. 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]:

  • , i.e. [θ i ] ∈ [θ i -∆θ, θ i +∆θ], are θ 1 = 45 • and θ 2 = 120 • .

Table 2 .

 2

	1: Center and radius of rings used in figure 2.10
	c 1	c 2	c 3

  Characterization of X using S defined in Equation (2.56). Circle c 1 is drawn in white, c 2 in red and c 3 in cyan. parts of the search space inside or outside the solution set X associated with S. The algorithm is given in Algorithm (2.1). Step 1 initializes a list L containing all boxes to be studied. Step 2 takes one box [x] in L. At Step 3, the separator S is then called to contract [x] into two boxes [x in ] and [x out ]. Step 4 stores ¬S in ([x]

	8
	4
	0
	-4
	-8
	-8	-4	0	4	8
	Figure 2.10:

Table 3 .

 3 1: Range and bearing measurements of the three landmarks

	landmarks	1	2	3
	m 1	6	-2	-3
	m 2	12	-5	10
	[y 1 ](m)	[11, 12]	[8, 10]	[5, 7]
	[y 2 ]( • )	[14, 100] [-147, -75] [63, 150]

Table 3

 3 .2 shows bounds used to quantify the errors on sensor readings. The map M is composed of 280 landmarks, and 90 observations have been done during the whole mission.

	data description	uncertainty	unit
	ψ(0) initial heading	[-10, 10]	deg
	v(t) linear speed	[-0.05, 0.05] ×2	m.s -1
	ω(t) angular velocity	[-0.001, 0.001]	deg.s -1
	ρ(t) range	ρ(t) • [-0.01, 0.01] m
	ψ(t) azimuth	[-1, 1]	deg

Table 3 .

 3 2: Uncertainties on data used for the application.

		400						#1
								#3
								#5
								#6
	w([p])(m)	200 300						#7 #8 #9 #10 #11
		100
								t(s)
		100	200	300	400	500	600	700

Table 3 .

 3 3: Time needed to contract the whole trajectory for each iteration. #min (resp. #max) is the smallest (resp. greatest) number of elements of [m](t i ) among all measurements. #ok denotes the number of good associations. The computing time is nearly constant.

	# time(s) #min #max #ok
	1	2.61	202	230	0
	2	2.65	27	81	0
	3	2.99	22	75	0
	4	2.77	15	73	0
	5	2.62	11	71	0
	6	2.48	4	71	0
	7	2.33	2	68	0
	8	2.02	1	29	10
	9	2.61	1	6	64
	10	2.44	1	3	82
	11	2.41	1	2	88
	12	2.40	1	1	90

4 Localization in an Unstructured Environment . . . . . 87 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . 89

	Chapter 4
	Shape-based Localization

Contents 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2 Shape Registration as a Projection . . . . . . . . . . . . 81 4.2.1 Projection of separators . . . . . . . . . . . . . . . . . . 81 4.2.2 Shape registration . . . . . . . . . . . . . . . . . . . . . 83 4.2.3 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3 Minkowsky Sum and Difference . . . . . . . . . . . . . . 85 4.3.1 Minkowski difference . . . . . . . . . . . . . . . . . . . . 85 4.3.2 Minkowski addition . . . . . . . . . . . . . . . . . . . . 86 4.

Table 6 . 1 :

 61 6.7: Illustration of the contraction between timest 1 and t 2 . The initial tube, in blue, is contracted into the red one. Shape registration. The two thick shapes, depicted on Figure

		[x](t 1 )	[x](t 2 )	[p 12 ]	w ([p 12 ])
	before	[-156.7, -154.3] ×[31.2, 36.6]	[-200, -174.7] ×[92.6, 107.7]	[-43.3, -20.3] ×[61.2, 71.1]	25.04
	after	[-156.7, -154.3] ×[31.2, 36.6]	[-192.1, -181.8] ×[93.8, 105.9]	[-37.8, -25.2] ×[61.2, 71.1]	16.02

  is the lower bound of the interval [a i ] + and ub([a i ] -) is the upper bound of [a i ] -. We thus get Algorithm A.3. The same type of algorithm applies to test

  1 , t 2 ] , ϕ (z, x(t 0 )) = 0 ∧ ψ (z, x(t 0 )) ≤ 0 which is not consistent with the previous proposition.Penumbra test. Given t 1 ,t 2 ≥ 0, two trajectories x a (•) , x b (•) ∈ [x] (t)and inclusion functions [ϕ],[ψ] for ϕ,ψ. Define the penumbra test T penumbra ([z]) as From Theorems B.3 and B.4, we get that T penumbra implies that neither T clarity nor T dark will able to conclude anything for any subbox of [z].

	
	  
	  

∀t, [ϕ] ([z] , x a (t)) = 0 ∨ [ψ] ([z] , x a (t)) > 0 [ϕ] [z] , x b (t 1 ) • [ϕ] [z] , x b (t 2 ) ≤ 0 ∧ ∀t ∈ [t 1 , t 2 ], [ψ] [z] , x b (t) ≤ 0. (B

.19)

Acknowledgments

The thin set X (which is not a ring and even not connected) drawn in red on the left belongs to the thick set [X] on the right. [X] contains all sets that enclose the red ring and which do not intersect the blue zone.

2006, [START_REF] Kreinovich | Interval Methods for Data Fitting under Uncertainty: A Probabilistic Treatment[END_REF] for the linear case or by allowing a given number of constraints to be unsatisfied [START_REF] Carbonnel | Q-intersection Algorithms for Constraint-Based Robust Parameter Estimation[END_REF]. A set inversion problem can be written as .28) where [Y] is a thick set and F is a set of functions, i.e, a subset of F(R n , R m ). For instance, F can be the set of all periodic functions, the set of positive functions, or a tube seen as an interval of functions from R to R m .

The set X is said to be a feasible solution if ∃f ∈ F,∃Y ∈ [Y] , X = f -1 (Y).

(5.29)

The set of all feasible solutions is not a thick set in general as illustrated by the following example. Solving a thick set-inversion problem will amount to finding the smallest thick set which encloses all feasible solutions for Equation (5.28).

Example 5.2. Let f (x) = x 2 1 + x 2 2 be a thin function from R 2 to R, and Y be a set such that [2,3] ⊂ Y ⊂ [START_REF]Density probability function of random variable x and[END_REF]4]. If X = f -1 (Y), we have: [START_REF]Density probability function of random variable x and[END_REF]4]) .

(5.30)

Now, all feasible X correspond to centered rings and it is clear that the inclusion (see Equation (5.30)) contains other types of sets as illustrated by Figure 5.5.

Set inversion theorem

The set inversion problem which involves an uncertain set and an uncertain function f is solved using the following theorem: of 2π. As illustrated by Figure 6.9, the scalar measurement y(t), assumed to belong to [y(t)] = [y -(t), y + (t)], can be cast into the interval shape [Z] = [Z -, R 2 ], with

which is a disk of center (0, 0), i.e. expressed in the robot frame. The impact area ∆Z is defined by

The quantities u 1 , u 2 and z are measured every 0.1 s with an error of 0.001 m.s -1 , 0.001 rad.s -1 , and 0.01 m, respectively. The initial state, taken as x = (0, 0) T , is assumed to be known. Fig. 6.10 provides a simulation of the robot moving inside an unknown map. This map is composed by segments but this is not required by the method. The shape of the map could be arbitrary and no parametric representation of the map is needed.

This problem is a particular case of the generic formulation of Equation (6.1). The translation invariant case can be expressed in terms of Minkowsky sums. It is a generalization of the localization problem of Chapter 4. It can be described by the following CN: which can be solved with Algorithm 6.3.

An illustration of the interval propagation method is depicted in Figure 6.12. The propagation process stops after 4 iterations in less than 40 seconds. The code of [Jaulin, 2011] provided by the author runs in around 90 seconds on the same computer.

Figure 6.12 shows the computed tube with the true path in black and corresponds to contractions of the inner approximation of the map. The true map, defined with segments, is superimposed to show the accuracy of the inner approximation of the map.

Figure 6.11 shows how Equation ( 3) of CN 6.24 is used to contract the position with respect to the map.

Chapter 7 Conclusion

This PhD thesis introduces new methods based on Interval Analysis, able to deal with constraints that involve uncertain sets. These types of constraints are mainly encountered in robotic applications in order to handle unstructured environments.

We proposed and illustrated a generic method that is able to solve the SLAM problem when measurements are shapes, i.e., subsets of R 2 of R 3 .

For that purpose, we introduced the shape registration and carving problem, where tools, algorithms, and applications have been presented along the chapters of this manuscript.

Summary of the Contributions

The key points of this thesis are now summarized. In Chapter 2, interval analysis is introduced. The main point of this chapter is the introduction of separators. Coupled with a Paver, which is a branching algorithm, separators allow to compute an inner and an outer approximation of a set defined by constraints. Their uses made it possible to characterize rigorously the pessimism introduced by intervalbased algorithms. This work was published in [START_REF] Jaulin | Introduction to the algebra of separators with application to path planning[END_REF] before the beginning of this thesis.

Now complex sets are defined by a combination of elementary ones, for which atomic separators are available. In Chapter 3, new theorems have been introduced in order to define a separator dedicated to the polar constraint, which is one of the important atomic constraints used in the context of localization. This constraint, associated with the change between Cartesian to Polar coordinates is often met in robotics, each time range and bearing measurements are involved. This work was published in [Desrochers and Jaulin, 2016a].

Based on the notion of projection of separators, the shape registration problem is presented in Chapter 4. this problem consists of estimating the smallest set of

A.1 Introduction

In this thesis, the set inversion problem X = f -1 (Y), where both f and Y are uncertain has been introduced. In the case where f depends on a parameter, a method has been given but it relies on projection algorithms that bisect along the parameter space.

In this chapter, we assume that f ∈

We also assume that the uncertain set Y belongs to an interval of sets, denoted by Y . We say that [f ] is a thick function and that Y is a thick set. As illustrated by the following example, existing interval methods can still be used to deal with this type of uncertainties, but they accumulate on a thick boundary which is called the penumbra. This accumulation makes classical interval methods inefficient, since they spend most of the computation time to test tiny boxes that are inside the penumbra.

Example A.1. Consider the set inversion problem

We assume that f is uncertain and that we only know that for all x

. Note that f is not necessarily a circular paraboloid, and may correspond to any weird function satisfying the enclosure condition. Since, for all x,

Using a classical interval arithmetic [Moore, 1979], we can easily test if a box

] is inside or outside the solution set X: Now, when dealing with practical applications, the penumbra often exists as for instance when we want to characterize the zone that has actually been explored by a robot, detailed in the Appendix B, in case of partial observability [START_REF] Grimson | Recognition and localization of overlapping parts from sparse data in two and three dimensions[END_REF] or in the shape carving problem introduced in Chapter 5. Characterizing the penumbra from inside will allow us to save computing time, but also to make the difference between the uncertainty due to the computation and that due to the initial uncertainties of the input parameters. The objective of this chapter is to extend set inversion to the thick case (where a penumbra exists) and to show how to conclude that a box is inside the penumbra without bisections Denote by (F(R n , R m ), ≤) the set of all functions from R n → R m equipped with the order relation ≤ defined as follows .6) the set F(R n , R m ) is a lattice where the meet and the join are defined by

and

For such a thick function [f ] , there exist two functions f -and, called the lower bound and the upper bound such that

is said to be thin and corresponds to a singleton ofF(R n , R m ), or equivalently to a classical function from R n to R m .

Remark A.1. The class of thick functions is not so restrictive. For instance, all set-valued functions of the form

where [a] is a box and where f is continuous with respect to a, are thick functions.

If now the box [a] is replaced by a disconnected set or when the function f is not scalar anymore, the function F (x) has no reason to be a thick function.

Thick set inversion problem

A thick set inversion problem can be written as

where Y is a thick set and [f ] is a thick function.

The set X is said to be a feasible solution if the subset bound is empty, i.e., a subset-supset representation may contain more intervals. For instance, if [b] = [4,5], we have 3,7], [6,8] . The corresponding lower and upper interval bounds are represented on FigureA.2 (right). Note that when the subset bound is not empty (as for the thick box a in red), no pessimism is added and both representations are equivalent.

Due to this pessimism, we will prefer to use a representation based on lowerupper bounds instead of the notation based on the subset-supset bounds. As already seen for thick sets (see (5.16)), set membership operations such as the union or the intersection can easily be extended to thick intervals. An extension for all classical operators of interval arithmetic is also valid. More precisely, if ∈ {+, -, •, ∩, , . . . }, we define

For instance, if a = [START_REF]Density probability function of random variable x and[END_REF]2], [4,5] and b = [3,7], [6,8] , we have 3,7], [4,5] + [6,8] = [4,9], [10, 13] a ∩ b = max( [START_REF]Density probability function of random variable x and[END_REF]2], [3,7]), min( [4,5], [6,8] [3,7]), max( [4,5], [6,8]) = [START_REF]Density probability function of random variable x and[END_REF]2], [6,8] .

An interpretation for these formulas is the following:

The following proposition shows how to compare, from a practical point of view, two thick intervals.

we have

Thus, the inclusion is true for all Note that using the subset-supset bounds, we could not reach this conclusion. Indeed, the subset-supset approximation of b is

The interval [b] = [6,6] is consistent with this inclusion and does not intersect

Thick boxes. Denote by IR n the set of all boxes of R n . A thick box x is a set of boxes of IR n which can be defined as

where [x -] and [x + ] are boxes of R n . The set of thick boxes of R n is denoted by IIR n . A thick box can be seen as an interval of boxes, i.e., an interval of intervals of R n . This is illustrated by

Since the two box bounds of [x -] and [x + ] are boxes of R n , we could decompose them as the Cartesian product of n intervals:

We define the ith component x i of the thick box [x -] , [x + ] as the thick interval

The following proposition will allow us to compare two thick boxes, with respect to the inclusion, from their interval components.

Proof. This proof is a direct consequence of Proposition A.1 and of the fact that

(A.23)

A.4 Thick Set Inversion

This section generalizes to set inversion algorithm [Jaulin and Walter, 1993c] to the thick case, as defined by Theorem 5.1.

A.4.1 Set inversion

Given a function f from R n to R m and a (thin) set Y ⊂ R n , solving the set inversion problem, denoted by X = f -1 (Y), is classically performed using an inclusion function [f ] : IR n → IR m of f [Moore, 1979], i.e., an interval function such that a

Most algorithms for set-inversion decompose R n into boxes [Moore, 1992][Jaulin andWalter, 1993c].

If it satisfies none of the previous tests, it is bisected until it becomes too small. A possible implementation for set inversion is given by Algorithm A.1 which is called Sivia (Set Inversion Via Interval Analysis) [Jaulin and Walter, 1993c].

When the algorithm terminates, we have [Jaulin and Walter, 1993c]

For the thick case, we have a thick function [f ] from R n to IR m and a thick set Y ∈ IP(R n ). We want to compute an approximation of the set of all sets Algorithm A.1 Set inversion algorithm: Sivia

perpendicularly to its largest side and push the two resulting boxes inL 6: end if 7: end while

We propose to compute an approximation of X by decomposing R n into three subsets: the clear zone X -, the penumbra X + \X -and the dark zone R n \X + . In our approach, a paver performs the decomposition of R n into boxes and a thick extension of an inclusion function is used to classify boxes.

A.4.2 Thick inclusion function

The function f :

tiny ε, log ε -1 is never large if we want the results in a reasonable time. It is thus difficult to get a complete experimental validation of our formula.

Enclosure

The algorithm computes guaranteed inner and outer approximations of the solution set of the thick set inversion problem. This is asserted by the following theorem. Theorem A.2. The algorithm returns an approximation of the thick set inversion problem X -, X + = [f ] -1 ( Y ) under the form of 3 subpavings (i.e., union of boxes): L clear , L penumbra , L dark . This approximation satisfies

Proof. To prove the Theorem, we need to show that, for a box [x], we have (see

Let us first prove (i). The left-hand side of (i) implies that

, and we show that x ∈ X -. Since

) and the previous formula implies.

The same reasoning applies to prove (ii). For (iii), assume that the left-hand side of (iii) is satisfied. Take one x ∈ [x] , the quantity f ([x]) becomes a singleton in IR n , i.e., a box of R n . We have

Thus, from (5.32), we get [x] ⊂ X + \X -.

Convergence

We now provide some convergence properties of our algorithm. We need first to define the convergence of a thick inclusion function f . This convergence can be interpreted as the continuity of the thick function f ([x]) around intervals [x] which are degenerated (i.e., the box [x] is a singleton).

Definition A.1. The thick inclusion function f ([x]) for [f ]([x]

) is said to be convergent if for all a ∈ R n , for all sequences of boxes [x] (k) and [y] (k), we have

where d H is the Hausforff distance between compact sets [Berge, 1963].

To make the link with morphological operations, the lower set X -corresponds to an erosion and X + to a dilation.

Test-case 2. Tolerable-United solution sets. Consider the interval linear system [START_REF] Kreinovich | Interval Methods for Data Fitting under Uncertainty: A Probabilistic Treatment[END_REF]]

The left-hand side corresponds to a thick function and the right-hand side corresponds to a thin set. The solution set X = X -, X + has for subset bound the tolerable solution set X -for supset bound the united solution set X + [START_REF] Goldsztejn | On the approximation of linear AE-solution sets[END_REF]] . Some techniques have been developed to approximate these sets [Shary, 1995] in the linear case. They are mainly based on the Kaucher interval arithmetic [Kaucher, 1980, Goldsztejn, 2005] and may be used to find boxes inside the penumbra. Now, these methods have mainly been developed to deal with linear interval problems and cannot be used to find boxes inside the penumbra for general nonlinear problems as for the two following test-cases. For this example, the thick set-inversion algorithm provides the paving of Figure A.9, in less than 0.4 sec, for ε = 0.01.

Test-case3. Parameter estimation.Consider the parametric model

where x = (x 1 , x 2 ) is the parameter vector and t ∈ R is the time. At time t i , we collect measurements y i with some interval uncertainties as written in Table A.2. Note that one of the main difficulties of this problem is that uncertainties exist on the independent variable (here the time) [Jaulin andWalter, 1999, Cerone, 1991].

In our formulation, the uncertainty of the t i is stored inside the model under the The set of all feasible parameter vectors is

If we define the thick function Test-case 4. Communication area.Consider p marks m (i) located at position (m 1 (i) , m 2 (i)) given by Table A.3 and a robot at the position x = (x 1 , x 2 ). i 1 2 3 4 m 1 (i) 1 ± 0.5 10 ± 0.5 10 ± 0.5 -2 ± 0.5 m 2 (i) 3 ± 0.5 -1 ± 0.5 6 ± 0.5 -5 ± 0.5

Table A.3: Location of the marks

The robot is able to communicate with the mark m (i) if its distance to the mark is smaller than 10m, i.e., if xm(i) < 10. The communication is not possible if the distance is larger than 20m. With a distance inside [10,20], the communication is uncertain. The set of all positions for the robot such that the robot can communicate with all marks is a thick set defined by

Our thick inversion algorithm provides in less than 0.3 seconds, the paving represented on Figure A.11. In the example, the clear and the dark boxes could have been obtained using existing interval algorithms. But these methods have to bisect

To make the link with the shape registration and carving problem, the thick set inversion algorithm allows to deal efficiently with uncertain transformations that can be modeled by thick functions. The translation of an uncertain vector, that will be used in applications of chapter 5, fit this definition.

Appendix B

Guaranteed Area Explored

B.1 Introduction

In Appendix A, the thick set inversion problem for which the uncertain function is described by a thick function, has been defined. A new algorithm has been provided to efficiently characterize the penumbra and was illustrated with some examples. In Chapter 5 a generic method, based on the projection algorithms, has also been introduced. Now, this chapter deals with a new method able to characterize the area explored by a robot taking into account its trajectory uncertainty. When the dimension of the observer space is smaller that those of the working space, previous methods cannot be used. The temporal evolution of the vehicle need to be correctly handled in by the set inversion algorithm.

Behind the practical case of the characterization of the area explored by a robot, the exposed method will allow us to compute a set that encloses the sensor data taking into account its geometry, its positioning uncertainty and the movement of its carrier.

It is organized as follows. Section B.2 defines the notion of explored zone and proposes a formalization of the problem. Section B.3 extends the formalism in order to take into account the uncertainty associated with the trajectory of the robot. Section B.4 gives the new algorithm which encloses the explored zone Z between two subpavings (union of boxes). An experiment involving the underwater robot Daurade is treated in Section B.5 in order to validate the feasibility of the approach. A conclusion is then given on Section B.6.

Note that, this chapter is extracted, almost with any changes, from [Desrochers and Jaulin, 2017b]. The original paper was written before the thick set formalism was clearly defined, but it helped us to better understand and define the notion of penumbra.

B.2 Problem Statement

We consider a robot moving inside an unknown environment in a dead reckoning manner, i.e., using the proprioceptive sensors only. As it is the case for most industrial underwater robots, the robot is equipped with some exteroceptive sensors (such as sonars or cameras) that are used exploration only and not for navigation.

After the mission, it has to find which part of the environment it has explored, taking into account the uncertainty on the localization. We assume that the interpretation of exploration sensors is not reliable enough to allow us using SLAM techniques such as in [Leonard et al., 1992] or [Frese, 2006]. When the mission is short and the quality of the proprioceptive sensors is good, occupancy map techniques [Elfes, 1987] can be used to mark all points that have been observed and a probability of being explored can be associated with each part of the space [Paull, 2013]. Now, in practice, due to the state noise, the prediction of the location of the robot is getting less accurate, and we have to take into account these uncertainties properly. We are in a typical situation where the uncertainty is combined with inaccuracy [START_REF] Dubois | Possibility Theory[END_REF]. The unknown variable of our problem is the explored set Z, and characterizing uncertain sets with classical probabilistic methods requires elaborated mathematical tools such as random sets [Molchanov, 2005a]. These random sets have already been used in the context of robot mapping (see [START_REF] Mullane | Random Finite Sets for Robot Mapping and SLAM[END_REF]) but the approach is limited to finite sets. In a set membership approach, an uncertain set Z can be bracketed by two sets Z -and Z + such that Z -⊂ Z ⊂ Z + . This representation is particularly adapted to represent uncertain maps [START_REF] Langerwisch | Building variable resolution occupancy maps assuming unknown but bounded sensor errors[END_REF]. Interval-type uncertainties can easily be propagated through nonlinear functions (see e.g., [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF]) or nonlinear state equations, as shown in [START_REF] Collins | The Reachand-Evolve Algorithm for Reachability Analysis of Nonlinear Dynamical Systems[END_REF] to

As motivated by the following examples, we may consider three different types of exploration.

1. Patch exploration [Drevelle et al., 2013]. The visible sets V (x) have a dimension q, as for Example 1. In this case, dim(Z)=q and it will be possible to compute an inner and an outer approximation of Z. This type of exploration exists when the sensor is a camera: a 3D zone is explored for each t.

2. Sweep exploration. The visible sets V (x) have a dimension q -1, as for Example 2. Thus, Z has a dimension of q. Again, it will be possible [START_REF] Olivier | General inner approximation of vector-valued functions[END_REF] to bracket Z from both inside and outside. This is the case when the robot is equipped with a side scan sonar, a multi-beam echo sounder or an airborne LiDAR.

3. Pen exploration. The dimension of V (x) is smaller than q -2, as for Example 3, and thus dim(Z) ≤q -1. It will neither be possible to get an inner approximation for Z nor to prove that a given point is actually inside Z. This is the case when the robot is equipped with a laser range-finder and takes one range measurement every second.

When the whole trajectory x (t) is known, the set Z is clearly defined. Different techniques could be thought in order to characterize the set Z. For instance, we could use subpaving-based methods [START_REF] Jaulin | Applied Interval Analysis[END_REF] or occupancy grid approaches. Now, in this paper, the trajectory of the robot is not known precisely and the set Z becomes uncertain which makes the problem much more difficult [Paull, 2013].

B.3 Uncertain Explored Zones

When the trajectory x (•) is uncertain, the explored zone Z cannot be approximated with an arbitrary accuracy. In a probabilistic context, we can associate to each z ∈ R q a probability of being explored. Now, even if this probability can be estimated using Monté-Carlo methods, due to the large (here infinite) dimension of the set of trajectories to be explored, the computational burden is high which makes probabilistic methods not so attractive. Denote by X (•), the set of all feasible trajectories, i.e., the trajectories x (•) consistent with (B.1) and with the initial condition x (0). It is important not to make the confusion between x (•) which is a trajectory and x (t) which is a vector of R n . We define the two following sets

The set Z -is called the certainly explored zone or the clear zone. It corresponds to the set of all points z of the environment that have certainly been seen by the the penumbra boxes which are neither dark nor clear and (d) the boxes which have not been classified yet. When a box cannot be classified, it is bisected into two non overlapping subboxes. To ensure the convergence of the algorithm, boxes which have a width smaller than a small given value are not bisected.

To classify boxes (clear, dark, penumbra), we need first to get a guaranteed enclosure of the set of feasible trajectories X (•). This can be done using interval guaranteed integration [START_REF] Wilczak | Cr-Lohner algorithm[END_REF], [Tucker, 2002], [Revol et al., 2005]. These techniques make it possible to compute a tube (i.e., a function [x] (•) mapping R into IR n , the set of boxes of R n ) which encloses the set of trajectories X (•). In what follows, we shall thus assume that the trajectory x (•) is inside a tube denoted by [x] (•). As a consequence, from (B.5), we have

B.4.1 Notion of logic

The presentation of the tests that will be used to classify the boxes of the workspace requires some notions on logic, involving quantifiers. Some useful rules are now presented.

(B.10)

B.4.2 Clarity test

The clarity test aims at proving that a box [z] has been explored whatever the uncertainties are. This test is based on the following theorem. Theorem B.1. Given an interval [t 1 , t 2 ] ⊂ R + and a tube [x] (•), we have :

Clarity test. Given t 1 ,t 2 ≥ 0 and inclusion functions [ϕ],[ψ] for ϕ,ψ, we define the clarity test as:

An illustration of this proposition is given by Figure B.3. We choose the pairs (t 1 , t 2 ) as in the figure, i.e., such that t 2 -t 1 are small and such that

Note that the test also concludes for the pair (t 5 , t 6 ) but fails for (t 3 , t 4 ). Since here, the test concludes for at least one pair, we get that the whole box [z] has been explored.

B.4.3 Dark test

The dark test aims at proving that no point in the box [z] has been explored whatever the uncertainties are. The dark test is based on the following theorem.

Proof. The proof is by contradiction, i.e., we show that assuming z ∈ Z + implies that left-hand side of (B.14) is false. We have

which corresponds to the negation of left-hand side of (B.14).

Dark test. Given two inclusion functions [ϕ],

[ψ] for ϕ, ψ. Define the dark test as

From Theorem B.2, we have

B.4.4 Penumbra test

Proving that a point z belongs to the penumbra ∆Z is difficult and not so useful.

The main interest to have an inclusion test for the penumbra is to limit the computation burden. Indeed, if [z] ⊂ ∆Z then we will never be able to prove that [z] ⊂ Z -or to prove that [z] ∩ Z + = ∅, and thus there is no need to bisect. Now, instead of proving that [z] ⊂ ∆Z, we propose here to test either the darkness test or the clarity test always fail for all subbox of [z]. In such a case, [z] will not be bisected. In such a situation, we shall say that [z] satisfies the penumbra test. Assessment of the covered area is usually done manually by an operator which looks at the sonar images. The proposed algorithm can be used to validate the mission plan or, at the end of the mission, to check the area to be explored has indeed been covered.

The horizontal kinematic model the robot is taken as

where (x 1 , x 2 ) corresponds to the 2D coordinates of the center of the robot expressed in an absolute inertial frame, ψ is the heading and v is the horizontal speed vector of the robot expressed in its own coordinate system. The aperture angle of the side scan sonar is taken as α = 80 • . This means that, with a flat seabed assumption, if the robot has an altitude a then the sonar is able to sense the part of the environment which is perpendicular to the robot and at a distance less than = a • sin α. The horizontal online estimated trajectory is depicted on Figure B.6. Note that due to the fact that the controller uses this trajectory for control, it looks perfect, which is not the case for the true trajectory (dotted line) directly measured by an USBL (Ultra-Short Base Line). Now, in our experiment, this true trajectory is neither used for control nor for estimation of the explored zone. It is only used for the validation of the results.

We assume that the initial position is x 1 = x 2 = 0. Taking the interval uncertainties into account and using a guaranteed integration of the state equation (B.20), we are able to compute a tube which contains the true trajectory. At the end of the mission the position error is around 17 meters. On this experiment, our method provides in less than 5 minutes an estimation of the explored zone as given by Figure B.7. This result is consistent with the true trajectory obtained by the USBL and by an observation of the sonar images.

Consider the same experiment, but, instead of having a sweep exploration, we assume that at each t the visible set V(t) is a disk of radius a • sin α (instead of a segment). In this case, the explored zone can still be obtained, but we do not need anymore to use the continuity of the trajectory. At each t ∈ {1, 2, . . . , m} where m = 2760, the visible set (i.e., the part of the bottom which is seen by the sonar of the robot) is a disk of radius 50 meters around the robot, the position of which is not exactly known. The explored zone X corresponds to the union of all patches that have been seen:

The complementary set of X is

where f (x) = (f 1 (x), . . . , f m (x)). The function f (x) is consistent with the intervals [a 1 ] (t) and [a 2 ] (t) containing the positions of the robot iff

B.6 Conclusion

This chapter has presented a new method to characterize the zone Z explored by a mobile robot. The method only uses the proprioceptive sensors for localization, which makes the approach reliable. The method is indeed not sensitive (i) to

Appendix C

Publications

Journal Articles

List of Tables