Andrea Cherubini
email: andrea.cherubini@lirmm.fr

François Chaumette
email: francois.chaumette@inria.fr

Visual Navigation of a Mobile Robot with Laser-based Collision Avoidance

Keywords: Vision-based Navigation, Visual Servoing, Collision Avoidance, Integration of vision with other sensors

In this paper, we propose and validate a framework for visual navigation with collision avoidance for a wheeled mobile robot. Visual navigation consists of following a path, represented as an ordered set of key images, which have been acquired by an on-board camera in a teaching phase. While following such path, the robot is able to avoid obstacles which were not present during teaching, and which are sensed by an on-board range scanner. Our control scheme guarantees that obstacle avoidance and navigation are achieved simultaneously. In fact, in the presence of obstacles, the camera pan angle is actuated to maintain scene visibility while the robot circumnavigates the obstacle. The risk of collision and the eventual avoiding behaviour are determined using a tentacle-based approach. The framework can also deal with unavoidable obstacles, which make the robot decelerate and eventually stop. Simulated and real experiments show that with our method, the vehicle can navigate along a visual path while avoiding collisions.

Introduction

A great amount of robotics research focuses on vehicle guidance, with the ultimate goal of automatically reproducing the tasks usually performed by human drivers [START_REF] Buehler | [END_REF], Zhang et al. 2008, Nunes et al. 2009, Broggi et al. 2010]. In many recent works, information from visual sensors is used for localization [Guerrero et al. 2008, Scaramuzza andSiegwart 2008] or for navigation [Bonin-Font et al. 2008, López-Nicolás et al. 2010]. In the case of autonomous navigation, an important task is obstacle avoidance, which consists of either generating a collision-free trajectory to the goal [Minguez et al. 2008], or decelerating to prevent collision when bypassing is impossible [Wada et al. 2009]. Most obstacle avoidance techniques, particularly those that use motion planning [Latombe 1991], rely on the knowledge of a global and accurate map of the environment and obstacles.

Instead of utilizing such a global model of the environment, which would infringe the perception to action paradigm [Sciavicco and Siciliano 2000], we propose a framework for obstacle avoidance with simultaneous execution of a visual servoing task [Chaumette and Hutchinson 2006]. Visual servoing is a well known method that uses vision data directly in the control loop, and that has been applied on mobile robots in many works [Mariottini et al. 2007, Becerra et al. 2011, López-Nicolás and Sagüés 2011, Allibert et al. 2008]. For example, in [Mariottini et al. 2007] the epipolar geometry is exploited to drive a nonholonomic robot to a desired configuration. A similar approach is presented in [START_REF] Becerra | [END_REF],

where the singularities are dealt with more efficiently. The same authors achieve vision-based pose stabilization using a state observer in [López-Nicolás and Sagüés 2011]. Trajectory tracking is tackled in [Allibert et al. 2008] by integrating differential flatness and predictive control.

The visual task that we focus on is appearance-based navigation, which has been the target of our research in [START_REF] Šegvić | [END_REF], Cherubini et al. 2009, Diosi et al. 2011]. In the framework that we have developed in the past1 , the path is a topological graph, represented by a database of ordered key images. In contrast with other similar approaches, such as [START_REF] Booij | [END_REF]], our graph does not contain forks. Furthermore, as opposed to [START_REF] Royer | [END_REF], Goedemé et al. 2007, Zhang and Kleeman 2009, Fontanelli et al. 2009, Courbon et al. 2009], we do not use the robot pose for navigating along the path. Instead, our task is purely image-based (as in [START_REF] Becerra | [END_REF]), and it is divided into a series of subtasks, each consisting of driving the robot towards the next key image in the database. More importantly, to our knowledge, appearance-based navigation frameworks have never been extended to take into account obstacles.

Obstacle avoidance has been integrated in many model-based navigation schemes. In [Yan et al. 2003], a range finder and monocular vision enable navigation in an office environment. The desired trajectory is deformed to avoid sensed obstacles in [START_REF] Lamiraux | [END_REF]]. The authors of [Ohya et al. 2008] use a model-based vision system with retroactive position correction. Simultaneous obstacle avoidance and path following are presented in [Lapierre et al. 2007], where the geometry of the path (a curve on the ground) is perfectly known. In [START_REF] Lee | [END_REF], obstacles are circumnavigated while following a path; the radius of the obstacles (assumed cylindrical) is known a priori. In practice, all these methods are based on the environment 3D model, including, for example, walls and doors, or on the path geometry. In contrast, we propose a navigation scheme which does not require neither the environment nor the obstacle model.

One of the most common techniques for model-free obstacle avoidance is the potential field method, originally introduced in [Khatib 1985]. The gap between global path planning and real-time sensor-based control has been closed with the elastic band [Quinlan and Khatib 1993], a deformable collision-free path, whose initial shape is generated by a planner, and then deformed in real time, according to the sensed data.

Similarly, in [Bonnafous et al. 2001, Von Hundelshausen et al. 2008], a set of trajectories (arcs of circles or "tentacles") is evaluated for navigating. However, in [Bonnafous et al. 2001], a sophisticated probabilistic elevation map is used, and the selection of the optimal tentacle is based on its risk and interest, which both require accurate pose estimation. Similarly, in [Von Hundelshausen et al. 2008], the trajectory computation relies on GPS way points, hence -once more -on the robot pose.

Here, we focus on this problem: a wheeled vehicle, equipped with an actuated pinhole camera and with a forward-looking range scanner, must follow a visual path represented by key images, without colliding with the ground obstacles. The camera detects the features required for navigating, while the scanner senses the obstacles (in contrast with other works, such as [START_REF] Kato | [END_REF], only one sensor is used to detect the obstacles). In this sense, our work is similar to [START_REF] Folio | A redundancy-based scheme to perform safe vision-based tasks amidst obstacles[END_REF], where redundancy enables reactive obstacle avoidance, without requiring any 3D model. A robot is redundant when it has more DOFs than those required for the primary task; then, a secondary task can also be executed. In [START_REF] Folio | A redundancy-based scheme to perform safe vision-based tasks amidst obstacles[END_REF], the two tasks are respectively visual servoing and obstacle avoidance. However, there are various differences with that work. First, we show that the redundancy approach is not necessary, since we design the two tasks so that they are independent. Second, we can guarantee asymptotic stability of the visual task at all times, in the presence of non occluding obstacles. Moreover, our controller is compact, and the transitions between safe and unsafe contexts is operated only for obstacle avoidance, while in [START_REF] Folio | A redundancy-based scheme to perform safe vision-based tasks amidst obstacles[END_REF], three controllers are needed, and the transitions are more complex. This compactness leads to smoothness of the robot behaviour. Finally, in [START_REF] Folio | A redundancy-based scheme to perform safe vision-based tasks amidst obstacles[END_REF], a positioning task in indoor environments is considered, whereas we aim at continuous navigation on long outdoor paths.

Let us summarize the other major contributions of our work. An important contribution is that our approach is merely appearance-based, hence simple and flexible: the only information required is the database of key images, and no model of the environment or obstacles is necessary. Hence, there is no need for sensor data fusion nor planning, which can be computationally costly, and requires precise calibration of the camera/scanner pair. We guarantee that the robot will never collide in the case of static, detectable obstacles (in the worse cases, it will simply stop). We also prove that our control law is always well-defined, and that it does not present any local minima. To our knowledge, this is the first time that obstacle avoidance and visual navigation merged directly at the control level (without the need for sophisticated planning) are validated in real outdoor urban experiments.

The framework that we present here is inspired from the one designed and validated in our previous work [START_REF] Cherubini | [END_REF]]. However, many modifications have been applied, in order to adapt that controller to the real world. First, for obstacle avoidance, we have replaced classical potential fields with a new tentacle-based technique inspired from [Bonnafous et al. 2001] and [Von Hundelshausen et al. 2008], which is perfectly suitable for appearance-based tasks, such as visual navigation. In contrast with those works, our approach does not require the robot pose, and exploits the robot geometric and kinematic characteristics (this aspect will be detailed later in the paper). A detailed comparison between the potential field and the tentacle techniques is given in [START_REF] Cherubini | [END_REF]. In that work, we showed that with tentacles, smoother control inputs are generated, higher velocities can be applied, and only dangerous obstacles are taken into account. In summary, the new approach is more robust and efficient than its predecessor. A second modification with respect to [START_REF] Cherubini | [END_REF] concerns the design of the translational velocity, which has been changed to improve visual tracking and avoid undesired deceleration in the presence of non-dangerous obstacles. Another important contribution of the present work is that, in contrast with the tentacle-based approaches designed in [Von Hundelshausen et al. 2008] and [Bonnafous et al. 2001], our method does not require the robot pose. Finally, the present article reports experiments, which, for the first time in the field of visual navigation with obstacle avoidance, have been carried out in real-life, unpredictable urban environments.

The article is organized as follows. In Section 2, the characteristics of our problem (visual path following with simultaneous obstacle avoidance) are presented. The control law is presented in full details in Section 3, and a short discussion is carried out in Section 4. Simulated and real experimental results are presented respectively in Sections 5 and 6, and summarized in the conclusion.

2 Problem Definition

General Definitions

The reader is referred to Fig. 1. We define the robot frame F R (R, X, Y ) (R is the robot center of rotation), image frame F I (O, x, y) (O is the image center), and camera frame F C (C, X c , Y c , Z c ) (C is the optical center). The robot control inputs are:

u = (v, ω, φ) .
These are, respectively, the translational and angular velocities of the vehicle, and the camera pan angular velocity. We use the normalized perspective camera model:

x = X c Z c , y = Y c Z c .
We assume that the camera pan angle is bounded: |ϕ| ≤ π 2 , and that C belongs to the camera pan rotation axis, and to the robot sagittal plane (i.e., the plane orthogonal to the ground through X). We also assume that the path can be followed with continuous v (t) > 0. This ensures safety, since only obstacles in front of the robot can be detected by our range scanner.

Visual Path Following

The path that the robot must follow is represented as a database of ordered key images, such that successive pairs contain some common static visual features (points). First, the vehicle is manually driven along a taught path, with the camera pointing forward (ϕ = 0), and all the images are saved. Afterwards, a subset (database) of N key images I 1 , . . . , I N representing the path (Fig. 1, right) is selected. Then, during autonomous navigation, the current image, noted I, is compared with the next key image in the database,

I d ∈ {I 1 , .
. . , I N }, and a relative pose estimation between I and I d is used to check when the robot passes the pose where I d was acquired.

