S Farahbakhsh
email: sahar.farahbakhsh@ens-paris-saclay.fr

L Chamoin

M Poncelet
email: martin.poncelet@ens-paris-saclay.fr

Damage tracking using distributed optic fiber sensors in structures

Keywords: Damage detection, SHM, Inverse problems, Optic fiber sensors, OFDR, Model Selection

La détection précoce de l'endommagement et le suivi de sa progression contribuent à améliorer la sécurité et la durabilité des structures. La capacité à modéliser et simuler ce phénomène, en complément de mesures, facilite le contrôle du système. Différents types de capteurs ont été utilisés dans le cadre de la surveillance de l'état des structures. Parmi eux, les capteurs à fibre optique distribués offrent divers avantages qui en font une option pratique pour la surveillance en temps réel des structures. Ce travail consiste à utiliser des capteurs distribués de déformation, basés sur la rétrodiffusion de Rayleigh et la technologie OFDR, dans des échantillons qui supportent une charge jusqu'à leur point de rupture. Les données recueillies sont ensuite post-traitées dans le cadre de l'erreur en relation de comportement modifiée (mCRE) pour effectuer la mise à jour du modèle d'endommagement et estimer l'état du matériau. Plusieurs modèles d'endommagement sont présentés et comparés, et le modèle le plus adapté est aussi sélectionné au moyen de la mCRE.

Introduction

Damage detection, tracking, and modeling are of great importance in today's complex structures. Early detection of damage makes it possible to observe damage progression and prevent critical issues. The ability to model this phenomenon leads to better control of the structure, preventing or delaying its failure and increasing its safety and durability. Distributed Optic Fiber Sensors (DOFS) are a proper and efficient tool for structural health monitoring and damage detection. Recent advances in optic fiber technology now permit temperature and strain sensing with high spatial resolution, leading to real-time and high accuracy measurements inside the material. These sensors are immune to electromagnetic fields and are simple to install on structures or embed in structural material, which brings them a wide variety of use [START_REF] Ding | Distributed optical fiber sensors based on optical frequency domain reflectometry : A review[END_REF]. A common type of DOFS employs Optical Frequency Domain Reflectometry (OFDR) and Rayleigh backscattering to measure strain or temperature alterations along the fiber [START_REF] Bao | Recent advancements in rayleigh scattering-based distributed fiber sensors[END_REF]. In this work, OFDR sensors are implanted in two different material sets to monitor their state and damage advancement through strain measurements. Damage models are presented for both materials. The acquired measurements are then compared to damage model outputs, for model updating purposes, by means of the modified Constitutive Relation Error (mCRE) concept. Since material behavior and enforced loading are generally the sole source of uncertainty in the models, mCRE formulates a cost function based on the constitutive relation while taking the measurement noise into account. Furthermore, mCRE enables determining the consistence of the model and measurements. It is thus a robust tool for parameter identification and remains valid even when measurements are corrupted, due to a strong mechanical regularization [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF]. For each damage model presented in this work, parameters are adjusted by solving an inverse problem that minimizes the mCRE. In addition, the accuracy of each damage model is assessed and discussed with regards to experimental measurements, through the modeling error part of the mCRE functional, so that model selection can be performed. In the next section of this paper, the principle behind optic fiber sensing is briefly explained. The concept of CRE and mCRE is presented in section 3, and section 4 consists of details about the experiments, including the sample preparations. Model updating and parameter identification are shortly discussed in the final section.

Optic fiber sensing

There are various technologies of optic fiber sensing, each with its advantages and drawbacks which makes it convenient for some cases and improper for others. In the case of this study, DOFS based on Rayleigh scattering and OFDR technique have been used for strain sensing. These sensors offer an acceptable sensing range (up to 100 meters), providing real time measurements with a spatial resolution of 2.6 millimeters and a sensing resolution of 1 µϵ.

Rayleigh scattering

Rayleigh scattering corresponds to the interference of light with small inhomogeneities in the medium, that is, inhomogeneities much smaller in size than the wavelength of the traveling light. In optic fibers, these inhomogeneities are the inevitable faults and fluctuations of the refractive index in the fiber core resulting from the manufacturing process. Although considered a disadvantage in optical fiber telecommunications, backscattering is the essential principle behind optic fiber sensors. Due to the Rayleigh scattering, the incident light interferes with a particle. The forward-propagating light loses energy ; thus, attenuation occurs. Since scattered light propagates in all directions, a fraction of scattered light is necessarily backscattered toward the input of the fiber (Figure 1), where the backscattered light is captured for sensing purposes. Rayleigh scattering does not change the frequency of the light, a matter which helps develop the concept of OFDR sensors.

OFDR sensors

The working principle of these sensors is presented schematically in Figure 2. OFDR sensors use interferometry and a frequency-modulated optical wave [START_REF] Bao | Recent advancements in rayleigh scattering-based distributed fiber sensors[END_REF]. A continuous signal with a linearly swept frequency is produced using a tunable laser source (TLS). The incident light is split into two identical beams, one is sent to a reference arm and the other to the fiber under test (FUT). The backscattered light from these two arms is captured and coupled via a coupler to produce a single signal which is further analyzed by an interferometer. The time delay between the backscattered light from the FUT and the reference light produces a unique beat frequency for each scattering point of the FUT. Moreover, the time delay of the backscattered light corresponds to the light flight time to the scattering point and back. Therefore, each beat frequency can reveal the position of a scattering point on the FUT. These beat frequencies can be obtained by applying a Fourier transform (FFT) to the coupled signal and converting it to the frequency domain. Performing this procedure on the fiber in its reference state, i.e., before applying any temperature or strain changes, provides a unique fingerprint for each fiber. The fingerprint signal, or fiber Rayleigh backscattering signature (RBS), remains the same unless the fiber undergoes an environmental change. Any alteration of temperature or strain causes a slight shift in the position of the scattering points, which in turn causes a change in the fiber signature. The same procedure is performed on the fiber while experiencing temperature or strain changes. The signature and measurement signals are transferred to spatial domain by performing another FFT, and are both divided into small windows ∆x, which indicates the spatial resolution of the method. Each divided part of these two signals, that are captured by the same window ∆x in spatial domain, corresponds to the same location on the fiber. For every window ∆x, the two divided signals are again transferred into frequency domain via an inverse Fourier transform. Subsequently, two spectra are resulted from each window, which are then cross-correlated. As a result of this cross correlation, a frequency shift appears in the windows where temperature or strain changes have occurred. This shift of frequency is linearly related to the temperature and strain alterations, and the position of the window containing the frequency shift reveals the measurement location [START_REF] Bao | Recent progress in distributed fiber optic sensors[END_REF]. Consequently, this process leads to mapping the temperature or strain measurements throughout the fiber.

Concept of mCRE

The modified Constitutive Relation Error (mCRE), as is proposed by its name, is developed based on the concept of constitutive relation error (CRE). The CRE functional was primarily introduced for error evaluation in the context of Finite Element method, by presenting an estimation of the discretization error [START_REF] Chamoin | The Constitutive Relation Error method : a general verification tool[END_REF]. The idea of CRE energy-based functional is to provide the model error that originates from the model specified constitutive relation, as it is generally the least reliable equation defining the problem compared to two other sets of equations, namely equilibrium and compatibility equations. In more recent applications, CRE was used as a cost function in inverse problems for parameter identification or model updating [START_REF] Kohn | A variational method for parameter identification[END_REF]. In this context, the observations are strongly imposed in the admissibility fields. The CRE functional is minimized in a two-step iterative manner by first finding an admissible displacement-stress couple and then updating the model parameters accordingly. The CRE functional and the minimization problem are presented for a linear elastic case in equations ( 1) and (2), respectively.

