Maarten V De Hoop
email: mdehoop@math.purdue.edu

Gunther Uhlmann
email: gunther@math.washington.edu

Andr Ás Vasy
email: andras@math.stanford.edu

Herwig Wendt
email: herwig.wendt@irit.fr

MULTISCALE DISCRETE APPROXIMATIONS OF FOURIER INTEGRAL OPERATORS ASSOCIATED WITH CANONICAL TRANSFORMATIONS AND CAUSTICS *

Keywords: Fourier integral operators, caustics, multiscale computations, wave packets, dyadic parabolic decomposition, operator compression, reflection seismology

We develop an algorithm for the computation of general Fourier integral operators associated with canonical graphs. The algorithm is based on dyadic parabolic decomposition using wave packets and enables the discrete approximate evaluation of the action of such operators on data in the presence of caustics. The procedure consists of constructing a universal operator representation through the introduction of locally singularity-resolving diffeomorphisms, thus enabling the application of wave packet-driven computation, and of constructing the associated pseudodifferential joint-partition of unity on the canonical graphs. We apply the method to a parametrix of the wave equation in the vicinity of a cusp singularity.

1. Introduction. In this paper, we develop an algorithm for applying Fourier integral operators associated with canonical graphs using wave packets. To arrive at such an algorithm, we construct a universal oscillatory integral representation of the kernels of these Fourier integral operators by introducing singularity-resolving diffeomorphisms where caustics occur. The universal representation is of a form such that the algorithm based on the dyadic parabolic decomposition of phase space previously developed by the authors applies [START_REF] Andersson | Multiscale discrete approximation of Fourier integral operators[END_REF]. We refer to [START_REF] Bradie | Fast numerical computations of oscillatory integrals related to acoustic scattering, I[END_REF][START_REF] Candès | Curvelets and Fourier integral operators[END_REF][START_REF] Candès | Fast computation of Fourier integral operators[END_REF][START_REF] Candès | A fast butterfly algorithm for the computation of Fourier integral operators[END_REF] for related computational methods aiming at the evaluation of the action of Fourier integral operators.

The algorithm comprises a geometrical component, bringing the local representations in universal form, and a wave packet component, which yields the application of the local operators. Here, we develop the geometrical component, which consists of the following steps. First we determine the location of caustics on the canonical relation of the Fourier integral operator. For each point on a caustic we determine the associated specific rank deficiency and construct an appropriate diffeomorphism, resolving the caustic in open neighborhoods of this point. We determine the (local) phase function of the composition of the Fourier integral operator and the inverse of the diffeomorphism in terms of universal coordinates and detect the largest set on which it is defined. We evaluate the preimage of this set on the canonical relation. We continue this procedure until the caustic is covered with overlapping sets, associated with diffeomorphisms for the corresponding rank deficiencies. Then we repeat the steps for each caustic and arrive at a collection of open sets covering the canonical relation.

In the special case of Fourier integral operators corresponding to parametrices of evolution equations, for isotropic media, an alternative approach for obtaining solutions in the vicinity of caustics, based on a redecomposition strategy following a multiproduct representation of the propagator, has been proposed previously [START_REF] Andersson | Multiscale discrete approximation of Fourier integral operators[END_REF][START_REF] Kumano-Go | Multi-products of phase functions for Fourier integral operators with an applications[END_REF][START_REF] Le Rousseau | Fourier-integral-operator approximation of solutions to first-order hyperbolic pseudodifferential equations I: Convergence in Sobolev spaces[END_REF]. Unlike multiproduct representations, our construction does not involve a subdivision of the evolution parameter and yields a single-step computation. Moreover, it is valid for the general class of Fourier integral operators associated with canonical graphs, allowing for anisotropy.

The complexity of the algorithm for general Fourier integral operators as compared to the noncaustic case arises from switching, in the sets covering a small neighborhood of the caustics, from a global to a local algorithm using a pseudodifferential partition of unity.

As an application we present the computation of a parametrix of the wave equation in a heterogeneous, isotropic setting for long-time stepping in the presence of caustics.

Curvelets, wave packets. We briefly discuss the (co)frame of curvelets and wave packets [START_REF] Candès | Fast discrete curvelet transforms[END_REF][START_REF] Duchkov | Discrete almost symmetric wave packets and multi-scale representation of (seismic) waves[END_REF][START_REF] Smith | A parametrix construction for wave equations with C 1,1 coefficients[END_REF]. Let u ∈ L 2 (R n ) and consider its Fourier transform, û(ξ) = u(x) exp[-i x, ξ ] dx. One begins with an overlapping covering of the positive ξ 1 axis (ξ ′ = ξ 1 ) by boxes of the form (1.1)

B k = ξ ′ k - L ′ k 2 , ξ ′ k + L ′ k 2 × - L ′′ k 2 , L ′′ k 2 n-1
, where the centers ξ ′ k , as well as the side lengths L ′ k and L ′′ k , satisfy the parabolic scaling condition

ξ ′ k ∼ 2 k , L ′ k ∼ 2 k , L ′′ k ∼ 2 k/2
as k → ∞. Next, for each k ≥ 1, let ν vary over a set of ∼ 2 k(n-1)/2 uniformly distributed unit vectors. Let Θ ν,k denote a choice of rotation matrix which maps ν to e 1 and B ν,k = Θ -1 ν,k B k . In the (co-)frame construction, one encounters two sequences of smooth functions on R n , χν,k and βν,k , each supported in B ν,k , so that they form a co-partition of unity, χ0 (ξ) β0 (ξ) + k≥1 ν χν,k (ξ) βν,k (ξ) = 1, and satisfy the estimates

| ν, ∂ ξ j ∂ α ξ χν,k (ξ)| + | ν, ∂ ξ j ∂ α ξ βν,k (ξ)| ≤ C j,α 2 -k(j+|α|/2) . One then forms ψν,k (ξ) = ρ -1/2 k βν,k (ξ), φν,k (ξ) = ρ -1/2 k χν,k (ξ), with ρ k = vol(B k ), satisfying the estimates (1.2) ∀N : |ϕ ν,k (x)| |ψ ν,k (x)| ≤ C N 2 k(n+1)/4 ( 2 k | ν, x | + 2 k/2 x ) -N .
To obtain a (co)frame, one introduces the integer lattice,

X m := (m 1 , . . . , m n ) ∈ Z n ; the dilation matrix, D k = 1 2π L ′ k 01×n-1 0n-1×1 L ′′ k In-1 ; det D k = (2π) -n ρ k ; and points x ν,k m = Θ -1 ν,k D -1 k X m .
The frame elements (k ≥ 1) are then defined in the Fourier domain as φγ (ξ) = φν,k (ξ) exp[-i x ν,k m , ξ ], γ = (m, ν, k), and similarly for ψγ (ξ). The function ϕ ν,k is referred to as a wave packet. One obtains the transform pair

(1.3) u γ = u(x)ψ γ (x) dx, u(x) = γ u γ ϕ γ (x).
2. Fourier integral operators and caustics. We consider Fourier integral operators, F , associated with canonical graphs. We allow the formation of caustics.

2.1. Oscillatory integrals, local coordinates. Let (y, x Ii , ξ Ji ) be local coordinates on the canonical relation, Λ say, of F , and S i a corresponding generating function: If, at a point on Λ, (dy, dx I ) are linearly independent and dx J vanishes, then (dy, dx I , dξ J ) are coordinates on Λ nearby, I ∪ J = {1, . . . , n}, I ∩ J = {∅}, and one can parameterize Λ as X J (y, x I , ξ J )x J , ξ J , where x J = X J (y, x I , ξ J ) locally on Λ (cf. [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators[END_REF]Thm. 21.2.18]). The fact that a (possibly empty) set I exists follows from the canonical graph property, i.e., that (y, η) are local coordinates and dy linearly independent. Then

x Ji = ∂S i ∂ξ Ji , ξ Ii = - ∂S i ∂x Ii , η = ∂S i ∂y .
(2.1)

The coordinates are defined in a standard way on (overlapping) open sets

O i in Λ, that is, (y, x Ii , ξ Ji ) → r(y, x Ii , ξ Ji ) is defined as a diffeomorphism on O i ; let i = 1, . . . , N .
The corresponding partition of unity is written as

(2.2) N i=1 Γ i (r) = 1, r ∈ Λ.
In local coordinates, we introduce