For key image selection, as well as visual point detection and tracking, we use the algorithm presented in [START_REF] Royer | [END_REF]]. The output of this algorithm, which is used by our controller, is the set of points visible both in I and I d . Then, navigation consists of driving the robot forward, while I is driven towards

I d .
We maximize similarity between I and I d using only the abscissa x of the centroid of the points matched on I and I d . When I d has been passed, the next image in the set becomes the desired one, and so on, until I N is reached.

Obstacle Representation

Along with the visual path following problem, we consider obstacles which are on the path, but not in the database, and sensed by the range scanner in a plane parallel to the ground. We use the occupancy grid in in the grid, we use past readings, which are progressively displaced using odometry.

We use, along with the set of all occupied grid cells:

O = {c 1 , . . . , c n } ,
a set of drivable paths (tentacles). Each tentacle j is a semicircle that starts in R, is tangent to X, and is characterized by its curvature (i.e., inverse radius) κ j , which belongs to K, a uniformly sampled set:

κ j ∈ K = {-κ M , . . . , 0, . . . , κ M } .
The maximum desired curvature κ M > 0, must be feasible considering the robot kinematics. Since, as we will show, our tentacles are used both for perception and motion execution, a compromise between computational cost and control accuracy must be reached to tune the size of K, i.e., its sampling interval.

Indeed, a large set is costly since, as we show later, various collision variables must be calculated on each tentacle. On the other hand, extending the set enhances motion precision, since more alternative tentacles can be selected for navigation. In the simulations and experiments, we used 21 tentacles. In Fig. 2(b-e), the straight and the sharpest counterclockwise (κ = κ M ) tentacle are dashed. When a total of 3 tentacles is used, these correspond respectively to j = 2 and j = 3.

Each tentacle j is characterized by three classification areas (collision, dangerous central, and dangerous external), which are obtained by rigidly displacing, along the tentacle, three rectangular boxes, with increasing size. The boxes are all overestimated with respect to the real robot dimensions. For each tentacle j, the sets of cells belonging to the three classification areas (shown in Fig. 2) are noted C j , D j and E j . Cells belonging to the dangerous central set, are not considered in the dangerous external set as well, so that

D j E j = ∅.
The sets O, C, D and E are used to calculate the variables required in the control law defined in Section 3.1: in particular, the largest classification areas D and E are used to select the safest tentacle and its risk, while the thinnest one C determines the eventual necessary deceleration.

In summary, as we mentioned in Sect 1, our tentacles exploit the robot geometric and kinematic characteristics. Specifically, the robot geometry (i.e., the vehicle encumbrance) defines the three classification areas C, D and E, hence the cell potential danger, while the robot kinematics (i.e., the maximum desired curvature, κ M ) define the bounds on the set of tentacles K. Both aspects give useful information on possible collisions with obstacles ahead of the robot, which will be exploited, as we will show in Sect. 3, to choose the best tentacle and to eventually slow down or stop the robot.

Task Specifications

Let us recall the Jacobian paradigm which relates a robot kinematic control inputs with the desired task. We name s ∈ IR m the task vector, and u ∈ IR m the control inputs. The task dynamics are related to the control inputs by:

ṡ = Ju, ( 1
)
where J is the task jacobian of size m × m. In this work, m = 3, and the desired specifications are:

1. orienting the camera in order to drive the abscissa of the feature centroid x to its value at the next key image in the database x d , 2. making the vehicle progress forward along the path (except if obstacles are unavoidable),

3. avoiding collision with the obstacles, while remaining near the 3D taught path.

The required task evolution can be written:

ṡ * = ṡd -Λ s -s d , (2)
with s d and ṡd indicating the desired values of the task, and of its first derivative, and

Λ = diag (λ 1 . . . λ m )
a positive definite diagonal gain matrix.

Since we assume that the visual features are static, the first specification on camera orientation can be expressed by:

ẋ * = -λ x x -x d , (3)
with λ x a positive scalar gain. This guarantees that the abscissa of the centroid of the points converges exponentially to its value at the next key image x d , with null velocity there ( ẋd = 0). The dynamics of this task can be related to the robot control inputs by:

ẋ = J x u = j v j ω j φ u, (4)
where j v , j ω and j φ are the components of the centroid abscissa Jacobian J x related to each of the three robot control inputs. Their form will be determined in Section 3.2.

The two other specifications (vehicle progression with collision avoidance) are related to the danger represented by the obstacles present in the environment. If it is possible, the obstacles should be circumnavigated. Otherwise, the vehicle should stop to avoid collision. To determine the best behaviour, we assess the danger at time t with a situation risk function H : IR * + → [0, 1], that will be fully defined in Section 3.3.

• In the safe context (H = 0), no dangerous obstacles are detected on the robot path. In this case, it is desirable that the robot acts as in the teaching phase, i.e., following the taught path with the camera looking forward. If the current pan angle ϕ is non-null, which is typically the case when the robot has just avoided an obstacle, an exponential decrease of ϕ is specified. Moreover, the translational velocity v must be reduced in the presence of sharp turns, to ease the visual tracking of quickly moving features in the image; we specify this using a function v s that will be detailed in Section 3.4.

In summary, the specifications in the safe context are:

   ẋ = -λ x x -x d v = v s φ = -λ ϕ ϕ , (5)
with λ ϕ a positive scalar gain. The corresponding current and desired task dynamics are:

ṡs =   ẋ v φ   , ṡ * s =   -λ x x -x d v s -λ ϕ ϕ   . (6)
Using (4) we can derive the Jacobian relating ṡs and u:

ṡs = J s u =   j v j ω j φ 1 0 0 0 0 1   u. (7)
Note that matrix J s is invertible if j ω = 0, and we will see in Section 3.2 that this condition is indeed ensured.

• In the unsafe context (H = 1), dangerous obstacles are detected. The robot should circumnavigate them by following the best tentacle (selected considering both the visual and avoidance tasks as we will see in Section 3.3). This heading variation drives the robot away from the 3D taught path. Correspondingly, the camera pan angle must be actuated to maintain visibility of the database features, i.e., to guarantee (3). The translational velocity must be reduced for safety reasons (i.e., to avoid collisions); we specify this using a function v u , that will be defined in Section 3.5. In summary, the specifications in the unsafe context are:

   ẋ = -λ x x -x d v = v u ω = κ b v u , (8)
where κ b is the best tentacle curvature, so that the translational and angular velocities guarantee that the robot precisely follows it, since: ω/v = κ b . The current and desired task dynamics corresponding to (8) are:

ṡu =   ẋ v ω   , ṡ * u =   -λ x x -x d v u κ b v u   . (9)
Using (4) we can derive the Jacobian relating ṡu and u:

ṡu = J u u =   j v j ω j φ 1 0 0 0 1 0   u. ( 10
)
Matrix J u is invertible if j φ = 0, and we will see in Section 3.2 that this condition is also ensured.

• In intermediate contexts (0 < H < 1), the robot should navigate between the taught path, and the best tentacle. The transition between these two extremes will be driven by situation risk function H.

3 Control Scheme

General Scheme

An intuitive choice for controlling (1) in order to fulfill the desired tasks s u and s s would be:

u = J -1 ṡ * (11)
with:

s = Hs u + (1 -H) s s
and therefore (considering Ḣ = 0):

J = HJ u + (1 -H) J s =   j v j ω j φ 1 0 0 0 H 1 -H  
In fact, away from singularities of J, controller (11) leads to the linear system:

ṡ -ṡd = -Λ s -s d
for which, as desired, s d , ṡd are exponentially stable equilibria, for any value of H ∈ [0, 1] (since Λ is a positive definite diagonal matrix). Note that replacing (2) in ( 11), leads to the well known controller for following trajectory s d = s d (t), given in [Chaumette and Hutchinson 2007]:

u = -ΛJ -1 (s -s d ) + J -1 ṡd .
However, the choice of controller (11), is not appropriate for our application, since J is singular whenever:

Hj φ + (H -1) j ω = 0 (12)
This condition, which depends on visual variables (j φ and j ω ) as well as on an obstacle variable (H), can occur in practice.

Instead, we propose the following control law to guarantee a smooth transition between the inputs:

u = HJ -1 u ṡ * u + (1 -H) J -1 s ṡ * s (13)
Replacing this equation in ( 7) and ( 10), guarantees that controller (13) leads to convergence to the desired tasks ( 5) and ( 8):

ṡs = ṡ * s if H = 0 ṡu = ṡ * u if H = 1
and that, in these cases, the desired states are globally asymptotically stable for the closed loop system.

In Section 4, we will show that global asymptotic stability of the visual task is also guaranteed in the intermediate cases (0 < H < 1).

In the following, we will define the variables introduced in Section 2.4. We will show how to derive the centroid abscissa Jacobian J x , the situation risk function H, the best tentacle along with its curvature κ b , and the translational velocities in the safe and unsafe context (respectively v s and v u ).

Jacobian of the Centroid Abscissa

We will hereby derive the components of J x introduced in (4). Let us define: v = (v c , ω c ) the camera velocity, expressed in F C . Since we have assumed that the features are static, the dynamics of x can be related to v by:

ẋ = L x v
where L x is the interaction matrix of x [Chaumette and Hutchinson 2006]. In the case of a point of depth Z c , it is given by [Chaumette and Hutchinson 2006]:

L x = -1 Zc 0 x Zc xy -1 -x 2 y (14)
In theory, since we consider the centroid and not a physical point, we should not use ( 14) for the interaction matrix, but the exact and more complex form given in [START_REF] Tahri | [END_REF]. However, using ( 14) provides a sufficiently accurate approximation [Cherubini et al. 2009]. It also has the strong advantage that it is not necessary to estimate the depth of all points, using techniques such as those described in [Davison et al. 2007, De Luca et al. 2008, Durand et al. 2010]. Only an approximation of Z c , i.e., one scalar, is sufficient. In practice, we set a constant fixed value. This strategy has proved successful for visual navigation in [Cherubini et al. 2009].

For the robot model that we are considering, the camera velocity v can be expressed in function of u by using the geometric model:

v = C T R u with: C T R =         sin ϕ -X C cos ϕ 0 0 0 0 cos ϕ X C sin ϕ 0 0 0 0 0 -1 -1 0 0 0        
In this matrix, X C is the abscissa of the optical center C in the robot frame F R . This parameter is specific of the robot platform. Since C belongs to the robot sagittal plane, and since the robot is constrained on the ground plane, this is the only coordinate of C in F R required for visual servoing.

Then, multiplying L x by C T R , we obtain the components of J x :

