José Raul Azinheira

Patrick Rives
email: patrick.rives@sophia.inria.fr

IMAGE-BASED VISUAL SERVOING FOR VANISHING FEATURES AND GROUND LINES TRACKING: APPLICATION TO A UAV AUTOMATIC LANDING

Keywords: Image-based visual servoing, Unmanned aerial vehicle (UAV), auto-landing, dynamic modeling, robust control

A unified modeling framework is defined and an image-based visual servoing scheme is introduced to evaluate the feasibility of a vision-based guidance for a UAV model during survey missions and for automatic landing. Relevant image features are selected to allow for the proposed objective, namely decoupling rotation and translation in the image, and allowing to respect the natural separation between the lateral and longitudinal motions of the aircraft. The controller design is adapted to include the aircraft dynamics and the vision output in an linear optimal control approach. Simulation results are presented in order to validate the proposed solution and allow for a critical evaluation.

NOMENCLATURE

∆ l , ∆ r , ∆ H left, right and horizon lines (δ a , δ e , δ r , δ T ) aileron, elevator, rudder and throttle deflections Ω = [p, q, r] T angular velocity in body frame (φ, θ, ψ) Roll, pitch and yaw Euler angles The main objective of the present work is the search for a feedback control strategy applied to an unmanned air vehicle (UAV) using the image provided by an airborne camera. The application to the automatic landing of an aircraft is the more demanding flight phase, both in terms of required accuracy and necessary robustness.

A specific characteristic of the landing (partly justifying the previous statement) is the fact that the relevant position variables are relative to a local ground frame and not to an absolute frame. Such a characteristic is also true for a large class of survey missions where the control objective is defined in terms of a relative positioning w.r.t. structures lying on the ground, like roads, rivers or a landing strip.

The landing problem may thus be considered as the extreme case of a need for perfect path tracking in a ground relative frame, for which a vision sensor (an airborne camera) may be regarded as a sensor with adequate characteristics, useful for the precise positioning relatively to the ground.

In terms of control strategy, two options are usually available:

• a first solution, referred as the position-based visual servoing (PBVS) or reconstruction approach, consists in considering the image output as another sensor specifically used to help for the landing phase in order to get a better estimate of the vehicle position, as compared to a desired path. The aircraft position w.r.t the ground is estimated from the image using vanishing points or lines or by computing an homography under the assumption that the ground plane is flat [START_REF] Rathinam | Vision based following of locally linear structures using an unmanned aerial vehicle[END_REF], [START_REF] Saripalli | Visually-guided of an unmanned aerial vehicle[END_REF], (Templeton et al. 2007).The problem of control is then stated as a trajectory tracking in the cartesian space.

• a second solution, so called image-based visual servoing (IBVS) approach, consists in switching from the flight control scheme to a landing scheme where the vehicle motion is only referenced to the image as compared to a reference image to be tracked. Until now, only a few IBVS approches have been applied to the control of aerial unmanned vehicles (see, for example, the USC Autonomous Flying Vehicle Project [START_REF] Mejias | A visual servoing approach for tracking features in urban areas using an autonomous helicopter[END_REF], (Silveira et al 2003).

The position-based visual servoing approach is a more conservative solution, probably easier to implement:

the flight control scheme is almost unchanged, only switching to an image enhanced estimation of the position.

On the other hand the reconstruction relies on a good calibration of the camera and the effect of the image noise or inaccuracy must be secured.

Some papers already touched the problem of vision based flight control, like (Dickmanns 1994), (Schell and[START_REF] Dickmanns | The 4D-approach to dynamic machine vision[END_REF], [START_REF] Chatterji | Vision based position and attitude determination for aircraft night landing[END_REF] or [START_REF] Kimmett | Vision based navigation system for autonomous aerial refueling[END_REF]. The first and the second papers The last paper presents the simulation results of a reconstruction control scheme for in-flight refueling.

The image-based visual servoing approach is more ambitious and a cautious analysis is required before it may be tested in real flight conditions. It must be secured that the control in the image does not excite the unstable or marginally stable dynamics of the vehicle. It is however a solution that is not so dependent on the camera calibration, as long as an acceptable reference image is provided for the visual tracking [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], [START_REF] Rives | Sensor-based tasks: From the specification to the control aspects[END_REF]. In contrast with the usual image-based control in robotics, in our particular case the target system may not be assumed as a pure velocity integrator, the aircraft is underactuated and the dynamics and couplings between the axes may not be neglected and are to be included into the controller design. Another approach to overpass this problem is presented by [START_REF] Hamel | Visual servoing of an under-actuated dynamic rigid-body system: An image-based approach[END_REF], with a dynamical visual servo to stabilize a scale helicopter over a landing pad and which treats image measurements as unit vectors defining a mapping from the image plane to the unit sphere which allows a decoupling between translations and rotations.

The approach presented here was originally introduced in [START_REF] Rives | Linear structures following by an airship using vanishing point and horizon line in a visual servoing scheme[END_REF] for the tracking of ground features with an airship. The scheme is adapted in the present paper in order to allow the automatic landing of an unmanned aircraft. A similar approach was recently presented in [START_REF] Bourquardez | Visual servoing of an airplane for alignment with respect to a runaway[END_REF] but dealing with the alignment phase and not considering the wind disturbance.

The present paper is structured as follows. In a first section, we briefly present the dynamic model of the aircraft, to be used for control design. The second section is devoted to the image modeling, introduces the visual servoing approach and the control design. The third section presents simulation results, which allow to evaluate the behaviour of the controled UAV and validate the approach. Some brief conclusions are finally drawn in the last section.

DYNAMIC MODELING OF THE AIRCRAFT

The basic equations describing the aircraft motion in atmospheric flight may be found in the bibliography on Flight Mechanics, by instance in [START_REF] Mclean | Automatic Flight Control Systems[END_REF], (Stevens and Lewis 1992), [START_REF] Phillips | Mechanics of Flight[END_REF]).

The present section only gathers the basic notations and equations further used in the simulation setup or for the control design phase.

Frames and Notations

The aircraft trajectory and the landing road coordinates are given in the earth frame (or NED, for North-East-Down), with the center by instance on the road axis (Figure 1) -for simplicity and without loss of generality,the road is chosen as aligned with the North axis. The local frame, linked to the aircraft (ABC, aircraft body centered), where the aircraft velocity V = [u, v, w] T is given, is centered at its center of gravity, u directed towards the aircraft nose, v towards the right wing and w downwards. The angular velocity is also expressed in the local frame: Ω = [p, q, r] T . The airborne camera is rigidly attached to the aircraft, with its optical axis aligned with the aircraft longitudinal axis. The change from earth frame to local frame is defined by the transformation matrix S, which may as usual be stated in terms of the Euler angles (φ, θ, ψ), respectively roll, pitch and yaw angles.

Flight Mechanics and Dynamic Modeling

The mechanics of atmospheric flight is generally deduced from the application of Newton's second law, considering the motion of the aircraft in the earth frame, assumed as an inertial frame, under the influence of forces and torques due to gravity, aerodynamics and propulsion.

If the aircraft motion is described as usual in the body frame by its inertial and angular velocity, and by the Euler angles, the dynamic system may be put in a synthetic form as the following state space equation:

Ẋ = f (X, U, D) (1)
where:

• Ẋ = dX dt is the usual notation for time derivation

• X = [u, v, w, p, q, r, φ, θ, ψ] T is the state vector,

• U = [δ a , δ e , δ r , δ T ] T is the input vector, with respectively the aileron, elevator and rudder deflections and the throttle input,

• the disturbance vector is usually representing the wind velocity, with its six components

D = [u w , v w , w w , p w , q w , r w ] T
The state equation (Eq.1) is established for a nearly constant air density and the aircraft motion only depends on its velocity and attitude (angular position) and the aircraft cartesian position P = [N, E, D] T in the earth frame does not appear. Since this position is to be controlled, three position states are added, as integrators of the aircraft inertial velocity V, considered in the earth frame:

Ṗ = S T V (2)

Linearized Models

The equations non-linearity, their complexity and a certain level of approximation in the aircraft models have justified the search for simplified versions and, as a first step, it is common to linearize the equations for small perturbations around an equilibrium flight. This equilibrium or trim flight is frequently taken as a horizontal straight leveled flight, with no wind.

Under these conditions, the equations are written as functions of the perturbations in the state vector x, in the input vector u or in the disturbance vector d, resulting in two differential matrix equations describing the dynamics of two independent (decoupled) motions:

ẋv = A v x v + B v u v + E v d v ẋh = A h x h + B h u h + E h d h (3)
where:

• the index v is associated with the longitudinal motion, in the vertical plane, and

• the index h is associated with the lateral/directional motion, mostly in the horizontal plane.

The longitudinal state vector is x v = [u, w, q, θ] T and the input vector is u v = [δ e , δ T ] T (where all the variables are changes from the trim value). In the lateral case, the state vector is x h = [v, p, r, φ] T and the input vector

u h = [δ a , δ r ] T .
The linear models described by Eq. 3 depend on the trim point chosen for the linearization: namely they are function of the airspeed V o and the altitude h o . The validity of these linear models is obviously limited to a surrounding near the trim point, namely by the linearization of the angles, and by the validity of the aerodynamic laws, within the flight envelope, which state that the airspeed must remain greater than the stall speed and lower than a maximum allowed airspeed. The inputs are also limited, in value or rate, both for the surface deflections and for the thrust input.

Experimental Setup: the ARMOR X7 UAV

The modelling process described above has been applied to the experimental UAV ARMOR X7 currently under development at IST in Lisbon (Figure 2). The half-scale flying model in use in our early experiments has 3 m wing span and 18 kg nominal weight. Its cruise airspeed is about 18 m/s.

The fully non-linear dynamic model (Eq. 1) was defined in order to enable a realistic simulation, including the effect of wind and atmospheric turbulence, as well as the ground effect when the aircraft is near touchdown.

A Matlab/Simulink simulation platform was developed to test control solutions and evaluate strategies for the desired autonomous or semi-autonomous operation [START_REF] Costa | Modelação da aeronave robotizada armor x7[END_REF]). The linearized model (Eq. 3) is used in the control design phase, as by instance in [START_REF] Azinheira | Longitudinal auto-landing for a low cost uav[END_REF].