(2.3) Γi (y, x Ii , ξ Ji ) = Γ i (r(y, x Ii , ξ Ji )). Then (F ϕ γ )(y) = N i=1 (F i ϕ γ )(y) with (F i ϕ γ )(y) = Γi (y, x Ii , ξ Ji )a i (y, x Ii , ξ Ji ) exp[i(S i (y, x Ii , ξ Ji ) -ξ Ji , x Ji φ(y,x,ξJ i ) ] ϕ γ (x) dxdξ Ji .
(2.4)

The amplitude a i (y, x Ii , ξ Ji ) is complex and accounts for the KMAH (Keller-Maslov-Arnold-Hörmander; cf. [START_REF] Červený | Seismic Ray Theory[END_REF]) index. We let Σ φ denote the stationary point set (in θ) of φ = φ(y, x, θ). The amplitude can be identified with a half-density on Λ. One defines the 2n-form d φ on Σ φ ,

d φ ∧ d ∂φ ∂θ 1 ∧ • • • ∧ d ∂φ ∂θ N = dy 1 ∧ • • • ∧ dy n ∧ dx 1 ∧ • • • ∧ dx n ∧ dθ 1 ∧ • • • ∧ dθ N .
In the above, we choose λ = (y, x I , ∂φ ∂xJ ) as local coordinates on Λ, while θ = ξ J . Then we get

d φ = |∆ φ | -1 |dλ 1 ∧ • • • ∧ dλ 2n |, ∆ φ = ∂ 2 φ ∂xJ ∂xJ ∂ 2 φ ∂ξJ ∂xJ ∂ 2 φ ∂xJ ∂ξJ ∂ 2 φ ∂ξJ ∂ξJ = -1;
λ is identified with (y, x I , ξ J ). The corresponding half-density equals

|∆ φ | -1/2 |dλ 1 ∧ • • • ∧ dλ 2n | 1/2 .
Densities on a submanifold of the cotangent bundle are associated with the determinant bundle of the cotangent bundle. Let a 0 i denote the leading-order homogeneous part of a i . The principal symbol of the Fourier integral operator then defines a halfdensity,

a 0 i d 1/2 φ .
That is, for a change of local coordinates, if the transformation rule for forms of maximal degree is the multiplication by a Jacobian , then the transformation rule for a half-density is the multiplication by || 1/2 . In our case, of canonical graphs, we can dispose of the description in terms of half-densities and restrict to zero-density amplitudes on Λ.

2.2. Propagator. The typical case of a Fourier integral operator associated with a canonical graph is the parametrix for an evolution equation [START_REF] Duchkov | Extended isochron rays in shot-geophone (map) migration[END_REF][START_REF] Duchkov | Evolution-equation approach to seismic image, and data, continuation[END_REF],

(2.5) [∂ t + iP (t, x, D x )]u(t, x) = 0, u(t 0 , x) = ϕ γ (x),
on a domain X ⊂ R n and a time interval [t 0 , T ], where P (t, x, D x ) is a pseudodifferential operator with symbol in S 1 1,0 ; we let p denote the principal symbol of P . For every (x, ξ) ∈ T * X\{0}, the integral curves (y(x, ξ; t, t 0 ), η(x, ξ; t, t 0 )) of (2.6)

dy dt = ∂p(t, y, η) ∂η , dη dt = - ∂p(t, y, η) ∂y ,
with initial conditions y(x, ξ; t 0 , t 0 ) = x and η(x, ξ; t 0 , t 0 ) = ξ, define the transformation, χ, from (x, ξ) to (y, η), which generates the canonical relation of the parameterix of (2.5), for a given time t = T ; that is, (y(x, ξ), η(x, ξ)) = (y(x, ξ; T, t 0 ), η(x, ξ; T, t 0 )). The perturbations of (y, η) with respect to initial conditions (x, ξ) are collected in a propagator matrix,

(2.7) Π(x, ξ; t, t 0 ) = W 1 W 2 W 3 W 4 = ∂ x y ∂ ξ y ∂ x η ∂ ξ η ,
which is the solution to the 2n × 2n system of differential equations

(2.8) dΠ dt (x, ξ; t, t 0 ) =     ∂ 2 p ∂η∂y (t, y, η) ∂ 2 p ∂η∂η (t, y, η) - ∂ 2 p ∂y∂y (t, y, η) - ∂ 2 p ∂y∂η (t, y, η)     Π(x, ξ; t, t 0 ),
known as the Hamilton-Jacobi equations, supplemented with the initial conditions [START_REF] Červený | Seismic Ray Theory[END_REF][START_REF] Červený | A note on dynamic ray tracing in ray-centered coordinates in anisotropic inhomogeneous media[END_REF] (2.9) Π(x, ξ; t 0 , t 0 ) = I 0 0 I .

Away from caustics the generating function of Λ is S = S(y, ξ) (I i = ∅), which satisfies

∂ 2 S ∂y∂ξ (y, ξ) = ∂x ∂y ξ = W -1 1 , (2.10) ∂ 2 S ∂ξ 2 (y, ξ) = ∂x ∂ξ y = - ∂x ∂y ξ ∂y ∂ξ y = -W -1 1 W 2 , (2.11) ∂ 2 S ∂y 2 (y, ξ) = ∂η ∂y ξ = ∂η ∂x ξ ∂x ∂y ξ = W 3 W -1 1 , (2.12)
where x = x(y, ξ; t 0 , T ) denotes the backward solution to (2.6) with initial time T , evaluated at t 0 . The leading-order amplitude is then (2.13) a(y, ξ/|ξ|) = 1/ det W 1 (x(y, ξ/|ξ|; t 0 , T ), ξ/|ξ|; T, t 0 ), reflecting that a is homogeneous of degree 0 in ξ.

In the vicinity of caustics, we need to choose different coordinates. Admissible coordinates are directly related to the possible rank deficiency of W 1 : One determines the null space of the matrix W 1 and rotates the coordinates such that the null space is spanned by the columns indexed by the set I i . Then (y, x Ii , ξ Ji ) form local coordinates on the canonical relation Λ, as in the previous subsection, and O i is given by the set for which the columns indexed by I i span the null space of W 1 .

3. Singularity-resolving diffeomorphisms. We consider the matrix W 1 (x(y, ξ; t 0 , T ), ξ; T, t 0 ) for given (T, t 0 ) at y 0 = y(x 0 , ξ 0 ; T, t 0 ) and ξ = ξ 0 and determine its rank. Suppose it does not have full rank at this point. We construct a diffeomorphism which removes this rank deficiency in a neighborhood of r 0 = (y 0 , η 0 ; x 0 , ξ 0 ) ∈ Λ, where η 0 = η(x 0 , ξ 0 ; T, t 0 ).

To be specific, we rotate coordinates such that ξ 0 = (1, 0, . . . , 0) (upon normalization). Let us assume that the row associated with the coordinate x 2 generates the rank deficiency. (There could be more than one row / coordinate.) We then introduce the diffeomorphism

Q : x → x = x 1 - α 2 (x 2 -(x 0 ) 2 ) 2 , x 2 , . . . , x n ;
to preserve the symplectic form, we map

ξ → ξ = (ξ 1 , ξ 2 + α(x 2 -(x 0 ) 2 ) ξ 1 , ξ 3 , . . . , ξ n ), yielding a canonical transformation C Q : (x, ξ) → (x, ξ). We note that C Q (x 0 , ξ 0 ) = (x 0 , ξ 0 ). The canonical transformation, C -1 Q , associated with Q -1 is given by x → x = x1 + α 2 (x 2 -x 0,2 ) 2 , x2 , . . . , xn , ξ → ξ = ( ξ1 , ξ2 -α(x 2 -x 0,2 ) ξ1 , . . . , ξn ).
We introduce the pull back,

Q * u(x) = u(Q -1 (x)) = u(x 1 + α 2 (x 2 -(x 0 ) 2 ) 2 , x2 , . . . , xn ). 3.1. Fourier integral representations of Q and Q -1 .
The diffeomorphism Q can be written in the form of an invertible Fourier integral operator with unit amplitude and canonical relation given as the graph of C Q . To see this, we write

(Q * u)(x) = u(X(x)), ((Q -1 ) * ũ)(x) = ũ( X(x)). That is, X = Q -1 and X = Q. The diffeomorphisms Q and Q -1 define the Fourier integral operators with oscillatory integral kernels, (3.1) A Q (x, x) = e -i ξ,x-X(x) dξ, A Q -1 (x, x) = e -i ξ,x-X(x) d ξ.
The generating functions are

S Q (x, ξ) = ξ, X(x) , S Q -1 (x, ξ) = ξ, X(x) ,
respectively. The canonical relations are the graphs of C Q and C Q -1 and are given by

Λ Q = {(x = X -1 (x), ξ, ∂ xX | x=X -1 (x) ; x, ξ)}, Λ Q -1 = {(x = X-1 (x), ξ, ∂ x X | x= X-1 (x) ; x, ξ)}.
The Hessians yield a unit amplitude:

det ∂ 2 ξ, X(x) ∂ x∂ξ = 1, det ∂ 2 ξ, X(x) ∂x∂ ξ = 1.
Substituting the particular diffeomorphism, we obtain

∂ x X| x= X-1 (x) =      1 -α(x 2 -x 0,2 ) 0 • • • 0 1 0 • • • 0 0 1 • • • . . . . . . . . . . . .      , ξ, ∂ x X | x= X-1 (x) =    ξ1 ξ2 -α(x 2 -x 0,2 ) ξ1 . . .    .
The corresponding propagator matrices are hence given by

Π Q = ∂ x ∂x ∂ x ∂ξ ∂ ξ ∂x ∂ ξ ∂ξ =               1 -α(x 2 -x 0,2 ) 0 • • • 0 0 0 • • • 0 1 0 • • • 0 0 0 • • • 0 0 1 • • • 0 0 0 • • • . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 • • • 1 0 0 • • • 0 αξ 1 0 • • • α(x 2 -x 0,2 ) 1 0 • • • 0 0 0 • • • 0 0 1 • • • . . . . . . . . . . . . . . . . . . . . . . . .               , (3.2)
Π -1 Q = ∂x ∂ x ∂x ∂ ξ ∂ξ ∂ x ∂ξ ∂ ξ =               1 α(x 2 -x 0,2 ) 0 • • • 0 0 0 • • • 0 1 0 • • • 0 0 0 • • • 0 0 1 • • • 0 0 0 • • • . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 • • • 1 0 0 • • • 0 -α ξ1 0 • • • -α(x 2 -x 0,2 ) 1 0 • • • 0 0 0 • • • 0 0 1 • • • . . . . . . . . . . . . . . . . . . . . . . . .               , (3.3)
which are easily verified to be symplectic matrices. In the more general case, each coordinate x j generating a rank deficiency yields additional nonzero entry pairs in the above propagator matrices. Fig. 1. Illustration of canonical relations χ (top) and χ (bottom) of operators F and F associated with a half-wave equation: (bi-)characteristics ("rays") in y for initial conditions (x 2 = x 2,0 , ξ = ξ 0 ) and (x 2 = x2,0 , ξ = ξ0 ), respectively, for evolution through a low velocity lens (see section 5). The black circles on the left indicate the conjugate points corresponding to the initial conditions.

Operator composition. It follows that the composition (x, ξ)

C -1 Q → (x, ξ) χ → (y, η)
generates the graph of a canonical transformation, χ say, which can be parametrized by (y, ξ) locally on an open neighborhood of (y 0 , ξ(x 0 , ξ 0 )) (see Figure 1 for an illustration). We denote the corresponding generating function by Š = Š(y, ξ). We can compose F with Q -1 as Fourier integral operators: F = F Q -1 . The canonical relation of F is the graph of χ. In summary, see the following diagram:

(x, ξ) Q: CQ F : χ / / (y, η) (x, ξ) F =F Q -1 : χ 8 8 r r r r r r r r r r r r r r r r r r r r r r r r r r r Q -1 : C Q -1 S S
For each given type of rank deficiency (here, in x 2 ) and each (x 0 , ξ 0 ) within this class, there is an open set O (x0,ξ0) on which the coordinates (I, J) are valid (see Figure 2). These sets form an open cover, and we obtain a family of diffeomorphisms parametrized by (x 0 , ξ 0 ); there exists a locally finite subcover, and we need only a discrete set to resolve the rank deficiencies everywhere. We index these by j = 1, . . . , N i and construct a set of diffeomorphisms, {Q ij } Ni j=1 , which locally resolve the rank deficiency leading to coordinates (y, x Ii , ξ Ji ). We write

(y, x Ii , ξ Ji ) κij -→ (y, ξ) ↑ ↓ r ↑ ↓ r Λ ∋ r = (y, η; x, ξ) CQ ij -→ (y, η; x, ξ) = ř ∈ Λij .
We write Ǒi for the image of O i under the diffeomorphism on the level of Lagrangians. Let the matrix

Γ ij (r) = 1, while Γij (y, x Ii , ξ Ji ) = Γ ij (r(y, x Ii , ξ Ji )), j = 1, . . . , N i . Then (F ϕ γ )(y) = N i=1 Ni j=1 (F ij ϕ γ )(y) with (3.5) (F ij ϕ γ )(y) = Γij (y, x Ii , ξ Ji )a i (y, x Ii , ξ Ji ) exp[i(S i (y, x Ii , ξ Ji )-ξ Ji , x Ji )] ϕ γ (x)dxdξ Ji .
Inserting the diffeomorphisms, we obtain

(3.6) (F ij ϕ γ )(y) = Ǎij (y, ξ) exp[i Šij (y, ξ)] Q * ij ϕ γ ( ξ) d ξ.
The amplitude Ǎij (y, ξ) and phase function Šij (y, ξ) -ξ, x are obtained by composing F ij with Q -1 ij as Fourier integral operators and changing phase variables. It is possible to treat this composition from a semigroup point of view. Then, to leading order, we get

(3.7) Ǎij (y, ξ) = Γ ij (y, ξ) ǎij (y, ξ),
where

(3.8) Γ ij (y, ξ) = Γij (ř(y, ξ)), in which (3.9) Γij (ř(r)) = Γ ij (r).
Moreover, ǎij (y, ξ) can be obtained as follows. If Π is the propagator matrix of the perturbations of χ, then the propagator matrix of the perturbations of χ is given by Πij

= Π Π -1 Qij . Then (3.10) ǎij (y, ξ) = 1/ det ∂ Šij (y, ξ) ∂y∂ ξ -1
, where det

∂ Šij (y, ξ) ∂y∂ ξ
-1 is obtained as the determinant of the upper-left subblock of Πij . To accommodate a common notation, we set

Q ij = I (N i = 1) if I i = ∅ and write Q i .
In the further analysis, we omit the subscripts i and ij where appropriate.

Expansion of the cutoff functions. To numerically evaluate (3.6) in reasonable time, we will use separated (in y and ξ) representations of Ǎij (y, ξ) and Šij (y, ξ) [START_REF] Andersson | Multiscale discrete approximation of Fourier integral operators[END_REF][START_REF] Candès | Fast computation of Fourier integral operators[END_REF][START_REF] Candès | A fast butterfly algorithm for the computation of Fourier integral operators[END_REF]. Such representations can be obtained by restricting the integration over ξ to domains following a dyadic parabolic decomposition. Here, these will be given by the boxes B ν,k ( ξ) following the redecomposition of

Q * ϕ γ into wave packets ϕ γ , Q * ϕ γ = γ u γ ϕ γ .
The key novelty is constructing a separated representation of the partition functions.

Consider our oscillatory integral in (y, ξ) including the cutoff Γ (y, ξ). Γ (y, ξ) is homogeneous of degree zero in ξ and is a classical smooth symbol (of order 0). We "subdivide" the integration over ξ. A possible procedure involves obtaining a (lowrank) separated representation of Γ (y, ξ) on the support of each relevant box in ξ [3, 5, 4], (3.11) Γ (y, ξ) =

J ν,k β=1 Γβ 1 (y) Γβ 2 ( ξ), ξ ∈ B ν,k .
(Basically, this can be obtained using spherical harmonics in view of the fact that the ξ is implicitly limited to an annulus.) One can view this also as windowing the directions of ξ into subsets (cones) using Γβ 2 ( ξ) and then constructing Γβ 1 (y) according to the smallest admissible set in y for the β-range of directions.

The oscillatory integral becomes

(3.12) (F ϕ γ )(y) = ν,k J ν,k β=1 Γβ 1 (y) ǎ(y, ν) exp[i Š(y, ξ)] Γβ 2 ( ξ) | χν,k ( ξ)| 2 Q * ϕ γ ( ξ)d ξ.
One can view Γβ 2 ( ξ) χν,k ( ξ) as a subdivision of the box B ν,k . We know that |J ν,k | → 1 as k → ∞ since the cone of directions in B ν,k shrinks as √ k. Hence, for large k this does not involve any action.

The procedure allows a subdivision for coarse scales, as long as the scaling is not affected for large k. If the subdivision is too "coarse," then parts of the integration will be lost.

Computation.

We describe an algorithm for applying Fourier integral operators in the above constructed universal oscillatory integral representation. The global hierarchy of operations is given by the following steps:

Q ij ; (c) construction of separated representation for Γij . 2. Evaluation of diffeomorphisms (Q * ij ϕ γ )(x). 3. Evaluation of actions of ( Fij (Q * ij ϕ γ ))(y).
Step 3 requires the evaluation of the action of Fourier integral operators associated with canonical graphs in microlocal standard focal coordinates. The choice of discretization and algorithm for step 3 determines how computations in steps 1 and 2 are to be organized. Here, we perform computations in the almost symmetric wave packet transform domain. We make use of the "box-algorithm" computation of the action of Fourier integral operators associated with canonical graphs in microlocal standard focal coordinates (y, ξ) [START_REF] Andersson | Multiscale discrete approximation of Fourier integral operators[END_REF]. The box algorithm is based on the discretization and approximation, to accuracy O(2 -k/2 ), of the action of Fij on a wave packet ϕ j,ν,k (x),

(4.1) ( Fij ϕ γ )(y) ≈ Ǎ(y, ν) R r=1 α (r) ν,k (y) ξ∈B ν,k e i T ν,k (y),ξ | χν,k ( ξ)| 2 θ(r) ν,k ( ξ).
The procedure relies on truncated Taylor series expansions of Šij (y, ξ) and Ǎ(y, ξ) near the microlocal support of ϕ γ , along the ν = ξ′ /| ξ′ | axis and in the ξ′′ directions perpendicular to the radial ν = ξ′ direction. Here, T ν,k (y) is the backwards solution

x(y) = T ν,k (y) = ∂ Šij (y, ν) ∂ ξ , and α (r)
ν,k (y) and ϑ (r) ν,k ( ξ) are functions realizing, on B ν,k , a separated tensor-product representation of the slowly oscillating kernel appearing in the second-order expansion term of Šij ,

(4.2) exp i 1 2 ξ′ ξ′′ , ∂ 2 Šij ∂ ξ′′2 (y, ν) ξ′′ B ν,k ( ξ) ≈ R r=1 α (r) ν,k (y) θ(r) ν,k ( ξ),
constructed from prolate spheroidal wave functions [START_REF] Beylkin | Wave propagation using bases for band-limited functions[END_REF][START_REF] Slepian | Prolate spheroidal wave functions, Fourier analysis and uncertainty-IV: Extensions to many dimensions, generalized prolate spheroidal wave functions[END_REF][START_REF] Slepian | On the symmetrized Kronecker power of a matrix and extensions of Mehler's formula for Hermite polynomials[END_REF][START_REF] Slepian | Prolate spheroidal wave functions, Fourier analysis and uncertainty-V: The discrete case[END_REF][START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF]. The number R of expansion terms is controlled by the prescribed accuracy ε of the tensor product representation. For a detailed description of the box-algorithm and its implementation, we refer the reader to [START_REF] Andersson | Multiscale discrete approximation of Fourier integral operators[END_REF].

Based on this tensor product representation, it is possible to group computations and to evaluate the action of Fij in step 3 for all data wave packets of the same frequency box B ν,k at once instead of for each ϕ γ individually. Consequently, steps 1 and 2 will also be organized in terms of frequency boxes B ν,k . We write u ν,k (x) = m u m,ν,k ϕ m,ν,k (x) for the data portion corresponding to a frequency box B ν,k . Starting from data u(x), each of the following steps are repeated for any frequency box B ν,k of interest.

4.1. Preparation step. We begin with determining the sets O i for the box B ν,k . To this end, we compute the integral curves (y(x, ξ), η(x, ξ)) and their perturbations with respect to initial conditions (x, ξ) and monitor the null space of the matrix ∂y ∂x following section 2.2. For parametrices of evolution equations, this involves solving the system (2.6) and evaluating the propagator matrices Π(x, ξ) by solving system (2.8). We evaluate the system of differential equations (2.8) in Fermi (or ray-centered) coordinates, in which the potential rank deficiencies of the upper left subblock W 1 (x, ξ) appear explicitly as zero entries in the corresponding row(s) and column(s) [START_REF] Červený | Seismic Ray Theory[END_REF]. The submanifolds Σ (x,ξ) on which W 1 is singular separate and define the sets O i . These computations on Λ F are performed by discretizing the set of orientations ν = ξ/|ξ| covering the frequency box B ν,k with resolution δ ν and the set in x for which u ν,k (x) has nonzero energy with resolution δ x .

Then, for each set O i , we detect Ǔij (and consequently Ǒij ) in a similar way, as the set on which the upper-left subblock W1,ij of Πij = Π Π -1 Qij has full rank. Here Π -1 Qij (x, ξ) is given by (3.3). The operators

Q ij are chosen such that { Ǒij } Ni j=1
overlappingly cover the set of singularities. For fixed α i , this induces a discrete set {x j 0 } Ni j=1 .

We then proceed with the construction of the partition of unity. Since the partition functions enter the computation as pseudodifferential cutoffs in the construction of the amplitude (cf. (3.7)), requiring the backwards solutions x(y, ξ) (compare (2.10)-(2.13)), we perform our numerical construction in coordinates (x, ξ). We obtain Γij (y, ξ) upon substituting y = y(x, ξ) implied by the canonical relation χij . For the construction of the partition functions Γij , we choose double-exponential cutoffs of the form exp(-exp(d(x, ξ))), mimicking a C ∞ 0 cutoff, with appropriate normalization and truncated to precision ε. Here d(x, ξ) is a function measuring the distance of the point (x, ξ) from the boundary ∂ Ūij of the set Ūij (x, ξ). The partition of unity is then formed by weighting Γij (x, ξ) on the overlaps of the sets Ūij (x, ξ) such that ij Γij (ř(x, ξ)) = 1.

Finally, we construct the separated representations of Γij (cf. (3.12)) in (x, ξ) coordinates by windowing the directions of ξ into subsets using Γβ 2 ( ξ), realizing a subdivision into ξ cones. This subdivision is performed for each frequency box B ν,k .

Evaluation of diffeomorphisms.

We evaluate each of the operators Q ij in the Fourier domain. This choice is guided by the property m u m,ν,k φm,ν,k (ξ) = û(ξ) βν,k (ξ) χν,k (ξ) of the discrete almost symmetric wave packet transform [START_REF] Duchkov | Discrete almost symmetric wave packets and multi-scale representation of (seismic) waves[END_REF], which enables the fast evaluation of the Fourier transform of the data at a set of frequency points ξ ν,k l limited to the box B ν,k . We obtain (Q * ij ϕ γ )(x) at once for all ϕ γ (x) belonging to the frequency box B ν,k by evaluation of their adjoint unequally spaced FFT [START_REF] Dutt | Fast Fourier transforms for nonequispaced data[END_REF][START_REF] Dutt | Fast Fourier transforms for nonequispaced data II[END_REF],

F US * ξ → x , at points x(x) = Q -1 ij (x), ǔij ν,k (x) = m u m,ν,k (Q * ij ϕ m,ν,k )(x) = F US * ξ=ξ ν,k l → x(x) û(ξ) βν,k (ξ) χν,k (ξ) .
In preparation for the evaluation of ( Fij ǔν,k )(y), we compute the discrete almost symmetric wave packet transform of the pull back ǔij ν,k (x), yielding its wave packet coefficients u ij ,גν,k . (3.6)) by evaluation of ( Fij ǔij ν,k )(y) using the box algorithm (cf. (4.1)). Note that numerically significant coefficients u ij ,גν,k of the pull back ǔij ν,k (x) are contained in a small set of boxes B ν,k neighboring the direction ν = ξ 0 /|ξ 0 |. We further subdivide each of these boxes according to the separated representation of Γij . Then, we apply the box algorithm to each subdivision, indexed by triples (β, ν, k), β = 1, . . . , J ν,k . Here, the Taylor series expansion of the generating function Šij (y, ξ) is constructed about the central ξ direction within the support of Γ β 2 ( ξ) χν,k ( ξ), accounting for the induced subdivision of the box B ν,k . Note that subdividing into ξ cones results in a reduction of the range of ξ orientations in each element (β, ν, k) of the subdivision, as compared to the ξ range contained in B ν,k . This reduces the number R of expansion terms in (4.1) and effectively counter-balances the increase by a factor J ν,k , evoked by the separated representation of Γij , of the number of times the box-algorithm has to be applied.