j v = -sin ϕ+x cos ϕ Zc j ω = X C (cos ϕ+x sin ϕ) Zc + 1 + x 2 j φ = 1 + x 2 . ( 15
)
From ( 15) it is clear that j φ ≥ 1 ∀x ∈ IR; hence, J u is never singular (see ( 10)). Furthermore, it is possible to ensure that j ω = 0, so that J s is also invertible (see ( 7)). In fact, in (15) we can guarantee that j ω = 0, by setting Z c > X C 2 in the J x components. Indeed, condition j ω = 0 is equivalent to:

X C (cos ϕ + x sin ϕ) Z c + 1 + x 2 = 0 (16)
Since |ϕ| ≤ π 2 : cos ϕ + x sin ϕ ≥ -x, ∀x ∈ IR. Hence, a sufficient condition for ( 16) is:

x 2 - X C Z c x + 1 > 0 which occurs ∀x ∈ IR when X C Zc < 2.
In practice, this condition can be guaranteed, since X C is an invariant characteristic of the robot platform, and Z c is a tunable control parameter, which can be set to a value greater than X C 2 . Besides, the value of X C on most robots platforms is usually smaller than 1 m, which is much less than the scene depth in outdoor environments. In [Cherubini et al. 2009], we have shown that overestimating Z c does not jeopardize navigation.

On the other hand, we can infer from (15) that the singularity of controller ( 11), expressed by ( 12) can occur frequently. For example, whenever Z c is large, yielding j φ ≈ j ω , and concurrently H ≈ 0.5, J becomes singular. This confirms the great interest in choosing control scheme (13), which is always well-defined if Z c > X C 2 .

Situation Risk Function and Best Tentacle

To derive the situation risk function H used in (13), we first calculate a candidate risk function

H j ∈ [0, 1]
for each tentacle, as will be explained below. Each H j is derived from the risk distance of all occupied cells in the dangerous areas.

This distance is denoted ∆ ij ≥ 0 for each c i ∈ O (D j E j ). For occupied cells in the central set When all risk distances on tentacle j are calculated, we compute ∆ j as their minimum:

D j , ∆ ij
∆ j = inf c i ∈(O∩D j )∪ Ēj ∆ ij . If (O D j ) Ēj ≡ ∅, ∆ j = ∞.
In the example of Fig. 2, ∆ 2 = ∆ 12 and ∆ 3 = ∆ 33 . Obviously, overestimating the bounding box sizes leads to more conservative ∆ j .

We then use ∆ j and two hand tuned thresholds ∆ d and ∆ s (0 < ∆ d < ∆ s ), to design the tentacle risk function:

H j =        0 if ∆ j ≥ ∆ s 1 2 1 + tanh 1 ∆ j -∆ d + 1 ∆ j -∆s if ∆ d < ∆ j < ∆ s 1 if ∆ j ≤ ∆ d . (17)
Note that H j smoothly varies from 0, when the dangerous cells associated to tentacle j (if any) are far, to 1, when they are near. If H j = 0, the tentacle is tagged as clear. In practice, threshold ∆ s must be set to the risk distance for which the context ceases to be safe (H becomes greater than 0), so that the robot starts to leave the taught -occupied -path. On the other hand, ∆ d must be tuned as the risk distance for which the context becomes unsafe (H = 1), so that the robot follows the best tentacle to circumnavigate the obstacle.

In our work, we used the values ∆ s = 6 m and ∆ d = 4.5 m. The risk function H j corresponding to these values is plotted in Fig. 3.

The H j s of all tentacles are then compared, in order to determine H in (13). Initially, we calculate the path curvature κ = ω/v ∈ IR that the robot would follow if there were no obstacles. Replacing H = 0 in (13), it is:

κ = λ x x d -x -j v v s + λ ϕ j φϕ /j ω v s ,
which is always well-defined, since j ω = 0 and we have set v s > 0. We obviously constrain κ to the interval

of feasible curvatures [-κ M , κ M ].
Then, we derive the two neighbours of κ among all the existing tentacle curvatures:

κ n , κ nn ∈ K such that κ ∈ [κ n , κ nn ) .
Let κ n be the nearest one, i.e., the curvature of the tentacle that best approximates the safe path2 . We denote it as the visual task tentacle. The situation risk function H v of that tentacle is then obtained by linear interpolation of its neighbours:

H v = (H nn -H n ) κ + H n κ nn -H nn κ n κ nn -κ n . ( 18
)
In practice, H v measures the risk on the visual path, by considering only obstacles on the visual task tentacle and on its neighbour tentacle. In particular, for the context to be safe (i.e., in order to follow the taught path and realize the desired safe task in ( 6)), it is sufficient that the neighbour tentacles are clear (H n = H nn = 0).

This way, obstacles on the sides do not deviate the robot away from the taught path.

Let us now detail our strategy for determining the best tentacle curvature κ b for navigation. This strategy is illustrated by the four examples in Fig. 4, where 5 tentacles are used. In the figure, the dangerous cells ambiguities are again solved first with the distance from κ n , then from κ nn .

Translational Velocity in the Safe Context

We will hereby define the translational velocity in the safe context v s . When the feature motion in the image is fast, the visual tracker is less effective, and the translational velocity should be reduced. This is typically the case at sharp robot turns, and when the camera pan angle ϕ is strong (since the robot is far from the taught 3D path). Hence, we define v s as:

v s (ω, ϕ) = v m + v M -v m 4 σ (19)
with function σ defined as:

σ : IR × -π 2 , π 2 → [0, 4] (ω, ϕ) → [1+tanh (π-k ω |ω|)] [1+tanh (π-k ϕ |ϕ|)] .
v M = 1 ms -1 , k ω = 13, k ϕ = 3.

Translational Velocity in the Unsafe Context

The unsafe translational velocity v u must adapt to the potential danger; it is derived from the obstacles on the best tentacle, defined in Section 3.3. In fact, v u is derived from the collision distance δ b , which is a conservative approximation of the maximum distance that the robot can travel along the best tentacle without colliding. Since the thinner box contains the robot, if R follows the best tentacle, collisions can only occur in occupied cells in C b . In fact, the collision with cell c i will occur at the distance, denoted δ ib ≥ 0, that the thinner box would cover along the best tentacle, before touching the center of c i . The derivation of δ ib is illustrated in Fig. 2 for four occupied cells.

Then, we define δ b as the minimum among the collision distances of all occupied cells in C b : The translational velocity must be designed accordingly. Let δ d and δ s be two hand tuned thresholds such that 0 < δ d < δ s . If the probable collision is far enough (δ b ≥ δ s ), the translational velocity can be maintained at the safe value defined in (19). Instead, if the dangerous occupied cell is near (δ b ≤ δ d ), the robot should stop. To comply with the boundary conditions v u (δ d ) = 0 and v u (δ s ) = v s , in the intermediate situations we apply a constant deceleration:

δ b = inf c i ∈O∩C b δ ib .
a = v 2 s /2(δ d -δ s ) < 0.
Since the distance required for braking at velocity v u (δ b ) is:

δ b -δ d = -v 2 u /2a,
the general expression of the unsafe translational velocity becomes:

v u (δ b ) =    v s if δ b ≥ δ s v s δ b -δ d /δ s -δ d if δ d < δ b < δ s 0 if δ b ≤ δ d , (20)
in order to decelerate as the collision distance δ b decreases. In practice, threshold δ d will be chosen as the distance to collision at which the robot should stop. Instead, threshold δ s must be defined according to the maximum applicable deceleration (noted a M < 0), in order to brake before reaching distance δ d , even when the safe velocity v s is at its maximum v M :

δ s > δ d - 2a M v 2 M .
In our work, we used δ d = 2.7 m and δ s = 5 m, as shown on the right of Fig. 5, where we have plotted v u in function of δ b for three values of v s : 0.4, 0.7, and 1 ms -1 .

Discussion

In this section, we will instantiate and comment our control scheme for visual navigation with obstacle avoidance. Using all the variables defined above, we can explicitly write our controller (13) for visual navigation with obstacle avoidance:

       v = (1 -H) v s + Hv u ω = (1 -H) λx(x d -x)-jvvs+λϕj φϕ jω + Hκ b v u φ = H λx(x d -x)-(jv+jωκ b )vu j φ -(1 -H) λ ϕ ϕ . ( 21
)
This control law has the following interesting properties.

1. In the safe context (H = 0), (21) becomes:

     v = v s ω = λx(x d -x)-jvvs+λϕj φϕ jω φ = -λ ϕ ϕ . ( 22
)
In Sect. 3.1, we proved that this controller guarantees global asymptotic stability of the safe task ṡ * s . As in [Cherubini et al. 2009] and [Diosi et al. 2011], where obstacles were not considered, the image error is controlled only by ω, which also compensates the centroid displacements due to v and to φ through the image jacobian components (15), to fulfill the visual task (3). The two remaining specifications in (5), instead, are achieved by inputs v and φ: the translational velocity is regulated to improve tracking according to (19), while the camera is driven forward, to ϕ = 0. Note that,

to obtain H = 0 with the tentacle approach, it is sufficient that the neighbour tentacles are clear (H n = H nn = 0), whereas in the potential field approach used in [START_REF] Cherubini | [END_REF],

even a single occupied cell would generate H > 0. Thus, one advantage of the new approach is that only obstacles on the visual path are taken into account.

2. In the unsafe context (H = 1), (21) becomes:

     v = v u ω = κ b v u φ = λx(x d -x)-(jv+jωκ b )vu j φ . ( 23
)
In Sect. 3.1, we proved that this controller guarantees global asymptotic stability of the unsafe task ṡ * u . In this case, the visual task ( 3) is executed by φ, while the two other specifications are ensured by the 2 other degrees of freedom: the translational velocity is reduced (and even zeroed to v = v u = 0 for very near obstacles such that δ b ≤ δ d ), while the angular velocity makes the robot follow the best tentacle (ω/v = κ b ). Note that, since no 3D positioning sensor (e.g., gps) is used, closing the loop on the best tentacle is not possible; however, even if the robot slips (e.g., due to a flat tire), at the following iterations tentacles with stronger curvature will be selected to drive it towards the desired path, and so on. Finally, the camera velocity φ in (23) compensates the robot rotation, to keep the features in view.

3. In intermediate situations (0 < H < 1), the robot navigates between the taught path, and the best path considering obstacles. The situation risk function H representing the danger on the taught path, drives the transition, but not the speed. In fact, note that, for all

H ∈ [0, 1], when δ b ≥ δ s : v = v s .
Hence, a high velocity can be applied if the path is clear up to δ s (e.g., when navigating behind another vehicle).