The simulation platform allows to handle simple models of the 3D scene. A ground road was defined, with a width of 5 m, and a simulated camera was introduced, which outputs the visual signals to be used in the image-based visual servoing.

IMAGE-BASED VISUAL SERVOING

In opposite to a 3D visual servoing method [START_REF] Dickmanns | The 4D-approach to dynamic machine vision[END_REF], [START_REF] Furst | A vision based navigation system for autonomous aircraft[END_REF], an imagebased visual servoing does not require an explicit 3D reconstruction of the scene. The basic idea is to assume that a task can be fully specified in terms of the desired configuration of a set of geometric features in the image. The task will be perfectly achieved when such a configuration is reached [START_REF] Samson | Robot control: the task function approach[END_REF]).

In terms of control, this can be formulated as a problem of regulation to zero of a certain output function directly defined in the image frame.

Let us consider the airborne camera C, which can be viewed as a mechanical system with several actuated degrees of freedom. The pose (position and orientation) of C is an element r of Ê 3 × ËÇ 3 , which is a six dimensional differential manifold. C interacts with its environment. We assume that the image given by C (see Figure 3) fully characterizes the relative position of C with respect to the NED frame attached to the scene. Moreover, let us consider that the information in the image may be modeled as a set of visual signals characterizing the geometric features which result from the projection onto the image of the 3D objects belonging to the scene. Each elementary signal s i (r) defines a differentiable mapping from Ê 3 × ËÇ 3 to Ê.

As shown in [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], the differential of this mapping (so-called interaction matrix) L T i relates the variation of the signal s i (r) observed in the image to the motion between the camera and the 3D target expressed by the camera velocity screw

T CT . ṡi = L i T T CT (4)
An analytical expression for the interaction matrix when the image features are general algebraic curves can be derived (for more details, see [START_REF] Chaumette | Visual Servoing, Chapter Classification and realization of the different vision-based tasks[END_REF], [START_REF] Rives | Sensor-based tasks: From the specification to the control aspects[END_REF]). In our peculiar case, we consider the set of geometric primitives in the 3D scene will be constituted by points and straight lines.

Modeling the Interaction Matrix

Let us assume hereafter that we use a pinhole camera model with a focal length equal to 1 and that both the points and lines in the 3D scene and their projection in the 2D normalized image plane are expressed in the camera frame F s .

Case of points:

Any point M with coordinates X = (X, Y, Z) projects onto the image plane as a point m with coordinates

x = (x, y) such that:

x = X/Z , y = Y /Z (5)
By differentiating Eq.( 5), it is obvious to compute the interaction matrix linking the 2D motion observed in the image to the camera motion in the 3D scene.

    ẋ ẏ     =     -1 Z 0 x Z xy -1 -x 2 y 0 -1 Z y Z 1 + y 2 -xy -x     T CT (6)

Case of straight lines:

A straight line in the 3D scene is here represented as the intersection of two planes described in the implicit form h(X, Q) = 0 such that :

h(X, Q) =        a 1 X + b 1 Y + c 1 Z = 0 a 2 X + b 2 Y + c 2 Z + d 2 = 0 (7)
with d 2 = 0 in order to exclude degenerated cases. In these equations, X = (X, Y, Z, 1) denotes the homogeneous coordinates, expressed in the camera frame, of the 3D points lying on the 3D line, and Q denotes a parameterization of the 3D lines manifold.

The equation of the 2D projected line in the image plane (see Figure 3) can also be written in an implicit form g(x, q) = 0 such that :

g(x, q) = x cos θ + y sin θ -ρ = 0 (8) with                cos θ = a 1 / a 2 1 + b 2 1 sin θ = b 1 / a 2 1 + b 2 1 ρ = -c 1 / a 2 1 + b 2 1
where x = (x, y, 1) denotes the homogeneous coordinates in the image of the 2D points lying on the 2D line, and q denotes a parameterization of the 2D lines manifold, hereafter the polar representation (ρ, θ).

A general form of the interaction matrix may then be obtained for each line:

L T θ = [ λ θ cos θ λ θ sin θ -λ θ ρ . . . -ρ cos θ -ρ sin θ -1 ] L T ρ = [ λ ρ cos θ λ ρ sin θ -λ ρ ρ . . . (1 + ρ 2 ) sin θ -(1 + ρ 2 ) cos θ 0 ] (9) with λ θ = (a 2 sin θ -b 2 cos θ)/d 2 and λ ρ = (a 2 ρ cos θ + b 2 ρ sin θ + c 2 )/d 2 .

Modeling the Visual Signals

A first step in order to include vision in the control loop is to define a reference scene, the image of which, as viewed from an airborne camera, would allow for a good vertical and lateral or attitude positioning, but leaving freedom enough to cope with the limitations of the vehicle dynamics.

As presented above (see Figure 1), let us consider a scene composed by a strip lying on the ground. We assume a smooth and limited curvature (piecewise linear), and the strip is parametrized by two parallel curves (i.e. the two sides of the road or river). These curves project onto the image as shown in Figure 3.

Let us define in the image the two tangents ∆ r and ∆ l to the right and left border lines of the road at a given image coordinate y T . ∆ r and ∆ l converge to a vanishing point x H which belongs to the horizon line ∆ H . The parameters of the lines ∆ r , ∆ l and ∆ H and the coordinates of the vanishing point x H in the image depend on the relative position and attitude of the camera frame w.r.t. the road.

In terms of a survey task for an aerial vehicle, following a road on the ground consists in keeping the longitudinal axis of the vehicle collinear with the road median axis and, simultaneously, keeping the height over the road constant. The longitudinal speed may be controlled separately and tunes the survey task.

These constraints between the reference ground frame and the camera frame have to be translated in terms

Z 2i Z 1i x 2i = ∞ 2 H 1 x 1i + 1 Z 1i c 1 ( 10
)
where

∞ 2 H 1 = K 2 R 1 K -1
is the homography of the plane at infinity and c 1 ∝ Kt is the epipole in the first image (K is the calibration matrix of the camera).

Let us consider now peculiar points in the scene which belong to the plane at infinity (in our case, the vanishing point and the horizon line), then Z 1i = Z 2i = ∞ and Eq. 10 yields:

∞ x 2i = K 2 R 1 K -1 ∞ x 1i (11)
From this equation, it appears that the observed motion in the image of the points which belong to the plane at infinity only depends on the rotation part of the camera displacement and are independent on the translations. Thus, controling such points in the image-based visual servoing allows a perfect decoupling of rotations and translations of the camera. To control the three rotations, we have chosen to use (see Figure 3)

as visual signals s 1 = [x H , θ H ] T ,
ṡ1 = 0 L T rot1 T CT (12)         ẋH ẏH θH         =         x H y H -(1 + x 2 H ) y H (1 + y 2 H ) -x H y H x H -ρ H cθ H -ρ H sθ H -1                 ω x ω y ω z        
where cθ = cos θ, sθ = sin θ and (ω x , ω y , ω z ) is the camera angular velocity.

Control of translations :

As was said at the beginning of this section, in a path tracking task, we consider that the longitudinal speed V x of the aerial vehicle does not require to be controlled from the vision task (the longitudinal speed is indeed to be regulated according to the vehicle requirements expressed in terms of airspeed or angle of attack). So, we want to select visual signals adequate to control the lateral and vertical translations. That can be done by relating the projection of the border lines observed in the image to the lateral and vertical position of the camera computed at the desired reference attitude (i.e. aligned with the road) (Figure 4).

Assuming the width L of the road is known and the optical axis of the camera is aligned with the road, we can compute the angles of the lines ∆ r , ∆ l , depending on the altitude h and the lateral position error e of the camera:

tan (θ r ) = L+2e 2h tan (θ l ) = -L+2e 2h ( 13
)
From this equation, a good choice of visual signals for controling the altitude and the lateral position of the camera will be:

s 2 =         t m t d         =         (tan θr+tan θl) 2 (tan θr-tan θl) 2         =         e h L 2h         (14)
The two lines are lying on the ground plane at a finite distance, thus the corresponding interaction matrix will depend both on rotations and translations. Derivating the tangent :

s = tan θ ⇒ ṡ = 1 cos 2 θ θ (15)
and using the general form (Eq. 9), yields an analytical form of the interaction matrix :

ṡ2 = L T tran2 L T rot2 T CT (16)
Finally, combining Eqs. 12 and 16, we obtain the global interaction matrix L T , which is a lower triangular matrix with good decoupling properties :

ṡ = L T T CT ṡ =     ṡ1 ṡ2     =     0 L T rot1 L T tran2 L T rot2         V CT Ω CT     (17)
Computing the desired visual signals s * = [s * 1 s * 2 ] T : Since a control in the image is used, the path reference to be tracked by the controller is converted into an image reference to be compared with the visual output and then the error is used by the controller. So, we need to express the Cartesian path tracking task in terms of image features trajectory. Let us consider a road following task at a constant altitude h * , centered w.r.t the road (e = 0) and a constant longitudinal speed V * x . We want to be aligned along the median axis of the road, which means in the image:

θ * = -θ * r = θ * l ,
x * H = 0 and θ * H = 0. If, for simplicity, we assume also that the camera optical axis z c is horizontal and aligned with the road axis, then y * H = 0, ρ H = ρ l = ρ r = 0, and the equation of the ground plane is such that (a = 0, b = 1, c = 0, d = h * ). Using these values, the desired visual signal is s * = (0, 0, 0, 0, θ * ) T and the interaction matrix L T s=s * computed at the desired position is :

L T s=s * =                 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 0 0 0 1 1 h * 0 0 0 0 -1 0 L 2(h * ) 2 0 0 0 0                 (18)
Let us note that only one coupling between the translation and the rotation remains, in the control of the translation along the x-axis of the camera (y-axis of the aircraft). We can also verify that the translation along the optical z-axis of the camera (x-axis of the aircraft) is not controlled by the visual output s.

Controller Design

The image-based auto-pilot is implemented according to the Block Diagram of Figure 5, with the use of an airspeed sensor and a camera as only feedback sensors, and with a controller regulating the airspeed and tracking the image reference constructed as above.

The idea of a visual control for an unstable platform of 12th order is challenging but implies a great concern with the robustness of the solution. As a first tentative, a solution is searched using optimal control, based on the linearized model of the vehicle dynamics and image output (Eqs. 3 and 17), looking for a pure gain applied on the measured output error:

u = D (y * -y) (19)
where the output y = [u, x H , y H , θ H , t d , t m ] T includes the longitudinal speed u and the visual output s.

The non-linear dynamics of the air vehicle (Eq. 1) is first linearized as was described above, around a trim equilibrium state corresponding to a stabilized leveled flight at a constant airspeed (V o ) and at a constant altitude (h o ) above the road and aligned with the road axis, the longitudinal and lateral variables of Eq. 3 are joined in a single state, and the deterministic case is considered:

ẋ = Ax + Bu (20)
The state x = v T , p T T of this model includes the 6D change in velocity v = [u, v, w, p, q, r] T and 6D change in position p = [n, e, d, φ, θ, ψ] T of the aircraft frame with respect to the ground frame. The input (u) includes the 3 control surface deflections (δ a , δ e , δ r ) and the engine thrust change (δ T ).

The visual output is also linearized, for the same trim condition, using the Jacobian of the image function (or interaction matrix Eq. 18) and including the change from aircraft frame to image frame (S a ):

s = L T (h o )S a p (21)
which, together with the longitudinal speed change (u), gives the output equation:

y = C (h o ) x =     C u 0 0 L T S a         v p     (22)
where C u = [1, 0, 0, 0, 0, 0] extracts the longitudinal velocity component. This output corresponds to the use of only two sensors: the camera, and an airspeed sensor as an approximate measure of the longitudinal speed

error u ≃ V t -V o .
For a fixed airspeed and altitude, the optimal state feedback gain of the LTI system is obtained with the standard Matlab LQR function, corresponding to the minimization of a cost function weighting the output error and control needs:

J ∞ 0 y T Qy + u T Ru dt (23)
through the definition of the appropriate weighting matrices Q and R.

The state feedback gain K is finally transformed into an output error feedback gain using the pseudo-inverse of the output matrix:

D = KC † (24)
PI controller :

In order to reduce the static error appearing in the previous solution with a constant wind disturbance1 , an error integrator y i was added in the previous procedure and an augmented output z is considered:

    ẋ ẏi     =     A 0 C 0         x y i     +     B 0     u z =     C 0 0 I         x y i     (25)
The resulting system was discretized and an optimal output error feedback gain was obtained, yielding the discrete controller:

x c k+1 = A c x c k + B c (z * k -z k ) u k = C c x c k + D c (z * k -z k ) (26)
Sliding Gain :

Since there is a change in the linearized system as altitude is changing, and namely because the Jacobian of the image and then C(h o ) is dramatically increasing when the vehicle is near touchdown, the interaction matrix is computed at each sample time according to the current desired altitude and the applied feedback gain is updated:

u k = C c k x c k + D c k (z * k -z k ) (27)
The weighting matrices are also corrected as altitude reduces in order to integrate the strict constraints near touchdown.

Obviously, the visual control is only valid in flight and is switched off once the aircraft has landed and is rolling on the ground.

SIMULATION RESULTS

In order to allow for a fair comparison, a similar setup was used for all simulations, in agreement with the Block Diagram presented above. The control was implemented with a 100 ms (10 Hz) sampling rate.

In agreement with the decoupling properties, the lateral and longitudinal behaviours are initially checked separately. The lateral tracking of the ground road is first presented, with a simulation at constant altitude; a pure longitudinal landing case is then analysed, and a realistic landing is finally considered, with a lateral initial error and wind disturbances.

For all the figures, the relevant parameters during visual control are in solid lines and in dotted lines before visual control or after touchdown. When they exist, the references are presented in dashed lines.

Lateral Tracking

The lateral behaviour of the visual tracking is analysed in a simulation with a 20 m constant altitude reference, a constant 16 m/s reference airspeed, and the simulation includes the correction of an initial alignment error of 18 m to the right of the road, and then the tracking of a road S turn (see Figure 6). Two sub-cases were considered: a nominal case without wind disturbance, and a wind case, with an intensity of 5 m/s, blowing with 15 • from the left of the path.

The horizontal path presented in Figure 6 clearly shows a smooth alignment and then a very good tracking of the road axis; the motion is fast and well damped. In the nominal case (left), the aircraft is perfectly aligned on the road axis and only the corners of the S shape are a little smoothed. On the other hand, the right figure shows the influence of the wind disturbance, introducing a static error on the horizontal tracking (in this case, an arrow was drawn to indicate the wind direction).

In Figure 7 are shown the time evolutions of horizontal and vertical tracking errors, as well as the corresponding visual signals t m and t d . The two sub-cases are very similar, the difference only revealing the static error, with 2.8 m in the lateral and 0.18 m in the vertical tracking.

The static error on the visual signal x H (see the windy case in Figure 8) is due to the characteristics of the aircraft dynamic model when it is flying with no side-slip (natural coordinated flight), expressing that an equilibrium flight above the road with a side wind component necessarily has a heading offset. The other visual signals are well regulated to zero but, as it is clear in Figure 7, the signal t m is also offset (with a value of 0.14), which corresponds to the lateral tracking error. A little coupling in the altitude is hardly visible in signal t d .

Longitudinal Landing

The following landing conditions were assumed:

• start at 700 m away from the touchdown point on the road; after a stabilization period, the landing control is switched on at 500 m from the desired touchdown, with alignment first and then descent;

• initial altitude at 20 m reference;

• initial airspeed at reference speed, equal to 16 m/s, which is an airspeed adequate for the approach phase, with a pitch attitude acceptable until touchdown (the model stall speed is slightly below 13 m/s).

Two altitude profiles were first considered:

• an usual linear descent at constant sinking speed, with a glide slope of 3 • , and flare for touchdown, with a final reduction of airspeed before ground contact;

• a cosine descent, varying continuously from the initial altitude to touchdown, with also an airspeed reduction for ground contact.

The simulation results comparing the nominal landings without wind and for the two altitude profiles are shown in Figure 9, with, from top to bottom, the altitude and airspeed curves, along with their reference profiles, the pitch angle and the two longitudinal inputs, thrust demand and elevator deflection.

The altitude profiles demonstrate a fair tracking of the reference, with a little lag at the start of descent, smoothly corrected before touchdown. The airspeed profiles exhibit the influence of the descent on the airspeed regulation, but the reference airspeed is still well responding for the deceleration before touchdown.

Both the flare in the altitude profiles and the airspeed curves seem to be more precise in the cosine case.

The pitch curves and the input curves clearly show the difference between the two profiles, with the cosine solution giving smoother curves but implying a steeper descent at mid altitude.

The relevant longitudinal visual signals, y H and t d , respectively associated with the pitch angle and altitude, are presented in Figure 10, along with their references. The characteristics are similar for the two descent profiles, maybe with a closer tracking in the linear case. The vanishing point coordinate exhibits a little overshoot at its maximum value, whereas the t d signal is fairly well tracked.

In terms of airplane automatic landing, the performance of this nominal case may be analysed through the impact vertical velocity (sinking speed) as presented in Figure 11, which is to be compared to the proposed UAVs regulation limit of 2 m/s (SBAC 1991). Both solutions are well inside the regulation limits, with an impact velocity near 0.2 m/s. The flare phase appears however very sudden for the linear case, whereas the cosine profile yields a very continuous curve till touchdown. According to these curves, and looking for a safer touchdown, the cosine profile was then selected for the more realistic landing simulations.

Realistic landing with wind

In order to have a first evaluation of the validity and robustness of the visual control scheme, a realistic windy landing simulation was run, again from a height of 20 m to touchdown, with a cosine descent and with an initial lateral error of 16 m to the right of the road. The wind conditions were defined with:

• a mean nose component with 5 m/s, with a 15 • angle to the left of the road (the wind intensity corresponds to 31% of the aircraft airspeed, regulated to 16 m/s, and is quite significant);

• plus an atmospheric turbulence component, simulated by a Dryden model, with an intensity of 3 which, in a scale from 0 to 7 m/s, corresponds to an intermediate gust case.