ε 2 CRE (û, ˆ σ; p) = 1 2 ∥ ˆ σ -K(p) ϵ(û) ∥ 2 K -1 (1)
p sol = argmin p∈P min (û, ˆ σ)∈(A + d ) ε 2 CRE (û, ˆ σ; p) (2)
In these formulations, p represents the parameters to be updated, and P is the parameter space, whereas û and ˆ σ show displacement and stress from kinematic and static admissible search spaces respectively. A + d is the ensemble of the admissibility space enriched by the observations, hence the plus sign. ϵ(û) is the strain corresponding to the displacement û, and K(p) stands for the constitutive relation of the problem whose parameters are to be updated. The advantage of using the CRE functional in inverse problems is that the defined cost function has a physical and mechanical sense. However, enforcing the measurement data directly into the admissible fields as additional constraints can decrease the accuracy of the optimal solution as the measurements always contain a level of noise, and might sometimes be corrupted. On the other hand, the modified version of the CRE functional, or the mCRE, solves this problem by separating the information into two sets. Reliable information, such as equilibrium equations, known boundary conditions, and sensor locations, is strongly imposed in the admissibility space. On the contrary, less reliable information, notably constitutive relations, noisy measurements, and unknown boundary conditions, is relaxed [START_REF] Ladevèze | Application of a posteriori error estimation for structural model updating[END_REF].

Linear mCRE

Considering the same linear elastic case when the observations are displacements in a subdomain of the structure, the mCRE functional and the following optimization problem are defined by :

ε 2 mCRE (û, ˆ σ; p) = ε 2 CRE (û, ˆ σ; p) + α 2 (d(û) -d obs ) T G -1 obs (d(û) -d obs ) (3)
p sol = argmin p∈P min (û, ˆ σ)∈(A - d ) ε 2 mCRE (û, ˆ σ; p) (4)
Here, ε CRE is the same functional defined in (1), and α is a scalar factor weighing the measurement part of the functional against the model error part, i.e. the term ε 2 CRE . d obs represents the observations and d(û) is the displacement û on the domain corresponding the measurements. G obs is a normalized scaling matrix usually diagonal, and A - d is the ensemble of the relaxed admissibility space, hence the negative sign. Similar to the minimization of ε 2 CRE , the optimization process of ε 2 mCRE is carried out in a two-step iterative manner. The first step of each iteration consists of calculating the optimal admissible fields (û (n) , ˆ σ (n) ) through a minimization process. In the second step the functional is minimized regarding the parameters p (n) while using the values obtained in the first step for the admissible fields (û

(n) , ˆ σ (n) ).
The superscript (n) indicates the iteration number. Equation [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF] shows that the given inverse problem can be observed from two different angles. The presented explanation is achieved from a mechanical point of view, where ε 2 mCRE is interpreted as a modified version of ε 2 CRE with more relaxed constraints to avoid strongly imposing the measurement error into the formulation. Furthermore, from an optimization perspective the ε 2 mCRE can be perceived as a least squares problem

[(d(û) -d obs ) T G -1 obs (d(û) -d obs )
] with a mechanical reqularization term ε 2 CRE and regularization parameter 1/α, improving the cost function convexity. In either approach the ε 2 mCRE formulation is in fact a compromise between model and measurement data, and the value of α demonstrates the confidence level in the information. When α → 0, the inverse problem is almost equivalent to the minimization of the ε 2 CRE with no regard to the measurement information, indicating high measurement uncertainty. Increasing α implies an increase in confidence towards the measurements. In the extreme case when α → ∞, the inverse problem turns into a classical least squares optimization, with high confidence in measurement data and low confidence in model definition. As α plays an important role in the quality of the optimization procedure, its value is set based on a priori knowledge of the problem, using Morozov principle in this case [START_REF] Morozov | The error principle in the solution of operational equations by the regularization method[END_REF].

Nonlinear mCRE

Extending this concept to a nonlinear framework with elasto-plastic behavior yields a more complex formulation of mCRE. This generalized case is presented through equations ( 5) to [START_REF] Morozov | The error principle in the solution of operational equations by the regularization method[END_REF], for a time interval [0, T ] and a structure volume Ω.

ε 2 CRE (ê e , êp , ŝ; p) = T 0 Ω η ψ (ê e , ŝ; p) + t 0 η φ ( ėp , ŝ; p)) (5)
η ψ (ê e , ŝ) = ψ(ê e ) + ψ * (ŝ) -ŝ.ê e ≥ 0 ; η φ ( ėp , ŝ) = φ( ėp ) + φ * (ŝ) -ŝ. ėp ≥ 0 (6)
ε 2 mCRE (ê e , êp , ŝ; p) = ε 2 CRE (ê e , êp , ŝ; p) + α 2 . 1 N t Nt nt=1 (d(û nt ) -d nt obs ) T G -1 obs (d(û nt ) -d nt obs ) (7)
p sol = argmin p∈P min (êe,êp,ŝ)∈(A - d ) ε 2 mCRE (ê e , êp , ŝ; p) (8)
In these equations, êe and êp are global flux variables regarding elastic and plastic strains. ŝ stands for global thermodynamic forces with respect to ( ˆ σ). ψ and φ denote free energy and dissipation potential, while the dual potentials are shown by ψ * and φ * . η ψ and η φ respectively represent residuals on state equations and evolution laws, a.b corresponds to duality product of a and b, and N t stands for number of time points. The optimization algorithm is similar to the linear elastic case, as shown in equation ( 8). Detailed discussion on the nonlinear mCRE can be found in references [START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF] and [START_REF] Marchand | Parameter identification and model updating in the context of nonlinear mechanical behaviors using a unified formulation of the modified constitutive relation error concept[END_REF].

Experiments

In a first attempt, a cement-based fiber-reinforced mortar is used for beam samples to ensure nearly homogenous, isotropic and nonfragile behavior. A single optic fiber sensor passes through the beam sample several times to provide measurements in different points of the beam section. A second type of sample is also prepared consisting in a fiber reinforced polymer plate. An optic fiber sensor is embedded inside the plate to monitor its damage during the tests. Both sets of samples undergo a quasi-static 4-point bending test (Figure 3). This arrangement provides an uniform pure-bending state in the beam sector between loading pins A and B. A few cycles of quasistatic loading and unloading are applied within the elastic limit of the material. In the final cycle, loading is increased until the rupture of the sample.

Concrete samples

The cement-based mortar is reinforced with polyamide and steel fibers to increase the ductility as well as the tensile strength of the concrete. The state of this material is isotropic and homogenous along the beam axis (z). However, as these samples undergo a 4-point bending test, the behavior of material along the height of the beam (y-axis) starts as homogenous, elastic and linear ; then turns heterogenous with the appearance and progression of cracks in the concrete.

The optic fiber sensor is crossed through both parts of the beam that are under tension and compression. Dimensions of the sample and placement of the fiber in the beam section are presented in Figure 3.

Composite samples

The second material set is a composite with carbon fiber reinforcements pre-impregnated by epoxy. This composite plate consists of 8 different layers with a quasi-isotropic stacking, leading to a homogenous characteristic at macroscale. However, it is considered heterogenous and anisotropic at microscale. The sensor is passed through the composite mid-layer. The chosen orientation of the layers and the placement of the sensor are shown in Figure 4. Since the fiber is positioned at the mid-layer of the composite, in a bending test it falls on the plate neutral axis where no strain is expected. To resolve this issue, the composite is as aforementioned glued on the bottom of a thicker plywood plate, and the assemblage is installed on the bending test machine. Therefore, the composite plate bears the maximum tension resulted from bending of the assemblage (Figure 4).

Model updating

Real-time data from the tests are exploited in the mCRE formulation for parameter updating in the introduced models. Measurements from each instant of the loading process can be used seperately in the formulation to provide information on the state of the sample and the progression of damage, if applicable, at that specific instant. The data obtained from loading and unloading cycles within the elastic limit of the material are used to identify and update undamaged model parameters, e.g. the correct initial Young modulus of the material. The resulting updated parameters are used in the definition of damage models for each material, and immediately after first signs of damage, those parameters characterizing the state of damage and its progression, are picked to be updated. In a first approach for each material, a simple isotropic damage model is presented as in [START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF] where the state of damage is declared by a local decrease in the Young modulus.