4. Control law (21) guarantees that obstacle avoidance has no effect on the visual task, which can be achieved for any H ∈ [0, 1]. Note that plugging the expressions of v, ω, and φ from ( 21) into the visual task equation:

ẋ = j v v + j ω ω + j φ φ (24)
yields (3). Therefore, desired state x d is globally asymptotically stable for the closed loop system, ∀ H ∈ [ 0, 1]. This is true even in the special case where v = 0. In fact, the robot stops and v becomes null, if and only if H = 1 and v u = 0, implying that ω = 0 and φ =

λx(x d -x) j φ
, which allows realization of the visual task. In summary, from a theoretical control viewpoint (i.e., without considering image processing nor field of view or joint limits constraints), this proves that if at least one point in I d is visible, the visual task of driving the centroid abscissa to x d will be achieved, even in the presence of unavoidable obstacles. This strategy is very useful for recovery: since the camera stays focused on the visual features, as soon as the path returns free, the robot can follow it again.

5. Controller ( 21) does not present local minima, i.e., non-desired state configurations for which u is null. In fact, when H < 1, u = 0 requires both v u and v s to be null, but this is impossible since from ( 19)), v s > v m > 0. Instead, when H = 1, it is clear from ( 23) that for u to be null it is suffcient that x d -x = 0 and v u = 0. This corresponds to null desired dynamics: ṡ * u = 0 (see ( 9)). This task is satisfied, since plugging u = 0 into (10), yields precisely ṡu = 0 = ṡ * u .

6. If we tune the depth Z c to infinity in (15), j v = 0, and j φ = j ω = 1 + x 2 . Thus, control law (21) becomes:

     v = (1 -H) v s + Hv u ω = (1 -H) λ x x d -x jω + (1 -H) λ ϕ ϕ + Hκ b v u φ = Hλ x x d -x jω -(1 -H) λ ϕ ϕ + Hκ b v u .
Note that, for small image error (x ≈ x d ), φ ≈ -ω. In practice, the robot rotation is compensated by the camera pan rotation, which is an expected behavior.

Simulations

In this section and in the following, we will detail the simulated and real experiments that were used to validate our approach. Simulations are in the video shown in Extension #1.

For simulations, we made use of Webots 4 , where we designed a car-like robot equipped with an actuated At first, no obstacle is present in the environment, and the robot is driven along a taught path. Then, up to 5 obstacles are located, near and on the taught path, and the robot must replay the visual path, while avoiding them. In addition, the obstacles may partially occlude the features. Although the sensors are noise-free, and feature matching is ideal, these simulations allow validation of controller (21).

By displacing the obstacles, we have designed the 6 scenarios shown in Fig. 6. For scenarios A, B and C, the robot has been previously driven along a 30 m straight path, and N = 8 key images have been acquired, whereas in scenarios D, E and F, the taught path is a closed loop of length 75 m and N = 20 key images, which is followed in the clockwise sense. In all scenarios, the robot is able to navigate without colliding, and this is done with the same parameters. The metrics used to assess the experiments are the image error with respect to the visual database x -x d (in pixels), averaged over the whole experiment and denoted ē, and the distance, at the end of the experiment, from the final 3D key pose ( , in cm). The first metric ē is useful to assess the controller accuracy in realizing the visual path following task. The latter metric is less relevant, since task (3) is defined in the image space, and not in the pose space.

In all six scenarios, path following has been achieved, and in some cases, navigation was completed using only 3 image points. Obviously, this is possible in simulations, since feature detection is ideal: in the real case, which includes noise, 3 points may be insufficient. Some portions of the replayed paths, corresponding to the obstacle locations, are far from the taught ones. However, these deviations would have been indispensable to avoid collisions, even with a pose-based control scheme. Let us detail the robot behaviour in the six scenarios:

• Scenario A: two walls, which were not present during teaching, are parallel to the path, and three boxes are placed in between. The first box is detected, and overtaken on the left. Then, the vehicle passes between the second box and the left wall, and then overtakes the third box on the right. Finally, the robot converges back to the path, and completes it. Although the walls occlude features on the sides, the experiment is successful, with ē = 5, and = 23.

• Scenario B: it is similar to A, except that there are no boxes, and that the left wall is shifted towards the right one, making the passage narrower towards the end. This makes the robot deviate in order to pass in the center of the passage. We obtain ē = 6, = 18.

• Scenario C: this scenario is designed to test the controller in the presence of unavoidable obstacles. Two walls forming a corner, are located on the path. This soon makes all tentacles unsafe:

∆ j ≤ ∆ d | ∀ j, yielding H = 1.
Besides, as the robot approaches the wall, the collision distance on the best tentacle δ b decreases, and eventually becomes smaller than δ d , to make v u = 0 and stop the robot (see ( 23)). Although the path is not completed (making metric irrelevant), the collision is avoided, and ē = 4 pixels. As proved in Section 4, convergence of the visual task (x = x d ) is achieved, in spite of u = 0. In particular, here, the centroid abscissa on the third key image in the database is reached.

• Scenario D: high walls are present on both sides of the path; this leads to important occlusions (less than 50% of the database features are visible), and to a consequent drift from the taught path. Nevertheless, the final key image is reached, without collisions, and with ē = 34, and = 142. Although this metric is higher than in the previous scenarios (since the path is longer and there are numerous 6). When the second obstacle is bypassed, the camera pan is reset to zero. The reduction of v at iteration 2800 is due only to the sharp turn (i.e., to the reduction of v s ), since the path is safe at this point. Then, the third obstacle triggers H, and is easily circumnavigated. From iteration 3700 onwards, the situation risk function is cancelled. Correspondingly, the variables are driven by (22). Note also that the camera angle ϕ is reset to 0 in less than 200 iterations, and remains null until the end of the experiment. The small spikes in the angular velocity ω, which appear throughout the experiment, correspond to the changes of database key images (except when they are provoked by the obstacles, as ), ω (in rad s -1 ), φ (in rad s -1 ), and ϕ (in rad).

discussed above).

The six simulations have been repeated by setting the feature depth Z c to infinity. For all six scenarios, the image accuracy, assessed with ē, is very near to the one obtained when Z c = 15 m. On the other hand, the pose accuracy, assessed with , is lower when Z c = ∞, as shown in Table 1. The difference is relevant on long paths (scenarios D, E and F). Although the navigation task is defined in the image space, these experiments show that tuning Z c , even coarsely, according to the environment, can contribute to the controller performance in the 3D space. This aspect had already emerged in part in our previous work [Cherubini et al. 2009].

Real Experiments

After the simulations, the framework has been ported on our CyCab vehicle, set in car-like mode (i.e., using the front wheels for steering), for real outdoor experimental validation. The robot is equipped with a coarsely calibrated 320 × 240 pixels 70 • field of view, B&W Marlin (F-131B) camera mounted on a TRACLabs Biclops Pan/Tilt head (the tilt angle is null, to keep the optical axis parallel to the ground), and with a 2-layer, 110 • scanning angle, laser SICK LD-MRS. A dedicated encoder on the TRACLabs head precisely measures the pan angle ϕ required in our control law (see ( 21)). The grid is built by projecting the laser readings from the 2 layers on the ground. Exactly the same configuration (i.e., same parameters, gains and grid size) tuned in Webots is used on the real robot. The centroid depth value that we used in simulations (Z c = 15 m) proved effective in all real experiments as well, although the scenarios were very variegate. This confirms, as shown in [Cherubini et al. 2009], that a very coarse approximation of the scene depth is sufficient to effectively tune Z c . The velocity (v M = 1, as in Webots) has been reduced due to the image processing rate (10 Hz), to limit the motion of features between successive images; the maximum speed attainable by the CyCab is 1.3 ms -1 anyway. Since camera (10 Hz) and laser (12.5 Hz) processing are not synchronized, they are implemented on two different threads, and the control input u derived from control law ( 21) is sent to the robot as soon as the visual information is available (10 Hz).

It is noteworthy to point out that the number of tentacles that must be processed, and correspondingly, the computational cost of the laser processing thread, increase with the context danger. For clarity, let us discuss two extreme cases: a safe and an occupied contexts. To verify that a context is safe (i.e., that H v = 0 in ( 17)), all the cells in the dangerous areas D E of only the two neighbour tentacles must be explored.

Instead, in a scenario where the grid is very occupied, all of the tentacles in K may need to be explored. In general, this second case will be more costly than the first. However, in practice, since only the minimum risk and collision distances (∆ j and δ b ) are required by our controller, exploration of a tentacle stops as soon as the nearest cell is found occupied, so that the tentacles are rarely entirely explored. The experiments showed that the computational cost of laser processing, using the chosen number of tentacles (i.e., 21, as mentioned in Sect. 2.3), was never a critical issue with respect to that of image processing.

First, we have assessed the performance of our control scheme when a very long obstacle is positioned perpendicularly on the taught path (denoted path A, and shown in Fig. 9). In Fig. 10, we have plotted the control inputs u, and the situation risk function H. The smooth trend of u at the beginning and end of the experiments is due to the acceleration saturation carried out at the CyCab low-level control. The obstacle is overtaken on the left, while the camera rotates right to maintain scene visibility (dotted black and dotted gray curves in Fig. 10). The robot is able to successfully reach the final key image and complete navigation, although it is driven over 5 meters away from the taught 3D path by the obstacle. In practice, soon after the obstacle is detected (i.e., after 5 s), tentacles with first positive (5 -16 s), and then negative (16 -25 s) curvature are selected. Since these tentacles are clear, v is reduced only for visual tracking, by (19) (solid black curve in Fig. 10). This is a major feature of the tentacle method, which considers only the actual collision risk of obstacles for reducing the translational velocity. After 25 s, the environment returns clear (H = 0), so the visual tentacle can be followed again, and the robot is driven back to the path. Then (38 -52 s) a small bush on the left triggers H and causes a small counterclockwise rotation along with a slight decrease in v. Then the context returns safe, and the visual path can be followed for the rest of the experiment. The translational velocity averaged over the experiment is 0.79 ms -1 , which is more than twice the speed reached in [START_REF] Cherubini | [END_REF].

After these results, we have run a series of experiments, on longer and more crowded paths (denoted B to E on Fig. 11) on our campus. All campus experiments here are also visible in the video shown in Extension #2. The Cycab was able to complete all paths (including 650 m long path E), while dealing with various natural and unpredictable obstacles, such as parked and driving cars and pedestrians. The translational velocity averaged over these experiment was 0.85 ms -1 . Again, by assessing the collision risk only along the visual path, non-dangerous obstacles (e.g., walls or cars parked on the sides) are not taken into account. This aspect is clear from Fig. 12 Another nice behaviour is shown in Fig. 12(center): if a stationing car is unavoidable, the robot decelerates and stops with (20), but, as soon as the car departs, the robot gradually accelerates (again with (20)), to resume navigation. In fact, as we mentioned in Section 4, when the best tentacle is clear up to distance δ s , a high velocity can be applied: v = v s , independently from the value of H. In the future, this feature of our approach could even be utilized for vehicle following.