The simulation results are presented in Figure 12. To some extent, the characteristics of this simulation are in agreement with the pure longitudinal and lateral cases, with the coupling from lateral to longitudinal during the initial alignment (say between 500 m and 400 m before the touchdown point), and then the influence from altitude to the lateral tracking, visible in the diminishing static error in the lateral track error (top-right curve). Globally the landing is well behaved and smooth, with a lateral error of 0.9 m at touchdown and an impact velocity near 0.2 m/s (the roll and yaw angles are respectively -0.2 • and -1.2 • ).

The visual signals are presented in Figure 13, showing again characteristics similar to the lateral and longitudinal cases. The influence of the altitude on the lateral tracking is here more visible, mostly on the t m signal which seems to go out of control: remember however that a trade-off has to be made for touchdown, allowing some lateral error in order to ensure the aircraft attitude is acceptable and permits that the undercarriage touches the ground safely, and then the aircraft starts to roll along the road (using a specific ground controller). The influence of the atmospheric turbulence is also more visible in the visual signals, namely in the angular signals x H , y H and θ H , but this influence remains very little.

CONCLUSIONS

In the work described in this paper, a linear structures following task by an aircraft is proposed using image-based visual servoing. This control approach has been analysed in order to set a first exploratory evaluation of a visual servoing technique applied to the automatic landing of an unmanned aircraft (UAV).

The classical image-based approach was adapted to the specific case:

• an adequate scene and image features were selected to allow for the proposed objective, namely decoupling rotation and translation in the image, and allowing to respect the natural separation between the lateral and longitudinal motions of the aircraft;

• the controller design was defined to include the aircraft dynamic characteristics and a sliding gain optimal control was chosen as a first robust solution.

The simulations used to analyse the close loop characteristics and the behavior of the control solution show quite a good performance, well in agreement with the specifications, and the visual-servoing scheme seems clearly able to land the aircraft in nominal and intermediate wind conditions. The conclusion is then that the idea looks feasible, and clearly justifies further studies to complete the validation and eventually implement such a visual servoing scheme on the real aircraft.

Aknowledgments

  (ρ, θ) polar representation of straight line in image plane A, B, C matrices of linear model K state feedback gain for controller L T interaction matrix P = [N, E, D] T aircraft position in earth frame (North, East, Down) Q, R matrices of LQR cost function r 6D pose of camera w.r.t. earth frame S transformation matrix from earth to body frame u, v, w] T inertial velocity in body frame (x H , y H ) coordinates of vanishing point in image x, u, y, z state, input and output vectors of linear model X U state and input vector in dynamic model x = (x, y) projected point coordinates in image plane X = (X, Y, Z) point coordinates in image frame 1. INTRODUCTION

  are representative of an impressive amount of works on vision-based flight control developed at UniBwM (Universitaet der Bundeswehr Munich) since 1985. The authors use full state reconstruction by recursive estimation and state feedback techniques for the design of longitudinal and lateral flight controllers. In complement to image features, inertial data are also used to deal with gusts and wind inputs. The third paper presents a pose reconstruction technique to be used as a pilot aid during landing, with no automatic control.

  of reference visual signals to be observed in the image. In order to cope with the natural decoupling of the aircraft dynamics (see section 2), it is interesting to choose these visual signals such that rotations and translations are, as much as possible, decoupled.Control of rotations :In the peculiar case of vision-based control applied to an aerial vehicle, we can take advantage of the projective properties of the sensor. O.D Faugeras et al.[START_REF] Faugeras | The geometry of multiple images[END_REF] and R.[START_REF] Hartley | Multiple View Geometry[END_REF] have shown that multiple views provided by a moving camera are related by projective constraints such that any point in the first image has a correspondence in the second image lying on a so-called epipolar line. The fundamental relationship between the i-th point in two different images is:

Fig. 3 .Fig. 4 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .-Fig. 10 .Fig. 12 .Fig. 13 .

 346789101213 Fig. 1. Frames and notations

  where x H = [x H , y H ] T are the coordinates of the vanishing point in the image, and θ H is the angle of the horizon line ∆ H = (ρ H , θ H ) expressed in polar coordinates. Using visual signals s 1 in Eq. 6 yields the analytical form of the interaction matrix :

In the case of a lateral wind component, the error cannot be completely cancelled because the reference output z * is not a stable trim solution.


Bruno Sareni
email: sareni@leei.enseeiht.fr

Jérémi Regnier
email: regnier@leei.enseeiht.fr

Xavier Roboam
email: roboam@leei.enseeiht.fr

Recombination and Self-Adaptation in Multi-objective Genetic Algorithms

This paper investigates the influence of recombination and selfadaptation in real-encoded Multi-Objective Genetic Algorithms (MOGAs). NSGA-II and SPEA2 are used as example to characterize the efficiency of MOGAs in relation to various recombination operators. The blend crossover, the simulated binary crossover and the breeder genetic crossover are compared for both MOGAs on multi-objective problems of the literature. Finally, a selfadaptive recombination scheme is proposed to improve the robustness of MO-GAs.

Introduction

In recent years, extensive research has been done in the field of Multi-Objective Genetic Algorithms (MOGAs) [START_REF] Corne | The Pareto Envelope-based Selection Algorithm for Multiobjective Optimization[END_REF], [START_REF] Costa | An Evolution Strategy for Multiobjective Optimization[END_REF], [START_REF] Deb | A fast-elitist non-dominated sorting genetic algorithm for multiobjective optimization: NSGA-II[END_REF], [START_REF] Zitzler | SPEA2: Improving the Strength Pareto Evolutionary Algorithm[END_REF]. However, most of works have been focused on selection, elitism and niching operators. Although the need of studying the influence of recombination and self-adaptation in MOGAs has been underlined [START_REF] Laumanns | Approximating the Pareto Set: Concepts, Diversity Issues, and Performance Assessment[END_REF], only a few contributions on these issues have been carried out [START_REF] Costa | An Evolution Strategy for Multiobjective Optimization[END_REF], [START_REF] Laumanns | Mutation Control and Convergence in Evolutionary Multi-Objective Optimization[END_REF]. In this paper, we investigate the efficiency of three crossover operators for real-encoded MOGAs and propose a self-adaptive recombination scheme which improves their robustness. The second versions of the Non-dominated Sorting Genetic Algorithm (NSGA-II) [START_REF] Deb | A fast-elitist non-dominated sorting genetic algorithm for multiobjective optimization: NSGA-II[END_REF] and of the Strength Pareto Genetic Algorithm (SPEA2) [START_REF] Zitzler | SPEA2: Improving the Strength Pareto Evolutionary Algorithm[END_REF] are used as example to characterize the efficiency of the studied recombination operators on test problems of the literature.

Elitist Multi-objective Genetic Algorithms

Since the mid-1990s, there has been a growing interest in solving multi-objective problems by Genetic Algorithms. In particular, elitist MOGAs based on Pareto approaches have become more and more popular because of their capabilities to approximate the set of optimal trade-offs in a single run [START_REF] Corne | The Pareto Envelope-based Selection Algorithm for Multiobjective Optimization[END_REF], [START_REF] Deb | A fast-elitist non-dominated sorting genetic algorithm for multiobjective optimization: NSGA-II[END_REF], [START_REF] Zitzler | SPEA2: Improving the Strength Pareto Evolutionary Algorithm[END_REF]. Elitist MOGAs use an external population, namely archive, which preserves non-dominated individuals in the population. At each generation, individuals (parents) selected from the archive (and/or from the population) following Pareto domination rules are crossed and mutated to create new individuals (children). The population of children and the archive are merged to assess the non-dominated set of the next generation. If the number of non-dominated individuals is higher than the size of the archive, a clustering method is used to preserve most representative solutions and eliminate others in order to keep a constant archive size. Note that niching is used in the selection scheme when individuals involved in a tournament have the same Pareto domination rank.

The second version of the Non-dominated Sorting Genetic Algorithm (NSGA-II) is based on the principles previously exposed. NSGA-II determines all successive fronts in the population (the best front corresponding to the non-dominated set). Moreover, a crowding distance is used to estimate the density of solutions surrounding each individual on a given front. In a tournament, if individuals belong to the same front, the selected one is that with the greater crowding distance. This niching index is also used in the clustering operator to uniformly distribute the individuals on the Pareto front. All details of the algorithm can be found in [START_REF] Deb | A fast-elitist non-dominated sorting genetic algorithm for multiobjective optimization: NSGA-II[END_REF].

The new version of the Strength Pareto Genetic Algorithm (SPEA2) is rather similar to the NSGA-II. It essentially differs in the clustering method used to merge non-dominated individuals in the archive and on the selection method based on a fitness assignment (called strength) related to Pareto ranking of individuals in the population [START_REF] Zitzler | SPEA2: Improving the Strength Pareto Evolutionary Algorithm[END_REF].

In this, study NSGA-II and SPEA2 are taken as reference to investigate the influence of recombination and self-adaptation in real-encoded MOGAs.

Recombination and Self-Adaptive Procedures

We examine various recombination and self-adaptive procedures for real-encoded MOGAs :

The Blend Crossover