F ij u ν,k )(y) = m u m,ν,k (F ij ϕ m,ν,k )(y) (cf.

Operator hierarchy. The operators

F ij for which Q ij = I, F (I)
ij say, are directly associated with the canonical relation Λ F and involve only computations on Λ F . In the algorithm, we reflect this physical hierarchy of the operators F ij in the construction of the partition of unity. First, we construct a partition of unity for these hierarchically higher operators. Then, we construct a joint partition of the remaining operators on the sets which are not covered by the sets for which Q ij = I.

Redecomposition. Starting from a single box B ν,k and applying Q ij to u ν,k (x), redecomposition of ǔij ν,k (x) results in a set of boxes B ν,k yielding numerically nonzero contribution to the solution. The number of boxes entering the computation is directly proportional to the computational cost of the algorithm. In applications, we therefore aim at keeping this number small and consider only a subset of boxes, yielding the most significant contributions. We choose this subset such that on an open neighborhood of (x 0 , ξ 0 ),

Q -1 ij Q ij ≈ I
to precision ǫ. We can estimate the energy loss induced by the restriction to subsets of B ν,k and renormalize the solution. We illustrate the impact of choices of subsets containing different numbers of boxes on the numerical accuracy of the diffeomorphic identity in Figure 3. Furthermore, the redecomposition of ǔij ν,k (x) yields in general, under the action of Q -1 ij , ξ-values outside the set B ν,k , ξ(x, ξ) ⊃ B ν,k . We monitor ξ(x, ξ) and do not consider their contribution in our computation if | χν,k (ξ(x, ξ))| is below the threshold ε.

The free parameters of the procedure are summarized in Table 4.1.

ϕγ(x) x2 x1 Q * ϕγ( x) x2 x1 Q -1 * γ u γ φ γ( x) x2 x1 γ u γ φ γ( x) x2 x1 Q -1 * γ u γ φ γ( x) x2 x1 γ u γ φ γ( x) x2 x1 Q -1 * γ u γ φ γ( x) x2 x1 γ u γ φ γ( x) x2 x1
Fig. 3. Illustration of diffeomorphism Q, Q -1 and redecomposition for a wave packet ϕγ(x) at frequency scale k = 2. Top row: ϕγ (x) (left) and pull back Q * ϕγ (x) (right). Rows 2 to 4: Redecompositions ν,k u γ φγ (x) of Q * ϕγ (x) using 3, 7, and 9 boxes B ν,k , respectively (right column), and the corresponding image Q * -1 ν,k u γ φγ (x) under the action of Q -1 (left column). An insufficient number of boxes alters the amplitudes and the minimum phase property of wave packets. Increasing the number of boxes yields satisfactory results in an open neighborhood of (x 0 , ξ 0 ).

Numerical example.

We numerically illustrate our algorithm for the evaluation of the action of Fourier integral operators associated with evolution equations. We consider wave evolution under the half-wave equation, that is, the initial value problem (2.5) with symbol

P (x, ξ) = c(x) 2 ||ξ|| 2 ,
in n = 2 dimensions. Here c(x) stands for the medium velocity.

Heterogeneous, isotropic model. We choose a heterogeneous velocity model,

c(x) = c 0 + κ exp(-|x -x 0 | 2 /σ 2 ),
containing a low velocity lens, with parameters c 0 = 2 km/s, κ = -0.4 km/s, σ = 3 km, and x 0 = (0, 14) km. As the initial data, we choose horizontal wave packets at frequency scale k = 2 and k = 3, respectively, in the vicinity of the point x ′ = (0, 5) km. set t 0 = 0 and fix the evolution time to T = 7 s. With this choice of parameters, most of the energy of the solution is concentrated near a cusp-type caustic. We illustrate the induced sets O i and the joint partition of unity Γ i in Figures 4 and5 Operator factorization. We partition the Lagrangian Λ into three sets O i , i = {1, 2, 3}. The sets i = {1, 3} are separated by the caustic. For these sets, we can choose coordinates (y, ξ); hence Q i = I. The set i = 2 contains the caustic. For illustration purposes, in the factorization F ij of F i for i 2, we choose to compute the operator j = 1, which resolves the singularity in an open neighborhood of the point indicated by a black dot on the Lagrangian plotted in Figure 6. This neighborhood contains the cusp of the caustic. Furthermore, we limit our separated representation to one term, J ν,k = 1. (For the corresponding admissible sets and partition functions, see Figure 7(left column).) We restrict the computation of F ij for the initial data at frequency scale k = 2 (k = 3) to 9 (11) boxes B ν,k neighboring the ν direction.

Results. In Figure 8, we plot the contributions of the different components in the factorization of the propagator acting on a single horizontal wave packet at frequency scale k = 2, and compare the results to a time domain finite difference computation. The support of the wave packet within the joint admissible set of the chosen factorization is mostly covered by the set O ij , such that most of its energy is contributed by the operator F ij for which Q ij = I.

We observe that, in the joint admissible set, our algorithm has effectively removed the singularity. We note that the phase of the operator computation matches the phase of the finite difference reference. This includes the KMAH index, which is best observed for operator F 3 , which exclusively contributes to the region beyond the caustic (cf. Figure 8(top left)). Furthermore, note that the amplitude obtained by our algorithm is slightly weaker than the true amplitude. This is consistent with the observations and discussion following Figure 3 and results from the energy leakage induced by restricting the number of boxes in the redecomposition step following the application of Q. We can compensate and renormalize the amplitude by monitoring the energy loss resulting from the restriction. (In Figure 8, we have not renormalized the amplitudes.) Finally, we note that our algorithm yields the correct result in an open neighborhood in the vicinity of the tip of the caustic, for which we have designed the operator F ij . Consistent with this fact, the operator F ij is ineffective for yielding the image of the entire wave packet which, at this low frequency scale, has support extending beyond the admissible set of the operator factors we compute.

These observations are further illustrated in Figure 9(left column), where we plot the contributions of the different components in the factorization of the propagator acting on horizontal wave packets, at higher frequency scale k = 3, centered at four different locations in the vicinity of the caustic tip. Results of a time domain finite difference reference computation are plotted in Figure 9(center column). With these initial data, we explore the open neighborhood about the point for which the operator composition with Q ij resolves the singularity. Indeed, at this frequency scale, we can obtain the image of an entire wave packet with only one operator factor F ij (cf. Figure 9(second row)). For the wave packet located slightly further above the tip of the caustic (top row), we observe a phase artifact in the region of overlap of O i=1 and O ij , which can be explained as follows: The restriction of the separated representation for F i to only one term implies that the computation of the geometry (bicharacteristics) for the entire box B ν,k is exclusively based on one single direction ν. This results in inaccuracies in regions close to the caustics where slight perturbations in ξ yield large variations in y. Furthermore, as discussed above, wave packets exploring the regions beyond the tip of the caustics eventually start to leave the admissible set for F ij (third and bottom line).

We note that, both for removing the phase artifact of F i close to the caustic and for enlarging the admissible set, it is necessary to increase the number of terms J ν,k in the separated representation (3.11) (compare with Figure 7). This is illustrated in Figure 9(right column), where the joint contributions of operators F i and F ij with J ν,k = 11 terms in the separated representation are plotted. Here, the expansion functions Γβ 2 ( ξ) are constructed as cones in ξ with a squared cosine cutoff window. For practical reasons and illustration purposes, only one single frequency box B ν,k has been used in the computation (cf. Figure 3). While this restriction to only one frequency box affects the amplitudes and the phases in the tails of the wave packet, the separated representation remains nonetheless effective in resolving the issues observed above: the admissible set is extended beyond the caustic, and the inaccuracies in the regions of overlap of sets O i and O ij as well as in the regions close to the caustics are considerably reduced. This covering with open sets induces a pseudodifferential partition of unity. Then, for each term of this partition, we apply a factorization of the associated operators using diffeomorphisms reflecting the rank deficiency and resolving the singularity in the set. This factorization admits a parametrization of the canonical graph in universal (y, ξ) coordinate pairs and enables the application of our previously developed box algorithm, following the dyadic parabolic decomposition of phase space, for numerical computations. Hence, our algorithm enables the discrete wave packet-based computation of the action of Fourier integral operators globally, including in the vicinity of caustics. This wave packet description is valid on the entire canonical relation. It can now be used with procedures aiming at the iterative refinement of approximate solutions, and it can drive the construction of weak solutions via Volterra kernels [START_REF] Andersson | A multi-scale approach to hyperbolic evolution equations with limited smoothness[END_REF][START_REF] De Hoop | Regularity and multi-scale discretization of the solution construction of hyperbolic evolution equations with limited smoothness[END_REF].

An alternative approach for obtaining solutions in the vicinity of caustics has been proposed previously [START_REF] Andersson | Multiscale discrete approximation of Fourier integral operators[END_REF][START_REF] Kumano-Go | Multi-products of phase functions for Fourier integral operators with an applications[END_REF][START_REF] Le Rousseau | Fourier-integral-operator approximation of solutions to first-order hyperbolic pseudodifferential equations I: Convergence in Sobolev spaces[END_REF] for the special case of Fourier integral operators corresponding to parametrices of evolution equations for isotropic media. It consists of a redecomposition strategy following a multiproduct representation of the propagator. Here, we avoid the redecompositions and operator compositions following the discretization of the evolution parameter, reminiscent of a stepping procedure. What is more, our construction is not restricted to parametrices of evolution equations, but is valid for the general class of Fourier integral operators associated with canonical graphs, allowing for anisotropy. The cost of the algorithm resides in the construction and application of the separated representation of the pseudodifferential partition of unity.

∂ 2

 2 Šij ∂y∂ ξ in the above be nonsingular on the open set Ǔij , and introduce

Fig. 2 .

 2 Fig. 2. Caustic surfaces Ξ(y, ξ) (dark gray) and Ξ(y, ξ) (light gray) of Λ and Λ corresponding to propagation through a low velocity lens (cf. section 5): The singular regions of Λ and Λ do not intersect.

1 .

 1 Preparation step. Preparation of universal oscillatory integral representation: (a) determination of open sets with local coordinates I i , J i on canonical relation Λ, inducing Q ij ; (b) construction of cutoff functions Γij for the locally singularity-resolving diffeomorphisms

Fig. 4 .

 4 Fig. 4. Top left: Iso-amplitude surface of the partition functions Γi (x, ξ), i = 1, 3, associated with Q i = I: the joint admissible set O 1 ∪ O 3 comprises the exterior of the two sheets. Top right: Iso-amplitude surface of Γij (x(x, ξ), ξ(x, ξ)) for ξ 0 = π/2, x 2,0 = 0, and α = 1 (i = 2, j = 1): the admissible set O ij contains the region on the back of the sheet. Bottom: Boundaries ∂O i , i = 1, 3 (dashed curves), and ∂O ij , i = 2, j = 1 (solid curves), of the admissible domains: clearly, the joint admissible set O 1 ∪ O 3 ∪ O 21 covers Λ.

Fig. 5 .

 5 Fig. 5. Illustration of joint partition of unity for the partition functions and sets in Figure 4 for ξ 0 = 1.67 fixed: Slice of Γi (x, ξ = ξ 0 ) (left), the admissible set U ij and the associated partition function Γij (x(x, ξ), ξ(x, ξ) = ξ 0 ) (center), and the partition function Γij (x(x, ξ), ξ(x, ξ) = 1.67) for O ij realizing the partition of unity with Γi (x, ξ = ξ 0 ).

  .

Fig. 6 .

 6 Fig. 6. Projection Λ(y, η 2 ) of a slice ξ = ξ 0 of the canonical relation Λ associated with a half-wave equation in the vicinity of a caustic (red solid line; the blue dashed lines indicate the neighborhood of the singularity) caused by a low velocity lens. The white solid lines are connected to a regular grid in x by bicharacteristics. The black dot indicates the center of an open neighborhood of conjugate points (x 0 , ξ 0 ) → (y 0 , η 0 ) for which the projection onto standard microlocal focal coordinates (y, ξ) is not diffeomorphic.

Fig. 7 .

 7 Fig. 7. Illustration of admissible sets and expansion functions Γβ 1 (y(x)) in (3.11) with J ν,k = 1 (top row) and J ν,k = 5 (lower rows) for the partition functions in Figure 4. Left column: The partition functions Γi (y(x), ξ), i = 1, 3. Center column: The partition functions Γij (y(x, ξ), ξ(x, ξ)), i = 2, j = 1. Right column: The joint admissible sets induced by the partition of unity for the expansion functions plotted in the left and center columns.

Fig. 8 .

 8 Fig. 8. Illustration of operator action on a wave packet ϕγ (x) at frequency scale k = 2: Contribution of operators F i (i = 1, 3) associated with Q i = I (top left), contribution of operator F ij (i = 2, j = 1) with diffeomorphism parameters (ξ 0 = π/2, x 2,0 = 0, α = 1), resolving the singularity in the tip of the caustic (top right), and joint action of F i and F ij (bottom left). Time domain finite difference reference (bottom right). In the operator computation, we consider 9 boxes B ν,k and a separated representation with J ν,k = 1 term.

Fig. 9 .

 9 Fig. 9. Left column: Joint contribution of the operators F i and F ij acting on a wave packet ϕγ (x) at frequency scale k = 3 (compare with Figure 8(bottom left) for a wave packet at frequency scale k = 2). Center column: Time domain finite difference reference. Right column: The equivalent to the left column when using a separated representation with J ν,k = 11 terms. (Note that for computational reasons, only one box B ν,k has been used in the numerical evaluation of F i and F ij with J ν,k = 11 terms.)

Table 4 . 1

 41 Table of user-defined parameters. Free parameter of operators Q ij , inducing the discrete sets of diffeomorphisms {Q ij } Evaluation of the actions of F ij . At this stage, we are ready to evaluate the action (

	δν , δx	Discretization steps in computations on Λ F and Λij
	α i J ν,k	Number of expansion terms in separated representation of cutoff functions Γ(y, ξ)	N i j=1
	ǫ	Precision of approximate redecomposition of ǔij ν,k (x)
	ε	Accuracy of the tensor-product representation in the box algorithm
	4.3.


Gilles Tagne
email: gilles.tagne@hds.utc.fr

Reine Talj
email: reine.talj@hds.utc.fr

Ali Charara
email: ali.charara@hds.utc.fr

Passivity Analysis and Design of a Robust Nested Passivity-Based Controller for Trajectory Tracking of Autonomous Vehicles

This paper focuses on the lateral control of intelligent vehicles. The aim is to minimize the lateral displacement of the autonomous vehicle with respect to a given reference trajectory. The control input is the steering angle and the output is the lateral error displacement. After passivity analysis of the system to establish the properties of passivity between the inputs and outputs, we present design and validation of a robust lateral nested controller based on passivity, to ensure robust stability and good performances with respect to parametric variations and uncertainties encountered in driving applications. To validate the control strategy, the simulation of the closed-loop system on Matlab-Simulink has been made using the experimental data acquired on the vehicle DYNA of Heudiasyc laboratory (a Peugeot 308), according to several real driving scenarios. The validation shows robustness and good performances of the proposed control approach, and puts in evidence the improvement brought by the Nested PBC controller.

INTRODUCTION

The developments in the applications of ADAS (Advanced Driver Assistance Systems) and autonomous driving have been favored by technological advances in recent years. Three main steps are necessary to ensure an autonomous navigation: the perception and localization, the path planning and the vehicle control. The vehicle control can be divided into two tasks: longitudinal control and lateral control.

Lateral control consists on handling the vehicle using the steering wheel to follow the reference trajectory. In recent years, considerable researches have been made to provide lateral guidance of autonomous vehicles. In literature, many control strategies have been developed. Simple PID controller have been proposed in [START_REF] Broggi | The ARGO autonomous vehicle's vision and control systems[END_REF]). We also have nested controller in [START_REF] Marino | Nested PID steering control for lane keeping in autonomous vehicles[END_REF]). Moreover, other classical techniques have been used. We can cite H ∞ [START_REF] Hima | Trajectory Tracking for Highly Automated Passenger Vehicles[END_REF]), state feedback [START_REF] Rajamani | Vehicle dynamics and control[END_REF]), Lyapunov stability based control [START_REF] Benine-Neto | Dynamic controller for lane keeping and obstacle avoidance assistance system[END_REF]), fuzzy logic [START_REF] Naranjo | Lane-Change Fuzzy Control in Autonomous Vehicles for the Overtaking Maneuver[END_REF]), fuzzy Takagi-Sugeno LQ [START_REF] Soualmi | Fuzzy Takagi-Sugeno LQ controller for lateral control assistance of a vehicle[END_REF]), linear quadratic optimal predictive control [START_REF] Kim | Control strategy for high-speed autonomous driving in structured road[END_REF]) and many others. Model Predictive Control (MPC) appears to be well suited to the trajectory following [START_REF] Levinson | Towards Fully Autonomous Driving : Systems and Algorithms[END_REF]).

In [START_REF] Hingwe | Experimental evaluation of a chatter free sliding mode control for lateral control in AHS[END_REF], Tagne et al. [2013a]), Sliding Mode Control (SMC) has been applied to the lateral control. This robust control strategy is well suited to driving applications, given its robustness against uncertainties and its capacity to reject disturbances. However, its main drawback is the chattering. Given the implicit resemblance between the SMC and the Immersion and Invariance (I&I) principle, in (Tagne et al. [2013b]), a controller based on the I&I approach have been developed to improve SMC perfomances.

Considering the high nonlinearity of the vehicle on one hand, and the uncertainties and disturbances in automotive applications on the other hand, an important issue to be considered in the control design is the robustness. The controller should be able to reject the disturbances and deal with parameter uncertainties and variations. For example, in [START_REF] Levinson | Towards Fully Autonomous Driving : Systems and Algorithms[END_REF]), a recent presentation of Junior; Stanford's autonomous vehicle (the second at the DARPA Urban Challenge), is carried out for the purpose of ensuring robust autonomous driving.