An experiment with a crossing pedestrian is presented in Fig. 12(right). The pedestrian is considered irrelevant, until it enters the visual task tentacle (second image). Then, the clockwise tentacles are selected to pass between the person and the right side walk. When the path is clear again, the robot returns to the visual task tentacle, which is first counter-clockwise (fourth image) and then straight (fifth image).

In October 2011, as part of the final demonstration of the French ANR project CityVIP, we have validated our framework in a urban context, in the city of Clermont Ferrand. The experiments have taken place in the crowded neighbourhood of the central square Place de Jaude, shown in Fig. 13. For four entire days, our Cycab has navigated autonomously, continuously replaying a set of visual paths of up to 700 m each, amidst a variety of unpredictable obstacles, such as cars, pedestrians, bicycles and scooters. In Fig. 14, we show some significant snapshots of the experiments that were carried out in Clermont Ferrand. These include photos of the Cycab, as well as images acquired by the on-board camera during autonomous navigation.

These experiments are also visible in the video shown in Extension #3.

In Fig. 14(a-c), Cycab is moving in a busy street, crowded with pedestrians and vehicles. First, in Fig. 14(a), we show a typical behaviour adopted for avoiding a crossing pedestrian: here, Cycab brakes as a lady with black skirt crosses the street. The robot would either stop or circumnavigate the person, and in spite of the robustness of the image processing algorithms. which has been proved in [START_REF] Royer | [END_REF], in some cases (e.g., when the camera was overexposed to sunlight), the visual features required for navigation could not be detected. Future work in adapting the camera automatic shutter settings should solve this issue.

In the current version of our framework, moving obstacles are not specifically recognized and modelled.

Although, as the experiments show, we are capable of avoiding slowly moving obstacles (e.g., crossing pedestrians or baby pushchairs as in Fig. 14(g)), the main objective of our future work will be to directly tackle this problem within the control law, in order to avoid fast obstacles as well. This can be done, for example, by estimating the velocity of the detected objects, and then using it to predict their future position.

In our opinion, the main difficulty, in comparison with the case of static obstacles, will concern the accuracy and computation cost of this estimation process.

Overall, Cycab has navigated using an average of approximately sixty visual points on each image, and some paths have even been completed using less than 30 points. Along with all the cited technical aspects, the experiments highlighted the reactions of non-robotic persons to the use of autonomous ground vehicles in everyday life. Most passer-bys had not been informed of the experiments, and responded with curiosity, surprise, enthusiasm, and -rarely -slight apprehension. Since our method is purely sensor-based and pose-independent, it is perfectly suited for visual navigation.

Conclusions

Extensive outdoor experiments, even in urban environments, show that it can be applied in realistic and challenging situations including moving obstacles (e.g., cars and pedestrians). To our knowledge, this is the first time that outdoor visual navigation with obstacle avoidance is carried out in urban environments at approximately 1 ms -1 on over 500 m, using neither GPS nor maps. In the near future, we plan to explicitly take into account the velocity of moving obstacles within our controller, in order to avoid fast obstacles, which are currently hard to deal with. Perspective work also includes automatic prevention of the visual occlusions provoked by the obstacles.

Figure 1 :

 1 Figure 1: General definitions. Left: top view of the robot (rectangle), equipped with an actuated camera (triangle); the robot and camera frame (respectively, F R and F C ) are shown. Right: database of the key images, with current and next key images emphasized; the image frame F I is also shown, as well as the visual features (circles) and their centroid (cross).

Fig. 2

 2 Fig. 2(a): it is linked to F R , with cell sides parallel to X and Y . Its longitudinal and lateral extensions are limited (X m ≤ X ≤ X M and Y m ≤ Y ≤ Y M ), to ignore obstacles that are too far to jeopardize the robot. The size of the grid should increase with the robot velocity, to guarantee the sufficient time for obstacle avoidance. An appropriate choice for |X m | is the length of the robot, since obstacles behind cannot hit it as it advances. In this work, we use: X M = Y M = 10 m, X m = -2 m, Y m = -10 m. Any grid cell c centered at (X, Y ) is considered occupied if an obstacle has been sensed in c. The cells have size 0.2 × 0.2 m. For the cells entirely lying in the scanner area, only the current scanner reading is considered. For all other cells

Figure 2 :

 2 Figure 2: Obstacle models with 4 occupied cells c 1 , . . . , c 4 . (a) Occupancy grid, straight (b, c) and sharpest counterclockwise (d, e) tentacles (dashed). When a total of 3 tentacles is used, the straight and sharpest counterclockwise are characterized respectively by index j = 2 and j = 3. For these two tentacles, we have drawn: classification areas (collision C j , dangerous central D j , dangerous external E j ), corresponding boxes and delimiting arcs of circle, and cell risk and collision distances (∆ ij , δ ij ). For tentacle j = 3 in the bottom right, we have also drawn the tentacle center (cross) and the ray of cell c 1 , denoted Γ 13 .

  is the distance that the middle boundary box would cover along tentacle j before touching the cell center. For occupied cells in the external set, only a subset Ēj is taken into account: Ēj ⊆ O E j . This subset contains only cells which reduce the clearance in the tentacle normal direction. For each external occupied cell, we denote Γ ij the ray starting at the tentacle center and passing through c i . Cell c i is added to Ēj if and only if, in D j E j , there is at least an occupied cell crossed by Γ ij on the other side of the tentacle.In the example of Fig.2(e), O E 3 = {c 1 , c 3 , c 4 }, whereas Ē3 = {c 1 , c 3 }. Cell c 4 is not considered dangerous, since it is external, and does not have a counterpart on the other side of the tentacle. Then, for cells in Ēj , ∆ ij is the sum of two terms: the distance from the center of c i to its normal projection on the perimeter of the dangerous central area, and the distance that the middle boundary box would cover along tentacle j before reaching the normal projection. The derivation of ∆ ij is illustrated, in Fig.2, for four occupied cells. Note that for a given cell, ∆ ij may have different values (or even be undefined) according to the tentacle that is considered.

Figure 3 :

 3 Figure 3: Risk H j , in function of the tentacle risk distance ∆ j (m) when ∆ s = 6 m and ∆ d = 4.5 m.

  (i.e., for each tentacle j, the cells in D j E j ) associated to the visual task tentacle and to the best tentacle are respectively shown in light and dark gray. The occupied cells are shown in black. The best tentacle is derived from the tentacle risk functions defined just above. If H v = 0 (as in Fig.4(a)), the visual task tentacle can be followed: we set κ b = κ n , and we apply (13) with H = 0. Instead, if H v = 0, we seek a clear tentacle (H j = 0). First, to avoid abrupt control changes, we only search among the tentacles between the visual task one and the best one at the previous iteration 3 , noted κ pb , and mid-gray in the figure. If many clear ones are present, the nearest to the visual task tentacle is chosen, as in Fig.4(b). If none of the tentacles with curvature in [κ n , κ pb ] is clear, we search among the others. Again, the best tentacle will be the clear one that is closest to κ n and, in case of ambiguity, the one closest to κ nn . If a clear tentacle has been found (as in Fig.4(c)), we select it and set H = 0. Instead, if no tentacle in K is clear, the one with minimum H j calculated using (17) is chosen, and H is set equal to that H j . In the example of Fig.4(d), tentacle 1 is chosen and we set H = H 1 , since ∆ 1 = sup{∆ 1 , . . . , ∆ 5 }, hence H 1 = inf{H 1 , . . . , H 5 }. Eventual

Figure 4 :

 4 Figure 4: Strategy for selecting the best tentacle among 5 in four different scenarios. The cells associated to the visual task tentacle, to the previous best tentacle, and to the best tentacle are shown in increasingly dark gray; the corresponding tentacles are dashed, and the occupied cells c 1 and c 2 are shown in black. (a) Since it is clear, the visual task tentacle with curvature κ n is selected: κ b = κ n . (b) The clear tentacle with curvature in [κ n , κ pb ] nearest to κ n is chosen. (c) Since all tentacles with curvature in [κ n , κ pb ] are occupied, the clear one nearest to the visual task tentacle is chosen. (d) Since all tentacles are occupied, we select the one with smallest H j , hence, largest risk distance ∆ j (here, ∆ 1 )given the middle boundary box).

Function ( 19 Figure 5 :

 195 Figure 5: Left: safe translational velocity v s (ms -1 ) in function of ω (rads -1 ) and ϕ (rad). Right: unsafe translational velocity v u (ms -1 ) in function of δ b (m) for three different values of v s .

  If all cells in C b are free, δ b = ∞. In the example of Fig. 2, assuming the best tentacle is the straight one (b = 2), δ b = δ 12 . Again, oversizing the box leads to more conservative δ b .

Figure 6 :

 6 Figure 6: Six obstacle scenarios for replaying two taught paths (black) with the robot (rectangle): a straight segment (scenarios A to C), and a closed loop followed in the clockwise sense (D to F). Visual features are represented by the spheres, the occupancy grid by the rectangular area, and the replayed paths are drawn in gray.

Figure 7 :

 7 Figure 7: Scenario A. For each of the 5 relevant iterations we show (top to bottom): the robot overtaking the first obstacle, the next key image, the current image, and the occupancy grid.

Figure 8 :

 8 Figure 8: Evolution of relevant variables in scenario E; top to bottom: H, v (in m s -1 ), ω (in rad s -1 ), φ (in rad s -1 ), and ϕ (in rad).

Table 1 :Figure 9 :

 19 Figure 9: Scenario A (a long obstacle is avoided): taught (white) and replayed (black) paths.

  (left), where a stage of the navigation experiment on path E is illustrated. From top to bottom, we show: the next key image in the database I d , the current image I, and three consecutive occupancy grids processed at that point of

Figure 10 :

 10 Figure10: Control inputs in scenario A: v (solid black, in ms -1 ), ω (dotted black, rads -1 ) φ (dotted gray, rads -1 ) and H (solid gray).

Figure 11 :

 11 Figure 11: Map of the four navigation paths B, C, D, E.

Figure 12 :

 12 Figure 12: Validation with: irrelevant obstacles (left), traffic (center) and a moving pedestrian (right). The visual task tentacle and best tentacle (when different) are respectively indicated as VTT and BT in the occupancy grids.

A

  framework with simultaneous obstacle avoidance and outdoor vision-based navigation, without any 3D model or path planning has been presented. It merges a novel, reactive, tentacle-based technique with visual servoing, to guarantee path following, obstacle bypassing, and collision avoidance by deceleration.

Figure 14 :

 14 Figure 14: Snapshots of the urban experiments. (a) Avoiding a crossing pedestrian. (b-c) Navigating close to moving cars and to a scooter, respectively. (d-e) Avoiding a stationing police patrol and a stationing vehicle, respectively. (f) Navigating with different light conditions, using the same taught database. (g) Avoiding a pedestrian with a baby pushchair.

See: www.irisa.fr/lagadic/demo/demo-cycab-vis-navigation/vis-navigation.

Without loss of generality, we consider that intervals are defined even when the first endpoint is greater than the second, e.g., [κn, κnn) should be read (κnn, κn] if κn > κnn.

At the first iteration, we set κ pb = κn.

Acknowledgment

This work has been supported by ANR (French National Agency) CityVIP project under grant ANR-07 TSFA-013-01.

occlusions), it is still reasonably low.