E d = E(1 -d) (9)
Here, E d is the reduced Young modulus, indicating damage has occured. d is the damage parameter to be updated and E is the initial Young modulus of the material, before it is damaged. The value of E is determined by parameter updating during the elastic phase of the loading cycles, as previously explained. Moreover, a more sophisticated and orthotropic damage model is presented where three damage parameters are used to better estimate the damage state of the structure. As for the progression of damage, several models based on different evolution equations are implemented in addition to state equations. Identification and updating the parameters in these models lead to the prediction of the damage state, or in other words, defining the damage progression. The attained functionals and determined damage localizations are compared to perform model selection.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 101002857).

Figure 1 -

 1 Figure 1 -Rayleigh scattering in the fiber core

Figure 2 -

 2 Figure 2 -(a) Working concept of OFDR sensors ; (b) Dividing the signal to ∆x windows, and crosscorrelation (from [1])

Figure 3 -Figure 4 -

 34 Figure 3 -(a) and (b) Dimensions of concrete sample and positioning of embedded fiber ; (c) Schematic set-up of a 4-point bending test


Rogelio Lozano

Sergio Salazar

Donovan Flores

Iván González-Hernández

PVTOL global stabilisation using a nested saturation control

This paper presents how the nested saturation control introduced in Teel (1992) can be modified to stabilise the model of the Planar Vertical Take-Off and Landing (PVTOL) system. A particular choice of the amplitudes of the nested saturations controller is given such that global stability of the closed loop system is obtained. The proof is given in detail. Numerical simulations illustrate the performance of the proposed control algorithm.

Introduction

In the last decades, the PVTOL (Planar Vertical Take-Off and Landing) problem has attracted many researchers. The interestismotivatedby thefactthatmanydifferentconfigurations of aerial systems are currently used in numerous applications. The PVTOL is a challenging nonlinear system mainly because it is an under-actuated mechanical system, i.e. the number of controlinputs islessthanthenumber ofdegrees-of-freedom. Therefore, finding an appropriate feedback control strategy will help achieving a high level of performance.

One of the first c ontrollers p roposed f or t he P VTOL was introduced by [START_REF] Hauser | Nonlinear control design for slightly non-minimum phase systems: application to V/STOL aircraft[END_REF]. They showed that applying an exact input-output linearisation technique to the flight control of the VTOL leads to a system with unstable internal dynamics. They proposed a solution with an approximate input-output linearisation procedure developed for slightly non-minimum phase nonlinear systems.

The beginning of the 90s was also marked by the interest on stabilising linear systems with bounded inputs. The nested saturation control technique introduced by [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF]was proposed to stabilise n integratorsincascadewith aboundedinput. However, it should be pointed out that it is not straightforward to apply such nested saturation approach to the PVTOL since it is a nonlinear system. [START_REF] Lozano | Global stabilization of the PVTOL: real-time application to a mini-aircraft[END_REF])andGarciaetal.(2006) proposed to use the nested saturation control to stabilise the PVTOL. The stability was ensured only when the orientation angle theta and its derivative were initially sufficiently sm all. Th e id ea to use the nested saturation technique surged from the fact that when the altitude is controlled by nonlinear compensation, then the resulting subsystem from the torque input τ to the output x (the horizontal displacement) reduces to four integrators in cascade when θ is small, i.e. when tan θ ≈ θ.

Theclosedloopsystemobtainedwhenusingthecontroller proposed in [START_REF] Lozano | Global stabilization of the PVTOL: real-time application to a mini-aircraft[END_REF] was only locally stable. Extending the proof to establish global stability of the closed loop has beenhamperedby thepossibilityofcrossingthesingularpoint cos θ = 0. The search for a global stability proof is motivated by the fact the controller based on nested saturations is relatively simple and is currently used in practice with good performance. [START_REF] Nicotra | Nested saturation control of an uav carrying a suspended load[END_REF] have also applied the nested saturation control to an UAV for the transportation of suspended loads. They proved stability of the equilibrium for small sway angles of the suspended load.

An almost-globally stable controller was proposed by [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF] using a Lyapunov approach. The results showed robustness with respect to aerodynamic drag disturbances.

A hierarchical controller was proposed for the PVTOL in [START_REF] Abdessameud | Global trajectory tracking control of VTOL-UAVs without linear velocity measurements[END_REF] for trajectory tracking in the absence of linear velocity measurements. [START_REF] Aguilar | Stabilization of the PVTOL aircraft based on a sliding mode and a saturation function[END_REF] presented a controller based on a combination of a PD controller and a sliding mode controller to stabilise both the horizontal and angular variables to the desired rest position. The closed loop system was proved to be locally stable provided that the initial angular rate and angular position belong to a given region where no singularity occurs. [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF] proposed a hierarchical control strategy for a miniature VTOL vehicle to track a desired trajectory globally with respect to the initial position and attitude. Their control strategy overcome the topological constraint on the orientation by using this orientation angle as a virtual input to stabilise the aircraft position. [START_REF] Tran | Adaptive control of a quadrotor aerial vehicle with input constraints and uncertain parameters[END_REF] proposed also a control algorithm to obtain a globally stable controller for UAVs. Furthermore, their control strategy is adaptive and can deal also with unknown system parameters. [START_REF] Zavala-Río | Global stabilization of a PVTOL aircraft model with bounded inputs[END_REF] propose a global stabilising control design for PVTOL aircraft, with bounded inputs. This approach is based on the use of non-linear combinations of linear saturation functions bounding the thrust input and the rolling moment to arbitrary saturation limits. Their proof of stability has been extended to the case of small lateral force coupling in [START_REF] López-Araujo | Global stabilisation of the PVTOL aircraft with lateral force coupling and bounded inputs[END_REF]. A global stabilising control law is proposed for a non-minimum phase under-actuated PVTOL aircraft by [START_REF] Poulin | Minimum energy oriented global stabilizing control of the PVTOL aircraft[END_REF]. Their approach is based on receding horizon control and takes into account the positiveness of the thrust and its saturation to any arbitrary level.

Both control schemes [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF][START_REF] Tran | Adaptive control of a quadrotor aerial vehicle with input constraints and uncertain parameters[END_REF]) used the orientation angle as a virtual control input in such a way that the PVTOL model appears as fully actuated. Such virtual control is based on smooth saturation functions of the position and velocity. The first and second derivatives of the virtual control are required so that the actual orientation angle converges to the required virtual control. These derivatives are needed even if the objective is to stabilise the system at hover. A drawback of this technique is that the expressions for the first and second derivatives of the orientation virtual control input are very complex. The approach based on virtual control forces requires the thrust input to be strictly positive, since for u = 0 the virtual control input for θ is undefined. To avoid this singularity, the virtual force for the altitude is chosen as F = g -σ (z, ż) > 0 where g is the gravitational acceleration and σ is a smooth saturation function. Therefore the thrust control input u will satisfy 0 < u < 2g (neglecting the virtual control input for x). The upperbound on the trust input u may not be a problem in many situations, but this constraint on u could limit the agility of those aerial vehicle which require to perform aggressive manoeuvres. These control schemes require measurement of the acceleration signal which is often noisy.