Passivity is a concept that can be used in several areas of science and explains some physical phenomena. The theory of passivity is a framework for analyzing systems and designing controllers using a description of the inputsoutputs based on energy considerations. The main idea is that a very large number of physical systems have some properties of passivity between the inputs and outputs. We study the passivity of the system to analyze the frequency behavior, and determine the passive outputs (to easily control the system). The interest of passivity results from the feedback interconnection of passive systems. We can therefore seek the passivity of the system as a way to impose robust stability and good performances. This is particularly relevant in applications of lateral control of vehicles where the controller must ensure robust stability and good performances with respect to variations of speed, curvature, coefficient of friction of the road, wind and uncertainties on parameters.

The design of the I&I controller in (Tagne et al. [2013b]) allows to prove a very strong stability criterion of the closed-loop system. This result led us to study the properties of passivity of the different input-output maps of the system to design robust controllers. The main purpose of this paper is to establish the properties of passivity between the inputs and outputs. Finally, we develop a Passivity-Based Controller (PBC). To design the controller, we consider that the vehicle is equipped with sensors and/or observers to measure the yaw rate, the lateral error and its derivative. To validate the proposed approach, tests were made with real data acquired on the vehicle DYNA, on the tracks and circuits of CERAM1 . The simulation results show the performance and robustness of the proposed approach. This paper is organized as follows. Section 2 presents the dynamical models of the vehicle used for control design and the control problem definition. And then the passivity analysis is presented in Section 3. The control strategy is developed in Section 4. Section 5 presents the simulation results. Finally, we conclude in Section 6, with some remarks and future work directions.