• Scenario E: two obstacles are located on the path, and two other are near the path. The first obstacle is overtaken on the left, before avoiding the second one, also on the left. Then, the robot converges to the path and avoids the third obstacle on the right, before reaching the final key image. We obtain ē = 33, and = 74. The experiment shows one of the advantages of our tentacle-based approach: lateral data in the grid is ignored (considering the fourth obstacle, would have made the robot curve away from the path).

• Scenario F: here, the controller is assessed in a situation where classical obstacle avoidance strategies (e.g. potential fields) often fail because of local minima. In fact, when detected, the first obstacle is centred on the X axis and orthogonal to it. This may induces an ambiguity, since occupancies are symmetric with respect to X. However, the visual features distribution, and consequent visual task tentacle κ n drive the robot to the right of the obstacle, which is thus avoided. We have repeated this experiment with 10 randomly generated visual feature distributions, and in all cases the robot avoided the obstacle. The scenario involves four more obstacles, two of which are circumnavigated externally, and two on the inside. Here, ē = 29, and = 75.

In Fig. 7, we show 5 stages of scenario A. In this figure, as well as later in Fig. 12, in Fig. 14, the segments linking the current and next key image points are drawn in the current image. In the occupancy grid, the dangerous cell sets associated to the visual task tentacle and to the best tentacle (when different) are respectively shown in gray and black, and two black segments indicate the scanner amplitude. Only cells that can activate H (i.e., cells at distance ∆ < ∆ s ) have been drawn. At the beginning of the experiment (iteration 10), the visual features are driving the robot towards the first obstacle. When it is near enough, the obstacle triggers H (iteration 200), forcing the robot away from the path towards the best tentacle, while the camera rotates clockwise to maintain feature visibility (iterations 200, 290 and 405). Finally (iteration 630), the controller drives the robot back to the path, and the camera to the forward-looking direction.

Further details can be obtained by studying some relevant variables. We focus on scenario E, for which we have plotted in Fig. 8 the values of H, v, ω, φ, and ϕ during navigation. The curve of H shows when the obstacles, respectively the first (at iterations 0 -1000), second (1100 -1700) and third (3200 -3700), have intervened in control law (21). As aforementioned, the fourth obstacle does not trigger H, since it is too far on the side to jeopardize the robot. Let us now discuss the trend of the five curves. Since the beginning, the Appendix A: Index to Multimedia Extensions

The multimedia extensions (3 videos) to this article (see Table 2) are at www.ijrr.org. In all three videos, the segments linking the current and next key image points are drawn in light green in the current image. In the occupancy grid, the dangerous cell sets associated to the visual task tentacle and to the best tentacle (when different) are respectively shown in red and blue, and two black segments indicate the scanner amplitude.

Only cells that can activate H (i.e., cells at distance ∆ < ∆ s ) have been drawn.


C Albea
email: cal-bea@laas.fr

F Gordillo
email: gordillo@esi.us.es

† F Gordillo

Estimating the attraction domain for the boost inverter

This work presents an approach for estimating the domain of attraction for polynomial systems with state and control-signal constraints including saturation. In many problems, it is possible to derive global stability properties for such systems, neglecting constraints. Consideration of the constraints usually makes the problem much more involved. In this paper, the stability analysis performed for the unconstrained case is used for the problem as a whole. For application of the method, there are powerful computational tools that can be employed in cases of polynomial systems. The technique is not only valid for the analysis of equilibrium points but also for other attractors, such as limit cycles.

As examples, the domain of attraction for given control laws is estimated for a nonlinear DC-DC boost converter as well as for a boost inverter.

Introduction

Many control designs do not take into account system constraints such as actuator saturation and forbidden regions in the state space, in spite of the fact that these kinds of constraints are present in almost any control application. Nonlinear control books contain an abundance of stability analysis examples in which constraints are not present [START_REF] Hassan | Nonlinear Systems[END_REF][START_REF] Vidyasagar | Nonlinear systems analysis[END_REF]. The behaviour of the closed-loop constrained system is not necessarily greatly degraded with respect to the case of unconstrained analysis, making this approach very successful. Nevertheless, in some situations the global stability character may be lost, resulting in local stability of the desired operating behaviour within a bounded region of attraction [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF][START_REF] Albea | Control of the Boost DC-AC Converter by Energy Shaping[END_REF]. Estimation of this region is desirable but represents a difficult problem.

There exist many published methods for estimating the region of attraction (see, for example [START_REF] Genesio | On the estimation of asymptotic stability regions: State of the art and new proposals[END_REF][START_REF] Hassan | Nonlinear Systems[END_REF] and the references therein). One example of this kind of method is based on Lyapunov theory, in which closed Lyapunov-function level surfaces are utilized to determine approximate sort of (conservative) estimations for the region of attraction [START_REF] Hassan | Nonlinear Systems[END_REF]. These methods often employ polynomial systems [START_REF] Levin | An analytical method for estimation the domain of attraction for polynomial differential equations[END_REF][START_REF] Romanchuk | Computing regions of attraction with polytopes: planar case[END_REF][START_REF] Tibken | Computation of subsets of the domain of attraction for polynomial systems[END_REF][START_REF] Topcu | Help on sos [ask the experts[END_REF]. There exist powerful mathematical tools that can be used in the computation of the maximum acceptable level for polynomial systems [START_REF] Huber | A polyhedral method for solving sparse polynomial systems[END_REF][START_REF] Parrilo | Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization[END_REF][START_REF] Verschelde | Algorithm 795: Phcpack: a general-purpose solver for polynomial systems by homotopy continuation[END_REF][START_REF] Chesi | Solving quadratic distance problems: an LMI-based approach[END_REF]. Some of these tools could be further developed for application to non-polynomial systems as well [START_REF] Chesi | Estimating the domain of attraction for non-polynomial systems via LMI optimizations[END_REF]. Application of these methods would imply to seek for Lyapunov functions in order to be able to deal with the constraints. The search for a Lyapunov function by means of the numerical estimation method may be seen as an advantage as the user would not be required to propose a Lyapunov function. In cases like these, however, the computational methods are faced with a problem that is very difficult to tackle alone and would encounter difficulty in resolving problems in moderately complex systems.

Furthermore, saturation-like functions, which are one of the most common nonlinearities in practice, usually fall outside of the scope of these techniques. Saturations-like functions are typically only considered in the case of linear systems [START_REF] Álamo | Estimation of the domain of attraction for saturated discrete-time systems[END_REF][START_REF] Blanchini | Any domain of attraction for a linear constrained system is a tracking domain of attraction[END_REF][START_REF] Te Pare | Synthesizing stability regions for systems with saturating actuators[END_REF][START_REF] Pittet | Stability regions for linear systems with saturating controls viacircle and Popov criteria[END_REF], or in specific cases, as for instance in a port-controlled Hamiltonian system [START_REF] Wei | Estimate of domain of attraction for a class of port-controlled hamiltonian systems subject to both actuator saturation and disturbances[END_REF].

This paper presents a simple idea that solves the problem of estimation of attraction domain estimation for polynomial non-linear systems (among others) with saturation-like constraints and state constraints. The term saturationlike constraints is used for non-linear functions γ(u) that appear in the system model and they become the identity, γ(u) = u, in certain regions of the state space that present the desired behaviour, (those regions are referred to regions in which such constraints are not active). Functions of this sort include typical control signal saturation as well as others, such as rate limiters, for example. Other constraints on the state variables can be considered as well, with an eye to being able to exploit the unconstrained global stability analysis and use the result to obtain a (conservative) estimation of the region of attraction for the constrained case. Employing the Lyapunov function, which is an integral part of this analysis, eliminates the need to search for a Lyapunov function to estimate the domain of attraction, thus simplifying the problem. On the other hand, the estimated attraction region is included in the domain where the saturation-like constraints are not active and, therefore, the method introduces a new source of conservatism, which could prove to be useful for solving problems when all other methods fail.

To demonstrate the effectiveness of this technique, two electronic application examples are presented. These applications use boost topology in order to obtain oscillating current in such a way that the output voltage can be higher than the input voltage. In the first application, the DC-DC boost converter, physical constraints prevent the output signal from passing through zero, and alternating current is therefore unattainable [START_REF] Biel | Sliding mode control of a boost-buck converter for ac signaltracking task[END_REF][START_REF] Fossas-Colet | Asymptotic tracking in DC-to-DC nonlinear power converters[END_REF][START_REF] Sira-Ramírez | A Dynamic Pole Assignment Approach to Stabilization and Tracking in DC-to-DC Full Bridge Power Converters[END_REF], leading to a configuration change in the second application in which the circuit is duplicated, thereby yielding an inverter [START_REF] Ro Cáceres | A boost DC-AC converter: analysis, design, and experimentation[END_REF]. The control laws analyzed are presented in [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF][START_REF] Gordillo | Autonomous oscillation generation in electronic converters[END_REF] for the first case and [START_REF] Albea | Control of the Boost DC-AC Converter by Energy Shaping[END_REF] for the second case. Application of these control laws does not achieve global stability due to two reasons: firstly, the ideal control signal cannot be implemented globally due to control signal saturation; and secondly, the circuit imposes physical constraints on certain state variables: the capacitor voltages, for example, cannot be negative. These reasons motivate the computation of an estimated attraction domain, in such a way that in the starting phase, the system must be driven into a point inside this region. The application of the method produces good, albeit conservative, results. However, certain difficulties should be mentioned here: 1) the system equations are quite involved and can even present non-polynomial terms, namely rational functions; 2) the desired attractor is not an equilibrium point but a limit cycle. Both complications make the use of any other analysis method a formidable task.