From two parent solutions p 1 (i) and p 2 (i), the blend crossover (BLX-α) creates one child c(i) as follows [START_REF] Eshelman | Real-coded genetic algorithms and interval schemata[END_REF] :

)) ( ) ( ( ) ( ) ( 1 2 1 i p i p i p i c - + = β (1)
where β is a random variable in the interval ] 1 , [ α α + -, i denoting the index related to the object variable of the child and parents solutions. If α is set to zero, this crossover creates a random solution inside the range defined by the parents similarly to the arithmetical crossover [START_REF] Michalewicz | Genetic Algorithms + Data Structures = Evolution Programs[END_REF]. Eshelman and Schaffer have reported that BLX-0.5 (with 5 . 0 = α

) performs better than BLX with any other α value in a number of test problems.

The Simulated Binary Crossover

The simulated binary crossover (SBX) simulates the working principle of the single point crossover operator on binary strings [START_REF] Deb | Self-Adaptation in Real-Parameter Genetic Algorithms with Simulated Binary Crossover[END_REF]. From two parent solutions p 1 (i) and p 2 (i), it creates two children c 1 (i) and c 2 (i) as follows :

® + + - = - + + = ] ) ( ) 1 ( ) ( ) 1 ( [ 5 . 0 ) ( )] ( ) 1 ( ) ( ) 1 ( [ 5 . 0 ) ( 2 1 2 2 1 1 i p i p i c i p i p i c β β β β (2)
with a spread factor β defined by (3), ( )

°°°°® 1 • © § - < = + + otherwise ) 1 ( 2 1 5 . 0 if 2 1 1 1 1 η η β u u u ( 3
)
where u is a random variable in the interval [0,1] and η is a nonnegative real number that characterizes the distribution of the children in relation to their parents. A large value of η gives a higher probability for creating children near parents. Acting alone and without any mutation operator, SBX presents interesting properties of selfadaptation similarly to Evolution Strategies [START_REF] Deb | Self-Adaptation in Real-Parameter Genetic Algorithms with Simulated Binary Crossover[END_REF].

The Breeder Genetic Crossover

From two parent solutions p 1 (i) and p 2 (i), the BGX crossover creates one child c(i) as follows [START_REF] Schlierkamp-Voosen | Strategy Adaptation by Competing Subpopulations. Parallel Problem Solving from Nature 3 -PPSN III[END_REF]:

δ i i p i p i p i p i p i c ∆ - - ± = ) ( ) ( )) ( ) ( ( ) ( ) ( 1 2 1 2 1 (4)
where ∆ i is normally set to 0.5 times the domain of definition of the object variable i and the metric denotes the Euclidean distance in the object variable space. δ is computed from a distribution that favors small values:

ku - = 2 δ ( 5
)
where u is a random variable in the interval [0,1], the precision constant k being typically set to 16. Note that in [START_REF] Schlierkamp-Voosen | Strategy Adaptation by Competing Subpopulations. Parallel Problem Solving from Nature 3 -PPSN III[END_REF], the child was placed more often in the direction to the best parent (p 1 (i) being the parent with the better fitness) and the minus sign in (4) was chosen with probability 0.9. In this work, the choice of p 1 (i) and the sign in (4) are made with a probability 0.5. In Fig. 1, the probability density function per child for the previous investigated crossover operators is depicted. The corresponding parents p 1 and p 2 are marked with a full circle. Note these three crossover operators are rather complementary since BGX with k=16 essentially reinforces the accuracy in the neighborhood of the parents whereas SBX-0.5 favors global exploration.

Self-Adaptive Mutation

We also investigate a self-adaptive mutation operator used in standard Evolution Strategies [START_REF] Th | Evolutionary algorithms in Theory and Practice[END_REF], [START_REF] Deb | Self-Adaptation in Real-Parameter Genetic Algorithms with Simulated Binary Crossover[END_REF]. Each individual is characterized by its object variables x(i) and associated standard deviation σ(i). Children (c(i),σ c (i)) are created from their parents (p(i),σ p (i)) using the following rule :

)) 1 , 0 ( ) 1 , 0 ( exp( ) ( ) ( i p c N N i i τ τ σ σ + ′ = (6) ) 1 , 0 ( ) ( ) ( ) ( i c N i i p i c σ + = (7)
where N(0,1) denotes a normally distributed random number with mean 0 and standard deviation 1. N i (0,1) indicates that the random number is generated anew for each value of i. The factors τ and τ ′are commonly set to (2m 1/2 ) -1/2 and (2m) -1/2 where m denotes the number of object variables [START_REF] Th | Evolutionary algorithms in Theory and Practice[END_REF].

Self-Adaptive Recombination

As it is not possible to a priori know which crossover operator will be the most efficient on a specific problem, we propose a self-adaptive scheme similar to that of Spears for binary encoded GAs [START_REF] Spears | Adapting crossover in evolutionary algorithms[END_REF]. It consists in associating in the chromosome of individuals an additional gene (X-gene) that codes the type of crossover to apply during the recombination. When recombining two parents, the operating crossover is randomly chosen from the X-gene of the parents. Using this procedure, the MOGA will favor the crossover that produces the best children through the selection operator.

To avoid premature convergence to a particular type of crossover, the X-gene also undergoes mutation.

Test Problems

We consider three multi-objective problems of the literature [START_REF] Zitzler | Comparison of multiobjective evolutionary algorithms : Empirical results[END_REF], [START_REF] Schaffer | Multiple objective optimization with vector evaluated genetic algorithms[END_REF] displayed in Table 1. ZDT4 is a multimodal continuous problem, which contains 21 9 local Pareto fronts. The global Pareto front is obtained with g=1 and is convex. ZTD6 has a nonuniformly distributed search space with solutions non-uniformly distributed along the Pareto front (the front is biased for solutions for which f 1 (x 1 ) is close to one). The Pareto front is obtained with g=1 and is non-convex. SCH is a generalization of the Schaffer's problem. A large variable space domain and a convex Pareto front characterize it.

¦ = - + = = ≤ ≤ - 1 • © § - = ≤ ≤ =
1 • © § + = = ≤ ≤ - = - - = ¦ = i i i x g x g f g x f x x x f π SCH ) 2 ( 40
¦ ¦ = = - = = ≤ ≤ = i i i i i x x f i x x x f

Performance Criteria

Three performance criteria are used to assess the efficiency of MOGAs :

Average deviation to the theoretical Pareto-optimal front. The average distance ε of the non-dominated set to the Pareto-optimal front is computed as follows [START_REF] Zitzler | Comparison of multiobjective evolutionary algorithms : Empirical results[END_REF],

{ } ¦ ∈ ∈ - = F a F a a a F * * * min 1 ε (8)
where F (respectively F*) denotes the non-dominated set in the final population (respectively the theoretical Pareto-optimal front), a and a* belonging to each subset. The metric in [START_REF] Michalewicz | Genetic Algorithms + Data Structures = Evolution Programs[END_REF] is the Euclidean distance computed in the objective space.

Spread. We define the spread min ε as the average minimum distance of the nondominated set to the Pareto-optimal solutions that minimizes each objective independently.

{ } ¦ = ∈ - = n i i F a a a n 1 * min min min 1 ε (9)
where n is the number of objectives and * min i a represents the theoretical solution of the Pareto-optimal front that minimizes the i-th objective. Spread characterizes the ability of MOGAs to detect boundary solutions of the Pareto-optimal front. Spacing. Spacing ∆ is a measure based on consecutive distances among the solutions of the non-dominated set [START_REF] Deb | A fast-elitist non-dominated sorting genetic algorithm for multiobjective optimization: NSGA-II[END_REF]. It assesses the ability of the MOGAs to distribute its population uniformly along the Pareto-optimal front.

¦

- = - - = ∆ 1 1 1 1 F i i d d F ( 10
)
where d i is the Euclidean distance between two consecutive solutions of the nondominated set, d being the average of these distances. A value of zero for this metric indicates all the non-dominated solutions found are equidistantly spaced. Unlike the definition of ∆ in [START_REF] Deb | A fast-elitist non-dominated sorting genetic algorithm for multiobjective optimization: NSGA-II[END_REF], we do not include in the non-dominated set the boundary solutions of the theoretical Pareto-optimal front to take into account the spread (spread is independently evaluated by ( 9)).

Note that these criteria are complementary to assess the efficiency of MOGAs as indicated by Fig. 2 which illustrates various situations of the non-dominated set in relation to the theoretical Pareto-optimal front.

Tests Results

We successively compare the efficiency of NSGA-II and SPEA2 on the previous test problems using each crossover operators and self-adaptive procedures of sect. 3. All tests are made with the same number of objective function evaluations. NSGA-II and SPEA2 are run for 200 generations with a population size of 100. The archive size is also set to 100 and the crossover probability is 1. Both MOGAs use the BGA mutation operator [START_REF] Schlierkamp-Voosen | Strategy Adaptation by Competing Subpopulations. Parallel Problem Solving from Nature 3 -PPSN III[END_REF] with a mutation rate of 1/m (where m is the number of variables). The self-adaptive recombination scheme employs the BGX with k=16, the SBX with η =1 and the BLX-0.5. The X-gene undergoes mutation with a probability of 5%.

NSGA-II and SPEA2 with self-adaptive mutations operate with initial standard deviations set to 1/10 times the domain of definition of the object variables and without any crossover operator. For all investigated tests, 100 runs are made with random populations to take into account the stochastic nature of the MOGA. An average statistic is taken from the final population for the performance criteria.

Influence of the Recombination Operators

We present in Tables 2-4 the values of the performance criteria on the investigated problems for the NSGA-II and the SPEA2 in relation to each recombination and selfadaptive procedures. From these results, we propose a sort of the recombination and self-adaptive operators on these test problems for each MOGA. Note that this sort was established by giving priority to the average deviation ε .

As also underlined in a earlier study [START_REF] Costa | An Evolution Strategy for Multiobjective Optimization[END_REF], it can been seen that MOGAs without crossover operator perform poorly even if they use self-adaptive mutations. Table 2 shows that NSGA-II and SPEA2 with BLX-0.5 give the best results on ZDT4. BLX-0.5 associated with the SPEA2 works well on ZDT6 but performs poorly when it is coupled with the NSGA-II. This behavior will be also noted for the self-adaptive recombination scheme in the following section. NSGA-II with the BGX crossover works well on ZDT6 but performs extremely poorly on SCH (convergence was not achieved after 200 generations; only one non-dominated individual was found in the final population in all runs). Except with the SBX crossover, NSGA-II and SPEA2 fail to spread correctly their population on SCH. Therefore, SBX is ranked at the top of the sort for this problem despite a slightly lowest quality for ε and ∆.

f 1 f 2 0 min ≈ ε 0 ≈ ε 0 ≈ ∆ f 1 f 2 0 min ≠ ε 0 ≈ ε 0 ≈ ∆ f 1 f 2 0 min ≈ ε 0 ≈ ε 0 ≠ ∆ f 1 f 2 0 min ≠ ε 0 ≈ ε 0 ≠ ∆ f 1 f 2 0 min ≠ ε 0 ≠ ε 0 ≈ ∆ f 1 f 2 0 min ≠ ε 0 ≠ ε 0 ≠ ∆
Finally, it can be seen that best results are always obtained by simple crossover acting alone. However, the efficiency of each crossover operator clearly depends on the characteristics of the test problem. Therefore, using simultaneously multiple crossover operators through a self-adaptive recombination scheme tends to improve the robustness of the MOGAs. We verify this property since the self-adaptive scheme performs extremely well in all cases (the ranking efficiency always equals 2).

Analysis of the Self-Adaptive Recombination Scheme

To understand the mechanism of the self-adaptive recombination scheme, we plot in Fig. 3-5 the origin of children for each investigated test problem as a function of the generation number. Results are extended to 400 generations and averaged on 100 runs. We also indicate the threshold from which the clustering is operating.

It can be seen that the self-adaptive recombination scheme is able to direct the MOGA towards the crossover operator which performs the best on a given test problem. Therefore, the BLX-0.5 which has the best global exploration properties, is rapidly favored on ZDT4 to avoid misleading local Pareto-optimal fronts (see Fig. 3). Because of the large search space, this crossover is also preferred on SCH at the beginning of the search but both MOGAs finally switch towards the SBX crossover to better distribute individuals on boundary points (see Fig. 5 and Table 3). Curiously, the behavior of the self-adaptive recombination scheme on ZDT6 is radically different for the SPEA2 and the NSGA-II. These results are in accordance with those of sect. 5.1 (see Table 3) which show that SPEA2 performs better with SBX and NSGA-II better with BGX. However, both MOGA need to use the BGX crossover to avoid misleading attractors of high f 2 values in the neighbourhood of the boundary point of the Pareto front defined by f 1 * ≈0.28. Note that the steady state behavior of the archive characterized by clustering operations does not necessary correspond to the steady state operations of the selfadaptive recombination scheme. Recombination rates of each crossover operator can evolve continuously during the search as shown in Fig 345.

Finally, we examine in Table 5-7 the influence of the X-gene mutation rate on the efficiency of MOGAs with the studied self-adaptive recombination scheme. Results show that performance criteria are not very sensitive to this factor. However, low mutation rates have to be preferred to exploit benefit of self-adaptation through the selection procedure.

Fig. 1 .

 1 Fig. 1. Probability distribution of children solutions related to each crossover operators

Fig. 2 .

 2 Fig. 2. Illustration of the performance criteria used in this study on various situations. The nondominated set is symbolized with full circles and the theoretical Pareto optimal front is represented by a continuous line. 0 ≈ ε (resp. 0 ≠ ε ) indicates a good (resp. a bad) average deviation, 0 min ≈ ε (resp. 0 min ≠ ε ) a good (resp. a bad) spread and ∆≈0 (resp. ∆≠0) a good (resp. a

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Origin of children in the self-adaptive recombination scheme. Results on ZDT4 (average on 100 runs)

Table 1 .

 1 Test problems used in this study (minimization of both objectives)

	Problem	Characteristics
	ZDT4

Table 2 .

 2 Performance criteria on problem ZDT4. Results are averaged on 100 runs with random initial populations. Best results are indicated in bold types and margin errors with 95% confidence are given in brackets

	MOGA	Operator	Av. dev. ε	Spread min ε	Spacing ∆	Ranking
		BGX (k=16)	2.231 [0.185] 1.769 [0.148] 0.112 [0.011] 4 (poor)
		SBX (η=1)	1.961 [0.170] 1.656 [0.138] 0.033 [0.009] 3 (good)
	NSGA-II	BLX-0.5	1.483 [0.151] 1.275 [0.122] 0.019 [0.006] 1 (excellent)
		Self-Ad. Mut. 10.41 [1.000] 9.439 [0.939] 0.017 [0.001] 5 (very poor)
		Self-Ad. Rec. 1.688 [0.164] 1.437 [0.136] 0.021 [0.008] 2 (very good)
		BGX (k=16)	5.343 [0.436] 4.367 [0.380] 0.301 [0.041] 4 (poor)
		SBX (η=1)	2.233 [0.221] 1.894 [0.163] 0.028 [0.008] 3 (good)
	SPEA2	BLX-0.5	1.527 [0.148] 1.373 [0.122] 0.021 [0.005] 1 (excellent)
		Self-Ad. Mut. 10.55 [1.090] 9.525 [1.015] 0.012 [0.002] 5 (very poor)
		Self-Ad. Rec. 1.616 [0.161] 1.427 [0.133] 0.024 [0.006] 2 (very good)

Table 3 .

 3 Performance criteria on problem ZDT6. Results are averaged on 100 runs with random initial populations. Best results are indicated in bold types and margin errors with 95% confidence are given in brackets

	MOGA	Operator	Av. dev. ε	Spread min ε	Spacing ∆	Ranking
		BGX (k=16) 0.000 [0.000] 0.000 [0.000] 0.006 [0.000] 1 (excellent)
		SBX (η=1)	0.185 [0.018] 0.001 [0.000] 0.202 [0.022] 4 (poor)
	NSGA-II	BLX-0.5	0.081 [0.044] 0.000 [0000] 0.096 [0.054] 3 (good)
		Self-Ad. Mut. 1.866 [0.114] 0.442 [0.077] 0.543 [0.063] 5 (very poor)
		Self-Ad. Rec. 0.014 [0.006] 0.000 [0.000] 0.024 [0.010] 2 (very good)
		BGX (k=16)	0.068 [0.008] 0.000 [0.000] 0.071 [0.006] 3 (good)
		SBX (η=1)	0.186 [0.016] 0.000 [0.000] 0.141 [0.015] 4 (poor)
	SPEA2	BLX-0.5	0.059 [0.005] 0.000 [0.000] 0.055 [0.005] 1 (very good)
		Self-Ad. Mut. 1.662 [0.126] 0.554 [0.090] 0.484 [0.049] 5 (very poor)
		Self-Ad. Rec. 0.061 [0.006] 0.000 [0.000] 0.056 [0.007] 2 (very good)

Table 4 .

 4 Performance criteria on problem SCH. Results are averaged on 100 runs with random initial populations. Best results are indicated in bold types and margin errors with 95% confidence are given in brackets

	MOGA	Operator	Av. dev. ε	Spread min ε	Spacing ∆	Ranking
		BGX (k=16)	no convergence achieved in 200 generations	5 (very poor)
		SBX (η=1) 0.007 [0.000] 0.086 [0.006]	0.006 [0.000] 1 (excellent)
	NSGA-II	BLX-0.5	0.004 [0.000] 0.209 [0.008]	0.004 [0.000] 3 (good)
		Self-Ad. Mut. 0.091 [0.006] 0.195 [0.059]	0.008 [0.001] 4 (poor)
		Self-Ad. Rec. 0.004 [0.000] 0.118 [0.009]	0.005 [0.000] 2 (very good)
		BGX (k=16)	no convergence achieved in 200 generations	5 (very poor)
		SBX (η=1) 0.006 [0.000] 0.100 [0.007]	0.005 [0.000] 1 (very good)
	SPEA2	BLX-0.5	0.009 [0.001] 0.289 [0.010]	0.002 [0.000] 3 (good)
		Self-Ad. Mut. 0.302 [0.306] 0.445 [0.299]	0.008 [0.002] 4 (poor)
		Self-Ad. Rec. 0.006 [0.000] 0.129 [0.009]	0.004 [0.000] 2 (very good)

Conclusion

In this paper, we have shown that the efficiency of real-encoded MOGAs strongly depends on the crossover operators used to explore new solutions of the search space. Therefore, we have proposed a self-adaptive recombination scheme based on three complementary crossover operators to reduce the sensitivity to the crossover proce-dure and improve MOGA robustness. First results are promising and further investigations on difficult constrained problems should confirm the advantage of using this technique in MOGAs.


Luca Marchetti

Daniele Pucci

Pascal Morin

Autonomous Shopping Cart Platform for People with Mobility Impairments

Providing a platform able to interact with a specific user is a challenging problem for assistance technologies. Among the many platforms accomplishing this task, we address the problem of designing an autonomous shopping cart. We assume that the shopping cart is set-up on a unicycle-like robot endowed with two sensors: an RGB-D camera and a planar laser range finder. To combine the information from these two sensors, a data fusion algorithm has been developed using a particle filter, augmented with a k-clustering step to extract person estimations. The problem of stabilizing the robot's position at a fixed distance from the user has been solved through classical control design. Results on a real mobile platform verify the effectiveness of the approach here proposed.

I. INTRODUCTION

Assistive technologies focus their efforts on providing reliable solutions to help people in the everyday life. One of the key components of an assistive system is the ability to actively follow a user, a task well exemplified by a mobile robot that follows the user. An autonomous shopping cart is a simple application that provides a good test-bed for a whole class of problem: a robotic butler that helps on carrying heavy objects; a robotic lift that has to follow a companion to accomplish a coordinated task; an automatic walking aid that should support elderly people and so on.

Platform of such kind should be able to detect and recognize the user, among other people, and be capable of following continuously the same user. The environment should be modelled in such a way the robot can avoid obstacles, and pursue the user at the same time. We focused our attention on developing methodologies to accomplish a safe following of user's trajectory, while maintaining a certain degree of freedom on the reference position of the robot w.r.t. the user.

The involved scientific challenges can be summarised in two aspects. On one hand, a module must be developed in order to estimate the user position while identifying other people in the environment. The selection of the user should be effective in such a way the continuous following of the user will not be confused by the presence of other subjects. This estimation must achieved in a cluttered and noisy environment. On the other hand, one has to provide a reliable and feasible way to control the mobile platform, respecting the peculiarity of human motion and exploiting the robotic aid. Combining these two challenges represents the main contribution of our work. People Position Estimation is a wellknown problem addressed in many scenarios, e.g. videosurveillance [START_REF] Fuentes | People tracking in surveillance applications[END_REF] or activity recognition [START_REF] Bodor | Visionbased human tracking and activity recognition[END_REF]. In the field of object recognition, the human body represents probably the most challenging one. The complexity of shape, as well as the multiplicity of configuration it can achieve require complex sensors to be captured. Color cameras are amongst the most effective sensors, even if the information are limited to the image plane [START_REF] Mckenna | Tracking groups of people[END_REF]. Recently the evolution of technology and availability of relatively cheap RGB-D sensors, capable of perceiving 3D structures, opened the possibility to extend the range image-based recognition [START_REF] Shotton | Real-Time Human Pose Recognition in Parts from Single Depth Images[END_REF]. It is then reasonable to choose such sensors to capture people positions [START_REF] Spinello | People detection in RGB-D data[END_REF].

These devices, however, usually have a limited field-ofview. While it is completely reasonable to use multiple camera to augment the virtual field-of-view [START_REF] Luber | People tracking in rgb-d data with on-line boosted target models[END_REF], other aspects of the application guided us on choosing a different solution to cope with this problem. In the case of a robot following a person, it is important to take care of obstacles that can limit the motion of the robot. To address this problem, usually a laser range finder is employed to map the surrounding [START_REF] Censi | An ICP variant using a point-to-line metric[END_REF]. Its high precision on a 2D plane is an effective way to detect obstacles for the robotic motion. Moreover, the laser information can also be used to detect and track the legs of several people [START_REF] Arras | Efficient people tracking in laser range data using a multi-hypothesis leg-tracker with adaptive occlusion probabilities[END_REF]. Thanks to the large field of view and range, and the sensor resolution, the user position can be estimated with high accuracy.

By exploiting the strengths of both sensors, an improved position estimation can be obtained. The result is an accurate person position estimation that can be given as input to a person following control module. While the human being is able to move along any direction, a wheeled robotic platform is usually subjected to kinematic constraints that limit the range of feasible trajectories [START_REF] Morin | Handbook of Robotics[END_REF]. Hence, to achieve human following by the robot, additional maneuvers may be necessary when the human trajectory does not satisfy the aforementioned constraints. This problem has been largely studied in the last decades and one can rely on existing techniques to address the control problem [START_REF] Morin | Handbook of Robotics[END_REF].

In this paper, we describe how a data fusion algorithm, that combines information from 2D and 3D sensors, can be effectively coupled with a trajectory stabilization method, able to drive the robot to solve the person following challenge.

II. NOTATION

We consider the class of unicycle-like robots sketched in Figure 1. The following notation is used. Let I = {O; ı 0 ,  0 } be a fixed inertial frame with respect to (w.r.t.) which the robot's absolute pose is measured. The point M is the middle point of the wheel's axis, and B = {M ; ı, } is a frame attached to the robot. The vector ı is perpendicular to the wheel's axis. The vector of coordinates of M in the basis of the fixed frame I is denoted as x m = (x m1 , x m2 ) T . Therefore, OM = x m1 ı 0 +x m2  0 . The robot's orientation is characterized by the angle θ m between ı 0 and ı. The rotation matrix of an angle θ m in the plane is R(θ m ). With {e 1 , e 2 } we denote the canonical basis in R 2 . In view of this notation, the kinematic model of the robot writes [START_REF] Morin | Handbook of Robotics[END_REF] ẋm = vR(θ m )e 1 θm = ω.

(

with v the robot's rolling velocity and ω its rotational velocity considered as kinematic control inputs. The position of the user is represented by a reference point p R . The vector of coordinates of p R in the basis of the fixed frame I is denoted as x r = (x r1 , x r2 ) T . Therefore, Op R = x r1 ı 0 + x r2  0 . The vector of coordinates associated with the linear velocity of p R w.r.t. I is denoted as ẋr .

III. PEOPLE POSITION ESTIMATION

People detection is one of the main component needed to have a reliable autonomous shopping cart. In Section I, we highlighted the advantages of using two different sensors to achieve the position estimation. Next, we describe the data fusion architecture pointing out the main characteristics. The functional blocks are described in Figure 2.

The architecture is two-tiered. The lowest level has a person estimation method for each sensors. The highest level combines the two estimations to obtain a more reliable position estimation. a) Laser-based position estimation: To detect people using the laser range finder, we use an implementation of the Kalman Filter leg tracker, described by Arras et al in [START_REF] Arras | Efficient people tracking in laser range data using a multi-hypothesis leg-tracker with adaptive occlusion probabilities[END_REF]. Therefore, we follow the notation presented in that paper to briefly describe the functioning of this estimator.

The KF-based multi-hypotheses tracker describes a leg track as p i L = (p x , p y , v p x , v p y ), with p x and p y position on the laser plane, and v p x and v p y the components of the velocity. The state prediction of the leg filter uses a constantvelocity model. The observations of legs are detected using an Adaboost algorithm [START_REF] Freund | A decision-theoretic generalization of on-line learning and an application to boosting[END_REF] that classifies the segments found in the laser scan according to a set of features. Using these features, several weak classifiers are used to separate leg candidates. The combination of all the classifiers generates the leg observations. The training procedure and how to obtain a valid set of informative threshold to classify the segments is explained in [START_REF] Arras | Using boosted features for detection of people in 2d range scans[END_REF]. The tracker labels assign each new measurement to existing leg tracks or creates new tracks. At any instant, a leg track can be detected (if measurements are assigned during the last observation phase) or deleted (if measurements are not assigned). New tracks are labeled as new track and a false alarm label is used when the measurements are mistakenly detected as track.

People tracks are extracted from leg tracks using the following heuristic:

• people have two legs;

• legs are close to each other; • legs move in similar direction;

• legs have a higher probability of occluding each other, than being occluded by other people's legs or objects. This model is implemented in such a way it takes into account the possible occlusions, thus avoiding deletion of track if legs are occluded for a short period of time. Other considerations about leg track labeling and probability association are described in details in [START_REF] Arras | Efficient people tracking in laser range data using a multi-hypothesis leg-tracker with adaptive occlusion probabilities[END_REF] and will be omitted here.

b) Kinect-based person estimation: The availability of a ROS-integrated library for user skeleton detection simplifies the problem of detecting human shapes using the RGB-D data from the Kinect sensor. As for OpenNI library version 1.3, the user needed to perform a peculiar calibration procedure (the so-called ψ pose) to be detected. This limited the usability of the bundle software for multiple people detection. As for version 1.5, however, the users can be detected using a User Generator that does not require any calibration at all. The output of this module is a set of people position estimation P K and can be extremely noisy, because the people's bodies are recognized applying statistical methods. This required to develop appropriate filtering procedure to establish correct person following. c) Data-fusion for person position estimation: We decided then to use both information, from the laser and the RGB-D camera, to provide more reliable information to the control layer. The camera can be really accurate on estimating the complexity of human body, while the laser provides a larger field of view and a better precision on distance estimation.

To combine the advantages of both sensors, we designed a particle filter with clustering. Each particle represents a possible position estimation as p i = (p x , p y , v p x , v p y ). The posterior density is approximated by:

P(X|z) ≈ |P | i=1 w i δ(X -x i ), (2)
where the X is the current state of the probability density function and w i is a weight associated to the sample x i ∈ X.

The interested reader can find more explanation about this representation in [START_REF] Arulampalam | A tutorial on particle filters for on-line non-linear/non-gaussian bayesian tracking[END_REF]. The algorithm is described below.

p i t ∼ π t (p i |P, v i t-1 ) 2 calculate weight: w i t ∝ L L ( p i t ) × L K ( p i t ) 3 resample: {p i t , w i t } |P | i=1 = resample({ p i t , w i t } |P | i=1 ) 4
get position estimation clusters: C = KClusterise(P )

At time t, the particle filter algorithm requires a proposal distribution (π t ) from which it draws samples during the prediction step. We use the previous set of particles evolved using a constant velocity model. The update step uses information from the laser and the RGB-D sensors. A sensor model calculate the likelihood of each particle to belongs to the set of laser measurements P L or camera measurements P K . The likelihood is evaluated as:

L L = prob(P L |p i t )prob( p i t |p i t-1 ) π( p i t | P , v i t-1 ) , (3)
L K = prob(P K |p i t )prob( p i t |p i t-1 ) π( p i t | P , v i t-1 )
.

The estimated posterior represents the distribution of people over the sensor's space. This posterior is usually multimodal, given the noisy nature of the RGB-D camera and ambiguity on laser estimation. Therefore, a clustering phase is necessary to extract all the possible tracks. A track, or a person estimation, will be the Gaussian approximation (mean and variance) of a single cluster. A selection procedure, not described here, selects the best candidate and assigns it to the control module.

A. K-Clustering

We implemented a k-clustering based technique [START_REF] Jain | Algorithms for Clustering Data[END_REF], described in Algorithm 2. KClusterise tries to detect up to N k clusters. Therefore, it is not a free-cluster algorithm, and this could potentially lead to a limitation. However, for the purpose of the presented applications, this is not a critical problem. First, the algorithm tries to find N k points (with N k ≥ 1) that are equally spaced out (Algorithm 3). Adding particles in line 6 is done using a priority queue principle, considering the distance. At the end of procedure, K will contain the most distant N K points, and they will be used as centroids.

Successively (Algorithm 4), for each point p, it calculates the Euclidean distance δ between p and clusters k ∈ K. Let δ min be the minimum distance between p and a cluster k.

Two cases are possible: if δ min is less than the threshold distance δ f ar (discussed in Section III-B), p will be put in cluster k. Otherwise, p will be put into the outliers set O.

When the clusterisation phase is finished, the points in the outliers set O are evaluated. Each point in O will be put in the nearest cluster by ignoring the threshold distance.

B. Threshold Distance Function

One of the major problems in clustering techniques is to find a threshold distance to approximate the correct number of clusters. There are two possible ways: a fixed value or a variable one. The first choice can be computationally efficient, but it is not flexible w.r.t. environment's changes. We adopted the second strategy, by using a dynamic threshold function. For each position estimation in P L , P K , we evaluate the average distance to other estimations as:

δ f ar = |P L | i=1,j=2 p i L , p j L + |P K | i=1,j=2 p i K , p j K |P L | + |P K | , i = j. (5)
The idea behind this function is to keep the position estimation as far as possible to each other. Using the average of pairwise distance helps us to obtain well-balanced cluster while keeping them separated.

IV. CONTROL DESIGN

Achieving a reliable person following requires correct control laws to smoothly let the robot follow a trajectory constrained by the person's motion. Our objective is to describe the desired position of the person w.r.t. the robot, then minimize the distance between this desired position and the actual position. From Figure 1 the position of the desired (or follower) point F is given by x f = x m + R(θ m )P f m , where P f m is the vector of coordinates of M F = OF -OM . Recall that the position of the person, is given by p R . This is the results of person position estimation described in Section III.

Let x be the position error in the fixed frame, defined by:

x = x f -x r , (6)
and w.r.t. the mobile frame by:

p = R(θ m ) T x. (7)
Note that here x r = (p x , p y ). Therefore the error dynamics w.r.t. the mobile frame writes:

ṗ = -ωS p + M u -R(θ m ) T ẋr (t), (8)
with

M = 1 -P f m 2 0 P f m1 , S = 0 -1 1 0 , u = (v, ω) T
and ẋr = (v p x , v p y ) T . Relying on the results in [START_REF] Morin | Handbook of Robotics[END_REF], one deduces the following lemma.

Lemma 1: Assume that P f m1 = 0, so that det(M ) = 0. Apply the control input

u = -M -1 K p -R(θ m ) T ẋr (t) , (9)
where

K = k 1 0 0 k 2 , K > 0. Then, ṗ = -ωS p -K p, (10)
so that P = 0 is a globally asymptotically stable equilibrium point for the closed-loop system.

The stability analysis follows by verifying that V (p) = p2 is a strict Lyapunov function for the closed-loop system. When P f m1 = 0, the control law ( 9) is not well defined and other solutions to the tracking problem must be considered. This problem will be addressed in a forthcoming paper.

V. RESULTS

The first implementations of the mobile platform have been conducted on an industrial wheeled robot. The robot (cfr. Figure 3) has a Kinect mounted on a pan/tilt unit, which allows a visual servoing procedure to follow a user moving all around the robot. The OpenNI1 library is used to detect and estimate the raw positions of all people in the field-of-view of the camera. A laser is mounted on the front side of the robot and captures scans of the environment. This scans are used to estimate the odometry of the robot (using the Canonical Scan Matcher 2 ). The onboard computer manages the visual servoing routine, controlling the pan/tilt unit. The attached laptop, instead, processes data from Kinect and laser, to model the pose of people and evaluate control inputs.

The software architecture has been developed using the ROS3 framework. The actual architecture, illustrated in Figure 4 is composed of three layers: the robot interface, the modeling core, and the behaviour component.

The robot interface provides the abstraction layer between the actual platform and the software components. It is easy to adapt this interface for several platforms and keep unchanged the higher levels. The software architecture and modules developed for the autonomous shopping cart robot.

The modeling layer detects the user within the sights of the sensors and evaluates the control inputs. The user position estimation is evaluated as described in Section III, while the control inputs are the results of control laws developed in Section IV.

The last layer, the behaviour component, handles the safety measures to ensure the robots avoid close obstacles (detected by the laser), starts and stops services on request, enables initialization procedures and so on. In particular, it selects the first user to be followed, among possible candidates.

A. Experimental Results

In this section, we focus on the practical experiments obtained using the robot. We present three different configuration, considering the behaviour and estimation performed by the robot when the follower point x f was placed in two different positions: A = (1.5, 0) and B = (1, 1). The positions are depicted in Figure V-A. Supplement material and high quality versions of the picture presented here can be found at http://goo.gl/NFnjg. 1) Test A: user in front of the robot: Figure 6 presents the outcome of the data fusion procedure. The data-fusion trajectory represents the trajectory of the person passed to the control module. Despite the high number of hypotheses, the data fusion algorithm is able to consistently track the movement of the person. In Figure 7 is presented the trajectory of the following point w.r.t. the person estimation: the control inputs are consistent with the estimation. The movement of the robot resulted much more smoothed than the estimation. This effect is due to the different frequency at which the control commands are sent to the motor (much lower than the estimation rate).

2) Test B: user at 45 • in front of the robot: In Figure 8 are presented the same results as for test A, showing the effectiveness of the data fusion algorithm. In Figure 9 are presented the trajectory of the following point while at 45 • w.r.t. the robot. The results are consistent with the test A.

VI. CONCLUSIONS

This article presented a reliable solution to the problem of person following. A data fusion algorithm has been presented to reliably detect and estimate the position of multiple people. Two kinds of sensors have been exploited to accomplish the position estimation: a laser range finder and a RGB-D camera. Based on this estimation, a feedback control law has been designed to track the user. To demonstrate the effectiveness of the proposed architecture, results have been presented using a real mobile platform. As future work, we are investigating how to address the singularity on the control design. A different control law needs to be developed to avoid the singularity, in order to have a completely arbitrary position for the following point x f .

Fig. 1 .

 1 Fig. 1. Unicycle-like robotic platform for autonomous shopping cart.

Fig. 2 .

 2 Fig. 2. The data fusion architecture for person position estimation.

Algorithm 1 :

 1 People Position Estimation Data: P L := {position estimation from laser} P K := {position estimation from Kinect} P := {particle set} for p i ∈ P do 1 draw particles:

Algorithm 2 : 2 Assign particle to cluster: 3 O 4 Algorithm 3 : 4 if isFar = true then 5

 2234345 KClusterise Data: P : particle set K: cluster set, of maximum size N K O: outliers set Initialise cluster set: K ← ∅ 1 Find cluster: K = F indCluster(P ) = Classif yP articles(K, P ) Redistribute outliers: SpreadOutliers(O, K) FindCluster // Find equally spaced out centroids for p ∈ P do 1 isFar = true 2 for k ∈ K and isFar do 3 isFar = ( p, k > δ f ar ) Add particle p as centroid: K ← p 6 return K 7

Algorithm 4 : 1 added = false 2 for k ∈ K and not added do 3 δ 6 added = true 7 if added = false then 8 9 return O 10 Algorithm 5 : 5 k candidate = k 6 δ min = δ 7

 41236789105567 ClassifyParticles // Classify particles for p ∈ P do min = p, k 4 if δ min < δ f ar then 5 Add particle p to cluster k: k ← p Add particle p to outliers set: O ← p SpreadOutliers // Redistribute outlier particles for p ∈ O do 1 δ min = ∞ 2 for k ∈ K do 3 δ = p, k 4 if δ < δ min then Classify particle p: k candidate ← p 8

Fig. 3 .

 3 Fig. 3. The robot used to run the experiments.

Fig. 4 .

 4 Fig. 4.The software architecture and modules developed for the autonomous shopping cart robot.

Fig. 5 .

 5 Fig. 5. The position of the follower point x f used during tests.

Fig. 6 .

 6 Fig. 6. Results for the data fusion algorithm during test A: x f = (1.5, 0).

Fig. 7 .

 7 Fig. 7. Trajectory of robot given the control input w.r.t. the position estimation for test A x f = (1.5, 0).

Fig. 8 .

 8 Fig. 8. Results for the data fusion algorithm during test B: x f = (1, 1).

1

  INRIA Sophia Antipolis Meditérranée, 2004, Route des Lucioles -BP 93 Sophia Antipolis Cedex 06902 <name>.<surname>@inria.fr

	2 ISIR-UPMC,	4,	place	Jussieu	75005	Paris

<surname>@isir.upmc.fr

http://www.openni.org

http://ros.org/wiki/csm

http://www.ros.org

ACKNOWLEDGEMENTS

We would thank the Large Scale Initiative Action PAL (Personally Assisted Living), that funded this ongoing research, and the people working at INRIA Sophia Antipolis Meditérranée that help on pursuing this work.