José Vilca

Lounis Adouane

Youcef Mezouar

A novel safe and flexible control strategy based on target reaching for the navigation of urban vehicles

Keywords: Mobile robot navigation, Control architecture, Target reaching and following, Lyapunov stability

This paper presents a complete framework for reactive and flexible autonomous vehicle navigation. A human driver reactively guides a vehicle to the final destination while performing a smooth trajectory and respecting the road boundaries. The objective of this paper is to achieve similar behavior in an unmanned ground vehicle to reach a static or dynamic target location. This is achieved by using a flexible control law based on a novel definition of control variables and Lyapunov synthesis. Furthermore, a target assignment strategy to enable vehicle navigation through successive waypoints in the environment is presented. An elementary waypoint selection method is also presented to perform safe and smooth trajectories. The asymptotic stability of the proposed control strategy is proved. In addition, an accurate estimation of the maximum error boundary, according to the controller parameters, is given. With this indicator, the vehicle navigation will be safe within a certain boundaries. Simulations and experiments are performed in different cases to demonstrate the flexibility, reliability and efficiency of the control strategy. Our proposal is compared with different navigation methods from the literature such as those based on trajectory following.

Introduction

Autonomous vehicle navigation is a complex problem of major interest to the research community. Systems capable of performing efficient and robust autonomous navigation are unquestionably useful in many robotic applications such as manufacturing technologies [START_REF] Sezen | Modeling automated guided vehicle systems in material handling, Otomatiklestirilmi Rehberli Arac Sistemlerinin Transport Tekniginde Modellemesi[END_REF], urban transportation [START_REF] Avanzini | On-line reference trajectory generation for manually convoying a platoon of automatic urban vehicles[END_REF], assistance to disabled or elderly people [START_REF] Martins | Assistive mobility devices focusing on smart walkers: Classification and review[END_REF] and surveillance [START_REF] Stoeter | Papanikolopoulos, A robot team for surveillance tasks: Design and architecture[END_REF]. Although much progress has been made, some specific technologies have to be improved for effective application in real environments. This paper particularly focuses on the problem of autonomous navigation of vehicles in an urban environment (cf. Fig. 1).

Different strategies for autonomous navigation have been proposed in the literature [START_REF] Lee | A unified framework of the automated lane centering/changing control for motion smoothness adaptation[END_REF], [START_REF] Adouane | Mobile robot navigation in cluttered environment using reactive elliptic trajectories[END_REF], [START_REF] Gu | On-road motion planning for autonomous vehicles[END_REF].

The most popular approaches are based on following a pre-defined reference trajectory [START_REF] Bonfè | Online trajectory generation for mobile robots with kinodynamic constraints and embedded control systems[END_REF], [START_REF] Kanayama | A stable tracking control method for an autonomous mobile robot[END_REF]. Most of the proposed control laws are dedicated to trajectory tracking (to track a time-parametrized reference) [START_REF] Blazic | Four-state trajectory-tracking control law for wheeled mobile robots[END_REF] and path following (to follow a path without explicit temporal references) [START_REF] Samson | Control of chained systems. application to path following and time-varying point-stabilization of mobile robots[END_REF]. These methods link the control to a reference trajectory which could be defined by a combination of path planning and trajectory generation techniques [START_REF] Lavalle | Planning Algorithms[END_REF].

Typically, to obtain the reference path to be followed by the robot, arc-lines, B-splines or polynomial equations are used over points [START_REF] Connors | Manipulating b-spline based paths for obstacle avoidance in autonomous ground vehicles[END_REF], [START_REF] Horst | Trajectory generation for an onroad autonomous vehicle[END_REF], [START_REF] Lee | A unified framework of the automated lane centering/changing control for motion smoothness adaptation[END_REF]. In [START_REF] Gu | On-road motion planning for autonomous vehicles[END_REF] a feasible path is obtained using a polynomial curvature spiral. In [START_REF] Bonfè | Online trajectory generation for mobile robots with kinodynamic constraints and embedded control systems[END_REF], the trajectory generation method provides a smooth path considering the kinodynamic constraints of the vehicle. In [START_REF] Courbon | Autonomous navigation of vehicles from a visual memory using a generic camera model[END_REF], straights line paths defined by the position and orientation of a single waypoint are considered. In this case, the orientation of the previous waypoint is not taken into account to simplify the implementation of the control law. Nevertheless, trajectory generation presents some drawbacks, such as the necessity of a specific planning method, the proof of guarantee of continuity between different segments of the trajectory and the complexity for replanning. A few works in the literature propose to use only specific set of way-points in the environment to lead the robot toward its final objective. In [START_REF] Aicardi | Closed loop steering of unicycle like vehicles via lyapunov techniques[END_REF], the authors propose a navigation strategy via assigned static points for a unicycle robot. This strategy does not allow accurate navigation since the kinematic constraints of the robot (maximum velocity and steering), the orientation error and the velocity profile of the robot when it reaches the assigned point are not considered. Harmonic Potential Field (HPF) is used to guide an unmanned aerial vehicle (UAV) to a global waypoint with a position and a direction of arrival in [START_REF] Masoud | A harmonic potential approach for simultaneous planning and control of a generic uav platform[END_REF]. The author proposes a virtual velocity field which allows to consider the UAV model. Each vector component of the field is treated as an intermediate waypoint with which the robot must comply with it in order to reach the global waypoint. Nonetheless, HPF requires a complex mathematical modeling for different shapes or dimension of the obstacles in the environment. In this paper, we propose a navigation strategy which avoids the pre-generation of any specific reference trajectory. Vehicle movements are obtained according to the proposed control law while considering vehicle kinematic constraints and sequential waypoints to reach (defined by its position, orientation and velocity). The vehicle can thus perform different movements between waypoints without the necessity of replanning any reference trajectory, and it can also add or change the location of the successive waypoints according to the environment configuration or to the task to achieve. Thus, this strategy allows flexible navigation while taking into account appropriate waypoints suitably placed in the environment.

Different control methods for trajectory tracking and path-following dedicated to wheeled mobile robots (unicycle, car-like robot, etc.) have been proposed in the literature [START_REF] Blazic | Four-state trajectory-tracking control law for wheeled mobile robots[END_REF], [START_REF] Fang | Trajectory tracking control of farm vehicles in presence of sliding[END_REF] and [START_REF] Pesterev | Stabilizing control for a wheeled robot following a curvilinear path[END_REF]. In [START_REF] Kanayama | A stable tracking control method for an autonomous mobile robot[END_REF], [START_REF] Blazic | Four-state trajectory-tracking control law for wheeled mobile robots[END_REF] and [START_REF] Vilca | Reactive navigation of a mobile robot using elliptic trajectories and effective online obstacle detection[END_REF], nonlinear control laws for trajectory tracking are synthesized for a unicycle robot using Lyapunov stability analysis. The Lyapunov functions used in these studies are based only on distance and orientation errors. A trajectory tracking control for a farm vehicle, incorporating sliding in the kinematic model, is proposed in [START_REF] Fang | Trajectory tracking control of farm vehicles in presence of sliding[END_REF]. For the path-following problem, a control law for a tricycle robot is proposed in [START_REF] Pesterev | Stabilizing control for a wheeled robot following a curvilinear path[END_REF] and [START_REF] Samson | Control of chained systems. application to path following and time-varying point-stabilization of mobile robots[END_REF]. They are based on feedback linearization and chained form representation [START_REF] Laumond | La Robotique Mobile[END_REF]. The path-following controller thus allows to make the lateral and longitudinal control of the vehicle independent along the reference trajectory. Furthermore, the path-following controller allows smoother convergence to the desired path than the trajectory tracking controller (designed for a time-parametrized trajectory) [START_REF]Springer Handbook of Robotics, Part E-34[END_REF]. The trajectory tracking controller allows to track the trajectory with a desired velocity profile, while the path-following controller acts only on the orientation to drive it along the path. Both, the pathfollowing and trajectory tracking controllers require the pose of the closest point to the trajectory (w.r.t. robot configuration) and/or the value of curvature at this point (cf. Fig. 15) at each sample time [START_REF] Laumond | La Robotique Mobile[END_REF] and [START_REF] Blazic | Four-state trajectory-tracking control law for wheeled mobile robots[END_REF]. Although there exist a multitude of techniques to compute these parameters, they can add an error in certain situations thereby influencing the control of the mobile robot negatively [START_REF]Springer Handbook of Robotics, Part E-34[END_REF] and [START_REF] Fang | Trajectory tracking control of farm vehicles in presence of sliding[END_REF]. In this paper, a control law based on a novel definition of control variables and Lyapunov function is proposed. The synthesized control law can perform either static or dynamic target reaching using only its current pose and velocity. Using dynamic target reaching, trajectory following can also be performed. The control law exhibits good flexibility properties and it could be adapted to different autonomous robotic applications such as multi-robot formations (cf. Subsection 5.2).

It is not always required to follow a fixed path with high fidelity, specifically in open or lowconstrained environments. We will demonstrate in this paper that only few waypoints will be sufficient to guarantee safe and flexible mobile robot navigation. A target assignment strategy is also proposed to perform autonomous navigation through pre-defined waypoints. We will also demonstrate that if we increase the number of these waypoints, the robot control performs as if we had applied common trajectory tracking control.

The rest of the paper is organized as follows: the following section presents the studied problem. Section 3 details the vehicle and target models, the proposed control law, its stability properties and setting the controller parameters. In Section 4 the proposed target assignment strategy and the waypoint selection approach are described. Simulation and experimental results are given in Section 5. A study with a few navigation methods from the literature is also provided. Finally, Section 6 provides a conclusion and perspectives for future studies.

Problem Statement

An important challenge in the field of autonomous vehicles consists of ensuring safe and flexible navigation in a structured environment (cf. Fig. 1 and2). In this work, safe navigation consists in not crossing over the road limits and bumping into obstacles while respecting the physical constraints of the vehicle. Flexible navigation consists in allowing different possible movements to achieve the task, while guaranteeing a smooth trajectory of the vehicle. The main idea of the proposed work is to guarantee both criteria simultaneously.

We consider the following scenario (cf. Fig. 2):

• The structured environment is a known road map where the roads have a specific width w R .

Figure 2: Nominal scenario with a road map and the task to achieve by the vehicle in its environment.

• The vehicle model (kinematic) is known.

• The vehicle starts at the initial pose P i and it has to reach the final pose P f (in certain conditions, P i = P f ).

As presented in section 1 and according to the presented scenario, a safe reference path in static environment can be obtained by different algorithms such as Voronoï diagram [START_REF] Latombe | Robot Motion Planning[END_REF], potential fields [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF] or others [START_REF] Lavalle | Planning Algorithms[END_REF]. In our case, specific key positions should be defined in the static environment, which we name waypoints. Their numbers and configurations in the environment are detailed in Subsection 4.2. Obviously, the waypoints can also be selected from a pre-defined trajectory if available. Consequently, the navigation problem is simplified to a waypoint tracking problem, i.e, the vehicle is guided by the waypoints (cf. Fig. 8) instead of following a specific fixed path. The vehicle has thus to reach each waypoint with a defined position, orientation and velocity while satisfying distance and orientation error limits (E dis and E angle respectively) to perform safe navigation (cf. Subsection 3.4).