The rest of this article is organized as follows: in Section 2 the problem statement as well as the proposed method for estimation of the attraction domain are presented. Section 3 will provide an overview of sums of squares optimization, the numerical tool that is employed in the application of the proposed technique. Section 4 is devoted to the practical cases. Conclusions are drawn in Section 5.

Estimation of the attraction domain

In every control system, the control signal is subject to physical constraints such as saturations, rate-limiters, etc. As for control designs, such constraints are typically disregarded and the resulting control law is applied to the actuator. In this way, if the designed control law is u d = α(x), where x is the state variable, the actual control signal is u = γ(u d ) = γ (α(x)), where γ(•) is a saturation-like function. This approach is valid when the actual expression for u is used in the stability analysis of the resultant system. It is however quite common to neglect the actuator constraints in the stability analysis so as to simplify the analysis. In fact, the local stability property is not usually affected by these constraints, since in a neighborhood D of the desired attractor they are not active, that is, γ(α(x)) = α(x) ∀x ∈ D. However, the resulting attraction domain may be affected by constraints γ. This study deals with the estimation of this attraction domain based on a stability analysis that neglects the constraints.

The analysis can also take into account state variable constraints in the following sense. Assume that there exists a "forbidden" region in the state space. This means that the system state must remain within the boundaries of a pre-specified admissible ("safe") region. The estimation of the domain of attraction should take into account these constraints.

Formally, the problem can be stated as follows:

Actual system Consider a control system type defined as such:

ẋ = f a (x, u), (1)
where x ∈ S ⊂ R n , u ∈ R nu . Function f a may include saturation-like functions. Furthermore, due to physical considerations, the state of the system must not go out of an admissible region T .

Unconstrained system Assume that an approximate model of the system is ẋ

= f (x, u), (2)
where function f : R n → R n and in D, f (x, u) ≡ f a (x, u). The reader should assume that this approximate model includes neither the saturation-like functions nor the state constraints.

Assume that a control law u = α(x) has been designed for the unconstrained model ( 2) for a given control objective.

Remark 1 The control objective is not necessarily the stabilization of an equilibrium point, but perhaps the stabilization of limit cycles, for instance, as seen in the examples in Section 4, can be considered.

Assumption 1 There is a widely known radially unbounded Lyapunov function V (x), in which a compact positively invariant set Ω, ∂V ∂x f (x, α(x)) ≤ 0. Let M be the largest invariant subset of the set for which V = 0 in Ω.

By the LaSalle invariance principle, Assumption 1 guarantees that the trajectories of the unconstrained model tend to M. It is implicitly assumed that this is the desired behaviour. Notice that if the original Lyapunov theorem is used to prove global stability, the previous assumption is also fulfilled.

Assumption 1 also guarantees local stability for system [START_REF] Álamo | Estimation of the domain of attraction for saturated discrete-time systems[END_REF]. The problem lies in the estimation of the domain of attraction.

The key is to ensure that the system state remains within the boundaries of the region where saturations are not active, thus introducing new constraints. A conservative estimation of the region of attraction can then be easily obtained. The idea of the method is schematically shown in Fig. 1. The advantage of the relative ease with which the estimation is obtained, however, is compromised by the fact that it may be far too conservative. Nevertheless, in many problems this simple idea may give satisfactory results. Assumption 2 Consider system (1) with control law u = α(x). Let be A = S ∩ T , that is, the intersection between the safe region and the region where the saturation-like functions are not active. It is assumed that this set can be estimated by a set of inequalities g(x) ≥ 0, where g : R n → R ng .

Now the problem can be transformed as follows:

Given a control system ẋ = f (x, u) with constraints in both the state variables and the control input g(x) ≥ 0, assume that a control law u = α(x) has been designed such that global stability is confirmed when no constraints are taken into account. The problem lies in estimating a region of attraction for the real system with constraints when this control law is applied. A (conservative) estimation for the attraction domain of the system with constraints is given by the following theorem:

Theorem 1 Under assumptions 1 and 2, assume that there exists a constant c > 0 such that in the set Ω c = {x : V (x) ≤ c}, the constraints g(x) ≥ 0 are satisfied. Then, all trajectories of the system with constraints starting at Ω c tend to M ∩ Ω c .

Proof Since in Ω c the constraints are satisfied, the results for the unconstrained system are valid in Ω c . Therefore, V ≤ 0 in Ω c and Ω c is positively invariant. Furthermore, since V (x) is radially unbounded Ω c is compact. The statement can be validated by applying LaSalle's invariance principle.

Remark 2 Since M ∩ Ω c ⊂ M, the theorem guarantees the desired asymptotic behaviour for the system with constraints.

Remark 3 As with other techniques for estimation of attraction domain, the present method is conservative. In this case the conservativeness is mainly due to two facts:

• The estimation of the region of attraction is restricted to domains in which V ≤ c.

• The method considers the saturation-like functions as hard constraints.

Nevertheless, there may be points where the saturations are active in the actual attraction domain.

Using Theorem 1, the problem is reduced to finding a value c > 0 such that g(x) ≥ 0 at the points where V (x) ≤ c. In order to use numerical tools for the determination of c, as will be seen in the next sections, the optimization problem can be stated as follows:

Problem 1 Maximize c subject to:

(V (x) -c) + p i (x)g i (x) -ε i ≥ 0 i = 1, . . . , N, (3)
where p i (x) are unknown semi-definite positive functions and ε i > 0; i = 1, . . . , N . The purpose of constraint ( 3) is to validate the Theorem 1 hypothesis. To observe this, notice that at the boundary of the set Ω c , V (x) = c and, thus, the above constraints are reduced to p i (x)g i (x) ≥ ε i > 0. As functions p i ≥ 0, the constraints g i (x) ≥ 0 are satisfied at the points on the boundary of Ω c . Furthermore, in the interior of this set, V (x) -c < 0 and, thus, these constraints are also satisfied. The p i functions lend even more degrees of freedom, thereby increasing problem feasibility. The small ε i constants are pre-specified and are necessary in order to avoid problems at the points where p i (x) = 0. The introduction of ε i parameters constitute a new source of conservatism.

Remark 4 Notice that the expression of V is not considered in the optimization Problem 1. This fact may greatly reduce the complexity of the numerical problem.

In this work, sum-of-squares optimization is used order to solve this problem. For this, a new assumption is needed.

Assumption 3 Functions f (x, u) and g(x) are polynomial.

Sum of squares optimization

Sum of squares optimization is an optimization technique based on the Sum Of Squares (SOS) decomposition for multivariate polynomials. A multivariate polynomial p(x) is said to be an SOS, if there exist polynomials f 1 (x), ..., f m (x), such that:

p(x) = m i=1 f 2 i (x)
and therefore, p(x) ≥ 0 [START_REF] Prajna | Sum of Squares Optimization Toolbox for MATLAB[END_REF].

An SOS program has the following form [START_REF] Prajna | Sum of Squares Optimization Toolbox for MATLAB[END_REF]:

Minimize the linear objective function

w T c,
where c is a vector formed from the (unknown) coefficients of:

• polynomials p i (x), for i = 1, 2, ..., N 1

• sum of squares p i (x), for i = N 1 + 1, ..., N 2 such that:

g 0,j (x) + N i=1 p i (x)g i,j (x) = 0 for j = 1, 2, , . . . , M 1 . g 0,j (x) + N i=1 p i (x)g i,j (x) are SOS, for j = M 1 + 1, . . . , M 2 .
where w is the linear objective function weighting coefficients vector, and g i,j (x) represent certain scalar constant coefficient polynomials.

Currently, SOS programs are solved by reformulating them as semi-definite programs (SDPs), which in turn are solved efficiently, e.g., using interior point methods. Several commercial as well as non-commercial software packages are available for solving SDPs. SOSTOOLS [START_REF] Prajna | Introducing SOSTOOLS: a general purpose sum of squares programming solver[END_REF] is a Matlab toolbox that performs this conversion automatically, calls the SDP solver, and converts the SDP solution back to that of the solution of the original problem.

The problem stated in the previous section can be addressed solving the following SOS problem: Problem 2 Maximize c subject to:

(V (x) -c) + p i (x)g i (x) -ε i are SOS; i = 1, . . . , N, (4)
where p i are unknown SOS polynomials. This problem is more restricted than that presented in the previous section. Nevertheless, any solution to SOS problem 2 is a solution to problem 1.

Remark 5 Assumption 3 can be relaxed since other types of functions, such as trigonometric functions [START_REF] Papachristodoulou | On the construction of Lyapunov functions using the sum of squares decomposition[END_REF] or rational functions (e.g., the application examples in the next section) can be considered.

Some practical cases: electronic converters

The method developed above is applied to two practical cases: a DC-DC boost converter and a boost inverter. The control laws used for these systems appear in [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF] and [START_REF] Albea | Control of the Boost DC-AC Converter by Energy Shaping[END_REF] respectively where the global stability analysis is performed neglecting saturations and forbidden regions.

Application 1: Oscillations in a DC-DC boost converter

A boost circuit is usually used as a DC-DC converter and is of particular interest due to the fact that it generates an output voltage larger than its input voltage. Nevertheless, in [START_REF] Fossas-Colet | Asymptotic tracking in DC-to-DC nonlinear power converters[END_REF][START_REF] Fossas | Iterative approximation of limit cycles for a class of Abel equations[END_REF] it is proposed to convert DC voltage into an oscillating one. However, alternating current cannot be generated with this circuit, as the output current cannot change its sign. In this section the attraction domain is analyzed when this converter is controlled with the law presented in [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF]. For completeness, the design is summarized here. The normalized model of the DC-DC boost converter is:

ẋ1 = -ux 2 + 1 (5) ẋ2 = ux 1 -ax 2 (6)
where x 1 is the normalized current and x 2 is the normalized capacitor voltage defining the system states and where a is a constant and u is the control input, which has a direct relationship with the average duty-cycle. System ( 5)-( 6) can therefore be assumed to be a continuous system [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF]. Under the assumption that a is a known constant, a nonlinear control law based on the Hamiltonian approach is proposed in [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF]. The design is based on the following change of coordinates:

ζ 1 = x 2 1 + x 2 2 2 ζ 2 = x 1 -ax 2 2 + ζ 20
The aim of the control design is to ensure that the following function:

Γ = ω 2 (ζ 1 -ζ 10 ) 2 + (ζ 2 -ζ 20 ) 2 -µ
tends to zero. This goal corresponds to an ellipse in the current-voltage plane and is assumed that correspond to the desired behaviour. Notice that this desired behaviour is a limit cycle. Based on this definition, the nonlinear control law proposed in [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF] can be expressed as follows:

u = 1 + 2a 2 x 2 2 + kΓ(ζ 2 -ζ 20 ) + ω 2 (ζ 1 -ζ 10 ) x 2 (1 + 2ax 1 ) (7)
where k > 0 defines the speed of the transient response.

In [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF] it is proven that, with this control law, for all initial conditions except the origin, the trajectories of the system tend to the curve Γ = 0. Nevertheless, there are several constraints in the state variable that make this analysis useless from the practical point of view. These constraints are of several types: C1. Constraints 0 ≤ u ≤ 1 make control law [START_REF] Chesi | Solving quadratic distance problems: an LMI-based approach[END_REF] unfeasible in the full state space. In practice, u is chosen equal to 1 when control ( 7) is greater than 1 and, u is chosen equal to zero when control ( 7) is negative. Therefore, this constraint is 'soft' in the sense that if the system arrives at a point where the constraints are violated, the analysis of [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF] is no longer valid for the system with constraints. The point may nonetheless still lie inside the attraction domain of the desired limit cycle. Therefore, these constraints are saturation-like.

C2. Capacitor voltage cannot be negative in this circuit, which implies x 2 ≥ 0. This is considered to be a 'hard' constraint, meaning that such a situation should always be avoided.

C3. Finally, the control law is not feasible when the denominator in [START_REF] Chesi | Solving quadratic distance problems: an LMI-based approach[END_REF] is zero. This constraint is actually contained in C1, since when the denominator is close to zero the control signal assumes large (positive or negative) values.

The objective of this work is to obtain a (possibly conservative) estimation for the region of attraction of the resultant system taking these physical constraints into account.

The stability proof in [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF] for the unconstrained (fictitious) case is based on LaSalle's invariance principle. The Lyapunov function used is:

V = Γ 2 2 . (8)
The constraints are (only constraints C1 and C2 are presented here; constraint C3 will be discussed later):

u(x) ≤ 1 (9)
u(x) ≥ 0 (10)

x 2 ≥ 0 (11)
The expression for u, which is given by ( 7) is not polynomial but rather a rational function. Nevertheless, by expressing it as a quotient u(x) = n(x)/d(x) polynomials, all the constraints can be formulated in standard form. It can then be assumed that polynomial d(x) does not vanish at any point along the objective curve Γ(x) = 0. Otherwise, control law [START_REF] Chesi | Solving quadratic distance problems: an LMI-based approach[END_REF] is not valid for this problem. Therefore, the sign of d(x) is constant along Γ(x) = 0 and, through continuity in the neighborhood of this curve. By means of numerical inspection, it can be verified that, for the circuit parameters given below, d(x) > 0 on Γ(x) = 0. With this consideration in mind, constraints ( 9)-( 11) can be written as polynomial constraints:

• d(x) -n(x) ≥ 0 • n ≥ 0 • x 2 ≥ 0
Therefore, the problem to be solved is as follows:

Minimize (-c) [START_REF] Huber | A polyhedral method for solving sparse polynomial systems[END_REF] subject to:

(V (x) -c) + p 1 (x)(d(x) -n(x)) -ε 1 ≥ 0 (V (x) -c) + p 2 (x)n(x) -ε 2 ≥ 0 (V (x) -c) + p 3 (x)x 2 -ε 3 ≥ 0 p 1 (x), p 2 (x), p 3 (x) are SOS
Notice that constraint C3 is considered in the previous set of constraints. Indeed, constraint [START_REF] Hassan | Nonlinear Systems[END_REF] implies n(x) ≥ ε 2 > 0 for V (x) ≤ c, while constraint (13) implies d(x) ≥ n(x) + ε 1 for V (x) ≤ c. This implies that d(x) > 0 in V (x) ≤ c.

The following analysis can be directly modified for the case when d(x) < 0.

Results

The values of the circuit parameters are taken from [START_REF] Pagano | Autonomous oscillation generation in the boost converter[END_REF], resulting in a value for parameter a equal to 0.9045. The desired output AC voltage of the boost circuit is: v * conv = 135 + 15 sin 2π50t Software SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF] was used as the SDP solver under SOSTOOLS. The values for parameters ε i are chosen equal to 10 -6 while the chosen order for the unknown polynomials p i is 6. The solution obtained is: c * = 208.5. This value is close to the optimal value as verified by numerical inspection and thus confirms that for the point x (1) = (6.9 1.1) , which corresponds to V (x (1) ) = 355, the control signal u takes a value equal to 1.0001.

Nevertheless, the result is probably conservative as was pointed out in Remark 3.

Figure 2 shows the resultant estimation of the attraction domain. Note, this attraction domain is presented in the voltage-current plane. Also shown in Fig. 2 is the smaller Lyapunov level curve which inscribes the forbidden point given before (x (1) ), as well as the desired limit cycle Γ = 0. ; Lyapunov level curve (dashed) corresponding to the a forbidden point x (1) ; and desired limit cycle (dotted).

Notice that the method proposed in this paper may be valid when the system presents complex fractional equations.

Application 2: boost inverter

A means of obtaining AC output using boost topology a double boost is proposed in [START_REF] Ro Cáceres | A boost DC-AC converter: analysis, design, and experimentation[END_REF]. In this section the attraction domain is analyzed when this converter is controlled with the law presented in [START_REF] Albea | Control of the Boost DC-AC Converter by Energy Shaping[END_REF]. The design is summarized here, with the normalized model of this system being:

ẋ1 = -u 1 x 2 + 1 ẋ2 = u 1 x 1 -a(x 2 -x 4 ) ẋ3 = -u 2 x 4 + 1 ẋ4 = u 2 x 3 + a(x 2 -x 4 )
The nature of the constants and variables here are similar to that in the case before: x 1 , x 3 are the normalized currents, x 2 , x 4 are the normalized capacitor voltages and u 1 , u 2 are the control inputs. Note that this system is of greater dimensions than the system described in example 1.

The control law designed in [START_REF] Albea | Control of the Boost DC-AC Converter by Energy Shaping[END_REF] for the unconstrained (fictitious) case employs the following Lyapunov function:

V = Γ 2 1 2 + Γ 2 2 2 with Γ 1 = ω 2 (ζ 1 -ζ 10 ) 2 + (ζ 2 -ζ 20 ) 2 -µ Γ 2 = ω 2 (ζ 3 -ζ 30 ) 2 + (ζ 4 -ζ 40 ) 2 -µ
In [START_REF] Albea | Control of the Boost DC-AC Converter by Energy Shaping[END_REF], it is proved that V < 0 by LaSalle's invariance principle.

The corresponding control laws are

u 1 = 1 + 2a 2 x 2 2 -3a 2 x 2 x 4 + a 2 x 2 4 + a 2 x 2 ẋ4 x 2 + 2ax 1 x 2 -ax 4 x 1 + k 1 Γ 1 (ζ 2 -ζ 20 ) + ω 2 (ζ 1 -ζ 10 ) x 2 + 2ax 1 x 2 -ax 4 x 1 u 2 = 1 + 2a 2 x 2 4 -3a 2 x 2 x 4 + a 2 x 2 2 + a 2 x 4 ẋ2 x 4 + 2ax 3 x 4 -ax 2 x 3 - k 2 Γ 2 (ζ 4 -ζ 40 ) -ω 2 (ζ 3 -ζ 30 ) x 4 + 2ax 3 x 4 -ax 2 x 3
This system has the same constraints as those of the previous application: C1, C2 and C3. It should be noted that in this case there are two control inputs (u i , i = 1, 2) and as before, two capacitor voltages (x i , i = 2, 4), and numerical inspection reveals the case in which d(x) > 0 is the one that should be analyzed. The optimization problem results Minimize (-c) [START_REF] Hassan | Nonlinear Systems[END_REF] subject to:

(V (x) -c) + p 1 (x)(n 1 (x) -d 1 (x)) -ε 1 ≥ 0 (V (x) -c) + p 2 (x)(n 2 (x) -d 2 (x)) -ε 2 ≥ 0 (V (

Results

The parameter values are taken from [START_REF] Albea | Control of the Boost DC-AC Converter by Energy Shaping[END_REF]. The desired output voltage is:

v * inv = 40 sin 2π50t The tuning parameters ε i are equal to 10 -6 and the degree of the unknown p i polynomials is assigned a value of 3.

The solution is obtained in approximately ten minutes on a PC (1.66 GHz Intel Core2): c * = 1.5. As a means of corroborating this result, a numerical inspection has found that for a value of x (2) = (0 -0.1 0.4 8.1) , which corresponds to V (x (2) ) = 2.38, the constraint x 2 ≥ 0 is violated.

Conclusions

A methodology for the estimation of the region of attraction that takes system physical constraints into consideration has been presented in this article. The kind of problems considered may present a high degree of difficulty due to the model and control law nonlinearities including saturation-like constraints. The method is based on the search for a Lyapunov level surface where the constraints are satisfied. This problem can be written as a sum of squares optimization problem, for which good numerical tools are available. This approach has general applicability to cases where the global stability of the unconstrained problem has been demonstrated (by means of Lyapunov methods) and in which extension to the constrained case is desired. This method is useful even when the degree and complexity of the equations are high. The closed-loop model and constraints must be polynomial or rational (there are cases, nevertheless, where SOS programming has been applied to trigonometric and other terms [START_REF] Papachristodoulou | On the construction of Lyapunov functions using the sum of squares decomposition[END_REF]). Conservativeness of the method has also been discussed.

The usefulness of the method has been shown with two electronic applications.

Figure 1 :

 1 Figure 1: Estimated attraction domain where the constraints are not active.

Figure 2 :

 2 Figure2: Example 1: Estimated attraction domain (solid); Lyapunov level curve (dashed) corresponding to the a forbidden point x(1) ; and desired limit cycle (dotted).

  x) -c) -p 3 (x)n 1 (x) -ε 3 ≥ 0 (V (x) -c) -p 4 (x)n 2 (x) -ε 4 ≥ 0 (V (x) -c) + p 5 (x)x 2 -ε 5 ≥ 0 (V (x) -c) + p 6 (x)x 4 -ε 6 ≥ 0 p 1 (x), p 2 (x), p 3 (x), p 4 (x), p 5 (x), p 6 (x) are SOS

Acknowledgments

The authors are grateful to the anonymous reviewers of this paper for their valuable comments and suggestions. This research has been partially supported by the MICINN-FEDER grant DPI2009-09961.