DYNAMIC MODEL OF VEHICLE AND CONTROL PROBLEM DEFINITION

Dynamic models of vehicle

In this work, we use two vehicle models. To design the controller, a simple and widely used dynamic bicycle model from [START_REF] Rajamani | Vehicle dynamics and control[END_REF]) is considered. This model is used to represent the lateral vehicle behavior and assumes that the vehicle is symmetrical, and tire's sideslip angles on the same axle are equal. The roll and pitch dynamics are neglected and angles are assumed to be small (steering, sideslip, yaw). With a linear tire force model we obtain a Linear Parameter Varying (LPV) model, where the longitudinal velocity V x is considered as a varying parameter. Dynamic equations in terms of slip angle and yaw rate of the bicycle model are given by:

β = - µ(C f + C r ) mV x β -1 + µ(L f C f -L r C r ) mV 2 x ψ + µC f mV x δ ψ = - µ(L f C f -L r C r ) I z β - µ(L 2 f C f + L 2 r C r ) I z V x ψ + µL f C f I z δ (1)
where β, ψ and δ represent respectively the sideslip angle, the yaw angle of the vehicle and the steering wheel angle (control input). Table 1 presents vehicle parameters and nomenclature.

To validate the proposed controller in simulation, we used a more representative model, namely, the 4-wheel model to represent the behavior of the vehicle and Dugoff's tire model for longitudinal and lateral forces.