The proposed navigation strategy can be extended easily in order to deal with dynamic environments, notably using the limit-cycle approach [START_REF] Adouane | Mobile robot navigation in cluttered environment using reactive elliptic trajectories[END_REF] and [START_REF] Vilca | Reactive navigation of a mobile robot using elliptic trajectories and effective online obstacle detection[END_REF]. This obstacle avoidance approach allows to modify locally the movement of the robot to avoid dangerous static or dynamic obstacles and to come back to its initial plan. Nevertheless, it is not the main focus of this paper. A single waypoint, defined as a set-point by its position, orientation and velocity, is called a target. The control law to reach this target is designed to generate smooth vehicle navigation. Moreover, if the successive waypoints are close to each other then the vehicle tends to perform trajectory tracking navigation (cf. Section 4 and Subsection 5.1.1). Fig. 3 gives an overview of the proposed navigation strategy. In this paper, waypoint selection (dashed blue block in Fig. 3) is simplified (cf. Subsection 4.2) to focus on navigation strategy (Control law and Target assignment blocks, cf. Section 3 and Subsection 4.1 respectively) .

Target-reaching control

This section gives the details of the proposed control law, the proof of its stability and the tuning of controller parameters.

Vehicle and target modeling

Our experimental vehicle (cf. Fig. 19) is devoted to urban transportation. The vehicle moves on asphalt with a low velocity (less than 3 m/s). Therefore, it appears quite natural to rely on a kinematic model, and to assume pure rolling and non-slipping at wheel ground contact. In such cases, vehicle modeling is commonly performed, for instance relying on a tricycle model as shown below (cf. Fig. 4):

   ẋ = v cos(θ) ẏ = v sin(θ) θ = v tan(γ)/l b (1)
where O G and O m are respectively the origin of the global and local reference frames, (x, y, θ) is the pose (configuration state) at vehicle point O m , γ is the orientation of the vehicle front wheels, v is the linear velocity at vehicle point O m . l b and w b are respectively the wheelbase and the track width of the vehicle (cf. Fig. 4). I cc is the instantaneous center of curvature of the vehicle trajectory. The radius of curvature r c is given by:

r c = l b / tan(γ) (2)
and c c = 1/r c is the curvature of the vehicle trajectory.

Let us consider a dynamic target modeled as a point with non-holonomic constraints (cf. Fig. 4). This model allows us to use the general model of robot motion and also to simplify the controller equations. Its kinematic characteristics are given by:

   ẋT = v T cos(θ T ) ẏT = v T sin(θ T ) θT = ω T (3)
where v T and ω T are respectively the linear and angular velocities of the target. The radius of curvature is computed by r c T = v T /ω T . An important consideration for target reaching is v T ≤ v max and r c T ≥ r cmin , where v max and r cmin are respectively the maximum linear velocity and the minimum radius of curvature of the vehicle, given by r cmin = l b / tan(γ max ). For static target reaching (point stabilization, i.e., to reach a specific point with a given orientation), ω T is considered equal to zero and v T is not necessarily equal to zero; v T is then considered as a desired velocity value for the vehicle when it reaches the desired target (x T , y T , θ T ).

Control law

Before presenting the control law, let us introduce the control variables of the system (cf. Fig. 5). The errors with respect to the local frame (X m Y m ) of the vehicle (e x , e y , e θ ) between the desired pose (x T , y T , θ T ) and the current vehicle pose (x, y, θ) are given by:

   ex = cos(θ)(xT -x) + sin(θ)(yT -y) ey = -sin(θ)(xT -x) + cos(θ)(yT -y) e θ = θT -θ (4)
A new error function e RT is added to the canonical error system (4) (cf. Fig. 5). Let us first define the distance d and the angle θ RT between the target and the vehicle position as (cf. Fig. 5):

d = (xT -x) 2 + (yT -y) 2
(5)

θRT = arctan ((yT -y)/(xT -x)) if d ≥ ξ θRT = θT if d < ξ ( 6
)
where ξ is a small positive value (ξ ≈ 0).

The error e RT is related to the vehicle position (x, y) with respect to the target orientation (cf. Fig. 5). It is defined as:

e RT = θ T -θ RT (7)
Furthermore, e RT can be written as a function of e x , e y and e θ as:

tan(e RT ) = tan(e θ -(θ RT -θ)) =
tan(e θ )e y /e x 1 + tan(e θ )e y /e x = e x tan(e θ )e y e x + tan(e θ )e y

where tan(θ RTθ) = e y /e x (cf. Fig. 5). Hence, e RT allows to consider an additional orientation error w.r.t. e x , e y and e θ , e.g., when e θ = 0 then e RT = -e y /e x . The stabilization of this error allows to decrease the lateral distance d l to zero (18) (cf. Fig. 5), and to always have the robot in the wake of the target. The derivatives of the errors ( 4) and ( 7) can be obtained using (1), ( 3), ( 5) and ( 6):

ėx = cos(θ)( ẋT -ẋ) + sin(θ)( ẏT -ẏ) -sin(θ)(xT -x) θ + cos(θ)(yT -y) θ = -v + ey θ + vT [cos(θT ) cos(θ) + sin(θT ) sin(θ)] = -v + eyv tan(γ)/l b + vT cos(e θ ) (9) ėy
= -sin(θ)( ẋT -ẋ) + cos(θ)( ẏT -ẏ) -cos(θ)(xT -x) θ -sin(θ)(yT -y) θ = -ex θ -vT cos(θT ) sin(θ) + vT sin(θT ) cos(θ) = -exv tan(γ)/l b + vT sin(e θ ) ( 10
) ėθ = θ T - θ = ωT -ω = vT rc T -v tan(γ)/l b (11) ėRT = θ T -θRT = vT rc T - d dt arctan yT -y xT -x = vT rc T -vT sin(θT )(xT -x) -cos(θT )(yT -y) d 2 - -v sin(θ)(xT -x) + v cos(θ)(yT -y) d 2 = vT rc T - vT ex sin(e θ ) d 2 + vT ey cos(e θ ) d 2 - eyv d 2 (12)
The control law to reach a target (static or dynamic) is obtained using the Lyapunov stability analysis framework. The desired linear velocity v and the front wheel orientation γ of the vehicle which lead the errors (e x , e y , e θ ) to converge to zero are chosen such that:

v = v T cos(e θ ) + v b (13) γ = arctan(l b c c ) ( 14
)
where v b and c c are given by:

v b =Kx [K d ex + K l d sin(eRT ) sin(e θ ) + Ko sin(e θ )cc] ( 15
) cc = 1 rc T cos(e θ ) + d 2 K l sin(eRT ) cos(eRT ) rc T Ko sin(e θ ) cos(e θ ) + K θ tan(e θ ) + K d ey -K l d sin(eRT ) cos(e θ ) Ko cos(e θ ) + KRT sin 2 (eRT ) sin(e θ ) cos(e θ ) ( 16
) K = (K d , K l , K o , K x , K θ , K RT
) is a vector of positive constants defined by the designer (cf. Subsection 3.4).

Stability analysis

In this section, the stability of the error system ( 9), ( 10), ( 11) and ( 12) is analyzed.

Assumption 1. The subsequent development is based on the assumption that the initial values of e RT and e θ satisfied:

e RT ∈ ]-π/2, π/2[ and e θ ∈ ]-π/2, π/2[ (17)
These conditions guarantee that the target is ahead of the vehicle with respect to its orientation.

Theorem 1. The control law given by ( 13) and ( 14) ensures that the differential system ( 9), ( 10), ( 11) and ( 12) is asymptotically stable according to Lyapunov-based analysis and if the Assumption 1 is satisfied [START_REF] Khalil | Nonlinear Systems[END_REF].

Proof. 1. The stability of the system is analyzed using the Lyapunov method [START_REF] Khalil | Nonlinear Systems[END_REF]. The proposed Lyapunov function V , given by [START_REF] Fang | Trajectory tracking control of farm vehicles in presence of sliding[END_REF], is a function of three parameters which depend on the distance d between the target and vehicle positions, the distance d l from the vehicle to the target line (line which passes through the target position with an orientation equal to the target orientation; this term is related to the Line of Flight and Sight of the target [START_REF] Siouris | Missile Guidance and Control Systems[END_REF]), and the orientation error e θ between the vehicle and the target (cf. Fig. 5).

The candidate Lyapunov function V is a positivedefinite function [START_REF] Khalil | Nonlinear Systems[END_REF] when considering [START_REF] Masoud | A harmonic potential approach for simultaneous planning and control of a generic uav platform[END_REF]. It is given by:

V = 1 2 K d d 2 + 1 2 K l d 2 l + K o [1 -cos(e θ )] = 1 2 K d d 2 + 1 2 K l d 2 sin 2 (e RT ) + K o [1 -cos(e θ )] (18)
It can be written with respect to e x , e y as follows:

V = 1 2 e 2 x + e 2 y [K d +K l sin 2 (e RT )]+K o [1-cos(e θ )] (19)
To guarantee the stability of the system, V must be negative-definite [START_REF] Khalil | Nonlinear Systems[END_REF]. By taking the derivate of (19) and using ( 9), ( 10), ( 11), ( 12), ( 13) and ( 14), V can be written:

V =(ex ėx + ey ėy)[Kd + K l sin 2 (eRT )] + K l d 2 sin(eRT ) cos(eRT ) ėRT + Ko sin(e θ ) ėθ =(ex [eyvcc -v + vT cos(e θ )] + ey [vT sin(e θ ) -exvcc])• [K d + K l sin 2 (eRT )] + K l d 2 sin(eRT ) cos(eRT )• vT rc T - vT ex sin(e θ ) d 2 + vT ey cos(e θ ) d 2 - eyv d 2 + Ko sin(e θ ) vT rc T -vcc (20)
Substituting ( 13) in ( 20)

V = [-exv b + vT ey sin(e θ )] [K d + K l sin 2 (eRT )] + K l sin(eRT ) cos(eRT ) d 2 vT rc T -vT ex sin(e θ ) -eyv b + Ko sin(e θ ) vT rc T -vT cos(e θ )cc -v b cc = ey(K d + K l sin 2 (eRT )) -exK l sin(eRT ) cos(eRT ) • vT sin(e θ ) + vT rc T d 2 K l sin(eRT ) cos(eRT ) + Ko sin(e θ ) -vT Ko sin(e θ ) cos(e θ )cc -v b ex(K d + K l sin 2 (eRT )) -v b [eyK l sin(eRT ) cos(eRT ) + Ko sin(e θ )cc] (21)
Using ( 8) in the first and last terms of ( 21) and factorizing the common terms, it holds that: -vT Ko sin(e θ ) cos(e θ )cc [START_REF]Springer Handbook of Robotics, Part E-34[END_REF] Finally, using ( 15) and ( 16) in [START_REF]Springer Handbook of Robotics, Part E-34[END_REF], we obtain:

V = -Kx[K d ex + K l d sin(eRT ) sin(e θ ) + Ko sin(e θ )cc] 2 -vT KoK θ sin 2 (e θ ) -vT KoKRT sin 2 (eRT ) ≤ 0 (23)
Equation ( 23) shows that the system is stable while the initial conditions [START_REF] Masoud | A harmonic potential approach for simultaneous planning and control of a generic uav platform[END_REF] are satisfied. To ensure the asymptotic stability of the proposed control law, V has to be a negative-definite function. Let us consider two cases, one where V = 0 with v T > ξ and another with v T = ξ, where ξ is a constant value (ξ ≈ 0). Firstly, when v T > ξ and using the initial assumption K > 0, it is straightforward to show that e x , e θ , e RT are equal to zero to satisfy [START_REF] Latombe | Robot Motion Planning[END_REF]; then according to [START_REF] Gu | On-road motion planning for autonomous vehicles[END_REF], ( 6) and [START_REF] Masoud | A harmonic potential approach for simultaneous planning and control of a generic uav platform[END_REF], d is equal to zero (e y = 0). Hence, V is equal to zero when v T > ξ, only if (e x , e y , e θ ) = (0, 0, 0).

Secondly, let us consider the case where v T = ξ. The initial assumption is identical. Therefore, the second and third terms of ( 23) are equal to zero when v T = ξ. Additionally, when v T = ξ (static case) then we can consider r c T → ∞ (cf. Subsection 3.1); consequently the first term of V is equal to zero when:

K d e x +
Using ( 8) in ( 25), we finally obtain:

K d d cos(eRT ) cos(e θ ) + KoK θ sin 2 (e θ ) cos(e θ ) + KoKRT sin 2 (eRT ) cos(e θ ) = 0 (26)
Equation ( 26) exhibits quadratic terms. Consequently, considering the initial conditions [START_REF] Masoud | A harmonic potential approach for simultaneous planning and control of a generic uav platform[END_REF], cos(e RT ) and cos(e θ ) are greater than zero. Therefore, all the terms of ( 26) are positive and they must be equal to zero (i.e., d, e θ , e RT = 0, and if d = 0 then e x , e y = 0). Hence, from [START_REF] Siouris | Missile Guidance and Control Systems[END_REF], V is equal to zero when v T = ξ and r c T → ∞, only if (e x , e y , e θ ) = (0, 0, 0).

In conclusion, if v T > ξ or v T = ξ, V is always strictly positive and V is always strictly negative while (e x , e y , e θ ) = (0, 0, 0). Therefore, the system is asymptotically stable while the initial vehicle conditions (17) are satisfied.

Safe target reaching

Synthesis of the proposed control law using a Lyapunov function enables us to confirm its asymptotic stability. Nonetheless, it does not allow us to obtain immediately the error values when the robot is in the immediate vicinity of the target to reach. The aim of this subsection is to determine a relation between the upper bound of the errors d and e θ , denoted E dis and E angle (cf. Fig. 5 and8) and the controller parameters K. Indeed, according to these dynamics and considering that the vehicle and target localization are always accurate, the values of K enable the vehicle constraints (maximum velocity, acceleration, steering, etc.) to be satisfied. Our analysis consists in determining the minimum d i (cf. Fig. 6) which allows to satisfy at the same time, the vehicle physical constraints and the errors (d and e θ ), when the vehicle reaches the target, which must be less or equal to the pre-defined bounds E dis and E angle respectively. Fig. 7 shows the global scheme that we will apply, discussed subsequently, to obtain the minimal initial distance d i while knowing K and e θo which will satisfy E dis and E angle .

The subsequent analysis considers a static target ( ẋT = ẏT = 0 and r c T → ∞) and a limit vehicle configuration, e RT ≈ 0 and e θ ≈ -pi/2, i.e., the vehicle has the maximum admissible orientation error with respect to the target and the convergence will be the slowest (cf. Fig. 6 and5). This analysis will allow reference values for the controller parameters K to be obtained in the case of a dynamic target, where the distance d between the vehicle and target change more slowly than in the static target case. The control law with the designed parameters K will have thus more time to converge the system errors to zero (4).

To simplify controller analysis, the errors in orientation e θ and distance d are dealt separately. Firstly, the orientation error is computed consider- ing enough initial distance d i and e RT ≈ 0. These considerations enable us to estimate the minimum time to attain effectively E angle (cf. Fig 6 and8). Using ( 13), ( 14) and ( 15) in [START_REF] Samson | Control of chained systems. application to path following and time-varying point-stabilization of mobile robots[END_REF], ėθ can be written as:

ėθ = vT rc T -v tan(γ)/l b (27) = - KxK d d cos(e θ ) + KxKoK θ sin 2 (e θ ) cos(e θ ) • K d ey Ko cos(e θ ) + K θ tan(e θ ) = - Kx Ko cos 2 (e θ ) K d d + KoK θ sin 2 (e θ ) • [K d d sin(e θ ) + KoK θ sin(e θ )] = - Kx(K d d + KoK θ ) Ko cos 2 (e θ ) K d d + KoK θ sin 2 (e θ ) sin(e θ ) (28)
To solve the differential equation ( 28), let us introduce the following notations A = K d d, B = K o K θ and C = (A/B + 1). Hence, the solution of (28) has the following form:

ln tan( e θ 2 ) C + cos(e θ ) C -cos(e θ ) C/2 e θ e θo = - KxAB Ko C 2 t t 0 ( 29
)
The objective is to compute the time t f necessary to obtain e θ value which is less than a given orientation error threshold E angle . Therefore, (29) can be written as:

e θ = 2f θ (t, K, e θo ) (30)
where:

f θ = tan e θo 2 (C + cos(e θo ))(C -1) (C -cos(e θo ))(C + 1) C/2 e -KxAB Ko C 2 t f (31)
Indeed, for a specific parameter value K, if the controller can reach the target (with e θ less than E angle for a time t f given according to (30)), then for the same K parameters, e θ must be less than E angle for a time t > t f . From [START_REF] Vilca | Adaptive leaderfollower formation in cluttered environment using dynamic target reconfiguration[END_REF], the orientation error e θ depends on the initial orientation error e θo and the controller parameters (K d , K x , K o , K θ ). Moreover, K d can be chosen as a function of the initial distance as K d = 1/d i , which allows us to obtain A ≤ 1. The controller parameters are designed to obtain a fast convergence rate, given by

t f = f -1 θ (E angle , K, e θo ) (32)
while taking into account vehicle constraints such as maximum vehicle velocity v max and minimum radius of curvature r cmin . Secondly, the time t f enables the maximum distance covered by the vehicle to be ascertained. With this aim, we can set e RT = 0 and e θ = 0 (straight line to the target). Hence, (15) can be written as

v b = K x K d d (33)
Introducing ( 9), ( 10), ( 13) and [START_REF] Vilca | Robust on-line obstacle detection using range data for reactive navigation[END_REF] in the derivative of the distance, we obtain: ḋ =(e x ėx + e y ėy )/d

= -e x v b = -K x K d d d di 1/d ∂d = t 0 -K x K d ∂t ⇒ d = d i e -KxK d t (34)
From (34), the convergence of the distance depends on K x and K d . Therefore, they are designed to obtain d i such that:

d i = E dis e KxK d t f (35)
while taking into account vehicle constraints v max and r cmin . Simulations (cf. Subsection 5.1.1) will validate the above approach. Moreover, (30) and (35) show the relations between initial configuration, controller parameters and error of reaching the target. Therefore, for certain initial configuration and defined error bounds, the controller parameters can be obtained.

Navigation Strategy Based on Sequential Target Assignment

The computation of a time-parameterized path while taking into account different vehicle constraints and environment characteristics is timeconsuming [START_REF] Choset | Principles of Robot Motion: Theory, Algorithms, and Implementation[END_REF], [START_REF] Kuwata | Motion planning for urban driving using rrt[END_REF]. Different algorithms that compute a safe path (without temporal reference) [START_REF] Avanzini | On-line reference trajectory generation for manually convoying a platoon of automatic urban vehicles[END_REF], [START_REF] Lavalle | Planning Algorithms[END_REF] require less computational time but provide trajectories which do not ensure the safe navigation of the vehicle. In the previous section, we showed that the proposed control law guarantees that the static or dynamic target will be reached. In the following, a strategy to use a finite set of targets (waypoints) to define the desired vehicle route is presented.

Sequential target assignment

The proposed strategy uses a sequence of sorted waypoints disposed in the environment. The aim of this sequence is to guarantee safe (cf. Subsection 3.4 and Subsection 4.2) and flexible navigation. An elementary method to select adequate waypoints (target set-points (x Tj , y Tj , θ Tj , v Tj )) to perform a safe vehicle navigation in structured environment is presented in Subsection 4.2. Each waypoint T j (x Tj , y Tj ) is defined by the following parameters: D j the euclidean distance between the last waypoint T j-1 (x Tj-1 , y Tj-1 ) and the current waypoint T j ; θ Tj is the orientation between T j and T j+1 (x Tj+1 , y Tj+1 ):

θ Tj = arctan (y Tj+1 -y Tj )/(x Tj+1 -x Tj ) (36)
The safety between the specified waypoints is guaranteed by imposing E dist and E angle (cf. Subsection 3.4). This allows to guide the vehicle when it reaches the target T j (cf. Fig. 8) in the appropriate conditions to reach the next target T j+1 while also guaranteeing safe navigation between T j and T j+1 .

The strategy to assign the target point is shown in Algorithm 1. The parameters of the control law (cf. Subsection 3.2) enable the vehicle to reach the next target point (cf. Section 5) while ensuring that the vehicle trajectory is always within the road limits (cf. Fig. 8). The error conditions (E d and E angle ) are used to allow switching to the next target, when the vehicle position enters a circle with a radius equal to E dis and center (x Tj , y Tj ). The current target is updated with the following waypoint in the list and the vehicle starts the movement to reach the new target. If the vehicle does not satisfy the error conditions then the perpendicular line L j to the line which connects T j and T j+1 is used to Switch from the current target Tj to the next sequential waypoint Tj+1 3: end if switch to the next target when the vehicle crosses this line (cf. Fig. 8).

It can be noticed that the proposed control law, ( 13) and ( 14), is well suited to the navigation task. It allows to reach each target with an assigned velocity ( 13) (which can be different from zero contrary to [START_REF] Aicardi | Closed loop steering of unicycle like vehicles via lyapunov techniques[END_REF]). The obtained robot movements become smoother and thus appropriate for public transportation.

Waypoint selection method

Waypoint selection consists in obtaining the minimum number of points (waypoints) on the road to be successively reached by the vehicle to perform safe navigation. These waypoints are selected considering a safe position on the road (as far as possible from the road limits) and the reliability of the obtained vehicle trajectory (smooth changes between the successive points).

In this paper, the proposed waypoint selection (cf. Algorithm 2) is maximally simplified to focus only on the navigation strategy, i.e., target-reaching control (cf. Section 3) and sequential target assignment (cf. Subsection 4.1). Future studies will focus on the optimization of the waypoints selection in any kind of environment. Therefore, to provide a complete framework to achieve the navigation strategy, in this paper we will consider in the subsequent discussion, the existence of a defined trajectory (infinite number of points); the aim of the method is to select an appropriate number of points (waypoints).

The reference path can be obtained by different algorithms [START_REF] Lavalle | Planning Algorithms[END_REF] or by using a recorded vehicle trajectory. Different criteria can be considered to obtain the minimum number of straight lines that closely fit the reference path. Criteria such as the euclidean or curvilinear distance, orientation or radius of curvature between waypoints can be used to fix the desired waypoints on the path. The discretized reference path r is composed of sorted po- sition r i = (x ri , y ri ) and its tangent orientation θ ri . The minimum number of straight line segments over the defined path is then computed while considering a constant threshold ∆ αmax for the orientation variation of the path ∆ α (cf. Algorithm 2).

Fig. 9 shows one vehicle trajectory and the obtained waypoints using Algorithm 2 with ∆ αmax = 5 • , 15 • and 30 • respectively. Obviously, the switch between waypoints is smoother with a small value of ∆ αmax .

Algorithm 2 Waypoint selection based on existing reference path

Require: Reference path r = (xr, y r ) and ∆α max ∈ R + Ensure: Set of waypoints Sp 1: Init j = 0, rw j = r0 (initial position of r) and θw j = θr 0 (tangent of the point along trajectory r) 2: for ri ∈ r (sorted set of trajectory points) do Set rw j = ri and θw j = θr i

7:

Add wj(rw j , θw j ) to Sp 8: end if 9: end for

Proposal validation

This section presents a set of experiments to demonstrate the efficiency of the control law for target reaching and autonomous navigation in a structured environment. Section 5.1 provides simulation results to show the validity of our proposal. Section 5.2 discusses experimental results applied to an urban electric vehicle.

Simulation results

In these simulations different aspects, such as the stability of the control law to reach a static target, the flexibility of the navigation strategy based on target assignment and the performance of the proposed control law compared with other approaches found in the literature, will be presented. The physical parameters of the urban vehicle VIPALAB (cf. Fig. 19) modeled using tricycle kinematics (1) were considered.

Target Reaching

The first simulation shows the performance (safety, smoothness and convergence) of the control law to reach a desired final configuration (pose and velocity). For each simulation, the vehicle starts at the same position but with different initial orientations.

This simulation validates the analysis presented in Subsection 3.4, where the minimum d i , obtained for a limit vehicle configuration e θ ≈ π/2, allows to satisfy the bound of the errors for other initial configuration. The desired final configuration is (x T , y T , θ T ) ≡ (15, 4, 0 • ) and v T = 1 m/s. The controller parameters K are designed considering E dist ≤ 0.1 m and E angle ≤ 5 • using (35) and [START_REF] Vilca | Adaptive leaderfollower formation in cluttered environment using dynamic target reconfiguration[END_REF]. Hence, it is designed such that K = (1/d i , 0.6, 10, 0.1, 0.3, 0.01), d i = 10.6 m is the minimum initial euclidean distance to the target. These parameters were chosen to obtain a safe and smooth trajectory, faster response (t f ≈ 10.5 s (32)) and velocity value within the limits of the vehicle (cf. Subsection 3.4), which are v max = 1.5 m/s and the minimum radius of curvature r cmin = 3.8 m. Fig. 10 shows that the convergence of the system depends on the initial orientation error. Fig. 11 shows that the system errors are bounded (30) (black line) and converge to zero (cf. Subsection -fÇ (t, K, eÇo) fÇ (t, K, eÇo ) 3.4). Furthermore, the Lyapunov function shows asymptotic stability (cf. Fig. 12). Fig. 12(a) shows the three terms of the Lyapunov function [START_REF] Fang | Trajectory tracking control of farm vehicles in presence of sliding[END_REF] where the first term is 0.5K d d 2 , second term is 0.5K l d 2 l and the third term is K o [1 -cos(e θ )]. These figures show that the vehicle satisfies the constraints (velocity, acceleration and steering) presented in Subsection 3.4.

Moreover, the case of successive target reaching through a set of waypoints is analyzed. Two set of waypoints selected from a reference trajectory are used, one set has a distance between waypoints equal to 2 m and the other equal to 4 m (cf. Fig. 13). Fig. 13 and 14 respectively, show the vehicle trajectories and lateral and angular errors w.r.t. the reference trajectory for two set of waypoints. It can be noted that the obtained vehicle trajectories are close enough from the reference trajectory; and as expected, the lateral and angular errors are smaller when the fixed distance between the waypoints decreases. Therefore, the proposed navigation strategy and control law permits to the vehicle to perform accurate trajectory tracking behavior if the waypoints are close enough.

Comparison study

In this simulation, two common approaches to follow a dynamic target by a tricycle are briefly presented and compared with the proposed control law. The comparison is focused on the position and orientation errors and the convergence time. In the sequel, the dynamic target is assumed to be another identical urban vehicle. These two strategies are briefly described below:

1. Approach based on a reference path: In [START_REF] Daviet | Platooning for small public urban vehicles[END_REF], a method for following a vehicle based on a Frenet frame was developed (cf. Fig. 15). It exploits the use of chained systems to separate the lateral and longitudinal control. Therefore, each controller can be designed independently. The lateral control is obtained using chained transformation (more details are given in [START_REF] Samson | Control of chained systems. application to path following and time-varying point-stabilization of mobile robots[END_REF] and [START_REF] Laumond | La Robotique Mobile[END_REF]). The longitudinal control consists in keeping a specific curvilinear distance d s between the target and the vehicle. One drawback of this approach is the dependency on a known reference path for the vehicle, i.e., if the vehicle follows a dynamic target then the target trajectory must be accurately known by the vehicle. 2. Approach based on a target model: In [START_REF] Laumond | La Robotique Mobile[END_REF] and [START_REF]Springer Handbook of Robotics, Part E-34[END_REF] a control law to track a reference vehicle (target) is proposed. A variable transformation to obtain the control is applied to the error system and commands. The control law is synthesized using a suitable Lyapunov function (more details can be found in [START_REF] Laumond | La Robotique Mobile[END_REF]).

The desired steering angle is computed by integration. Nevertheless, the control law considers a non-zero linear target velocity, i.e, if the target is static then the commands sent to the vehicle are zero. The approaches presented above are implemented in simulation and compared to our proposal. In order to do this, the main target (equivalent to a Leader robot in formation control [START_REF] Vilca | Adaptive leaderfollower formation in cluttered environment using dynamic target reconfiguration[END_REF]) tracks a sinusoidal trajectory and the followers must maintain a distance of 3 m w.r.t. this first robot, i.e., the secondary target to be reached is located at 3 m (curvilinear distance) from the main target (cf. Fig. 16).

Fig. 16 and 17 show the trajectories and the control output of each vehicle (leader and followers). It can be noticed that the proposed control law has a similar performance to the controller based on a reference trajectory (Frenet control).

Table 1 shows the convergence time to satisfy the error threshold in distance d T arget and orientation e θ T arget with regards to the target pose (cf. Fig. 5). The proposed control law has the smallest convergence time to satisfy simultaneously both threshold errors (maximum value between convergence time of d T arget and e θ T arget ). The difference with the Frenet control is equal to 1.81 s and with chained system is equal 2.3 s.

Table 2 shows the convergence time to satisfy the error threshold in distance y i and orientation θ with regards to the target trajectory (Frenet reference frame X f Y f (cf. Fig. 15)). The Frenet control has the smallest convergence time to satisfy simultaneously both threshold errors, since it is dedicated to follow the reference trajectory. However, there is only a small difference (0.28 s) compared to our proposal while our method uses only the current pose of the target (thus more flexible).

The proposed control law was not designed to take into account the reference trajectory, however the obtained results are very close to those designed for trajectory tracking. In addition to its accuracy, the proposed control law has more flexibility (cf.

Section 3) to perform the autonomous navigation of the vehicles. Indeed, we need only to know the current pose and dynamics of the target instead of all recorded trajectory.

Experimental results

The navigation strategy was also experimented with a pair of real urban vehicles to perform dynamic target reaching. The scenario was built to show different situations, such as: multi-vehicle navigation in formation, static and dynamic target reaching and obstacle avoidance situation.

Testbed and scenario

Navigation was performed in a structured environment named PAVIN (Plate-forme d'Auvergne pour Véhicules INtelligents) (cf. Fig. 18). The proposed strategy was implemented using VIPALAB urban vehicles (cf. Fig. 19). This vehicle carry different embedded proprioceptive and exteroceptive sensors such as odometers, gyrometer, steering angle sensor and an RTK-GPS (more details are given in [START_REF]The Institut Pascal Data Sets[END_REF]). In these experiments, each vehicle uses a combination of RTK-GPS and gyrometer to estimate its current position and orientation at a sample time of T s = 0.1 s (cf. Table 3). The vehicles have a range sensor (LIDAR) with a maximum detected range equal to 10 m. These sensors provide enough accurate data w.r.t. the vehicle dynamic. Indeed, in these experiments, the vehicles move at maximum velocity of 1.5 m/s due mainly to the relative short dimensions of the used urban platform (cf. Fig. 18). Moreover, the vehicles communicate by WI-FI, enabling the transmission of the leader's pose data.

Results analysis

Experiments were carried out to show the performance of the proposed control law and target assignment strategy using waypoint selection based on Algorithm 2 (with ∆ αmax = 15 • ) on an already defined reference trajectory. The Leader vehicle has to reach successively static waypoints. Moreover, the proposed control law was implemented in another vehicle (Follower) which takes the first vehicle (Leader) as dynamic target to track at a curvilinear distance equal to 5 m (behind the Leader). The tracking of the dynamic target allows to apply the proposed control law to multi-robot systems where the dynamic set-point is given by the leader and the desired geometric formation shape [START_REF] Benzerrouk | Obstacle avoidance controller generating attainable set-points for the navigation of multi-robot system[END_REF]. The configuration of the dynamic target is sent by the Leader to the Follower via WI-FI. This experiment can be found online 1 .

Furthermore, to exhibit the flexibility of the proposed navigation strategy, a scenario with the presence of an obstacle is presented (cf. Fig. 19). An obstacle is placed between the waypoints. As mentioned in Section 2, the proposed strategy can easily integrate the obstacle avoidance behavior (cf.

1 http://maccs.univ-bpclermont.fr/uploads/Profiles/VilcaJM/Navigation.avi Fig. 3). Therefore, the vehicle can perform different maneuvers between waypoints, in this case the obstacle avoidance without the use of any trajectory replanning method. The used obstacle avoidance method is based on limit-cycles as given in [START_REF] Adouane | Mobile robot navigation in cluttered environment using reactive elliptic trajectories[END_REF], [START_REF] Vilca | Reactive navigation of a mobile robot using elliptic trajectories and effective online obstacle detection[END_REF] and [START_REF] Benzerrouk | Obstacle avoidance controller generating attainable set-points for the navigation of multi-robot system[END_REF]. It was selected because it is a stable and robust method which could use only local information from range sensors. Let us briefly inspect this method (for more details see [START_REF] Adouane | Mobile robot navigation in cluttered environment using reactive elliptic trajectories[END_REF] and [START_REF] Vilca | Reactive navigation of a mobile robot using elliptic trajectories and effective online obstacle detection[END_REF]). A limit-cycle is a reactive safe trajectory which encloses the hinder obstacle. According to that, the vehicle avoids the obstacle while tracking the direction of the limit-cycle trajectories. The obstacle avoidance is activated as soon as the vehicle detects at least one obstacle which can hinder the future vehicle movements toward the current assigned waypoint [START_REF] Vilca | Reactive navigation of a mobile robot using elliptic trajectories and effective online obstacle detection[END_REF] and [START_REF] Vilca | Robust on-line obstacle detection using range data for reactive navigation[END_REF].

It can be seen in Fig. 20(a) that the Leader reaches accurately the successive static waypoints successfully and the Follower tracks accurately also the dynamic target (Leader). Moreover, the Follower trajectory using the proposed control law is close to the leader trajectory (cf. Fig. 20(a)). Fig. 20(b) focuses on the vehicles' trajectories when the obstacle avoidance is activated. The Leader detects the hinder obstacle between the waypoints and it applies the reactive limit-cycle method [START_REF] Vilca | Reactive navigation of a mobile robot using elliptic trajectories and effective online obstacle detection[END_REF], [START_REF] Vilca | Robust on-line obstacle detection using range data for reactive navigation[END_REF]. The Follower avoids also the obstacle since it tracks accurately the Leader trajectory. It can be noted that the proposed navigation strategy allows flexible and smooth movements between the waypoints and also to perform different behaviors, such as: obstacle avoidance, emergency stop or waypoint re-

Conclusion

This paper has presented a novel safe and flexible control strategy based on target reaching for the navigation of autonomous vehicles in structured environment. A control law was presented and synthesized using a suitable Lyapunov function, which takes into account the position, the angle between the robot and the target, and its orientation with respect to set-points. Moreover, it enables static and dynamic target reaching. The stability of the overall control architecture was proved using a suitable Lyapunov function based on a new set of variables. A sequential target assignment strategy to perform safe navigation was also proposed. It is based on target switching using appropriate reference frames, linked to the current selected waypoint and to the next one in the list. This target assignment strategy enables a smooth and flexible vehicle trajectory while satisfying an upper bound on distance and orientation errors (to not collide with the boundary of the environment). An elementary waypoint selection method was also presented to perform safe and smooth trajectories. Simulations and experiments using real urban vehicles show the efficiency and the flexibility of the proposed control strategy for the navigation of urban vehicles. Furthermore, different target-reaching methods from the literature were presented and their performances were compared with the proposed control law. These comparisons showed the interesting features of the proposed control law in term of stability and flexibility for different tasks.

In future works, a robust analysis will be developed to determine the sensor inaccuracy effects on our control law. Different methods to obtain the optimal number of waypoints on the road map will be developed. In addition, this strategy will be applied to high dynamic multi-vehicle system based on the combination of virtual structure and Leader-Follower approaches.

Highlights:

 Novel control strategy based on static/dynamic target reaching for autonomous navigation.  Control law is synthesized using a Lyapunov function based on a new set of variables.  Interesting properties of the control law in term of stability and flexibility.  Sequential target assignment strategy allows performing safe navigation.  Experiments using actual vehicles and several simulations show the advantages of the proposal.

Figure 1 :

 1 Figure 1: Autonomous navigation of an electric vehicle in an urban environment (Clermont-Ferrand, France).

Figure 3 :

 3 Figure 3: Overview of the proposed navigation strategy.

Figure 4 :

 4 Figure 4: Vehicle and target configuration in Global (X G Y G ) and Local (XmYm) reference frames.

Figure 5 :

 5 Figure 5: Control variables according to Lyapunov synthesis.

V

  =vT sin(e θ )[K d ey -K l d sin(eRT ) cos(e θ )] + vT rc T [d 2 K l sin(eRT ) cos(eRT ) + Ko sin(e θ )] v b [K d ex + K l d sin(eRT ) sin(e θ ) + Ko sin(e θ )cc]

  K l d sin(e RT ) sin(e θ ) + K o sin(e θ )c c = 0 (24) Replacing (16) with r c T → ∞ in (24), the following expression is obtained: 0 =K d e x + K l d sin(e RT ) sin(e θ ) + tan(e θ )[K d e y -K l d sin(e RT ) cos(e θ )] + K o sin(e θ ) K θ tan(e θ ) + K RT sin 2 (e RT ) sin(e θ ) cos(e θ ) =K d [e x + e y tan(e θ )] + K o K θ sin 2 (e θ ) cos(e θ ) + K o K RT sin 2 (e RT ) cos(e θ )

Figure 6 :

 6 Figure 6: Limit vehicle configuration for tuning controller parameters.

Figure 7 :

 7 Figure 7: Block diagram of the analysis.

Figure 8 :

 8 Figure 8: Description of waypoints and target assignment.

Algorithm 1

 1 Sequential target assignment Require: Vehicle pose, current target Tj and a set of N sorted waypoints Ensure: Switch between target set-points 1: if ( (d ≤ E dis and e θ ≤ E angle ) or (x T j ≥ 0) ) { x T j is the coordinate of the vehicle in the local Target frame XT j YT j (cf. Fig. 8) } then 2:

Figure 9 :

 9 Figure 9: Example of waypoint selection based on a reference path and Algorithm 2.

3 : 4 :

 34 Compute ∆α = |θriθw j | if ∆α ≥ ∆α max then

  Fig. 10 shows the trajectory of the vehicle for different initial configurations. The orientation errors are shown in Fig. 11. The Lyapunov function values are shown in Fig. 12.

Figure 10 :

 10 Figure 10: Trajectories of the vehicle for several initial orientations.

Figure 11 :Figure 12 :

 1112 Figure 11: Orientation errors (e θ ) for several initial vehicle orientations.

Figure 13 :

 13 Figure 13: Vehicle trajectories for different distances between waypoints.

Figure 14 :

 14 Figure 14: Errors w.r.t. reference trajectory for different distances between waypoints.

Figure 15 :

 15 Figure 15: Vehicle modeling in a Frenet frame.

Figure 16 :

 16 Figure16: Trajectory of the leader and followers (proposed control law, Frenet control[START_REF] Samson | Control of chained systems. application to path following and time-varying point-stabilization of mobile robots[END_REF] and Chained system[START_REF] Laumond | La Robotique Mobile[END_REF]).

Figure 17 :

 17 Figure 17: Control output.

Figure 18 :

 18 Figure 18: PAVIN experimental platform. (a) to (f) correspond to the locations shown in Fig. 19.

Figure 19 :

 19 Figure 19: Some images from the performed experiment.

  (a) Vehicle trajectories. (b) Obstacle avoidance.

Figure 20 :

 20 Figure 20: Vehicle trajectories obtained using GPS and a set of waypoints positioned in the environment using Algorithm 2 (∆α max = 15 • ).

  assignment.

Fig. 21

 21 Fig.21shows the velocity and steering angle of the vehicles. These actual values have been filtered, during the experimentation, using an Extended Kalman Filter (EKF) to reduce the sensor noise. Fig.22and 23 show the Lyapunov function values which highlight that each vehicle is stable and it converges to each static waypoint for the Leader and to the dynamic target for the Follower. Therefore, smooth, flexible and safe trajectories for the vehicles were obtained.

Figure 21 :

 21 Figure 21: Control output (real experiment).

Figure 22 :

 22 Figure 22: Lyapunov function of the leader (based on static waypoint reaching).

Figure 23 :

 23 Figure 23: Lyapunov function of the follower (based on dynamic target tracking).

Table 1 :

 1 Comparison with dynamic target pose (Fig. 5) Time [s] to always keep: d T arget < 0.15 m |e θ T arget | < 5 •

	Proposed control Frenet control Chained control	13.17 s 14.98 s 15.47 s	4.24 s 3.92 s 6.62 s

Table 2 :

 2 Comparison of the errors defined according to Frenet reference frame (Fig. 15) Time [s] to always keep: |y i | < 0.15 m | θ| < 5 •

	Proposed control Frenet control	3.33 s 2.93 s	4.14 s 3.86 s

Chained control 10.22 s 10.57 s

Table 3 :

 3 VIPALAB platform (cf. Fig.19)

Elements Description

Chassis (l, w, h)= (1.96, 1.30, 2.11) m Computer Intel Core i7, CPU:1.73 GHz RAM:8Go OS(32bits):Ubuntu12.04 RTK-GPS NacTechGPS, accuracy: 2 cm framerate: 10 Hz Gyrometer Xsens MTi, accuracy: 0.2 • /s framerate: 2 KHz Proprioceptive sensors Steering angle, resolution: 0.02 • framerate: 50 Hz Linear speed, resolution: 0.1 m/s framerate: 50 Hz

$ Supported by the French National Research Agency (ANR) through the SafePlatoon project.


M Aharrouche

The Atlas Liquid Argon Calorimeter: Integration, Installation and Commissioning

The construction and the installation of the ATLAS liquid argon calorimeters are now completed. In this paper we will briefly review the design of the calorimeters and describe the integration and installation steps. The result of the calorimeters' commissioning on the surface and in the cavern, and the results of the combined test beam in 2004 are also shown.

Introduction

ATLAS (A Toroidal LHC Apparatus) is one of the two general-purpose detectors which will operate at the LHC. Its conception is very constrained by the complexity of the hadronic interactions and the high center of mass energy. The optimisation of the detector is driven by the requirement to detect new particles, in particular the Higgs boson in its decay to two photons (H → γγ). ATLAS consists of three main sub-detectors: the inner tracking detector, the muon spectrometer and the calorimeter bathing in the magnetic field of two magnetic toroidal and solenoidal systems. The calorimetry technique in ATLAS is a sampling calorimetry with two detectors types: liquid argon and scintillating tiles.

The construction of the Liquid Argon (LAr) calorimeter system (see Fig. 1) was completed since mid-2004. The evaluation of the production quality and the validation of the performance needed multiple tests using beam at CERN; the most recent one is the combined run in 2004 where a full slice of ATLAS has been tested. Calorimeter modules were constructed in several institutes, assembled into wheels at CERN in 2003-2004 and integrated into four cryostats: two half-barrels and two end-cap cryostats. End 2005, all calorimeter sub-detectors were installed in the AT-LAS cavern. Since autumn 2006, a cosmic run based on the EM barrel together with some hadronic tile calorimeter modules is on-going.

In this paper we describe briefly the ATLAS LAr calorimeter and present some results from the combined test beam. The status of the LAr calorimeter integration is presented. Finally, we show some results of the calorimeter commissioning on the surface and in the cavern.

The ATLAS Liquid Argon Calorimeter

The liquid argon has been chosen for ATLAS because of its intrinsic linear behaviour, its stability of the response in time and its radiation tolerance. The ATLAS LAr calorimeter [START_REF] Aharrouche | ATLAS Liquid Argon Calorimeter Technical Design Report[END_REF][START_REF] Aubert | [END_REF] system consists of an electromagnetic barrel calorimeter and two end-cap cryostats with electromagnetic (EMEC), hadronic (HEC) and forward (FCal) calorimeters.

The electromagnetic calorimeter is an accordion-shaped lead/liquid argon detector comprising a barrel and two endcaps. The accordion geometry allows a very good hermeticity (without cracks) and minimizes inductances in the signal path.

The barrel, covering a pseudo-rapidity range of |η| < 1.475, shares its cryostat with the superconducting solenoid, the calorimeter being behind the solenoid. Each wheel consists of 16 modules of 1024 lead absorbers interleaved with read-out electrodes. Each end-cap, covering the pseudo-rapidity range 1.4 < |η| < 3.2, consists of 8 modules (see Fig. 2 (right)). To correct for the energy lost in the dead material in front of the calorimeter, both endcaps and barrel are completed with pre-sampler detectors, covering the |η| < 1.8 range. These pre-samplers are thin layers of argon equipped with read-out electrodes, but without absorber. Figure 2 (left) shows the segmentation in the (η,φ) plane and along the longitudinal EM shower development in the barrel.

The design of the EM calorimeter was driven by the search of the Higgs, decaying to two photons and to four electrons, and the aim to reach a precision on its mass of 1%. This requires a sampling term of 10%/ √ E associated with a constant term of better than 0.7%.

The HEC is structured in two wheels: the front HEC1 and rear HEC2 with 32 modules each, covering the pseudorapidity range 1.5< |η| <3.2. It consists of copper-plate absorbers of 25 mm (50 mm) thickness for the HEC1 (HEC2) interleaved with electrodes. The read-out technique (Electrostatic Transformer) of the HEC splits the gap between two absorbers into four sub-gaps, each with a width of 1.85 mm (see Fig. 3). The read-out granularity is ∆η×∆φ= 0.1×0.1 for 1.5< |η| <2.5 and 0.2×0.2 for 2.5< |η| <3.2. The HEC is designed to measure the jets in the forward direction.

The FCal shares the same cryostat with the EMEC and HEC and covers the pseudo-rapidity range 3.2< |η| <4.9. It consists of three consecutive modules along the beam line. The modules consist of a tungsten matrix in FCal2 and FCal3 (copper matrix in FCal1) housing cylindrical electrodes (see Fig. 4), consisting of copper rods. The LAr gaps are very small: 0.250, 0.375 and 0.500 mm for FCal1, FCal2 and FCal3, respectively. The read-out granularity is ∆η×∆φ ≈ 0.2×0.2. The design of the FCal was constrained by the high radiation level in the very forward region.

Combined Test Beam

During the 2004 summer a full slice of the ATLAS central detector has been exposed to beams of electrons, photons, pions, muons and protons in the energy range 1 GeV≤ E ≤ 350 GeV at CERN in the H8 beam line. This slice contains one module of the electromagnetic barrel. The electron or photon reconstruction in the calorimeter consists of first forming a cluster around the electron or photon shower, and then summing the cells' energies within the cluster. However, the shower is not fully contained in the cluster due to the energy lost in the dead material, longitudinal and lateral leakage of the shower outside the cluster. In order to recover the energy lost and the energy not contained in the cluster, while obtaining good results on energy resolution and linearity, the following calibration ansatz [3,4] (equation 3) was developed on data obtained from a Geant4 simulation of the beam test setup. The formula was applied to the real data taken during the combined run. It parameterizes the dependence of the energy lost as a function of the energy deposit in each layer of the calorimeter:

E particle = offset + W 0 E 0 + W 01 E 0 E 1 + +λ 3 i=1 E i + W 3 E 3 (1)
where E i is the energy deposit in the layer i of the EM calorimeter and W i , of f set and λ are the calibration weights.

The calorimeter performances were evaluated by applying this calibration scheme to the electron data, and excellent results were obtained: linearity of better than 0.2%, energy resolution with the sampling term of 10.6% GeV 0.5 EM Barrel 49(0.04) 1(0.01) 31(0.03) 8.5(1.9) 0(0) Table 1 Results of the calorimeter tests at LAr temperature, performed in the surface test area. A is the number of bad channels, B is the number of bad calibration lines, C is the number of dead signal channels, D is the number of HV corrected channels and E is the number of dead HV sectors. In each column, the fraction in %, relative to the total number of the same objects in the full calorimeter part, is indicated in parentheses. and a local constant term of 0.43% (see Fig. 5 and Fig. 6). These performances meet the ATLAS requirements.

LAr Calorimeter Integration

After their production [START_REF] Aubert | [END_REF] in the institutes of the collaboration between 2001 and 2004, the LAr calorimeter modules were delivered to CERN where they have been assembled into wheels in clean rooms. Then, each calorimeter wheel was inserted into a cryostat (see Fig. 7). Each step was followed by control measurements to check the functionality and connectivity of all the channels. Assembly and integration of the complete LAr calorimeter was finalized from 2002 to 2004.

Calorimeter Tests at LAr temperature

A full test of the calorimeter cryostats at LAr temperature (87 K) was done during a total period of 24 weeks in the surface area, where the calorimeters' integration took place. The important tests consisted, first, of the checking of connectivity and integrity of all read-out channels by pulsing them via the calibration lines. To ensure the integrity of all HV channels, a test was performed at room temperature and after filling with LAr: all channels were kept at nominal voltages during several weeks to check the stability. The table 1 shows the test results of the 190304 read-out channels, the 14592 calibration lines and the 4248 HV channels.

Installation and Commissioning

All LAr calorimeters' cryostats have already been installed, the barrel cryostat was lowered to the experimental area in October 2004 (see Fig. 7), the end-cap cryostats in December 2005. In November 2005 the barrel was moved to its final position around the interaction point. The installation and connection of the front-end electronics was finished for the barrel and the end-caps. The HV system for the barrel and the end caps has been completed and is operational. The back-end electronics read-out is already functional for the barrel and is partially functional for the end-caps.

The commissioning of the calorimeters is organized in three phases. The phases 1 and 2 are technical, they consist of the installation and testing of the read-out channels and the integration of the calorimeter data acquisition with the ATLAS TDAQ system. The phase 3 commissioning involves the cold tests of the detector and the checking of the performance in noise studies and calibration runs. The data taking with cosmic rays is the main part of the phase 3.

The pedestals, random and coherent noise, the autocorrelation matrix, gains and pulse shapes of all channels are being monitored regularly. Two run periods of data taking with cosmics were done in August and October 2006 with 15% of the EM Barrel modules and some tile calorimeter modules to trigger on cosmic muons. The aim of these runs was to check the uniformity of the cell response at a 1% level. The first analyses of those events have started and give already valuable results. Figure 8 shows the uniformity response of the EM barrel calorimeter, along η, to the cosmic muons events.

Conclusions

The ATLAS LAr calorimeter has been constructed in modules, then assembled into wheels, integrated into four cryostats and installed in the ATLAS cavern during the last five years. The most recent combined test beam, in a configuration very close to the final ATLAS experiment, shows good results meeting the ATLAS requirements. The cold commissioning in the surface hall demonstrated the excellent condition of those calorimeters. Electronics is now almost fully installed and is being integrated in the ATLAS global data acquisition system. The ongoing cosmic muon runs will allow us to prepare the detector to tackle the first LHC data as efficiently as possible.

Fig. 1 .

 1 Fig. 1. The ATLAS Liquid Argon Calorimeter. 1: the EM barrel, 2: the EM end cap, 3: the hadronic end cap, 4: the forward calorimeter.

Fig. 3 .

 3 Fig. 3. Schematic view of the HEC module (left), and the gap between the two absorber plates.

Fig. 4 .

 4 Fig. 4. The face of the electromagnetic FCal module (right).

Fig. 2 .Fig. 5 .Fig. 6 .

 256 Fig. 2. Left: projective view of the EM module showing the longitudinal and lateral segmentation. Right: schematic view of the EM end cap wheel.

Fig. 7 .

 7 Fig. 7. Top left: 25 modules of the HEC wheel. Top middle: EM half-barrel in its cryostat. Bottom left: The HEC integration in its cryostat(behind the EM end cap). Bottom middle: insertion of the FCal in the cryostat. Right: lowering of the EM barrel cryostat to the ATLAS cavern.

Fig. 8 .

 8 Fig. 8. Uniformity response to muons versus η compared to the cell depth.

Preprint submitted to Elsevier


Francisco Chinesta

Antonio Falcó

Mariano González

MODEL REDUCTION METHODS IN OPTION PRICING

Keywords: Model Reduction, Proper Orthogonal Decomposition, Finite Difference Schemes, Crank-Nicolson Scheme

In this work we introduce the Proper Orthogonal Decomposition (POD) approach to the valuation of contingent claims for one-dimensional price models.

First, we present the POD in the context of an abstract Hilbert space and we give an application for the numerical pricing of Double Barrier Options. In a finite dimension setting, we show the model reduction method for Finite Difference schemes of implicit type. In particular, we construct the reduced version of the Crank-Nicolson scheme and some numerical examples are given.

Introduction

The aim of this paper is to introduce the Proper Orthogonal Decomposition (POD) for model reduction in the framework of Partial Differential Equations arising in Option Pricing Models. The Proper Orthogonal Decomposition is a powerful and elegant method for deriving low order models of dynamical systems. The POD provides a basis for the modal decomposition of an ensamble of functions, composed of theoretical, experimental o computer data. Its properties suggest that it is the preferred basis to use in various applications.

It was successfully used in different fields including signal analysis and pattern recognition (see for example [START_REF] Fukunaga | Introduction to Statistical Recognition[END_REF]), fluid dynamics (see [START_REF] Berkooz | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF] and [START_REF] Sirovich | Turbulence and the dynamics of coherent structures, parts Quart[END_REF]) and more recently in control theory (see [START_REF] Atwell | Reduced Order Controllers for Spatially Distributed Systems via Proper Orthogonal Decomposition[END_REF]). The process is also known as Principal Component Analysis (PCA) (see [START_REF] Fukunaga | Introduction to Statistical Recognition[END_REF] ) or the Karhunen-Loève expansion the latter name comes from two individual papers by Karhunen [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlinchkeitsrechnung[END_REF] and Loève [START_REF] Loève | Functions aleatoires de second ordre[END_REF].

The most striking feature of the POD is its optimality. It provides the most efficient way of capturing the dominant components of an infinite-dimensional process with only a finite number of "modes" and often there are surprinsigly few. The main idea is to find a set of ordered orthonormal basis vector in a subspace associated to a set of data. This set of data can be composed of theoretical, experimental, or computed data. The optimal ordering of the basis is such that the first basis element best represents the data in the input collection, the second basis element is the next best, and so forth. For the purposes of model reduction, one type of input collection is a set of time snapshots, where each snapshot contains spatial data obtained from a numerical simulation at a fixed time. These snapshots are chosen so that the basis reflects the system dynamics.

In the present work we examine the feasibility and efficiency of the POD method for pricing numerically an European type derivative. To this end we introduce the reduced finite difference scheme of implicit type.

The paper is organized as follows. In Section 2 we introduce the POD in the context of a general Hilbert space and we give an example of application to the valuation of Double Barrier Options by using a set of data composed by data of theoretical type. In Section 3, by using the POD in the context of a finite Hilbert space, we introduce the reduced finite difference scheme of implicit type. We applied it by using the reduced version of the Crank-Nicolson scheme to study, at the numerical level, the valuation of the European Put in the Black-Scholes model. Finally, in 4 we give some concluding remarks.

Model Reduction with Proper Orthogonal Decomposition

The tools of POD and Galerkin method provides a systematic way for producing reduced-order models from data. In this section, we give an overview of these methods, first we introduce the proper orthogonal decomposition in the context of an abstract Hilbert space. The abstract setting is used due to the fact that, in this paper we will consider different inner products, and keeping the exposition general allows one to see precisely where the dependence on the inner products lies. Next, we introduce the POD in the context of pricing Double Barrier Options

The Proper Orthogonal Decomposition

Let H be a real Hilbert space with inner product ( and we refer to V snap as ensamble consisting of the snapshots, at least one of which is assumed to be non-zero. Assume that dim V snap = p and let {Ψ 1 , Ψ 2 , . . . , Ψ p } denote an orthonormal basis of V snap . Then each member of the ensamble can be expressed as

y j = p k=1 (y j , Ψ k ) H Ψ k (1)
for j = 1, 2, . . . , n. The proper orthogonal decomposition consists in choosing the orthonormal basis such that for every l ∈ {1, 2, . . . , p} the mean square error between the elements of V snap and the corresponding l-th partial sum of ( 1) is minimized on averaged, that is,

min {Ψ 1 ,Ψ 2 ,...,Ψ l } 1 n n j=1 y j -l k=1 (y j , Ψ k ) H Ψ k H subject to (Ψ i , Ψ j ) H = δ i,j for 1 ≤ i ≤ l, 1 ≤ j ≤ i. (2)
A solution {Ψ 1 , Ψ 2 , . . . , Ψ l } to (2) is called a POD basis of rank l. The solution of ( 2) is characterized by the first order necessary conditions of optimality, which can be written as an eigenvalue problem. To see this let us define the bounded linear operator L n : R n -→ H by

L n [u 1 u 2 • • • u n ] T = n i=1 u i y i .
Then the adjoint L * n : H -→ R n is given by

L * n (v) = [(v, y 1 ) H • • • (v, y n ) H ] T . Let us next define R n = L n L * n and K n = L * n L n , given by R n (v) = n i=1 (v, y i ) H y i and (K n ) i,j = (y j , y i ) H ,
respectively. It is now easy to show that the solution to (2) is characterized by the following first order necessary condition of optimality

R n Ψ = λΨ,
where R n is a linear self-adjoint compact operator. Therefore, by the Hilbert-Schmidt theory, there exists a complete orthonormal basis {Ψ i } ∞ i=1 so that

R n Ψ i = λ i Ψ i ,
and λ 1 ≥ λ 2 ≥ • • • and λ i = 0 for i > k. Moreover, we have the following error formula

n j=1 y j - p i=1 (y j , Ψ i ) H Ψ i H = k i=p+1 λ i .
In practical applications the POD basis of rank l may be computed as follows:

Step 1 Solve the eigenvalue problem

K n w = λw
where (K n ) i,j = (y j , y i ) X and the nonnegative eigenvalues satisfy

λ 1 ≥ λ 2 ≥ • • • ≥ λ k > 0.
Step 2 Find

Ψ i = 1 √ λ i n j=1
w j i y j where w j i denotes the j-th component of the eigenfunction w i , for i = 1, 2, . . . , k.

The optimality properties of POD state among all linear combinations the one which correspond to POD is the best in the sense that it will capture the most information possible in the average sense. Since the eigenvalues can be used to find how close the low-dimensional subspace approximant is to the date set, one can seek the number of modes l such that such that the fraction

l i=1 λ i / n i=1 λ i
is close to one and yet l << n. The l-dimensional reduced basis subspace is defined as

V POD = span{Ψ 1 , Ψ 2 , . . . , Ψ l }.

A case study: the Double Barrier Option

Now we would to use the POD method to valuate an European style derivative. Assume that the payoff of this derivative at expiration time T, is given by a function f of a single state variable X, which follows a scalar diffusion process under the risk neutral measure:

dX t = b(X t )dt + a(X t )dW t , X 0 = x. (3)
Suppose that the derivative contract has a double barrier provision given by 0 < L < U < ∞, such that if at any time between the contract inception and expiration time either the lower barrier L or the upper one U is reached, the contract is canceled. It is well-known that the value u of the derivative security, solves the Partial Differential Equation (PDE):

-u t (x, t) + A(u(t, x)) = 0, ( 4
)
where A is the negative of the infinitesimal generator of X given by

A(φ(x)) = - 1 2 a 2 (x)φ xx (x) -b(x)φ x (x) + r(x)φ(x).
and with boundary conditions u(x, T ) = f (x) and u(L, t) = u(U, t) = 0. Our main goal is to construct an approximation û(x, t) of the true solution u(x, t), by using the Proper Orthogonal Decomposition.

Consider the speed density of the diffusion process X given by

M(x) = 2 a 2 (x)S(x) (5)
where

S(x) = exp - x • 2b(y) a 2 (y) , (6)
is the scale density of X. Let H = L 2 ([L, U ], M) be the Hilbert space of square integrable functions on (L, U ) with the speed density M and endowed with the inner product

(g(x), h(x)) H = U L g(x)h(x)M(x)dx, (7)
(see Proposition 1 in [START_REF] Davidov | Pricing Options on scalar diffusions: An eigenfunction expansion approach[END_REF]).

Let u(x, t) be a flow defined on [L, U ] × [0, T ]. For now, imagine a time average over an ensamble with members {u (k) (x)} = {u(x, t k )}, obtained from successive measurements during a single run for

0 = t 1 < t 2 < • • • < t n = T,
a grid in the trading interval [0, T ], represents the "snapshots" or a sampled data representation of the flow. Assume that we obtain as snapshots the following finite sequence of functions

u(x, t k ) = u (k) (x) = k j=1 c j e -µ j (T -t j ) φ j (x) (8)
for k = 1, 2, . . . , n, where

{φ j (x)} ∞ j=1 is a complete orthonormal basis of L 2 ([L, U ], M), of eigenvectors with associated eigenvalues 0 < µ 1 < µ 2 < • • • < µ j < • • • obtained as solution of the regular Sturm-Liouville boundary value problem A(φ(x)) = µφ(x), φ(U ) = φ(L) = 0, (9)
and c j = (f (x), φ j (x)) H for j = 1, 2, . . . , n. Then, it follows that

K n = L Kn L T Kn where L Kn =      c 1 e -µ 1 (T -t 1 ) 0 0 . . . 0 c 1 e -µ 1 (T -t 2 ) c 2 e -µ 2 (T -t 2 ) 0 . . . 0 . . . . . . . . . . . . . . . c 1 e -µ 1 (T -tn) c 2 e -µ 2 (T -tn) c 3 e -µ 3 (T -tn) • • • c n e -µn(T -tn)      .
We denote by

w i = [w 1 i , w 2 i , . . . , w n i ]
T the i-th eigenvector of K n with associate eigenvalue equal to λ i for i = 1, 2, . . . , l and where l is choosed by means the condition

l i=1 λ i / n i=1 λ i > 0.99. Recall that A(φ i (x)) = µ i φ i (x)
for i = 1, 2, . . . . Thus the modal basis is given by

ψ i (x) = 1 √ µ i n j=1 w j i u (j) (x)(x), for i = 1, 2, . . . , l.
The approximate solution û of ( 4) can be expanded in terms of these POD basis functions as

û(x, t) = l i=1 a i (t)ψ i (x) (10)
and forms the Galerkin expansion. We derive the reduced order model by employing the Galerkin weighted residual discretization of the model ( 4). The model residual can be expressed as

R(û(x, t)) = -û t (x, t) + A(û(x, t))
Applying the Galerkin projection which enforces the residual to be orthogonal to each basis functions (R(û(x, t)), ψ i (x)) H = 0 (11) for i = 1, 2, . . . , l leads to the reduced order model of l ordinary differential equations for the amplitude coefficients

da i (t) dt = l j=1 a j (t) (A(ψ j (x)), ψ i (x)) H ( 12
) for i = 1, 2, . . . , l. Recall that (ψ i (x), ψ j (x)) H = δ i,j .
Then it is easy to see that

(A(ψ i (x)), ψ j (x)) H = w T i √ µ i Z w j √ µ j ,
where Z is an n × n matrix satisfying that

Z i,j = (u (i) (x), A(u (j) (x))) H .
Then it follows that Z = L√ µ L T √ µ . where

L√ µ =      c 1 √ µ 1 e -µ 1 (T -t 1 ) 0 . . . 0 c 1 √ µ 1 e -µ 1 (T -t 2 ) c 2 √ µ 2 e -µ 2 (T -t 2 ) . . . 0 . . . . . . . . . . . . c 1 √ µ 1 e -µ 1 (T -tn) c 2 √ µ 2 e -µ 2 (T -tn) • • • c n √ µ n e -µn(T -tn)      .
Moreover, if consider the residual at time T given by

R(û(x, T )) = f (x) - l j=1 a j (T )ψ j (x), (13)
and in a similar way we enforce the residual to be orthogonal to each basis functions

(R(û(x, T )), ψ i (x)) H = 0 (14)
for i = 1, 2, . . . , l. We obtain that

a i (T ) = (f (x), ψ i (x)) H (15)
for i = 1, 2, . . . , l.

Let the n × l matrix

W µ = w 1 √ µ 1 w 2 √ µ 2 • • • w l √ µ l and a(t) = a 1 (t) a 2 (t) • • • a l (t) T .
From (15), it follows that

a(T ) = 1 T L c W µ T . ( 16
)
where

L c =      c 2 1 e -µ 1 (T -t 1 ) 0 0 . . . 0 c 2 1 e -µ 2 (T -t 2 ) c 2 2 e -µ 2 (T -t 2 ) 0 . . . 0 . . . . . . . . . . . . . . . c 2 1 e -µ 1 (T -tn) c 2 2 e -µ 2 (T -tn) c 2 3 e -µ 3 (T -tn) • • • c 2 n e -µn(T -tn)      and 1 = [1 1 • • • 1] T ∈ R n .
Then the reduced model can be written as

d dt a(t) = W T µ ZW µ a(t), a(T ) = 1 T L c W µ T , (17)
with solution

a(t) = exp -W T µ ZW µ (T -t) 1 T L c W µ T . (18)
Thus, we can write

û(x, t) = a(t) T      ψ 1 (x) ψ 2 (x) . . . ψ l (x)      .

Numerical Results

We consider that X follows a Geometric Brownian Motion with b(x) = rx and a(x) = σx. Assume that r(x) = r and let

ν = 1 σ r - σ 2 2
Then, from [START_REF] Davidov | Pricing Options on scalar diffusions: An eigenfunction expansion approach[END_REF], we known that

φ j (x) = σ ln(U/L) x -2ν/σ sin πj ln(x/L) ln(U/L) (19)
for j = 1, 2, . . . , and x ∈ [L, U ], form the complete orthonormal basis in H obtained for the regular Sturm-Liouville boundary value problem [START_REF] Thomas | Numerical Partial Differential Equations: Finite Difference Methods[END_REF] and where

µ j = r + ν 2 2 + σ 2 π 2 j 2 2 ln 2 (U/L)
.

for j = 1, 2, . . . . We consider the payoff of an European Call f (x) = (x -K) + on [L, U ]. Then it follows from [START_REF] Davidov | Pricing Options on scalar diffusions: An eigenfunction expansion approach[END_REF], that

c j = L ν/σ ln(U/L) [Lγ j (ν + σ) -Kγ j (ν)] ,
where

γ j (a) = 2 ω 2 j + a 2 e aκ (ω j cos(ω j k) -a sin(ω j k)) -(-1) j ω j e aη , with η = 1 σ ln(U/L), ω j = jπ u , κ = 1 σ ln(K/L).
for j = 1, 2, . . . . Figure 1 plots call prices for L = 9, U = 12, r = 0.1, σ = 0.6, K = 10 T = 1/12 by using the representation of the call price by means 20 eigenfunctions against the call price computed by û(x, t) with n = 4.

We remark that we need to known previously a basis function obtained by solving a Sturm-Liouville problem, that allow to us to construct a solution of the problem under consideration. Thus, in practical application when no closed solution is available the above methodology doesn't runs. The main problem in real applications is that the price u(x, t) is known in a discrete form, that is, at the nodes of a spatial grid or mesh and for some times. In this way, in the next section we apply the model reduction to a numerical scheme by using the information of either a previous simulation or a empirical data set.

Reduced Finite Difference Schemes of implicit type in Option Pricing Models

In this section we would to introduce the Model Reduction in the context of Finite Difference schemes of implicit type. To this end consider the boundary value problem (4) where u(x, t) = g(x, t), for (x, t) ∈ [0, T ] × {0, a},

and u(x, T ) = f (x) for x ∈ [0, a].
Note that f (T ) = g(0, T ) and f (a) = g(a, T ). Here, f and g are given data. We assume that the problem under consideration has a unique solution u with certain smoothness that makes the following calculation meaningful.

To develop a finite difference method, we need to introduce grid points. Let N x and N t be integers, h x = a/N x , k t = T /N t and define the partition points

x j = jh x , j = 0, 1, . . . , N x , t m = T -mh t m = 0, 1, . . . , N t .
A point of the form (x j , t m ) is called grid point and we are interested in computing approximate solutions values at the grid points. We use the notation U m j for an Then, by using the boundary conditions we have

U 0 j = f (x j ) for 0 ≤ j ≤ N x ; U m 0 = g(0, t m ), and U m Nx = g(a, t m ), for 0 ≤ m ≤ N t . Now, set b j = b(x j ),
a j = a(x j ) and

r j = -r(x j )
for 0 ≤ j ≤ N x . From the Crank-Nicolson scheme we obtain that

1 k t U m j -U m+1 j + b 2 j 4 1 h 2 x U m+1 j+1 -2U m+1 j + U m+1 j-1 + b 2 j 4 1 h 2 x U m j+1 -2U m j + U m j-1 + a j 2 1 2h x U m+1 j+1 -U m+1 j-1 + a j 2 1 2h x U m j+1 -U m j-1 + r j 2 U m+1 j + r j 2 U m j = 0, for 1 ≤ j ≤ N x -1 and 0 ≤ m ≤ N t -1.
This can be written as

-A j U m+1 j-1 + (1 -B j ) U m+1 j -C j U m+1 j+1 = A j U m j-1 + (1 + B j ) U m j + C j U m j+1 where A j = 1 4 ν 1 b 2 j + 1 4 ν 2 a j , B j = -1 2 ν 1 b 2 j + 1 2 kr j , C j = 1 4 ν 1 b 2 j -1 4 ν 2 a j .
Written in matrix form the above expressions we have the difference equation

      -A 1 1 -B 1 -C 1 0 • • • 0 -A 2 1 -B 2 • • • • • • • • • 0 • • • • • 1 -B Nx-2 -C Nx-2 0 • • • 0 -A Nx-1 1 -B Nx-1 -C Nx-1                    U m+1 0 U m+1 1 • • • • U m+1 Nx-1 U m+1 Nx              =       A 1 1 + B 1 C 1 0 • • • 0 A 2 1 + B 2 • • • • • • • • • 0 • • • • • 1 + B Nx-2 C Nx-2 0 • • • 0 A Nx-1 1 + B Nx-1 C Nx-1                   U m 0 U m 1 • • • • U m Nx-1 U m Nx             .
Since the values U m 0 and U m Nx for 0 ≤ m ≤ N t are known, then the above system is equivalent to (20) where

U m =    U m 1 . . . U m Nx-1    , Q =       1 -B 1 -C 1 0 • • -A 2 1 -B 2 • • • • • • • 0 • • • 1 -B Nx-2 -C Nx-2 • • 0 -A Nx-1 1 -B Nx-1       , g m =        A 0 U m 0 + A 0 U m+1 0 0 . . . 0 C Nx U m Nx + C Nx U m+1 Nx        and f (U m ) =        (1 + B 1 ) U m 1 + C 1 U m 2 A 2 U m 1 + (1 + B 2 ) U m 2 + C 2 U m 3 . . . A Nx-2 U m Nx-3 + (1 + B Nx-2 ) U m Nx-2 + C Nx-2 U m Nx-1 A Nx-1 U m Nx-2 + (1 + B Nx-1 ) U m Nx-1       
.

We need to solve N t linear systems with a (N x -1)×(N x -1) coefficient matrix. We would to reduce the dimension of the associated linear system. To this end we construct an orthogonal basis of dimension l in order to restrict the map Q onto this subspace.

Proper Orthogonal Decomposition

Now, we introduce the POD method to be used to reduce the finite-difference scheme of implicit type (20). Let consider a grid in [0, a] × [0, T ] of points

{(x j , tm ) : j = 1, 2, . . . , N x -1, m = 0, 1, . . . , Nt },
here Nt is the snapshots number. We assume that in this case we can obtain the snapshots

u (m) s = [u(x 1 , tm ) u(x 2 , tm ) • • • u(x Nx-1 , tm )] T ∈ R Nx-1 .
for m = 0, 1, . . . , Nt , either by experimental data or by previous computer simulation of a similar problem. In this situation we consider the Hilbert space R Nx-1 with the usual inner product. It easy to see that the map R n is given by

R n (v) = Nt i=1 (v, u (i) s ) R Nx-1 u (i) s = u s u T s v where u s = u (1) s u (k) s • • • u ( Nt) s ∈ R (Nx-1)× Nt .
Thus, the l-modal basis, that we denote by {Ψ 1 , . . . , Ψ l }, can be computed solving the eigenvalue problem

u s u T s v = λv. (22)
This is called the direct method.

The normalized eigenvectors in (22) are the left singular vectors in the Singular Value Decomposition (SVD) of u s , that is, the eigenvectors Ψ j , for j = 1, 2, . . . , k are the columns of Ψ where

u s = ΨΣV T . (23)
Here Ψ and V have orthonormal columns and Σ is a square diagonal matrix. The diagonal elements of Σ are the non-zero singular values, that is, √ λ. The matrix representation Ψ for the POD basis vectors can also be computed by finding the right singular vectors of u s as described below. The eigenvalues and eigenvectors for

u T s u s v = λv (24)
correspond to the right singular vectors and singular values of (22), because

u T s = V ΣΨ T .
Since, ΣV T V Σ -1 = I, then, from (23), we obtain

Ψ = u s V Σ -1 . (25)
The latter method is called the snapshot method or the sample method. We note that if the data u

(i)
s are linearly independent and if Nt < N x -1 then it is convenient to use this methodology.

A reduced finite-difference scheme of implicit type

Now we construct a reduced finite-difference scheme by using the representative modes computed either by the direct method or by the snapshot method. To this end we consider the matrix of representative modes

Ψ * = [Ψ 1 Ψ 2 • • • Ψ l ]
where l << N x -1 is choosed by means the condition

l i=1 λ i / Nx-1 i=1 λ i > 0.99.
Then, we would to search an approximate solution Ûm+1 of the difference equation (20) satisfying that

Ûm+1 = Ψ * a m+1 = l j=1 a m+1 j Ψ j , (26)
where

a m+1 = [a m+1 1 a m+1 2 • • • a m+1 l
] T , for m = 0, 1, . . . , N t -1. Observe that (26) is equivalent to the Galerkin expansion (10). Moreover,

Q| ColΨ * = (Ψ * ) T QΨ * ,
where Col Ψ * = span{Ψ 1 , Ψ 2 , . . . , Ψ l }. Thus, a m+1 is the solution of the reduced model given by the difference equation

Q| ColΨ * a m+1 = (Ψ * ) T f ( Ûm ) + g m , Û0 = u 0 . ( 27
)
for m = 0, 1, . . . , N t -1. Solving (27) we obtain the matrix

a = [a 1 a 2 • • • a Nt ]
and finally we have that Û = Ψ * a.

Note that in the reduced version we solve N t linear systems with a l × l coefficient matrix.

Numerical Results

We consider a vanilla European put with strike K = 8 and maturity T = 1/3 and we assume that the volatility and interest rates are constant σ = 0.30 and r = 0.05. Moreover, we consider N x = 200 and N t = 400. Plotting the pointwise error for the Crank-Nicolson scheme Figure 2 shows that the error is concentrated around the singularity, that is t = 0 and x = K. Thus if we would to improve the accuracy the grid should be highly refined near t = 0 and x = K. However, finite difference methods of implicit type presented so far rely on uniform grids in the x-variable. Therefore, at this point, we keep in mind to use a reduced Crank-Nicolson scheme. First, we generate the snapshots by using the closed form solution for the European Put in a grid where x j = xj for j = 1, 2, . . . , N x -1 and Nt = 10. We obtained, in the numerical experiments implemented for different values of r and σ and Nt ≥ 10, only one basis function, that is, l = 1. In Figures 3 and4 we plot the pointwise errors produced at t = 0 from the reduced Crank-Nicolson scheme using the direct method or the snapshot method respectively, against the Crank-Nicolson scheme. As we can see in all cases the Reduced scheme has more accuracy near x = K. A similar result is obtained when we consider the snapshots generated by a previous Crank-Nicolson taking tm = 40m for m = 1, 2, . . . , 10; (see Figures 5 and6). Note, that due to the fact that the modal basis is composed by only one vector, we solve on each iteration, in the case of the Reduced Crank-Nicolson scheme, a one-dimensional linear equation. Now, we compute by means a Crank-Nicolson with parameter values r P OD = 0.10 and σ P OD = 0.3 a modal basis by using the snapshot method. Then we use this modal basis to run a Reduced Crank-Nicolson with parameter values r = 0.05 and σ = 0.4. In Figures 7 and8 we show as the accuracy near x = K is, in both cases, unexpected good. In Figure 9 we plot ||u Nt -U Nt || R Nx against ||u Nt -ÛNt || R Nx as functions of (r, σ), where ÛNx was computed by using the Reduced Crank-Nicolson scheme by means the modal basis with parameter values r P OD and σ P OD . From this we can conclude that it is possible to use a modal basis obtained from a previous computer simulation in order to valuate, by using the reduced finite difference scheme, options with a different parameter values.

Concluding Remarks

We have presented an effective method for the numerical solution of a Pricing Partial Differential Differential Equation of European type, due to the use of PODbased reduced order models applied to a Finite Difference scheme of implicit type. Moreover we observe from the numerical experiments that the efficiency of the reduced scheme, constructed from a previously computed basis, remains for a wide set of parameters values.

Figure 1 :

 1 Figure 1: The maps u(x, 0) = 20 j=1 c j e -µ j T φ j (x) and û(x, 0) with n = 4.

Figure 2 :

 2 Figure 2: The pointwise error produced by the Crank-Nicolson scheme for an European Put in [0, 16] × [0, 1/3], K = 8, σ = 0.30, and r = 0.05.

Figure 3 :

 3 Figure 3: The pointwise error at t = 0 produced by the Crank-Nicolson and Reduced Crank-Nicolson (direct method) schemes for an European Put in [0, 16] × [0, 1/3], K = 8, σ = 0.30, and r = 0.05.

Figure 4 :

 4 Figure 4: The pointwise error at t = 0 produced by the Crank-Nicolson and Reduced Crank-Nicolson (snapshot method) schemes for an European Put in [0, 16] × [0, 1/3], K = 8, σ = 0.30, and r = 0.05.

Figure 5 :

 5 Figure 5: The pointwise error at t = 0 produced by the Crank-Nicolson and Reduced Crank-Nicolson (direct method) schemes for an European Put in [0, 16] × [0, 1/3], K = 8, σ = 0.30, and r = 0.05.

Figure 6 :

 6 Figure 6: The pointwise error at t = 0 produced by the Crank-Nicolson and Reduced Crank-Nicolson (snapshot method) schemes for an European Put in [0, 16] × [0, 1/3], K = 8, σ = 0.30, and r = 0.05.

Figure 7 :

 7 Figure 7: The pointwise error at t = 0 produced by the Crank-Nicolson and Reduced Crank-Nicolson (snapshot method) schemes for an European Put in [0, 16] × [0, 1/3], K = 8.

Figure 8 :

 8 Figure 8: The pointwise error at t = 0 produced by the Crank-Nicolson and Reduced Crank-Nicolson (snapshot method) schemes for an European Put in [0, 16] × [0, 1/3], K = 8.

Figure 9 :

 9 Figure 9: The Crank-Nicolson Error function ||u Nt -U Nt || R Nx against the Reduced Crank-Nicolson Error function ||u Nt -ÛNt || R Nx and the cross section at r = 0.03. 26

approximation to u m j = u(x j , t m ). We remark that

for m = 0, 1, . . . , N t and j = 0, 1, . . . , N x are known values. In order to compute the approximate values for the solution of (4) at the grid points we introduce the following vectors. Let

and the corresponding the approximation vector

for m = 0, 1, . . . , N t . Now, assume that by using a typical finite difference scheme of implicit type for equation ( 4) we are to be able to construct the following two-level difference equation

for m = 0, 1, . . . , N t -1 (see for example [START_REF] Thomas | Numerical Partial Differential Equations: Finite Difference Methods[END_REF] and Chapter 3 of [START_REF] Achdou | Computational Methods for Option Pricing[END_REF]). The vector g m is usually constructed from values {u m j : j = 0, N x , m = 0, 1, . . . , N t } and f : R Nx-1 -→ R Nx-1 . Here the matrix Q may depend on k t and h x . Thus, by solving this difference equation we obtain the matrix

which represents the approximate solution to the matrix

An example: The Crank-Nicolson Scheme

Now, we apply the Crank-Nicolson scheme to the PDE (4), for (x, t) ∈ (0, a) × (0, T ) and where

In this case we introduce the following parameters