The present paper proposes a controller which is simpler compared to those based on virtual control proposed by [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF]; [START_REF] Tran | Adaptive control of a quadrotor aerial vehicle with input constraints and uncertain parameters[END_REF]. The control strategy is presented in a multistage constructive procedure as follows. The torque control input τ is chosen as a nested saturation controller which, after a transient time, ensures bounds for the orientation rate θ and angular position θ so that the orientation θ is bounded away from the singular point of 90 • . We exploit the fact that provided the control inputs u and τ are bounded, the PVTOL system cannot have a finite escape time. Boundedness of ẋ and x is established later after the altitude converges to its desired value. The control input u is based on the compensation of the nonlinearity cos θ and therefore can only be computed after θ is bounded away from the singular point θ = 90 • . The bounds of the saturation functions have been modified from those proposed by [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF] in order to guarantee that after an explicit finite time the orientation angle θ is constrained to evolve in an interval which does not contain the singular point θ = 90 • . Consequently, once θ is inside the region which is free of singularities, then the altitude can be stabilised by compensating the nonlinear term cos θ. We then resume the analysis of the nested saturation control of the torque input τ to study the boundedness and convergence of ẋ and x. It turns out that provided some conditions on the bounds of the saturation functions we can prove that the state of the closed loop system converges to the origin. The proposed controller does not require computing the first and second derivatives of a virtual force. The controller requires only the measurements of the system state and does not require measurement of linear acceleration nor angular acceleration. The presented controller can be interpreted as an extension of nested saturation controller for n integrators in cascade including nonlinearities as the tangent function. [START_REF] Aguilar | Robust trajectory-tracking control of a PVTOL under crosswind[END_REF] solve the trajectory-tracking control problem of the PVTOL aircraft under crosswind. The PVTOL is linearised by differentiating twice the system equations which is called dynamic extension. This method requires the obtained control input to be different from zero, which is often the case in practical applications. This method is based on disturbance observers designed using sliding mode techniques. [START_REF] Escobar | PVTOL control using feedback linearisation with dynamic extension[END_REF] focused on finding t he conditions for local asymptotic stability when using a control based on feedback linearisation of the PVTOL based on dynamic extension. The domain of attraction is obtained by using a Lyapunov approach. [START_REF] Hernández | Priority altitude PVTOL aircraft control via immersion and invariance[END_REF] use the immersion and invariance control technique to design a controller that globally stabilises the PVTOL aircraft system. The controller gives priority to the aircraft's altitude before controlling the lateral displacement resulting in less thrust effort and less roll moment. [START_REF] Aguilar | A simple approach to regulate a PVTOL system using matching conditions[END_REF] present an energy shaped approach to stabilise the PVTOL aircraft, in conjunction with an exact differentiator o bserver, t o e stimate t he n on-available velocities. Asymptotic convergence is proved by applying LaSalle's theorem, assuming that the system is initialised in the interval

(-π/2, π/2).
This paper is organised as follows : Section 2 presents the PVTOL model. Section 3 presents the nested saturation control algorithm for the torque control input and shows that the angular rate θ ˙ has an upperbound. Section 4 gives the upperbound on the orientation angle θ.Section5 presents the thrust control algorithm. Section 6 shows anupperboundforthe velocity x ˙.Section7 finally shows convergence of the position x to zero. The proposed algorithm is tested in numerical simulations in Section 8. Final remarks are given in the conclusion.

PVTOL model

ThemodelofthePVTOLisgivenbyLozanoetal.( 2004):

mẍ = u sin θ (1) mz = u cos θ -mg (2) θ = τ ( 3
)
where m is the mass, θ is the angle of the aircraft with respect to the horizontal line, g is the gravitational acceleration, x is the horizontal displacement, z is the vertical displacement. u and τ are the total thrust and torque respectively.

Upperbound on the orientation rate θ ẋ1

= x g x 2 = ẋ g x 3 = θ x 4 = θ (4)
Let us define the following variables:

ν 1 = x 3 + x 4 ν 2 = x 2 + 2x 3 + x 4 ν 3 = x 1 + 3x 2 + 3x 3 + x 4
(5)

In this section, we will present the nested saturation controller appliedtothex-θ subsystem. We will show that after some finite time the upperbound on θ ˙ is smaller than a given value.

The state of the x -θ subsystem is given as

consider the subsystem ẋ4 = τ (6) ν1 = x 4 + τ (7)
and the torque input given by the following nested saturation control law:

τ = -σ a (x 4 + σ b (ν 1 + σ c (ν 2 + σ d (ν 3 )))) (8)
where σ a is the saturation function:

σ a = ⎧ ⎪ ⎨ ⎪ ⎩ a x > a x -a ≤ x ≤ a -a x < -a (9)
Introducing ( 8) into ( 6)

ẋ4 = -σ a (x 4 + σ b (w 1 )) ( 10
)
with

w 1 = ν 1 + σ c (ν 2 + σ d (ν 3 )) (11) Assume that a ≥ 2b + 1 (12)
for some 1 > 0.

Appendix 1 gives the expression for a finite time t 1 such that for some 1 > 0

| θ| = |x 4 | ≤ b + 1 ∀ t ≥ t 1 (13)

Upperbound on the orientation angle θ

In this section, we will show that after some finite time θ will belong to an interval which does not include the singular point θ = 90 • . After t = t 1 and in view of ( 7), ( 8), ( 12) and ( 13), ( 7) reduces to

ν1 = -σ b (ν 1 + σ c (w 2 )) (14)
with

w 2 = ν 2 + σ d (ν 3 ) (15) Let b ≥ 2c + 2 (16)
for c > 0 and some 2 > 0. Appendix 2 gives the expression for a finite time t 2 such that

|θ + θ| = |ν 1 | ≤ c + 2 ∀ t ≥ t 2 (17) Let us choose c such that c + 2 < π n + 1 (18)
for some integer n = 3, 4, . . . to be determined. From ( 4) and ( 5) it follows:

θ = -θ + ν 1 (19)
Appendix 3 gives the expressions for a finite time t 3 such that

|θ| ≤ π n ∀ t ≥ t 3 (20) π π
The integer n = 3, 4, 5, ... will be chosen later so that (tan θθ) is small enough to guarantee the stability of the overall system.

Thrust control input u

In this section, we present the thrust control input u which is based on nonlinear compensation. Notice that after time t 3 the value of θ in ( 20) is far away from the singular points 2 and -2 .Wecanthenproposethe following thrust control input u:

u = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ [-2ż -(z -z d )]m + mg cos θ |θ | < π n mg π n ≤ |θ| ≤ π 3 0 |θ| > π 3 (21)
where z d is the constant desired altitude. The constant values for u above when |θ| > π n do not play a role in the convergence analysis. However, they can be used to improve the performance of the algorithm in the presence of disturbances as will be shown in the numerical simulations.

Introducing ( 58) into (55) gives

z + 2ż + (z -z d ) = 0 (22) or ż z = 0 1 -1 -2 z -z d ż (23)

Upperbound on the velocity xİ

n this section, it will be shown that the velocity x ˙ will converge to a bounded interval. Introducing ( 58) into ( 54) and ( 4) gives

ẍ1 = ẍ g = tan θ + e 1 ( 24
)
where

e 1 = tan θ [-2ż -(z -z d )] g ( 25
)
which can be rewritten as

ẍ1 = θ + e 1 + e 2 ( 26
)
where

e 2 = tan θ -θ (27)
Thus, from ( 54) and ( 57), the subsystem for x, ẋ, θ and θ is given by (see ( 4))

ẋ1 = x 2 ẋ2 = x 3 + e ẋ3 = x 4 ẋ4 = τ (28)
with e = e 1 + e 2 [START_REF] Hernández | Priority altitude PVTOL aircraft control via immersion and invariance[END_REF] From ( 5)

ν2 = ẋ2 + 2ẋ 3 + ẋ4 = x 3 + e + 2x 4 + τ = x 4 + ν 1 + e + τ ( 30
)
Taking into account ( 13) and ( 17), ( 8) reduces to

τ = -x 4 -ν 1 -σ c (ν 2 + σ d (ν 3 ))
Introducing the above into (30) leads to

ν2 = -σ c (ν 2 + σ d (ν 3 )) + e (31)
Choose d so that for some 3 > 0

c ≥ 3d + 3 ( 32
)
Appendix 4 gives the expression for a finite time t 4 such (see ( 25))

|e 1 | ≤ d 6 ∀ t ≥ t 4 (33)
The integer n ≥ 3 will be chosen such that for some k > 0 the following holds (see ( 20) and ( 27))

|e 2 | = | tan θ -θ| ≤ k|θ| ≤ k π n ≤ d 6 (34)
Introducing ( 33) and ( 34) into ( 29)