Control problem definition

The aim of the lateral control of intelligent vehicles is to minimize the lateral displacement of the autonomous vehicle with respect to a given reference trajectory. The lateral error dynamic at the center of gravity of the vehicle, with respect to a reference trajectory, is given by: ë

= a y -a y ref (2)
where a y and a y ref represent respectively the lateral acceleration of the vehicle, and the desired one on the reference trajectory. Assuming that the desired lateral acceleration of the vehicle can be written as a

y ref = V 2 x ρ
, where ρ is the curvature of the road, and given that a y = V x ( β + ψ), we have:

ë = V x ( β + ψ) -V x 2 ρ (3)
Replacing β by its expression in equation ( 1), the new system state variables become x = (β, ψ, ė, e) , corresponding to the sideslip angle, the yaw rate, the lateral error and its derivative. The new system has the following dynamics:

ẋ = Ax + B 1 δ + B 2 ρ (4)
where,

A =          - µ(C f + C r ) mV x -1 - µ(L f C f -L r C r ) mV 2 x 0 0 - µ(L f C f -L r C r ) I z - µ(L 2 f C f + L 2 r C r ) I z V x 0 0 - µ(C f + C r ) m - µ(L f C f -L r C r ) mV x 0 0 0 0 1 0          B 1 =         µC f mV x µL f C f I z µC f m 0         , B 2 =    0 0 -V 2 x 0   
The aim of the lateral control is to cancel the lateral error displacement. Then, for a given curvature ρ and longitudinal velocity V x , the desired behavior corresponds to ė1 = e 1 = 0. Hence, it is easy to prove that the desired equilibrium point is:

(β, ψ, ė, e) = (β , ψ , 0, 0) with β = (L r - L f mV 2 x µC r (L f + L r ) )ρ ψ = V x ρ (5)
At the equilibrium point, the control input is:

δ = L f C f -L r C r L f C f β + L 2 f C f + L 2 r C r L f C f V x ψ (6)
Hence, define the new error variables:

     β = β -β ˙ψ = ψ -ψ δ = δ -δ (7)
The error dynamics of the system (4) having the origin as equilibrium point ( β, ˙ψ, ė, e) = (0, 0, 0, 0) become:

ẋ = Ax + B 1 δ (8)
where, A and B 1 have been defined above (4).

PASSIVITY ANALYSIS

Theory of passivity is a framework for analyzing physical systems and designing controllers using a description of the input-output relationship based on energy considerations. The Kalman-Yakubovich-Popov (KYP) lemma is considered to be one of the pillars for control and systems theory.

To demonstrate the passivity of an output, this lemma is considered. It establishes an equivalence between the conditions in the frequency domain (positive realness of a system), in the time domain (an input-output relationship of the system), and in the state-space representation where conditions on the matrices describing the system should be verified. The Lemma is given as follows: Lemma 1. (KYP lemma) [START_REF] Lozano | Dissipative Systems Analysis and Control: Theory and Applications[END_REF] In the frequency domain

H(s) is stable [H(jω)] ≥ 0, ∀ω ∈ R (9)
In the time domain

A system with input u and output y where

u(t), y(t) ∈ R n is passive if there is a constant such that T 0 y T (t)u(t) ≥ υ (10)
for all function u and all T ≥ 0.

In the state-space representation

Let us consider a system described by the following statespace representation: ẋ = Ax + Bu and y = Cx + Du.

Where

x ∈ R n , u, y ∈ R m with n ≥ m. The transfer function H(s) = C T (sI -A) -1 B +D, with A ∈ R n×n , B ∈ R n×m , C ∈ R m×n , D ∈ R m×m is PR with H(s) ∈ R m×m ,
if and only if there exists matrices P > 0, P ∈ R n×n , L ∈ R n×m and W ∈ R m×m such that:

   P A + A T P = -LL T P B -C T = -LW D + D T = W T W (11)
Remark : In view of the above conditions, it is clear that unstable systems or non-minimum phase systems are not positive real. Assume that the road coefficient of friction µ = 1. Combining dynamic equations of the bicycle model in terms of slip angle and yaw rate (1) with the equation of lateral error dynamic (3), and after some calculations, one can find that the transfer function H 0 (s) between δ as input and ë as output is given by:

H 0 (s) = ë(s) δ(s) = as 2 + bs + c s 2 + ds + f (12)
where,

                             a = C f m b = L r C f C r (L f + L r ) mI z V x c = C f C r (L f + L r ) mI z d = (C f + C r ) mV x + (L 2 f C f + L 2 r C r ) I z V x f = C f C r (L f + L r ) 2 mI z V 2 x + (L r C r -L f C f ) I z (13)
The transfer function H 0 (s) has a relative degree equal to 0. Consequently, according to the KYP Lemma applied in the frequency domain, the proof of SSPR is tantamount to verifying :

[H 0 (jω)] ≥ ζ > 0, ∀ω ∈ (-∞, +∞), (14)
Whatever the uncertainties and variations encountered, the model parameters (C f , C r , L f , L r , m, I z ) are always positive. Thus, the coefficients a, b, c, d, f of H 0 (s) are always positive. Note that f is always positive even if the second fraction of f could become negative. According to the criterion of Routh-Hurwitz, zeros and poles of H 0 (s) are strictly stable, so H 0 (s) is minimum phase and Hurwitz. Setting s = jω, the frequency response function becomes

H 0 (jω) = -aω 2 + bjω + c -ω 2 + djω + f = R N 0 + jI N 0 R D + jI D ( 15
)
where R N 0 , I N 0 , R D and I D are the real and imaginary parts of the numerator and the denominator respectively, given by:

R N 0 = c -aω 2 , I N 0 = bω R D = f -ω 2 , I D = dω Hence, [H 0 (jω)] = R N 0 R D + I N 0 I D R 2 D + I 2 D (16) Knowing that R 2 D + I 2 D > 0, [H 0 (jω)] is positive if and only if R N 0 R D + I N 0 I D > 0. R N 0 R D + I N 0 I D = aω 4 + (bd -af -c)ω 2 + cf (17)
Given the order of magnitude of the parameters mentioned above and whatever is V x > 0, one can find that:

a > 0 cf > 0 bd -af -c is not always positive. ( 18
)
Given that (bd -af -c) is not always positive, let x = ω 2 , hence

R N 0 R D + I N 0 I D = f (x) = ax 2 + (bd -af -c)x + cf (19)
The derivative of this function is zero for

x 0 = bd -af -c -2a So, the unique extremum of the function is f (x 0 ) = (bd -af -c) 2 -4a + cf > 0 (20) Moreover, given that f (0) = cf > 0 and f (∞) → a > 0, then f (x) = R N 0 R D + I N 0 I D ≥ ζ 0 > 0 (21) Given that lim ω→∞ [H 0 (jω)] = a > 0 (22) and lim ω→0 [H 0 (jω)] = c f > 0 (23)
one can prove that:

[H 0 (jω)] ≥ ζ > 0, ∀ω ∈ (-∞, +∞) (24)
hence, the transfer function H 0 (s) is Strongly Strictly Positive Real (SSPR), yielding to the desired result.

Passivity of the map δ → ė

Proposition 2. The map output δ → ė is passive (P).

The transfer function H 1 (s) of the output ė with respect to the input δ is given by:

H 1 (s) = 1 s H 0 (s) (25)
The transfer function H 1 (s) is a cascade connection of an integrator with the strongly strictly positive real transfer function H 0 (s). Hence, H 1 (s) is positive real, yielding the passivity of the map δ → ė, the desired result. For more details see [START_REF] Lozano | Dissipative Systems Analysis and Control: Theory and Applications[END_REF].

Strict Passivity of the map δ → ˙ψ

Proposition 3. The map δ → ˙ψ is Strictly Passive (SP).

Proof.

The transfer function H 2 (s) of the yaw rate error ˙ψ relative to the input δ is given by:

H 2 (s) = ˙ψ(s) δ(s) = gs + h s 2 + ds + f (26)
where,

        
d and f defined in (13),

g = L f C f I z h = (C f + C r )(L f + L r ) mI z V x (27)
The coefficients d, f , g and h are always positive. Then, according to the Routh-Hurwitz criterion, poles of H 2 (s) are stable. The proof of Strict Passivity of the map δ → ˙ψ is equivalent to the proof of Strict Positive Realness (SPR) of the transfer function H 2 (s). Given that the transfer function H 2 (s) has a relative degree equal to 1, the SPR is established verifying the following:

(P1) [H 2 (jω)] > 0, ∀ω ∈ (-∞, +∞), (P2) lim ω→∞ ω 2 [H 2 (jω)] > 0.
Indeed, the frequency response function of H 2 (s) for s = jω is

H 2 (jω) = gjω + h -ω 2 + djω + f = R N 2 + jI N 2 R D + jI D ( 28
)
where R N 2 , I N 2 , R D and I D are the real and imaginary parts of the numerator and the denominator respectively, given by:

R N 2 = h, I N 2 = gω R D = f -ω 2 , I D = dω Hence, [H 2 (jω)] = R N 2 R D + I N 2 I D R 2 D + I 2 D (29) As above, R 2 D + I 2 D > 0, then [H 2 (jω)] is positive if and only if R N 2 R D + I N 2 I D > 0. R N 2 R D + I N 2 I D = (gd -h)ω 2 + f h (30)
Given the order of magnitude of parameters mentioned above and whatever is V x > 0, one can find that: Proof. The transfer function H 3 (s) of β, relative to the input δ is given by:

gd -h > 0 f h > 0 ⇒ [H 2 (jω)] > 0, ∀ω ∈ (-∞,
H 3 (s) = β(s) δ(s) = ks + l s 2 + ds + f (33)
where,

        
d and f defined in ( 13),

k = C f mV x l = L r C f C r (L f + L r ) mI z V 2 x - L f C f I z
The transfer function H 3 (s) has a relative degree equal to 1. The variables d, f and k are always positive, but l can be either positive or negative depending on the value of V x and other parameters. When l < 0, H 3 (s) has a negative zero, and is not positive real, what proves the proposition.

The zero of the transfer function H 3 (s) is stable if and only if l is positive, i.e.

- L f C f I z + L r C f C r (L f + L r ) mI z V 2 x > 0, hence, V x < L r C r (L f + L r ) L f m . ( 34
)
This fact can be understood as if the map δ → β is passive when the speed is limited as in (34), and this characteristic is stolen at high speeds. Using the parameters given in Table I, the speed limit of equation ( 34) corresponds to V x < 16.5m/s 60km/h.

Proposition 5. The map ˙ψ → ë is Strictly Passive and the map ˙ψ → ė is Passive.

Proof. The proof is established at the same manner of the previous subsections.

To summarize, considering the error system (8), the passivity maps can be resumed in the Fig. 1. Probably the main property used in PBC for passive systems is the fact that the feedback interconnection of passive systems is passive. So passivity is invariant under negative feedback interconnection. Therefore, passive systems can be decomposed into passive subsystems. Thus, in this methodology, the controller can be designed as a passive system. However, it will be also useful to know that interconnections not only preserve the passivity properties of the subsystems but, in certain cases, passivity can be strengthened.

Considering the Fig. 1 and knowing that the yaw rate dynamic is faster than the lateral error dynamic, we can decompose the lateral controller into two passive nested controllers. This is a major interest since e and ˙ψ are controlled simultaneously. We can therefore create two separate controllers for each own dynamic based on its characteristic of passivity. In this subsection, we develop a control strategy of an autonomous vehicle using two passive controllers to ensure guidance and stability (lateral error and yaw rate error) in two separate loops.

The outer controller minimizes lateral error. The inner controller minimizes the yaw rate error by providing the corresponding steering angle. The particularity of such a strategy is doubled: it helps to design a robust controller while preserving the passivity properties of the closedloop system to ensure good performance. If necessary, two controllers of different natures can be used, in order to benefit from the robustness of each. It also allows an independent control of the lateral error and yaw rate error. Consequently, the scheme in Fig. 2 illustrates the proposed control strategy.

Fig. 2. Nested PBC control strategy Proposition 6. Consider the diagram in Fig. 2, where the control is achieved by two PI controllers, the closed-loop system is stable and passive.

The proof is relatively simple, the corollary 4.1 in [START_REF] Lozano | Dissipative Systems Analysis and Control: Theory and Applications[END_REF]) allows to show it.