|e| ≤ |e 1 | + |e 2 | ≤ d 3 (35)
Using ( 31), (32) and the above, the expression for a finite time t 5 is given in Appendix 5 such that

|ν 2 | ≤ 2d + 3 ∀ t ≥ t 5 (36) Therefore (31) reduces to ν2 = -ν 2 -σ d (ν 3 ) + e (37)

Convergence of the position x

In this section, we will prove that the horizontal displacement x converges to zero. From (5) and using ( 28) and ( 8)

ν3 = ẋ1 + 3ẋ 2 + 3ẋ 3 + ẋ4 = x 2 + 3x 3 + 3e + 3x 4 + τ = x 4 + ν 1 + ν 2 + 3e + τ = -σ d (ν 3 ) + 3e (38)
From ( 35) and (38) a finite time t 6 is given in Appendix 6 such that

|ν 3 | ≤ d ∀ t ≥ t 6
From (4), ( 5) and ( 28) we get 14) and ( 15)) ( 40)

ẋ3 = -x 3 + ν 1 (39) ν1 = -ν 1 -ν 2 -ν 3 (see (
ν2 = -ν 2 -ν 3 + e (see (37)) (41) ν3 = -ν 3 + 3e (see (38)) (42)
Define

X T = x 3 ν 1 ν 2 ν 3 (43) b T = 0 0 1 3 (44)
and

A = ⎡ ⎢ ⎢ ⎣ -1 1 0 0 0 -1 -1 -1 0 0 -1 -1 0 0 0 -1 ⎤ ⎥ ⎥ ⎦ (45)
Therefore ( 39) through ( 42) can be rewritten as

Ẋ = AX + be (46)
The state space representation above has a system matrix A that is stable with 4 eigenvalues located at -1. Recall that e = e 1 + e 2 (see [START_REF] Hernández | Priority altitude PVTOL aircraft control via immersion and invariance[END_REF] where e 1 converges to zero, (see ( 22) and ( 25)). Furthermore |e 2 | ≤ k|θ | = k|x 3 | for a value of 'k' that has to be determined so that it satisfies (34) and is simultaneously small enough such that (46) remains stable.

The values for k and n such that ( 46) is a stable system could be obtained by using a Lyapunov function X T PX where P > 0 is the solution of the Lyapunov equation A T P + PA = -Q for some Q > 0. Nevertheless, for simplicity we will proceed as follows.

From ( 20) and ( 22) it follows that e 1 in (25) converges to zero. Since A in ( 46) is a stable matrix, the contribution of e 1 in X will also converge to zero. From (34) |e 2 | ≤ k|θ |. Thus, following the procedure developed in Appendix 3, we can prove using (42) that there exists a finite time t 7 such that

|ν 3 | ≤ 3k|θ | + ∀ t ≥ t 7 (47)
for some > 0.

Similarly, using (41), there exists a finite time t 8 such that

|ν 2 | ≤ 4k|θ | + 2 ∀ t ≥ t 8 (48)
Similarly, from (40) it follows that there exists t 9 such that

|ν 1 | ≤ 7k|θ | + 3 ∀ t ≥ t 9 (49)
Let us assume that k also satisfies

μ 2 = 1 -7k > 0 (50)
Finally from ( 4) and ( 39)

θ = -θ + ν 1 (51) Let V = 1 2 θ 2 then, using (49) V = θ θ = θ(-θ + ν 1 ) ≤ -θ 2 + 7kθ 2 + 3|θ| ≤ -μ 2 θ 2 + 2 μ|θ| 2 3 μ ≤ - 3 4 μ 2 θ 2 + 9 2 μ 2 ≤ - 3 2 μ 2 V + 9 2 μ 2 (52)
Since can be chosen arbitrarily small it follows from the above that V → 0 and therefore θ → 0. From ( 51) and ( 4), we conclude that x 3 → 0 and x 4 → 0. From ( 47)-( 49) it follows that ν 1 → 0, ν 2 → 0 and ν 3 → 0. From ( 5), x 2 → 0 and x 1 → 0.

From (4) it finally follows x → 0, ẋ → 0, θ → 0, θ → 0.

The constraints required for n and k are given in ( 34) and ( 50). Appendix 6 proves that for n = 9 and k = 0.0424 both constraints are verified. Notice that π 9 = 20 • . Finally let us summarise the constraints on the saturation upperbounds a, b, c and d. See ( 12), ( 16), ( 18), ( 32) and ( 34)

a ≥ 2b + 1 b ≥ 2c + 2 0.1π ≥ c + 2 ≥ 3d + 2 + 3 ≥ 0.085π + 2 + 3 (53)
Thechoice oftheparametera dependsonthe upperboundof thetorque controlinput τ . In practice the higher the value of a thefastertheangularrateθ ˙ will decrease after a disturbance and thereforetheaerialvehicle attitudewillbemorerobusttodisturbances. Nevertheless, when the parameter a is chosen too high, it is common that mechanical vibrations appear. The parameters b, c and d are determined by the choice of the parameter a. Concerningtheparameters i , they have been introduced so thatthecorrespondingstatevariablesaresmallerthanagiven bound in finite time. The higher the value of i the shorter the corresponding time period t i will be.

Numerical simulation results

In this section, we will present a comparison between the proposed control algorithm and the control strategy introduced in [START_REF] Lozano | Global stabilization of the PVTOL: real-time application to a mini-aircraft[END_REF]. We will also show the robustness properties of the proposed algorithm with respect to plant parameter changes. We will consider that the PVTOL is represented by the following equations:

mẍ = u sin θ (54) mz = u cos θ -mg (55) θ = θ 1 + w p (56) θ1 = kτ (57)
where w p is a perturbation given as follows: We have introduced a disturbance at time 15 s in the orientation rate θ ˙. As can be seen in Figure 2 the disturbance produces a much smaller displacement x in the case of the proposedalgorithmandthedisplacement x takesmuchlesstimeto come back to the origin. Figure 2(b) shows that the thrust input issmallerfortheproposedalgorithm. Theangulardisplacement in Figure 2(c) is just a little smaller for our algorithm. Figure 2(d) shows that the price to be paid for obtaining the above advantages with the proposed algorithm is to reduce the altitude of the PVTOL during a few seconds.

w p = 2.3 15s ≤ t ≤ 16s 0 otherwise (58)
Concerning the robustness of the proposed algorithm with respect to changes in the plant parameter, we have considered two different plant parameters: the mass m and the torque control gain k.

We have used the same initial conditions as before as well as the same disturbance in the angular orientation rate θ ˙. Figure 3 shows the performance of the proposed controller when the plant mass m increases or decreases 20%. As can be seen a change in the plant m producesanerrorinbetweenthealtitude and the desired altitude. The rest of the closed-loop system behaves with no significant changes.

Figure 4 shows the performance of the proposed controller when the torque gain k increases or decreases 20%. As can be seen when k increases the closed loop system behaves well. However, when k decreases the system requires a larger control input and more time to converge.

Conclusion

Thispaperaddressedtheproblemofobtainingasimplecontrol algorithm for the PVTOL. The nonlinearity of the PVTOL has been an obstacle for obtaining stable controllers. Since the PVTOL can be approximated by four integrators in cascade, the nested saturations controller proposed in [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF] has been currently used to control the PVTOL. However, no proof of stability of the closed loop has been obtained in the past. This paper has shown that a particular choice of the amplitudes of the saturations involved in the technique in [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF] allows for the proof of global stability of the closed loop. The performance of the proposed controller is shown in numerical simulations.

Figure 1 .

 1 Figure 1. PVTOL configuration.

Figure 2

 2 Figure 2 presents the comparison between the control algorithm proposedinthepresentpaperandthealgorithminLozano et al. (2004). Thisfigurepresents(a) t hedisplacementx, (b) the control input u, (c) the orientation angle theta and (d) the altitude z. The initial conditions have been chosen as z(0) = 0, θ(0) = 0, x(0) = 0.We have introduced a disturbance at time 15 s in the orientation rate θ ˙. As can be seen in Figure2the disturbance produces a much smaller displacement x in the case of the proposedalgorithmandthedisplacement x takesmuchlesstimeto come back to the origin. Figure2(b)shows that the thrust input issmallerfortheproposedalgorithm. Theangulardisplacement in Figure2(c) is just a little smaller for our algorithm. Figure2(d)shows that the price to be paid for obtaining the above advantages with the proposed algorithm is to reduce the altitude of the PVTOL during a few seconds.Concerning the robustness of the proposed algorithm with respect to changes in the plant parameter, we have considered two different plant parameters: the mass m and the torque control gain k.We have used the same initial conditions as before as well as the same disturbance in the angular orientation rate θ ˙. Figure3shows the performance of the proposed controller when the plant mass m increases or decreases 20%. As can be seen a change in the plant m producesanerrorinbetweenthealtitude and the desired altitude. The rest of the closed-loop system behaves with no significant changes.Figure4shows the performance of the proposed controller when the torque gain k increases or decreases 20%. As can be seen when k increases the closed loop system behaves well. However, when k decreases the system requires a larger control input and more time to converge.

Figure 2 .

 2 Figure 2. between the proposed control algorithm and the original Teel's strategy when a disturbance is introduced at time 15 s: (a) x-position, (b) u control input, (c) θ orientation angle and (d) z-position.

Figure 3 .

 3 Figure 3. Robustness of the proposed algorithm with respect to variations of the mass value m:(a) x-position, (b) u control input, (c) θ orientation angle and (d) z-position.

Figure 4 .

 4 Figure 4. Robustness of the proposed algorithm with respect to variations of the gain value k:(a) x-position, (b) u control input, (c) θ orientation angle, (d) z-position.

Appendices Appendix 1. Bound for x

Introducing (10) into above

Notice that when |x 4 | > b + 1 it follows that V ˙ 1 < 0, and V 1 will decrease. So if |x 4 |≤ b + 1 at some time, it will remain so from then on.

Appendix 2. Bound for

thus

Therefore

Introducing (A20)

From ( 18)

or

Let t 3 be such that the RHS of the above satisfies

Then

Appendix 4. Bound for e 1 in (25)

Consider ( 23) and define

then the Lyapunov equation

holds for

then using ( 23) and ( A32) -(A35)

with

(10) and ( 25) it follows that

From ( A34) and ( A35)

Introducing the above into (A44) and (A42) leads to (33), i.e.

then using (31) 32) and ( 35)

Integrating and proceeding as in appendices 1 and 2 we conclude that there existsafinitetime t 5 such that (36) is satisfied. Therefore k = 0.0424 verifies both constraints. Notice that tan θ -θ is a strictly increasing function and the result will hold for values of |θ| smaller than π 9 .


Naoya Ohnishi
email: ohnishi@graduate.chiba-u.jp

Yukiko Kenmochi
email: y.kenmochi@esiee.fr

Atsushi Imiya
email: imiya@faculty.chiba-u.jp

TEMPORAL MODEL BASED ROBOT NAVIGATION USING OPTICAL FLOW 1

Keywords:

In this paper, we aim to develop an algorithm for the navigation of an autonomous robot using optical flow observed through a vision system mounted on a mobile robot. Optical flow is used for generating a temporal-model of the robot motion. The temporalmodel enables us to detect a planar area on which robot moves by matching with optical flow, and estimate a ego-motion of the mobile robot. If we obtain the planar area and a ego-motion, we can navigate the mobile robot autonomously. We show experiments for the detection of the planar area and estimation of the ego-motion of the robot using the image sequence observed through the mobile robot.

...

INTRODUCTION

Model-based method for object recognition provides a modelbased method for obstacle detection of mobile robots. CADbased method detects obstacles by matching geons [START_REF] Biederman | Human image understanding: Recent research and a theory[END_REF] in a database [START_REF] Dickinson | Decoupling recog-nition and localization in cad-based vision[END_REF]. Environment-map-based method detects obstacles using a environmental map without obstacles. The first method requires to propose geons detect-base. The second method requires a learning stage for the construction of environmental map in which the robot works. In this paper, we propose a temporal-model-based navigation strategy for autonomous mobile robots. Temporal-model-based method detects the feasible region where robot can move using three successive images observed form mobile robot. First two images are used for the generation of a short term model, second two images are used for data. Matching these two features from successive three images. Therefore, temporal-model-based method dose not require neither model in a database nor a learning stage for the generation model. For the detection of obstacles or planar areas, there are many methods using vision systems [START_REF] Guilherme | Vision for mobile robot navigation: A survey[END_REF]. For example, the methods of objects recognition using edge detection or landmarks observation are dependent on the environment around a robot. These methods involve difficulties, such as the need to predetermine landmarks, when applied to general environments. If a robot captures an image sequence of moving objects, the optical flow [START_REF] Barron | Performance of optical flow techniques[END_REF] [START_REF] Horn | Determining optical flow[END_REF] [6] [START_REF] Nagel | An investigation of smoothness constraint for the estimation of displacement vector fields from image sequences[END_REF], which is the flow of movement in the scene, is obtained for the fundamental features to generate environment information around the mobile robot. Additionally, optical flow is considered to be fundamental information for obstacle detection in the context of biological data processing [START_REF] Mallot | Inverse perspective mapping simplifies optical flow computation and obstacle detection[END_REF]. Therefore, the use of optical flow is an appropriate method for the detection of obstacles from the viewpoint of the affinity between the robot and human beings.

The obstacle detection using optical flow is proposed in [START_REF] Enkelmann | Obstacle detection by evaluation of optical flow fields from image sequences[END_REF] [START_REF] Santos-Victor | Uncalibrated obstacle detection using normal flow[END_REF]. Enkelmann [START_REF] Enkelmann | Obstacle detection by evaluation of optical flow fields from image sequences[END_REF] proposed an obstacle-detection method using the model vectors from motion parameters. Santos-Victor and Sandini [START_REF] Santos-Victor | Uncalibrated obstacle detection using normal flow[END_REF] also proposed an obstacledetection algorithm for a mobile robot using the inverse projection of optical flow to ground floor. Those methods require the learning stage for the construction of the model of a robot motion and a environment. Figure 1 shows the algorithm for the model-based method. In this method in the learning stage, a model for the environment without obstacles is generated.

Our algorithm generates a short term model of the flow field using first two images from successive three images. Furthermore, second two images from successive three images generates features for obstacle detection. Therefore, matching the short term model and features, the algorithm detects a dominant plane as the feasible region on which the mobile robot moves. Figure 2 shows the algorithm for the temporal-model-based method. This method generates temporal-models.

DOMINANT PLANE DETECTION FROM IMAGE SEQUENCE

In this section, we develop an algorithm for the detection of the dominant plane from the image sequence observed by a camera mounted on a mobile robot. When the mobile robot moves over the ground plane, we obtain successive images which include a dominant plane area and obstacles. The dominant plane is adapted as the model of own navigation strategy. In own method, from a triplet of successive images, the model is estimated and renewed temporally.

First, we show that the corresponding points in a pair of successive images, which are a projection of the dominant plane in a space, are connected by an affine transformation. Therefore, the computed optical flow from the successive images describes the motion of the dominant plane and obstacles on the basis of camera displacement. Since the corresponding points in dominant plane motion are combined by an affine transformation, we can compute the affine coefficients in Eq.( 1) from the optical flow on the dominant plane. Once the affine coefficients are computed, we can estimate the dominant plane motion in the image from the affine coefficients. The dominant plane motion is described by a planar flow field, as shown in Fig. 3. The difference between the estimated planar flow field and the computed optical flow field enables us to detect the dominant plane area on the image by matching the flow vectors on the image planes.

Approximation of dominant plane motion by affine transformation

We show that the corresponding points in a pair of successive images, which are a projection of the dominant plane in a space, are connected by an affine transformation. Therefore, the corresponding points (u, v) and (u , v ) on the dominant plane are expressed by

u v = a b d e u v + c f , (1)
when the camera displacement is small. This means that homography between the two images of a planar surface is approximated by an affine transformation when the camera displacement is small. Setting H to be the 3 × 3 matrix [START_REF] Hartley | Multiple view geometry in computer vision[END_REF], the homography between the two images of a planar surface is expressed as

p = Hp , (2)
where p = (u, v, 1) and p = (u , v , 1) are the corresponding points on the two images. The matrix H is expressed as where K, R, t, and n are the projection matrix, the rotation matrix, the translation vector, and the plane normal of the planar surface, respectively. The matrixes K and K -1 are affine transformations since these matrixes are the projection matrix of the pinhole camera. Assuming that the camera displacement is small, the matrix R and the matrix tn are approximated by an affine transformation. This approximation is illustrated in Fig. 4. These geometrical and mathematical assumptions are valid when the camera is mounted on a robot moving over the dominant plane. Then above assumptions allow us to describe the relationship between (u , v ) and (u, v) as the affine transformation. Therefore, setting a,b,c,d,e,f to be affine coefficients, Eq.( 2) can be expressed as Eq.( 1). In the next section, we develop an algorithm for the estimation of these six parameters from a sequence of images.

H = K(R + tn )K -1 , (3)

Planar flow estimation

Using a pair of successive images from a sequence of images obtained from the camera during robot motion, we compute the optical flow field ( u, v). Since optical flow is the correspondence of dense points between an image pair, Eq.( 1) can be applied to each point on the dominant plane. Let (u, v) and (u , v ) be the points of the corresponding image coordinates in the two successive images. Point (u , v ) can be calculated by adding the flow vector ( u, v) to (u, v), if (u, v) and (u , v ) are a pair of corresponding points. Thus, by setting ( u, v) to be the optical flow at point (u, v) in the image coordinate system, as shown in Fig. 5, we have the relation If we obtain (u, v) and (u , v ), we can compute the affine coefficients in Eq.( 1). If three points are non-collinear, Eq.( 1) has a unique solution. Since the dominant plane occupies the largest domain in the image, we select three points randomly to compute the affine coefficients in the successive pair of images. After we compute the affine coefficients, using Eq.( 1) again, we can compute the motion of the image sequence of the dominant plane, that is, the collection of corresponding points

u v = u v + u v . (4)
û v = u v - u v (5)
in the image sequence can be regarded as motion flow of the dominant plane on the basis of camera displacement

t t t+∆t u v u' v' a=(u,v) a'=(u',v') o o' (u,v) . .
(u , v ) is (u + u, v + v) and planar flow (û, v) is (u -u, v -v).
as shown in Fig. 6. We call this flow planar flow (û, v). This planar flow is used as temporal-model for dominant plane detection.

Dominant plane detection

Next, we compute the dominant plane area using the estimated planar flow and the computed optical flow. Setting ε to be the tolerance of the difference between the optical flow vector and the planar flow vector, if

ut ut - ût ût < ε (6)
is satisfied, we accept point (u t , v t ) as a point in the dominant plane.

In the case that at least one point on the obstacle area in the image is selected, the estimated planar flow is no longer the dominant plane motion. Therefore, the detected dominant plane area is very small. Since the dominant plane occupies the largest domain in the image, in such cases, it becomes evident that the selection of points is incorrect.

In those cases, we consequently select another three points randomly. Figure 7 shows examples of each case.

Once we have detected the dominant plane in a certain frame of the image sequence, the planar flow of subsequent images can be estimated robustly using the leastsquares method, because dense optical flows are used for the estimation of affine coefficients. Assuming that the robot displacement is small, the dominant plane of the successive images changes negligibly. Therefore, using the optical flow on the estimated dominant plane in the previous image, we estimate the affine coefficients using the least-squares method, as shown in Fig. 8. Setting (u i , v i ) and (u i , v i ) (0 ≤ i ≤ n) to be corresponding points, the Incorrect case. Correct case.

Figure 7: Examples of random sampling. Bottom-left is incorrect case, since the point is selected on the obstacle area. Select another points randomly. Bottom-right is correct case. This planar flow is equal to optical flow in 50% or more of area.

mean-squared errors E u and E v associated with Eq.( 1) are

E u = n i=1 {u i -(au i + bv i + c)} 2 , (7)
E v = n i=1 {v i -(du i + ev i + f )} 2 , ( 8
)
where n is the number of points used for estimation. Therefore, we can compute affine coefficients which minimize errors E u and E v . This algorithm generates planar flow as temporal-model at subsequent frames

Procedure for dominant plane detection

Our algorithm is summarized as follows.

1. Compute optical flow ( u, v) from two successive images.

2. Compute affine coefficients in Eq.( 1) by random selection of three points.

3. Estimate planar flow (û, v) from affine coefficients.

4. Match the computed optical flow ( u, v) and estimated planar flow (û, v) using Eq.( 6).

5. Detect the dominant plane. If the dominant plane occupies less than half of the image area, then return to step(2).

Figure 9 shows the procedure of dominant plane detection from the image sequence.

Navigation using the dominant plane

We describe an algorithm for navigation of a mobile robot using the dominant plane. The percentages of the dominant plane in the image determine for the robot a strategy of a robot motion. If the percentage of the dominant plane in the image is greater than 80% and the dominant plane of the left side of the image is greater than that of right side, the robot moves forward. If the percentage of the dominant plane in the image is less than 80% and the dominant plane of the left side of the image is greater than that of right side, the robot rotates counterclockwise to avoid collision of obstacles. If the percentage of the dominant plane in the image is less than 80% and the dominant plane of the right side of the image is greater than that of left side, the robot rotates clockwise. These algorithm are summarized in Fig. 10.

EXPERIMENT

In this section, we evaluate the performance of our method for steps (2) and ( 5) of the algorithm listed in Fig. 9, since the estimation of affine coefficients and dominant plane detection are essential to the resolution of our problem. In Step (2), the parameters are estimated from the image sequence. In Step [START_REF] Horn | Determining optical flow[END_REF], the area of the dominant plane in the images is evaluated. These experiments are described in subsections 4.1 and 4.2, respectively. For the computation of optical flow, we use the Lucas-Kanade method with pyramids [START_REF] Bouguet | Pyramidal implementation of the lucas kanade feature tracker description of the algorithm[END_REF]. The tolerance for the matching of flow vectors in Eq.( 6) is set to be ε = 0.2, which was determined experimentally. To evaluate the accuracy of the dominant plane detection, we define the error ratio E(t). Setting D(u, v, t) and M (u, v, t) to be the dominant planes on the image coordinate system (u, v) at frame t detected from the optical flow and detected manually, respectively, and setting

D(u, v, t) = 0 for dominant plane 1 for other area , (9)
M (u, v, t) = 0 for dominant plane 1 for other area ,

we evaluate the ratio

E(t) = 100 u,v |M (u, v, t) -D(u, v, t)| u,v 1 (11)
as a function of time t.

The error rate of affine coefficients with uniform motion

First, we show the validity of the affine transformation in two successive images. From Eq.( 1), we define matrix A t as

A t =   a b c d e f 0 0 1   , ( 12
)
where t is the frame number in the image sequence. Assuming that the velocity of camera displacement is constant in the two successive of images, the relation A t+1 A t is satisfied. Therefore, we define the error of affine coefficients as

e(t) = |A t+1 -A t | 2 . (13)
In this experiment, we use the marbled-block image sequence [START_REF]Marbled-block sequence: recorded and first evaluated by otte and nagel[END_REF], since the velocity of camera displacement should be constant. Figure 11 shows the error of affine coefficients. Figs.12 and 13 show the image sequence of the detected dominant plane and the error ratio of the dominant plane, respectively. Figure 11 shows that |A t+1 -A t | 2 approaches zero with time. Therefore, the affine coefficients are constant when the camera displacement is constant. We show the validity of the affine transformation in two successive images.

Detection of dominant plane

We show next the validity of the random selection of three points for the computation of affine coefficients. In this experiment, we use image sequences of forward motion is the error of affine coefficients e(t), as given by Eq.( 13).

The horizontal axis is the frame number of the image sequence. and rotational motion, because the motion of a mobile robot mainly consists of these motions. For forward motion, we prepared image sequences in which the percentages of the dominant plane are 90%, 80% and 70%, since we assume that the dominant plane occupies more than half of an image. These percentages of the dominant plane were computed using M (u, v, t) in Eq. [START_REF] Santos-Victor | Uncalibrated obstacle detection using normal flow[END_REF]. For rotational motion, we use the rotating blocks image sequence [START_REF]Rotating blocks: Otago optical flow evaluation sequences[END_REF].

The result of dominant plane detection is shown in Fig. 14. The first row is the original image, and from the second row to the fifth row, we show image sequences of the dominant plane in the first, second, third and fourth frames, respectively. In the three images of the dominant plane, the white areas are the dominant planes and the black areas are the obstacles. Table 1 lists error ratios evaluated with Eq. [START_REF] Hartley | Multiple view geometry in computer vision[END_REF].

Dominant plane detection in the long sequence

We present the dominant plane detection in a long image sequence for the navigation and path planning of the mobile robot. Figure 15 and Fig. 16 show the results for the dominant plane detection and the error ratio, respectively. As evident from Fig. 16, even if the error ratio increases in a sequence, the error ratio becomes small with a subsequent frame. Therefore, our method detects the dominant plane without accumulation of the error ratio. Our algorithm runs on a 2.0 GHz Pentium4 for 320 × 240 images and takes 0.25 seconds per frame for the detection of the dominant plane.

Autonomous navigation of the mobile robot

Using the algorithm of Section 3.4., we experiment for the obstacle avoidance of the mobile robot. Figure 17 shows the result of the experiment. In the top row of Fig. 17, at the 16 frame, the robot rotates counterclockwise since the robot detects the obstacle in the right side of the image.

In the middle row of Fig. 17, at the 68 frame, the robot rotates clockwise since the robot detects the obstacle in the left side of the image. In the bottom row of Fig. 17, at the 100 frame, the robot moves forward since the percentage of the dominant plane in the image is greater than 80%.

Figure 18 shows the top view of the experimental environment. Figure 18 is generated by projecting dominant planes in the image plane onto the ground plane. Figure 18 is generated by projecting dominant planes in the image plane onto the ground plane. The trajectory of the mobile robot is also obtained by projecting the planar flow in the image plane onto the ground plane.

CONCLUSION

We developed an algorithm for the navigation of a mobile robot by temporal-model-based planar area detection using optical flow. The algorithm allows the positions of obstacles and the location of the robot to be detected using optical flow using optical flow computed from the image sequence observed through the camera mounted on the mobile robot. Furthermore, we showed that corresponding points on dominant planes in a pair of succes- The line started at the origin is the trajectory of the location of the mobile robot. We set the vertical axis to coincide with the direction of the robot at initial location, and the horizontal axis to be orthogonal the vertical axis.

sive images are combined by affine transformation. Using this idea, if we compute the affine coefficients which relate the corresponding points in two successive images, we can easily obtain a dense planar flow which expresses a camera motion. This property of the points in a dominant plane allows us to design an algorithm which enable the dominant plane to be detected by simple pattern matching of the flow vectors in a series of dominant planes.

Although model-based approach to dominant plane detection has been proposed in [START_REF] Enkelmann | Obstacle detection by evaluation of optical flow fields from image sequences[END_REF], our method is a nonmodel-based approach. In addition, our algorithm allows the dominant plane to be detected without camera calibration, since our algorithm uses short-term model of the robot motion. Results of experiments using real image sequences confirmed that the dominant plane can be detected accurately. These experiments allow the application of our method to the navigation and path planning of a mobile robot with a vision system.

Figure 1 :

 1 Figure 1: Model-based navigation using optical flow. The upper part is the learning stage for the construction of the model of a robot motion and a environment. During a navigation, obstacles are detected by matching optical flow of this model and observed ones.

Figure 2 :

 2 Figure 2: Temporal-model-based navigation using optical flow. This method detects obstacles without a model. The first two images from triplet of images are used for the construction of temporal-models.

Figure 3 :

 3 Figure 3: Planar flow of the image sequence. top: Example of camera displacement and the environment. leftbottom: Computed optical flow. right-bottom: Estimated planar flow. In a top-middle area, where exists the obstacle, the length of optical flow is bigger than planar flow.

2 XFigure 4 :

 24 Figure 4: Approximation of dominant plane motion by affine transformation. top: Homography is expressed by x = Hx. bottom: Affine transformation is expressed by x = Ax.

Figure 5 :

 5 Figure 5: Optical flow generated by affine transformation. Optical flow ( u, v) on the dominant plane is described as an affine transformation.

  v)=(u'-u,v'-v)

Figure 6 :

 6 Figure 6: Optical flow computed from corresponding points in two successive images. Optical flow is described as the correspondence of the same points in the two successive images. Therefore, point (u , v ) is (u + u, v + v) and planar flow (û, v) is (uu, vv).

  (u t ,v t ),(u t+1 ,v t+1 ) dominant plane

Figure 8 : 9 :

 89 Figure8: For the first frame, our algorithm uses a pair of images for the detection of the dominant plane. For subsequence frames, the dominant plane in the previous frame is used for the application of the least-squares method.

Figure 10 :

 10 Figure 10: Strategy of a mobile robot. left: Robot moves forward. middle: Robot rotates counterclockwise. right: Robot rotates clockwise.

Figure 11 :

 11 Figure11: Error of affine coefficients. The vertical axis is the error of affine coefficients e(t), as given by Eq.(13). The horizontal axis is the frame number of the image sequence.

Figure 12 :

 12 Figure 12: Detected dominant plane in frames 0, 10 and 20. The original image(top-left), and detected dominant plane in frames 0(top-right), 10(bottom-left) and 20(bottom-right). In the images of the dominant plane, the white areas are the dominant planes and the black areas are obstacles.

Figure 13 :

 13 Figure13: Error ratio of dominant plane detection using the image sequence shown in Fig.12. The vertical axis is the error ratio E(t) expressed by Eq.(11). The horizontal axis is the frame number of the image sequence.

Figure 14 and

 14 Figure 14 and Tab.1 show that the validity of regions detected as the dominant plane increases with time.

Figure 14 :

 14 Figure 14: Detected dominant plane. The first and second column are the original image and the dominant plane, respectively.

Figure 15 :

 15 Figure 15: Detected dominant plane at the 50, 100, 150 and 200 frames.

Figure 16 :

 16 Figure16: Error ratio E(t). The vertical axis is the error ratio E(t) in Eq.[START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. The horizontal axis is the frame number of the image sequence.

Figure 17 :

 17 Figure 17: The left column shows the images of the experimental environment, the middle column shows the images observed by the mobile robot, and the right column shows the dominant planes computed by our algorithm. The top, middle, and bottom rows are captured at frame 16, 68, and 100, respectively.

Figure 18 :

 18 Figure18: Environmental map computed from dominant plane. White and black area are planar area and obstacles, respectively. The line started at the origin is the trajectory of the location of the mobile robot. We set the vertical axis to coincide with the direction of the robot at initial location, and the horizontal axis to be orthogonal the vertical axis.

Table 1 :

 1 Error ratios of dominant plane detection in Fig.14. In the second row, F90, F80 and F70 are the image sequences of forward motion in which the percentages of the dominant plane in the image are 90%, 80% and 70%, respectively. Rotation gives the image sequence of rotational motion. This image sequence is shown in the first row of Fig.14.

			error ratio E(t) [%]
	frame	F90	F80	F70	Rotation
	0 23.22 24.01 20.33	19.37
	1 14.40 12.99 12.29	14.50
	2 10.91 11.83 10.45	11.09
	3 10.54 11.08 10.55	11.29