Synthesis of the outer controller:

This controller is designed to cancel the lateral error e.

The control input is the yaw rate error.

Subsystem 1 : ˙ψ -→ ė is passive (P), so any Input Strictly Passive (ISP) controller guarantees stability and passivity of the closed-loop system. In this paper, we choose a simple PI. We remind that a PI applied to the passive output ė is equivalent to a PD applied to the output e. So, the control input given by the outer controller is: ˙ψ

= -K D1 ė -K P 1 e (35)
where K D1 and K P 1 are positive gains.

Synthesis of the inner controller:

This controller is designed to cancel yaw rate error with respect to the reference given by the outer controller. The control input is the steering angle δ.

Subsystem 2 : δ -→ ˙ψ is strictly passive (SP), so any passive (P) controller guarantees stability and passivity of the closed-loop system. In this paper, we choose a PI. The control input given by the inner controller is:

δ = -K I2 ˙ψ -K P 2 ˙ψ ( 36
)
where K I2 and K P 2 are positive gains.

Given the strict passivity (SP) of the output ˙ψ for an input δ, the closed-loop system with a PI controller is stable and passive. Note that any passive controller, a simple proportional for example, would achieve the same result. The addition of integral action, well known to reject constant disturbances, also has the advantage that the controller can be implemented without the knowledge of δ which depends on the uncertain model parameters.

SIMULATION RESULTS

To validate our control law, the experimental data used were acquired by the vehicle DYNA (a Peugeat 308) on the CERAM test circuits. This vehicle is equipped with different sensors and observers providing all the dynamic variables of the vehicle and its traveled path.

The simulations in closed-loop are performed using as a reference, the real data with the full vehicle model.

For the control law, we used the following values of the gains: K D1 = 0.08, K P 1 = 10, K P 2 = 5 and K I2 = 1. To highlight the improvements brought by the PBC controller, a comparison with a previously developed Immersion and Invariance (I&I) controller was made (Tagne et al. [2013b]).

Test of the controller with nominal parameters

This test (Fig. 3) was made with the purpose of testing the ability of the controller to track the reference trajectory during normal driving with known nominal parameters. Note that in this scenario we have a large turn to test the ability of the controllers to ensure proper path tracking for hard maneuvers. In this scenario, the assumption of small angles is not respected (the steering angle is greater than 12 degrees during the turning). In spite of that, the controllers are able to follow the path with low errors. Both nonlinear controllers are robust to variations in longitudinal speed. So, this test show the good performance of the I&I and Nested PBC controller during normal driving and at varying speeds for known nominal parameters. We will evaluate the robustness of these controllers with respect to parametric uncertainties.

Robustness to vehicle parameter uncertainties

As mentioned previously, one of the major challenges for trajectory tracking is the robustness of the controller.

In this subsection, we evaluate the robustness of the controllers with respect to parametric uncertainties of the vehicle. It is important to note that the parametric uncertainties can be due to the fact that the parameters may vary, but are considered to be fixed for the command. This is the case of the mass of the vehicle for example.

It is difficult to estimate accurately the cornering stiffness of the tires. Moreover, this parameter varies greatly depending on the type of the road, the vertical load, camber, etc. It is therefore important to assess the robustness of the controllers over cornering stiffness and mass variations. Fig. 4 presents lateral errors with uncertainties of +/-10% on the value of the cornering stiffness. For uncertainties in the order of 10%, on the value of cornering stiffness or mass, Nested PBC controller is able to follow the path with similar errors (to those of nominal conditions). The errors of the I&I remain acceptable.

Several other tests were performed, including scenarios of driving at high speed as well as situations where there is a saturation of lateral forces. With the I&I controller, the error depends on the value of the parametric uncertainty.

Therefore with sensitivity study, we can estimate in advance the maximum error according to the knowledge of the parametric uncertainty. The passive controller (Nested PBC) is more robust to parameter uncertainties. Indeed, the inputs of the controller depend only on passive outputs ness (+/ -10%) (errors), which do not depend on uncertain parameters. However, the input of the I&I controller depends on these parameters. Furthermore, the PBC controller also has the advantage of using less measurements than the previous one. Indeed, with the I&I controller, the estimation of the sideslip angle is necessary. So, these results put in evidence the improvement brought the Nested PBC controller.

CONCLUSION

In this paper, after a detailed study of the passivity properties of the model, a controller based on these properties has been proposed to ensure robust trajectory following of autonomous vehicles. The design of the controller has been presented. It guarantees robust stability and the passivity of the closed-loop system.

Results have shown that a greater consideration of structural features of the model during the design of the controller, and the passivity properties of the closed-loop system, can significantly improve the robustness of the controller for autonomous driving application. To illustrate this improvement, a comparison has been made with the I&I controller developed previously. The validation has shown robustness and good performance of the proposed Nested PBC controller.

The study of passivity is done considering the vehicle speed as a varying parameter. This is a non-restrictive assumption because the lateral dynamic is relatively faster than longitudinal dynamic. In this paper, the study of passivity properties has been done considering the road coefficient of friction µ = 1. An analysis for variant µ will be studied by the authors. In addition, a study of robustness based on sensitivity analysis will be performed. We will also evaluate the robustness against sensor noise like in [START_REF] Althoff | Set-based computation of vehicle behaviors for the online verification of autonomous vehicles[END_REF]). Very soon, we will test these control laws on the laboratory robotized vehicle.

3. 1

 1 Strong Strict Passivity of the map δ → ë Proposition 1. The map δ → ë is Strongly Strictly Passive. The proof of strong strict passivity is established showing that the transfer function H 0 (s) of the map δ → ë is Strongly Strictly Positive Real (SSPR).

Fig. 1 .

 1 Fig. 1. Passivity maps of the system

   shows the longitudinal speed variations. Fig.3-(b) presents different curves: the reference path and the trajectories followed by the controlled vehicle with two different controllers; the lateral error and the yaw angle error. Both controllers ensure the reference tracking with low errors (the lateral displacement of the controlled vehicle remains less than 10cm in this test conditions).

Fig. 3 .

 3 Fig. 3. (a) Longitudinal speed, (b) Trajectories ; reference and closed-loop simulations

Fig. 4 .

 4 Fig. 4. Robustness against uncertainties of cornering stiffness (+/ -10%) (errors), which do not depend on uncertain parameters. However, the input of the I&I controller depends on these parameters. Furthermore, the PBC controller also has the advantage of using less measurements than the previous one. Indeed, with the I&I controller, the estimation of the sideslip angle is necessary. So, these results put in evidence the improvement brought the Nested PBC controller.

Table 1 .

 1 Vehicle Parameters and Nomenclature

	Vx	Longitudinal velocity	-	[m/s]
	β	Sideslip angle	-	[rad]
	ψ	Yaw rate	-	[rad/s]
	δ	Steering wheel angle	-	[rad]
	µ	Road friction coefficient	1	-
	m	Mass	1719	[kg]
	Iz	Yaw moment of inertia	3300	[kgm 2 ]
	L f	Front axle-COG distance	1.195	[m]
	Lr	Rear axle-COG distance	1.513	[m]
	C f	Cornering stiffness of the front tire	170550	[N/rad]
	Cr	Cornering stiffness of the rear tire	137844	[N/rad]

  Proposition 4. The map δ → β is not passive.

	Note that lim ω→∞	[H 2 (jω)] = 0. Then, the Strict Positive
	Realness (SPR) of H 2 (s) is not Strong, yielding the desired
	result of Strict Passivity of the map δ → ˙ψ.
	3.4 Characteristics of the sideslip angle error β,
		+∞), (31)
		verifying (P1). Furthermore,
		lim

ω→∞ ω 2 [H 2 (jω)] = gd -h > 0,

(32)

verifying (P2). Hence, H 2 (s) is a Strictly Positive Real transfer function.

This work was carried out in the framework of the Labex MS2T, which was funded by the French Government, through the program Investments for the future managed by the National Agency for Research (Reference ANR-11-IDEX-0004-02).

CERAM -"Centre d'Essais et de Recherche Automobile de Mortefontaine" is an automobile testing and research center located in France.

ACKNOWLEDGEMENTS

We would like to thank Rogelio Lozano for our exchange about the theory of passivity, which has been a great help for us.