Alex Bombrun
email: alex.bombrun@gmail.com

Jean-Baptiste Pomet
email: jean-baptiste.pomet@inria.fr.

THE AVERAGED CONTROL SYSTEM OF FAST OSCILLATING CONTROL SYSTEMS *

Keywords: Averaging, control systems, small control, optimal control, Finsler geometry AMS subject classifications. 34C29, 34H05, 49J15, 93B11, 93C15, 93C70, 53B40

For control systems that either have a fast explicit periodic dependence on time and bounded controls or have periodic solutions and small controls, we define an average control system that takes into account all possible variations of the control, and prove that its solutions approximate all solutions of the oscillating system as oscillations go faster.

The dimension of its velocity set is characterized geometrically. When it is maximum the average system defines a Finsler metric, not twice differentiable in general. For minimum time control, this average system allows one to give a rigorous proof that averaging the Hamiltonian given by the maximum principle is a valid approximation.

1. Introduction. We consider either a "fast-oscillating control system" (1):

ẋ = u 1 X 1 ( t ε , x) + • • • + u m X m ( t ε , x) , u ≤ 1 ,
where all X i 's are 2π-periodic with respect to t/ε, or a "Kepler control system" (46):

ξ = f 0 (ξ) + v 1 f 1 (ξ) + • • • + v m f m (ξ) , v ≤ ε
where all solutions of ξ = f 0 (ξ) are periodic.

Averaging techniques for conservative -periodic or not-ordinary differential equations (ODEs) date back at least to H. Poincaré; see [2, §52] or [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF] for recent expositions. Roughly speaking, on a fixed interval, the solutions of ẋ = F (t/ε, x) differ from those of ẋ = F (x) by a term of order ε, with F the average of F with respect to its first argument.

If u or v above is assigned to be a fixed function of state and time (or computed from additional state variables as in u = α(p, x), ṗ = g(p, x)), then these techniques for ODEs can be applied to give an approximation at first order with respect to small ε of the movement of the slow variables. Averaging is usually used in this way in control theory: in vibrational control [START_REF] Meerkov | Principle of vibrational control: theory and applications[END_REF], fast oscillating controls are designed and averaging techniques allows analysis and proof of stability; in the same way, it solves stability and path planning questions in control of mechanical systems, see for instance [START_REF] Bullo | Averaging and vibrational control of mechanical systems[END_REF]; in [12, §5], high frequency control is used to approach a non-flat system by a flat one; one may also mention many applications to control [START_REF] Murray | Nilpotent bases for a class of nonintegrable distributions with applications to trajectory generation for nonholonomic systems[END_REF][START_REF] Liu | Averaging theorems for highly oscillatory differential equations and iterated Lie brackets[END_REF][START_REF] Morin | Design of homogeneous time-varying stabilizing control laws for driftless controllable systems via oscillatory approximation of lie brackets in closed loop[END_REF] of the work [START_REF] Kurzweil | Iterated Lie brackets in limit processes in ordinary differential equations[END_REF] that mimics Lie brackets by highly oscillatory controls along the original vector fields. A common feature to these references is that the use of oscillations "creates" new independent controls used for design. The use of averaging in optimal control of oscillating systems [START_REF] Chaplais | Averaging and deterministic optimal control[END_REF][START_REF] Geffroy | Généralisation des techniques de moyennation en contrôle optimal -Application aux problèmes de rendez-vous orbitaux en poussée faible[END_REF][START_REF] Geffroy | Optimal low-thrust transfers with constraints-generalization of averaging technics[END_REF][START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF] is similar in spirit to the above, but closer to the framework of this paper because oscillations are present in the system instead of being introduced by the control. Very interesting results are obtained applying averaging to the Hamiltonian equations arising from Pontryagin Maximum principle. For instance, in [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF], the authors have studied in this way the problem of minimal energy transfer between two elliptic orbits; extremals are the same as those giving the geodesics of a Riemannian metric. Again, averaging introduces "new independent controls": Riemannian geodesics are minimizers of a problem where all velocity directions are allowed whereas the velocity set of the original system at each point had positive codimension. The same averaging computation may be applied to the Hamiltonian differential equation obtained for minimum time, but, since this differential equation is discontinuous, there is no theoretical justification for averaging in that case.

Our contribution is to introduce a different way of averaging that takes into account all possible variations of the control -hence the control strategy can be decided after performing averaging-and to prove that it has satisfying regularity properties and is a good first order approximation of the above systems as ε → 0. This gives, as a side result, a justification of the use of averaging for minimum time in [START_REF] Geffroy | Généralisation des techniques de moyennation en contrôle optimal -Application aux problèmes de rendez-vous orbitaux en poussée faible[END_REF][START_REF] Geffroy | Optimal low-thrust transfers with constraints-generalization of averaging technics[END_REF]. This procedure also "creates new independent control", i.e. increases the dimension of the velocity set, that we characterize in terms of the original vector fields. When this dimension is maximum, the average system defines a Finsler metric [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF] on the manifold, whose geodesics are the limits of minimum time trajectories for the original systems as ε → 0. This Finsler metric is in general not twice differentiable (hence it is not a Finsler metric in the sense of [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF], indeed); we however prove that, at least in the less degenerate case, the Hamiltonian system governing extremals, although it is not locally Lipschitz, generates a flow on the cotangent bundle. Low thrust planar orbit transfer belongs to this less degenerate case.

The average control system may be used for other purposes than optimal control, for instance [START_REF] Bombrun | Les Transferts Orbitaux à Faible Poussée : Optimalité et Feedback[END_REF] designs a Lyapunov function for feedback control in the average system and uses it for the oscillating systems; indeed the present work was developed out of comparing feedback control based on a priori chosen Lyapunov functions with minimum time control for low thrust orbital transfer.

Preliminary versions of this paper can be found in [START_REF] Bombrun | On the averaged control system[END_REF][START_REF] Bombrun | Les Transferts Orbitaux à Faible Poussée : Optimalité et Feedback[END_REF]. It is organized as follows: the construction and results are developed for "fast-oscillating control system" in §3 and then transferred in §4 to "Kepler control systems", and applied to minimum time orbit transfer in the planar 2-body problem in §5.

2. Notations and conventions. 2.1. M is a smooth connected manifold of dimension n; its tangent and cotangent bundles are denoted by TM and T * M . One may assume for simplicity

M = R n , TM = R n × R n , T * M = R n × (R n ) * , and, for x ∈ M , T x M = R n , T * x M = (R n ) * . For v ∈ T x M , p ∈ T *
x M (or any v, p taken in a vector space and its dual), we denote by p, v (rather than p(v)) their duality product.

2.2. If E is a subset of a vector space V , then E ⊥ is its annihilator, the vector subspace of its dual V * made of all p's such that p, v = 0 for all v in E.

2.3. We assume that M is endowed with an arbitrary Riemannian distance d. If M = R n , just choose the canonical Euclidean distance.

In local coordinates, . and (. |. ) stand for the canonical Euclidean norm and scalar product. On a compact coordinate chart, k 1 x -y ≤ d(x, y) ≤ k 2 x -y for some positive k 1 , k 2 (Lipschitz equivalence). We also denote operator norms by . . 2.4. S 1 is R/2πZ. For θ in S 1 (an angle), we denote by µ(θ) the unique real number in [0, 2π) such that µ(θ) ≡ θ mod 2π. For a real number s ∈ R, we denote the angle it represents by s mod 2π; it belongs to the quotient S 1 .

Maps S1 → E (arbitrary set) are identified with 2π-periodic maps R → E. For instance, if f is such a map S 1 → E and τ ∈ R, we write f (τ ) instead of f (τ mod 2π); the average of f is denoted by

1 2π 2π 0 f (θ)dθ, or 1 2π π -π f (θ)dθ, or 1 2π θ∈S 1 f (θ)dθ; one identifies L p (S 1 , R m ) with the subset of L p (R, R m ) made of 2π-periodic functions.
2.5. The Euclidean norm in R m or (R m ) * is denoted by . , and the ball of radius one centered at the origin by B m . We view an element of R m as m × 1 matrix (column) of real numbers and an element of (R m ) * as a 1 × m matrix (line); transposition, denoted . ⊤ , sends R m to (R m ) * and vice-versa.

3. Fast oscillating control systems. We call fast oscillating control system on M a family of non-autonomous systems, linear in the control u ∈ R m :

ẋ = G( t ε , x) u = m i=1 G i ( t ε , x) u i , u ≤ 1 (1)
indexed by a positive number ε. Each G i is a smooth "periodic time-varying" vector field:

G i ∈ C ∞ (S 1 ×M, TM ). An admissible control is a measurable u(.) : [0, T ] → B m
for some T > 0. For a given control u(.) and initial condition x(0), there is a unique solution x(.), defined either on [0, T ] or only on a maximal interval [0, T ′ ), T ′ < T . Remark 3.1. Apart from being a notation defined by the double equality in (1),

G(θ, x) defines a linear map R m → T x M that sends (u 1 , . . . , u m ) ⊤ to m i=1 G i ( t ε , x)u i . 3.1. Average control system of fast oscillating control systems. Define the map G : M × L ∞ ([0, 2π], R m ) → TM by G(x, U) = 1 2π 2π 0 G(θ, x) U(θ) dθ . (2)
It allows one to define, for all x ∈ M , the subset E(x) ⊂ T x M by

E(x) = G(x, U), U ∈ L ∞ ([0, 2π], R m ) , U ∞ ≤ 1 ⊂ T x M , (3)
and the average control system of (1) as follows 1 .

Definition 3.2. The average control system of (1) is the differential inclusion

ẋ ∈ E(x). (4)
A solution of ( 4) is an absolutely continuous x(.)

: [0, T ] → M such that ẋ(t) ∈ E(x(t))
for almost all t.

Proposition 3.3. For all x in M , E(x) is convex, compact and symmetric with respect to the origin.

Proof. It is closed, convex and symmetric because it is the image of the unit ball of L ∞ S 1 , R m by a linear map; it is compact because G(x, .) is bounded on S 1 .

Further characterizations of E(x) use the map

H : T * M → [0, +∞) defined by H(x, p) = 1 2π 2π 0 p, G(θ, x) dθ (5)
where the Euclidean norm is used according to §2.5 and, for each

(θ, x), p, G(θ, x) = ( p, G 1 (θ, x) , . . . , p, G m (θ, x) ) ∈ (R m ) * . ( 6
)
Proposition 3.4. For all (x, p) ∈ T * M , one has, with H defined in (5),

E(x) = v ∈ T x M , sup p∈T * x M H(x,p)≤1 p, v ≤ 1 , (7)
H(x, p) = sup v∈E(x) p, v = sup U∈L ∞ (S 1 ,R m ), U ∞ ≤1 p, G(x, U) = p, G(x, U * p,x ) , ( 8
)
with U * p,x ∈ L ∞ S 1 , R m defined by: U * p,x (θ) = 0 if p, G(θ, x) = 0 , p,G(θ,x) ⊤ p,G(θ,x) if p, G(θ, x) = 0 . (9)
Proof. The last equality in ( 8) is a straightforward maximization, the second one comes from the definition (3) of E(x) and a simple computation yields H(x, p) = p, G(x, U * p,x ) ; this proves [START_REF] Bullo | Averaging and vibrational control of mechanical systems[END_REF]. Being closed and convex, E(x) is the intersection of all its supporting half-spaces [START_REF] Schneider | Convex bodies: the Brunn-Minkowski theory[END_REF]Corollary 1.3.5]; according to [START_REF] Bullo | Averaging and vibrational control of mechanical systems[END_REF], this yields the following relation, equivalent to [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF]:

E(x) = p∈T * x M {v ∈ T x M , p, v ≤ H(x, p)}.
A convenient characterization of solutions of (4). According to Definition 3.2, a solution x(.) is such that, for almost all t, there is

U(t) ∈ L ∞ ([0, 2π], R m ) such that ẋ(t) = G(x(t), U(t)); the map (t, θ) → U(t)(θ) is measurable with respect to θ only.
It turns out that it may always be chosen jointly measurable with respect to (t, θ) according to the following "measurable selection" result:

Proposition 3.5. A map x : [0, T ] → R n is a solution of the differential inclu- sion (4) if and only if there exists u ∈ L ∞ ([0, T ] × S 1 , R m ), u ∞ ≤ 1 such that ẋ(t) = 1 2π 2π 0 G(θ, x(t)) u(t, θ) dθ (10)
for almost all t in [0, T ].

Proof. After possibly partitioning [0, T ] into intervals where ẋ(t) remains in the same coordinate chart, we work in coordinates and use a Euclidean norm when useful.

Sufficiency is clear: from Fubini theorem, θ → u(t, θ) is measurable for almost all t, hence x(.) is a solution of (4). Conversely, let x(.) be a solution of (4): ẋ(.) is measurable and, for almost all t, there exists ũt

∈ L ∞ (S 1 , R m ), ũt ∞ ≤ 1 such that ẋ(t) = G(x(t), ũt ) = 1 2π 2π 0 G(s 1 , x(t))ũ t (s 1 )ds 1 . ( 11
) Let φ : L ∞ [0, T ] × S 1 , R m → L 2 ([0, T ], R n ) be the linear map defined by φ(u)(t) = G(x(t), u(t, .)) = 1 2π 2π 0 G(s 1 , x(t)) u(t, s 1 ) ds 1
and I the image by φ of the unit ball of L ∞ ([0, T ]×S 1 , R m ). Since, by [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF], ẋ(.) is essentially bounded, it is in L 2 ([0, T ], R n ); since I is closed and convex in that Hilbert space, the distance from ẋ to I is reached for a unique element ξ ∈ I:

ξ = φ(ū) , ū ∈ L ∞ [0, T ] × S 1 , R m , ū L ∞ ≤ 1 .
Let us prove by contradiction that ξ = ẋ, i.e. ẋ(.) ∈ I; this will end the proof. If ẋ = ξ, one has, for all u in the unit ball of

L ∞ [0, T ] × S 1 , R m , ẋ -ξ φ(u) -φ(ū) L 2 ≤ 0 (12)
with equality only if φ(u) = φ(ū). Define u by u(t, s) = U * ( ẋ(t)-ξ(t)) ⊤ , x(t) (s) with U * p,x defined by [START_REF] Caillau | Coplanar control of a satellite around the Earth[END_REF]; clearly, u is in the unit ball of L ∞ [0, T ] × S 1 , R m , and, for all

(t, s) ∈ [0, T ] × S 1 and all u ∈ R m , u ≤ 1 ⇒ ẋ(t) -ξ(t) ⊤ G(s, x(t)) ( u(t, s) -u) ≥ 0 , (13)
hence ẋ(t) -ξ(t) ⊤ G(s 1 , x(t)) ( u(t, s 1 ) -ū(t, s 1 )) is non-negative for almost all (t, s 1 )

and, since it is the integrand of the left-hand side of ( 12), it must be zero; hence ξ = φ(ū) = φ( u) and ξ(t) = G(x(t), u(t, .)) for almost all t.

In [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF], ũt satisfies ũt (s 1 ) ≤ 1 for almost all s 1 , hence, according to [START_REF] Geffroy | Généralisation des techniques de moyennation en contrôle optimal -Application aux problèmes de rendez-vous orbitaux en poussée faible[END_REF],

ẋ(t) -ξ(t) ⊤ G(s 1 , x(t)) ( u(t, s 1 ) -ũt (s 1 )) ≥ 0 . Since ẋ(t) = G(x(t), ũt ), ξ(t) = G(x(t), u(t, .
)), the integration with respect to the variable s 1 yields -ẋ(t) -ξ(t) 2 ≥ 0 for almost all t; this contradicts ẋ = ξ.

Remark 3.6. The differential inclusion (4) is equivalent to the "control system"

ẋ = G(x, U) , U ∈ L ∞ (S 1 , R m ) , U ∞ ≤ 1
where, by Proposition 3.5, admissible controls are maps t → U(t) such that u : (t, θ) → U(t)(θ) is measurable with respect to (t, θ). Since this "control" is infinite dimensional, and we could not find a representation of the type ẋ = f (x, v), v ∈ U ⊂ R r , r finite, we stay with the differential inclusion (4), with E(x) described by ( 8) and (5).

Convergence theorem.

The following result relates solutions of the fast oscillating systems as ε tends to zero to solutions of the average system. To our knowledge, this kind of theorem where the control is not chosen prior to averaging has never been stated in the literature.

Theorem 3.7 (Convergence for fast-oscillating control systems). 1. Let x 0 (.) : [0, T ] → M be an arbitrary solution of (4). There exist a family of measurable functions u ε (.) : [0, T ] → B m , indexed by ε > 0, and positive constants c, ε 0 , such that, calling x ε (.) the solution of (1) with control u = u ε (t) and initial condition x ε (0) = x 0 (0), one has:

x ε (.) is defined on [0, T ] for all ε smaller than ε 0 and converges to x 0 (.) as ε → 0, with an error of uniform order ε:

d (x ε (t), x 0 (t)) < c ε , t ∈ [0, T ] , 0 < ε < ε 0 . (14)
2. Let K be a compact subset of M , (ε n ) n∈N a decreasing sequence of positive real numbers converging to zero, and, for each n, x n (.) : [0, T ] → K a solution of (1) with ε = ε n and for some control

u = u n (t), u n (.) ∈ L ∞ ([0, T ], R m ), u n (.) ∞ ≤ 1.
Then the sequence x n (.) n∈N is compact for the topology of uniform convergence on [0, T ] and any accumulation point is a solution of the average system (4).

The statement is more complex than the one for ODEs, e.g. [2, §52.C], due to underdetermination (choice of control in (1), multi-valued right-hand side in (4)).

Informally, "1" states that any solution of the average system is the limit of solutions of fast oscillating systems with well chosen controls and "2" states that, conversely, any limit of solutions of oscillating systems, with arbitrary controls, is a solution of the average control system. There is an estimate on the error in "1" but not in "2" because some sequences may converge slower than others. Remark 3.8. One may consider systems that are affine instead of linear in the control by adding a drift vector field G 0 (t/ε, x) to (1). Then, in the average control system, E(x) is replaced by G 0 (x) + E(x), with G 0 (x) = 1 2π 2π 0 G 0 (θ, x)dθ. By a straightforward extension, convergence does hold for these systems too.

In the proof of Theorem 3.7, the following technical lemma is needed. Lemma 3.9. Let ε > 0 and a < b be real numbers and u : [a -2πε, b] × S 1 → R m be measurable. One has the following identity (see §2.4 for the notation µ(.)) :

θ∈S 1 a≤s≤b G(θ, x(s)) u(s, θ) dθ ds = θ∈S 1 a≤s≤b G( s ε , x(s)) u(s + εµ(θ), s ε ) dθ ds + ∆ ε (15) with ∆ ε = T a ε G( s ε , x(s + εµ(θ))) u(s + εµ(θ), s ε ) dθ ds - T b ε G( s ε , x(s + εµ(θ))) u(s + εµ(θ), s ε ) dθ ds + θ∈S 1 a≤s≤b G( s ε , x(s + εµ(θ))) -G( s ε , x(s)) u(s + εµ(θ), s ε
) dθ ds [START_REF] Kailath | Linear systems. Information and System Sciences[END_REF] and the set T a ε defined by T a ε = {(s, θ), θ ∈ S 1 , a -ε µ(θ) ≤ s ≤ a} and T b ε accordingly. Proof. Thanks to the change of variables θ = τ /ε mod 2π, s = τ + ε µ(φ), with µ(θ) as defined in §2.4, the left-hand side of ( 15) is equal to

φ∈S 1 a-εφ≤τ ≤b-εµ(φ) G( τ ε , x(τ + εµ(φ))) u(τ + εµ(φ), τ ε ) dτ dφ .
Keeping the names (s, θ) instead of (τ, φ), one gets [START_REF] Hartman | Ordinary differential equations[END_REF], the correcting term ∆ ε coming from the modified domain of integration and argument of x.

Proof of Theorem 3.7, point 1. Consider a solution x 0 : [0, T ] → M n of (4). By Proposition 3.5 there exists

u 0 ∈ L ∞ ([0, T ] × S 1 , R m ), u 0 ∞ ≤ 1 satisfying (10). For ε > 0, define u ε (.) ∈ L ∞ ([0, T ], R m ) by (see §2.4 for notations): u ε (t) = 1 2π 2π 0 u 0 (t + ε µ(θ), t ε ) dθ , (17)
where u 0 is prolonged by zero outside [0, T ]:

u 0 (t + ε µ(θ), t ε ) = 0 if t + ε µ(θ) > T .
Let us prove that this construction of u ε satisfies the two announced properties.

Step 1. Let us first assume that M is an open subset of R n and G is zero outside a compact subset of M . Then G(θ, x) is a n × m matrix for all (θ, x) and, denoting by . the Euclidean norm for vectors and the operator norm for matrices, there are global constants Lip G and sup G such that, for all

x, x ′ , θ in M × M × S 1 , G(θ, x) -G(θ, x ′ ) ≤ Lip G x -x ′ , G(θ, x) ≤ sup G . ( 18
)
Let b be a non-negative constant and consider, for each ε > 0, a solution x ε (.) of ( 1) with control u = u ε (t) and initial condition x ε (0) such that

x ε (0) -x 0 (0) ≤ b ε . (19)
In fact b = 0 in the theorem itself, but we need a nonzero b in step 2. By definition, expanding u ε (s) as in [START_REF] Kurzweil | Iterated Lie brackets in limit processes in ordinary differential equations[END_REF] and using Lemma 3.9, one has

x ε (t) = x ε (0) + 1 2π t 0 2π 0 G( s ε , x ε (s)) u 0 (s -ε µ(θ), s ε ) dθ ds , = x ε (0) + 1 2π t 0 2π 0 G(θ, x ε (s)) u 0 (s, θ) dθ ds -∆ ε (20)
with ∆ ε given by ( 16), that satisfies

∆ ε ≤ 4π 2 Lip G (1 + T sup G) ε because, in particular, u 0 ≤ 1, |εµ(θ)| < 2πε and G( s ε , x ε (s)) -G( s ε , x ε (s + ε µ(θ))) u 0 (s + ε µ(θ), s ε ) ≤ 2π Lip G sup G ε .
Using ( 19), ( 20) the bound on ∆ ε and the relation

x 0 (t) = x 0 (0) + 1 2π t 0 2π 0 G(θ, x 0 (s)) u 0 (s, θ) dθ ds , one gets x ε (t) -x 0 (t) ≤ b + 2π Lip G 1 + T sup G ε + Lip G t 0 x ε (s) -x 0 (s) ds for all t in [0, T ],
and finally, by Gronwall lemma,

x ε (t) -x 0 (t) ≤ b + 2π Lip G (1 + T sup G) e T Lip G ε (21)
for all t in [0, T ] and ε in [0, ε 0 ]. This proves the theorem if M is an open subset of R n and G is zero outside a compact subset, with an explicit constant c corresponding to the distance d defined from the Euclidean norm and with ε 0 = +∞.

Step 2. General case. Let x ε (.) be the solution of (1) with control u = u ε (t) defined in [START_REF] Kurzweil | Iterated Lie brackets in limit processes in ordinary differential equations[END_REF] from u 0 and with initial condition x ε (0) = x 0 (0); it is not necessarily defined on [0, T ] but may have a maximum interval of definition [0, T ε ) with T ε < T . Let T ∈ [0, T ] be the supremum of the set of numbers τ ∈ [0, T ] such that, for some ε 0 and some c, that may depend on τ , the solution x ε (.) is defined on [0, τ ] and satisfies d(x ε (t), x 0 (0)) < c ε for all t ∈ [0, τ ] and ε ∈ [0, ε 0 ]. Let us prove by contradiction that T = T . This will end the proof of Theorem 3.7, point 1.

Assume T < T , and let -O be a coordinate neighborhood of x 0 ( T ),

-α > 0 be such that 0 < T -α < T + α ≤ T and x 0 ([ T -α, T + α]) ⊂ O, -c > 0, ε 0 > 0 be such that d(x ε (t), x 0 (t)) < c ε for all t ∈ [0, T -α] and ε ∈ [0, ε 0 ]. Taking ε 0 possibly smaller, one also has x ε ( T -α) ∈ O for ε < ε 0 . Let K be a compact neighborhood of x 0 ([ T -α, T + α]) contained in O, K ′ a compact neighborhood of K contained in O, and ρ : M → [0, 1] a smooth map, zero outside K ′ and constant equal to 1 in K. Defining G ρ by G ρ (θ, x) = ρ(x)G(θ, x), let us apply Step 1 in coordinates in O, with G ρ instead of G and [ T -α, T + α] instead of [0, T ]. Call x ρ 0 (resp. x ρ ε , ε > 0) the solution of (10) (resp. of (1) with control u = u ε (t)), replacing G by G ρ , with initial condition x ρ ε ( T -α) = x ε ( T -α), ε ≥ 0. One clearly has, as in (19), x ρ ε ( T -α) -x ρ 0 ( T -α) < b ε with b deduced from c via Lipschitz equivalence of
the distance d and the Euclidean norm in coordinates (see §2.3); then Step 1 provides ε ′ 0 > 0 such that, by [START_REF] Murray | Nilpotent bases for a class of nonintegrable distributions with applications to trajectory generation for nonholonomic systems[END_REF], the inequality

x ρ ε (t) -x ρ 0 (t) ≤ b + 2π Lip G ρ (1 + 2α sup G ρ ) e 2α Lip Gρ ε (22) is valid for t ∈ [ T -α, T + α] and ε ∈ [0, ε ′ 0 ]. Possibly choosing a smaller ε ′ 0 , this implies that x ε ([ T -α, T + α]) ⊂ K for ε < ε ′ 0 ;
since G coincides with G ρ in K, the conclusion holds for x ε and G as well as for x ρ ε and G ρ if ε is no larger than ε ′ 0 . We have shown that, for all ε < ε ′ 0 , the solution x ε is defined on [0, T + α] and satisfies d(x ε (t) -x 0 (t)) ≤ c ′ ε for t in [0, T + α] where c ′ is larger than c and than a bound deduced from [START_REF] Pontryagin | Théorie mathématique des processus optimaux[END_REF] and from Lipschitz equivalence of d with the Euclidean distance. This contradicts the definition of T .

Proof of Theorem 3.7, point 2. Since G is bounded on S 1 × K (one may cover K with a finite number of coordinate charts and define this bound in coordinates), the maps x n (.) have a common Lipschitz constant and the sequence (x n (.)) is equicontinuous, hence compact by Ascoli-Arzela Theorem: one may extract a uniformly convergent sub-sequence. Still denoting by (x n (.)) n∈N such a converging sub-sequence and by x * (.) its (uniform) limit, we need to prove that this limit is a solution of (4).

Define, for each n,

u n : [0, T ] × S 1 → R m by u n (t, θ) = u n ( β n (t, θ) ) , (23)
where

u n (.) ∈ L ∞ ([0, T ], R m ) is associated to x n (.)
according to the assumption of the theorem and where the map

β n : [0, T ] × S 1 → R is defined by t -2πε n < β n (t, θ) ≤ t , β n (t, θ) ε n ≡ θ mod 2π . ( 24
)
Clearly

u n is in L ∞ ([0, T ] × S 1 , R m ) and u n ∞ ≤ 1.
Hence, after possibly extracting a sub-sequence, ( u n ) converges in the weak- * topology to some u * . Let us prove that, for almost all t ∈ [0, T ],

ẋ * (t) = 1 2π 2π 0 G(θ, x * (t)) u * (t, θ)dθ . ( 25
)
Let T ∈ [0, T ] be the supremum of the set of numbers τ ∈ [0, T ] such that this is true for almost all t in [0, τ ], and let us prove by contradiction that T = T . Assume T < T , and let O be a coordinate neighborhood of x 0 ( T ) and α be such that 0

< T -α < T + α ≤ T and x 0 ([ T -α, T + α]) ⊂ O. Uniform convergence implies x n ([ T -α, T +α]) ⊂ O for n large enough and then, in coordinates, for t ∈ [ T -α, T +α], x n (t) -x n ( T -α) = t T -α G( s ε n , x n (s)) u n (s) ds . (26)
From ( 24), one has β n (s + ε n θ, s εn ) = s, hence, from [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF], u n (s + ε n θ, s εn ) = u n (s) for all θ ∈ S 1 , s ∈ R; using this in Lemma 3.9, one has

1 2π T -α≤s≤t θ∈S 1 G(θ, x n (s)) u n (s, θ) dθ ds = 1 2π T -α≤s≤t θ∈S 1 G( s ε n , x n (s)) u n (s) dθ ds + ∆ εn .
Since the integral in the right-hand side -whose integrand does not depend on θis equal to the right-hand side of (26), one gets, using uniform convergence of x n to x * , weak convergence of u n to u * and convergence of ∆ εn to zero,

x * (t) -x * ( T -α) = 1 2π t T -α 2π 0 G(θ, x * (s)) u * (s, θ) dθ ds , for t in [ T -α, T + α],
and finally that (25) hold for almost all t in [0, T + α], thus contradicting the definition of T .

Dimension of the velocity set E(x).

Recall that, for a convex subset C of a linear space, containing the origin, its linear hull is the smallest linear subspace that contains C, the interior of C in its linear hull is always nonempty, and dim C is the dimension of this linear hull.

Viewing ∂ j G ∂θ j (θ, x) as a linear map R m → T x M (see Remark 3.1), and denoting by Σ a sum of linear subspaces of T x M , define the integer r(θ, x) by:

r(θ, x) = dim j∈N Range ∂ j G ∂θ j (θ, x) . ( 27
)
It is also the rank of the collection of vectors

∂ j Gi ∂θ j (θ, x) ∈ T x M , 1 ≤ i ≤ m, j ≥ 0.
In the following proposition, and it is the sole place where this property is used, "system (1) is real analytic with respect to θ" means that the vector fields G i are real analytic with respect to θ for fixed x (while being smooth with respect to (θ, x)).

Proposition 3.10.

1. The linear hull of E(x) satisfies the following two properties for all x in M , where the inclusion ( 29) is an equality if (1) is real analytic with respect to θ:

Linear hull E(x) = θ∈S 1 Range G(θ, x) , (28)
Linear hull E(x) ⊃ j∈N Range ∂ j G ∂θ j (θ, x) for all θ ∈ S 1 . ( 29
)
2. If r(θ, x) = n for at least one θ in S 1 , then E(x) has a nonempty interior in T x M , i.e. dim E(x) = n.

3. If the system (1) is real analytic with respect to θ, then r(θ, x) does not depend on θ and r(θ, x) = dim E(x).

Proof.

If p is in Range G(θ, x) ⊥ for all θ, then any v = G(x, U) in E(x) sat- isfies p, v = 0 because p, G(θ, x)U(θ) is identically zero on [0, 2π]. Conversely, let p be in E(x) ⊥ , and consider v = G(x, U * p,x ) ∈ E(x); then p, v = 0 implies p, G(θ, x) = 0, i.e. p ∈ Range G(θ, x)
⊥ for all θ. We have proved the identity

E(x) ⊥ = θ∈S 1 (Range G(θ, x)) ⊥ , hence (28). If p is in Range G(φ, x)
⊥ for all φ, differentiating p , G(φ, x) = 0 with respect to φ yields p , ∂ j G/∂φ j (x, φ) = 0, j ∈ N; we have proved, taking φ = θ, the inclusion

φ∈S 1 (Range G(φ, x)) ⊥ ⊂ j∈N Range ∂ j G ∂θ j (θ, x)
⊥ , hence (29). To prove the reverse inclusion in the real analytic case, fix θ ∈ S 1 and p ∈ j∈N

∂ j G ∂θ j (θ, x) ⊥
, and consider the real analytic mapping S 1 → (R m ) * , φ → p , G(φ, x) ; the assumption on p implies that this map vanishes for φ = θ, as well as its derivatives at all orders, hence it is identically zero: p ∈ φ∈S 1 (Range G(φ, x)) ⊥ . This ends the proof of Point 1. Point 2 is an easy consequence and Point 3 is classical.

3.4. Further properties in the full rank case. We now assume that the mapping G in (1) is such that the rank r(θ, x) defined by ( 27) is maximal: r(θ, x) = n for all x in M and θ in S 1 .

(30)

3.4.1. Controllability. Condition (30) is strongly related to controllability of the linear approximation of (1) around equilibria, i.e. around solutions where x is constant and u is identically zero. Indeed, the linear approximation of the timevarying nonlinear system (take ε = 1 in (1)):

ẋ = G(t, x)u (31)
around the equilibrium x = x 1 is the time-varying linear system ξ = G(t, x 1 )u; according to [16, p.614], it is "controllable with impulsive controls at any time" if and only if r(t, x 1 ) = n for all t. If this is true at all points x 1 then all end-point mappings are submersions around zero controls; we shall need the following more precise result:

Proposition 3.11. Assume that (30) holds.

1. For all x 1 ∈ M and T > 0, there exist a coordinate neighborhood W of x 1 (the ball B below refers to the Euclidean norm in these coordinates), positive constants α 0 , c 3 , and, for all y ∈ W, a smooth map χ y : B(y, α 0 ) → L ∞ ([0, T ], R m ) with Lipschitz constant c 3 , which is a right inverse of the end-point mapping of (31) on [0, T ] starting from y, i.e. for all y f ∈ B(y, α 0 ), the control χ y (y

f ) : [0, T ] → R m is such that the solution of ẋ = G(t, x)χ y (y f )(t), x(0) = y satisfies x(T ) = y f .
2. For all ε > 0, the system (1) is fully controllable, i.e. there exists, for any ε > 0 and any two point x 0 , x 1 in M , a time T and a measurable control u : [0, T ] → B m such that the solution of (1) with x(0) = x 0 satisfies x(T ) = x 1 .

Proof. Let E y : L ∞ ([0, T ], R m ) → M be the end-point mapping with starting point y. Condition (30) implies that the derivative of E x1 at the zero control has rank n; hence there exists an n-dimensional subspace V of L ∞ ([0, T ], R m ) such that the restriction of E x1 , and hence of E y for y close enough, to V is a local diffeomorphism at zero; the χ y 's are the local inverses of these local diffeomorphisms; they depend smoothly on y, hence the common α 0 and c 3 in Point 1.

This implies that the reachable set from any point at any positive or negative time contains a neighborhood of this point; a classical argument then tells us that the reachable set from a point x 0 is M , assumed to be connected, for it is both open (obvious) and closed (if x is in the closure of the reachable set, some points in the reachable set can be reached in negative time, hence x can be reached from x 0 ).

Let us now turn to the average system [START_REF] Bombrun | Les Transferts Orbitaux à Faible Poussée : Optimalité et Feedback[END_REF]. From H : T * M → [0, +∞) defined by (5), we define N : TM → [0, +∞] by

N (x, v) = max p∈T * x M, H(x,p)≤1 p, v . (32)
Proposition 3.12. Assume the rank condition (30).

1. For all x ∈ M , H(x, .) defines a norm on the cotangent space T * x M , its dual norm on the tangent space T x M is N (x, .), and E(x) is the unit ball for N (x, .), i.e.

E(x) = {v ∈ T x M, N (x, v) ≤ 1}.
2. System (4) is fully controllable, i.e. there exists, for any points x 0 , x 1 in M , a time T and a solution x(.) : [0, T ] → M of (4) such that x(0) = x 0 , x(T ) = x 1 .

Proof. From (5), H(x, p) = 0 implies p, G(θ, x) = 0 for all θ and, differentiating with respect to θ and using (30), this implies p = 0; this makes p → H(x, p) a norm, the other properties being straightforward. Hence N given by (32) is finite for any (x, v) and it is, by definition, the dual norm of H(x, .); E(x) is its unit ball by [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF] in Proposition 3.4. To prove Point 2, take a continuously differentiable curve γ : [0, 1] → M such that γ(0) = x 0 and γ(1) = x 1 and σ : [0, T ] → [0, 1] for some T > 0, differentiable, such that

t ≥ σ(t) 0 N (γ(s), dγ ds (s))ds
(N and H are obviously continuous), then t → x(t) = γ(σ(t)) is a solution of (4) such that x(0) = x 0 and x(T ) = x 1 .

3.4.2. On the differentiability of H. It is clear that H, given by ( 5), is as smooth as G on T * M \ Z with

Z = {(x, p) ∈ T * M, ∃θ ∈ S 1 , p, G(θ, x) = 0} . ( 33
)
Unfortunately, Z is not empty in general: it is generically a 2n -m + 1 dimensional submanifold of T * M . One however has the following result, valid also at these points.

Theorem 3.13. If condition (30) holds, H 2 is continuously differentiable.

It is stated for H 2 : (x, p) → H(x, p) 2 , because H itself, homogeneous of degree 1 with respect to p, cannot be differentiable on {p = 0}, that coincides with {H(x, p) = 0} by Proposition 3.12 item 1.

The map H fails in general to be twice differentiable on Z. We have the following estimate of the of the modulus of continuity of its first derivative, that even fails to be Lipschitz continuous. Its main consequence is Theorem 3.21.

Theorem 3.14. Assume that the rank condition (30) holds and that (i) for (x, p) ∈ T * M , p = 0, there is at most one θ ∈ S 1 such that p, G(θ, x) = 0, and p, ∂G ∂θ (θ, x) does not vanish at the same point, (ii) for all (θ, x) ∈ S 1 × M , one has rank G(θ, x) = m, then any point (x, p) has a constant c and a coordinate neighborhood in T * M such that for all X and Y in R 2n , coordinates of points in the neighborhood,

dH(Y ) -dH(X) ≤ c X -Y ln 1 X -Y . (34)
In the left-hand side, . stands for the operator norm in coordinates, see §2.3.

Remark 3.15 (Finsler geometry). If H 2 was an least twice continuously differentiable, with a positive definite Hessian with respect to p, so would be N 2 (see (32)), and it would define a (reversible) Finsler metric [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF] on M . The lack of differentiability calls for further developments.

Before proving these theorems, we state a more generic result, whose notations are totally independent from the rest of the paper. Its proof is in the appendix.

H(X) = 1 2π 2π 0 V(θ, X) dθ , Z = {X ∈ O d , ∃θ ∈ S 1 , V(θ, X) = 0} . (35) 1. Assume that, for all X in O d , the set {θ ∈ S 1 , V(θ, X) = 0} has measure zero in S 1 . ( 36
)
Then H is continuously differentiable and, for all X, dH(X).h = 1 2π

2π 0 ∂V ∂X (θ, X).h V(θ, X) V(θ, X) dθ . (37)
2. Let V satisfy the following assumptions:

(a) V(θ, X) = 0 for all (θ, X) ∈ S 1 × O d such that V(θ, X) = 0, (b) for any X ∈ O d
, there is at most one θ such that V(θ, X) = 0, (c) V and ∂V/∂θ do not vanish simultaneously,

   ( 38
)
and let X be in Z. There is a neighborhood U of X in O d , and a constant K > 0 such that, for all X, Y in U ,

dH(X) -dH(Y ) ≤ K X -Y ln 1 X -Y . ( 39
)
In the left-hand side of (39), . stands for the operator norm, see §2.3.

Proof of Theorem 3.13. This is a local property. We operate in coordinates and apply Proposition 3.16 (Point 1) with d = 2n, O d a neighborhood of a point where p = 0, X = (x, p) ∈ R 2n and V(θ, X) = p, G(θ, x) . The rank condition implies that derivatives of all orders of the map θ → V(θ, X) never vanish at the same point, so that its zeroes are isolated and the set {θ ∈ S 1 , V(θ, X) = 0} is finite and a fortiori has measure zero; hence H is continuously differentiable outside {p = 0}. Since 0 ≤ H(x, p) ≤ k p for some local constant k the derivative of H 2 is zero at all points (x, 0) and, since (37) implies that the norm of dH(x, p) at neighboring points where p = 0 is bounded, the derivative of H 2 at these points tends to zero as p → 0. H 2 is therefore continuously differentiable everywhere.

Proof of Theorem 3.14. Smoothness outside Z is obvious from the expression (5) of H; inequality (34) is a consequence of Proposition 3.16 (Point 2), applied with d = 2n, X = (x, p) ∈ R 2n , V(θ, X) = p, G(θ, x) and O d a neighborhood of a point of Z \ {p = 0}; it is clear that points (i) and (ii) imply the three conditions (38).

3.5. Application to the minimum time problem. Fix two points x 0 , x 1 in M and consider the time optimal problem associated to (1) for ε > 0:

(P ε ), ε > 0 : ẋ(t) = G(t/ε, x(t))u(t), u(t) ∈ B m , t ∈ [0, T ], x(0) = x 0 , x(T ) = x 1 min T , (40)
and the time optimal problem associated to the average system:

(P 0 ) : ẋ(t) ∈ E(x(t)), t ∈ [0, T ], x(0) = x 0 , x(T ) = x 1 min T . (41)
Call T ε (x 0 , x 1 ) the minimum time for (P ε ), ε > 0 and T 0 (x 0 , x 1 ) the one for (P 0 ); when no confusion arises, we write T ε and T 0 . Let us develop (40)-(41): concerning (40), T ε is the infimum of the set of T 's such that there is an admissible control u(.): [0, T ] → B m , and x(.): [0, T ] → M satisfying x(0) = x 0 , x(T ) = x 1 and ẋ(t) = G(t/ε, x(t))u(t) for almost all t; Proposition 3.11, point 2 implies that this set is nonempty, hence T ε is finite. Concerning (41), T 0 is the infimum of the set of T 's such that there is x(.) : [0, T ] → M satisfying x(0) = x 0 , x(T ) = x 1 and ẋ(t) ∈ E(x(t)) for almost all t, T 0 is finite from Proposition 3.12, point 2. A solution to (P ε ) (resp. to (P 0 )) is x(.), u(.) (resp. x(.)) as above with T = T ε (resp. T = T 0 ). In general, the minimum T ε or T 0 need not be reached, i.e. there need not be a solution.

Lemma 3.17. Assume the rank condition (30). 1. There is a neighborhood W of any x 1 and two constants α 0 > 0 and C 3 > 0 such that, for all y in W, T ε (y, x 1 ) ≤ 2πε + C 3 d(x 1 , y).

2. For any x 0 , x ′ 1 , x 1 in M , one has T ε (x 0 , x 1 ) ≤ T ε (x 0 , x ′ 1 ) + T ε (x ′ 1 , x 1 ) + 2πε. Proof. Apply Proposition 3.11, point 1 with T = 2π, using as a distance in W the Euclidean norm in some coordinates: for any two points y, y ′ in W such that y -y ′ ≤ α 0 , there is a control defined on [0, 2π], with L ∞ norm smaller than c 3 y -y ′ that brings y ′ at time 0 to y at time 2π for system (31); rescaling time and control by ε yields, if c 3 y -y ′ ≤ ε, a control with L ∞ norm less than 1 that brings y ′ at time 0 to y at time 2πε for system (1) and hence, by concatenating controls and using periodicity of G, for any positive integer k, a control with L ∞ norm less than 1 that brings y ′ at time 0 to y at time 2kπε for system (1) if c 3 y -y ′ ≤ kε. In other words, T ε (y ′ , y) ≤ 2π(ε + c 3 y -y ′ ). Take y ′ = x 1 and 2πc 3 /C 3 the ratio between the Euclidean norm and the distance d; this proves point 1. Point 2 follows from using periodicity of G and concatenating controls while inserting a zero control between time T ε (y ′ , y) and the next multiple of 2π. Theorem 3.18 (limit of minimum time). Assume the rank condition (30).

1. T ε is bounded as ε → 0 and lim sup ε→0 T ε ≤ T 0 .

2. If, for ε > 0 small enough, each (P ε ) has a solution x ε : [0, T ε ] → M and there exists a compact K ⊂ M such that x ε ([0, T ε ]) ⊂ K for all ε > 0 small enough, then all accumulation points of the compact family (x ε (.)) ε>0 in C 0 ([0, T 0 ], M ) are solutions of (P 0 ) and lim ε→0 T ε = T 0 .

Proof. Consider a minimizing sequence for problem (P 0 ), i.e. solutions x k : [0, T 0 + β k ] → M of the average system (4) with (β k ) a sequence of positive numbers that tends to zero and x k (0) = x 0 , x k (T 0 + β k ) = x 1 for all k. For each x k , there is, according to Theorem 3.7, a family ( x k ε (.) ) ε>0 such that each x k ε (.) is a solution of (1) with x k ε (0) = x 0 and d(x k ε (t),

x k (t)) ≤ c 1 ε for all t in [0, T 0 + β k ]. In particular d(x k ε (T 0 + β k ), x k (t) 1 ) ≤ c 1 ε. Now, from Lemma 3.17, T ε ( x k ε (T 0 + β k ) , x 1 ) ≤ (2π + c 1 C 3 )ε; hence, from the second point of that lemma (with x ′ 1 = x k ε (T 0 + β k ) ), one has T ε = T ε (x 0 , x 1 ) ≤ T 0 + β k + (4π + c 1 C 3 )ε and, letting k go to infinity, T ε ≤ T 0 + (4π + c 1 C 3 )ε; this implies Point 1. Let us turn to point 2.
Extend x ε on [0, T ], with T an upperbound of T ε , by taking x ε (t) = x 1 for t in [T ε , T ]. Any sequence (x ε k (.)) k∈N with lim ε k = 0 is compact in C 0 ([0, T ], M ): take a convergent subsequence such that T ε k also converges to some T * . The uniform limit goes through x 0 at time 0 and x 1 at time T * and is, by Theorem 3.7, a solution of the average system (4), hence T * ≥ T 0 by definition of T 0 . This, together with Point 1, implies Point 2 because T * can be any accumulation point of (T ε ) as ε → 0.

Let us now write the Pontryagin Maximum Principle [START_REF] Pontryagin | Théorie mathématique des processus optimaux[END_REF] both for (P ε ), ε > 0 and for (P 0 ) and see how they are related.

The extremals of problem (P ε ), ε > 0, are absolutely continuous maps t → (x(t), p(t)) solution to

ṗ = - ∂H ε ∂x , ẋ = ∂H ε ∂p with H ε (t, p, x) = p, G(t/ε, x) , (42)
whose right-hand side is discontinuous on S ε = {(x, p, t), p, G(t/ε, x) = 0} (the "switching surface"), where it is in fact not defined.

The extremals of (P 0 ) are absolutely continuous t → (x(t), p(t)) solution to

ṗ = - ∂H ∂x , ẋ = ∂H ∂p . ( 43
)
with H given by ( 5). The right-hand sides are continuous according to Theorem 3.13.

Theorem 3.19. If an absolutely continuous map t → x(t) defined on [0, T ] is a solution of (P ε ), ε > 0, (resp. of (P 0 )), then there exists t → p(t) defined on [0, T ] such that t → (p(t), x(t)) is an extremal of (P ε ), ε > 0 (resp. of (P 0 )).

Proof. Problem (P ε ), ε > 0 deals with a classical smooth control system; according to [START_REF] Pontryagin | Théorie mathématique des processus optimaux[END_REF][START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF], the pseudo-Hamiltonian is h(t, x, p, u) = p, G(t/ε, x) u ; an extremal is a curve on the co-tangent bundle solution, in local coordinates, of:

ṗ = -∂h ∂x (t, x, p, u * ) = -p, ∂G ∂x u * , ẋ = ∂h ∂p (t, x, p, u * ) = G u * , (44)
with u * (t) a control that maximizes the pseudo-Hamiltonian for almost all time; it is defined by

u * = p,G p,G
if p, G(t/ε, x) = 0; the maximized Hamiltonian H ε (t, p, x) = max u h(t, x, p, u) is the one in (42), and ( 44) is then the differential equation ( 42), whose right-hand side is discontinuous at points where p, G(t/ε, x) vanishes.

Let us now turn to (P 0 ). Since the set of admissible velocities is not a priori smooth with respect to the state variable we use a non-smooth version of the Pontryagin maximum principle for differential inclusions, that we recall for self-containedness: Theorem 9. 3) is indeed locally Lipschitz: in local coordinates , for x 1 , x 2 in R n , denoting by δ the Hausdorff distance between two sets, one has:

δ (E(x 1 ), E(x 2 )) = max sup v1∈E(x1) inf v2∈E(x2) v 1 -v 2 , sup v2∈E(x2) inf v1∈E(x1) v 1 -v 2 = max U1 ∞ ≤1 min U2 ∞≤1 G(x 1 , U 1 ) -G(x 2 , U 2 ) ≤ Lip G x 1 -x 2 .
According to [START_REF] Bullo | Averaging and vibrational control of mechanical systems[END_REF], the Hamiltonian H defined in the above quoted theorem coincides with the map H defined in [START_REF] Bombrun | On the averaged control system[END_REF].

Remark 3.20. This result and Theorem 3.18 have two interpretations: 1. They prove that the operations of averaging and computing the Hamiltonian for the minimum time problem commute. Indeed, the Hamiltonian H was obtained by applying the maximum principle to problem (41), i.e. minimum time for the average system (4), but it also the average of the one in (42) with respect to the fast variable. 2. They prove indirectly an averaging result for the minimum time control problem (41); the averaging techniques in [START_REF] Chaplais | Averaging and deterministic optimal control[END_REF] do not apply to minimum time for they require smoothness of the Hamiltonian, while averaging is used in [START_REF] Geffroy | Optimal low-thrust transfers with constraints-generalization of averaging technics[END_REF][START_REF] Geffroy | Généralisation des techniques de moyennation en contrôle optimal -Application aux problèmes de rendez-vous orbitaux en poussée faible[END_REF] for minimum time with only partial theoretical justifications but numerical evidence of efficiency.

Let us now focus on the differential equations (43) that govern the extremals of (P 0 ). It is of great importance to know whether it defines a Hamiltonian flow on T * M , i.e. whether solutions trough all initial conditions are unique or not. Its right-hand side is continuous because, from Theorem 3.13, H is continuously differentiable; this ensures existence of solutions. We saw that H is smooth (C ∞ ) on T * M \ Z (see (33)), hence solutions through points outside Z are always unique. The following result gives uniqueness of solutions even on Z in the less degenerated case possible.

Theorem 3.21 (Hamiltonian flow for (P 0 )). Assume that the rank condition (30) holds, as well as conditions (i) and (ii) in Theorem 3.14. Then the differential equation (43) has a unique solution from any initial condition.

Proof. For an autonomous ODE ż = f (z) in a finite dimensional space, where c ln ln(1/u). Proving existence of a flow for (43) in more general situations (weaker sufficient condition) than this theorem is an interesting program to be pursued. However, it turns out to be applicable to the control of orbit transfer with low thrust, see §5.

f satisfies f (z 1 ) -f (z 2 ) ≤ ω( z 1 -z 2 ) with ω : [0, +∞) → [0, +∞) non-
Point (ii) is very mild and only states that the control vector fields are linearly independent. Point (i) is more artificial: the fact that p, G(θ, x) = 0 has at most one solution θ has to be checked by hand, while the fact that ∂G/∂θ does not vanish at the same time is equivalent to the rank condition

rank G(θ, x), ∂G ∂θ (θ, x) = dim Range G(θ, x) + Range ∂G ∂θ (θ, x) = n . (45)
It is true for the Kepler problem and used in [START_REF] Caillau | Coplanar control of a satellite around the Earth[END_REF] to show that the discontinuities in (42) are always "π-singularities", i.e. the control u * switches to its opposite.

4. Kepler control systems. We call Kepler control system with small control a family of control system on S 1 × M of the form

(K ε ) θ = ω(θ, x) + g(θ, x) v ẋ = G(θ, x) v , v ≤ ε , (46)
where G and g can be viewed, with the same convention is in (1), as n × m and 1 × m matrices smoothly depending on (θ, x) and ω is a smooth function S 1 × M → R that remains larger than a strictly positive constant:

ω(θ, x) ≥ k ω > 0 ∀(θ, x) ∈ S 1 × M . (47)
In fact, this is an affine control system on

S 1 × M ξ = f 0 (ξ) + m i=1 v i f i (ξ) (48)
with ξ = (θ, x), f 0 = ω ∂ ∂θ and, for 1 ≤ i ≤ m, the smooth vector field f i is represented by the i th column of the matrix notations G and g. If, in (48), one only assumes that, all solutions of ξ = f 0 (ξ) are periodic, additional conditions are needed for the orbits to induce a nice foliation that splits the state manifold into a product M × S 1 .

4.1. Relation with fast oscillating systems. For a solution t → (θ(t), x(t)) of (K ε ) in (46), let Θ(t) be the cumulated angle i.e. Θ(.) is continuous [0, T ] → R with Θ(t) ≡ θ(t) mod 2π for all t and Θ(0) ∈ [0, 2π), and define a new "time"

λ = R(t) ∆ = ε (Θ(t) -Θ(0)) . (49)
Proof. We assume that M is R n , d the Euclidean distance and all vector fields have a common compact support, hence all maps share a global Lipschitz constant and a global bound; by "a constant", we mean a number that depends only on these bounds and Lipschitz constants. It is left to the reader to check that, as for the proof of Theorem 3.7, the present proof extends to M with any distance d described in §2.3.

Let τ → x 0 (τ ) be a solution of (53) on [0, T ]. Define P(.) by

P(τ ) = τ 0 ω (x 0 (t)) dt (56)
and x 0 (.) by x 0 (λ) = x 0 (P -1 (λ)). The latter is a solution on [0, P(T )] of

d x 0 dλ ∈ 1 ω( x 0 ) E( x 0 ) (57)
with E defined by ( 3) and (54). This is the average system (in the sense of Definition 3.2) of the fast oscillating control system

( Σ θ0,ε ) d x dλ = G(θ 0 + λ ε , x) u ω(θ 0 + λ ε , x) , u ≤ 1 . (58)
Theorem 3.7 (Point 1) yields a family of controls u ε such that the solutions x ε (.) of ( Σ θ0,ε ) with initial condition x 0 (0) and control u ε converge to x 0 (.) uniformly:

d( x ε (λ), x 0 (λ)) ≤ c ′ ε for all λ ∈ [0, P(T )] (59)
for some constant c ′ . For each ε, let x ε (.) be the solution of ( Σ θ0,ε ) -see (51)with same initial condition and same control. Since (51) can be re-written as

d x dλ = 1 - ε g(θ 0 + λ ε , x) u ω(θ 0 + λ ε , x) + ε g(θ 0 + λ ε , x) u G(θ 0 + λ ε , x) u ω(θ 0 + λ ε , x) , (60)
the norm of the difference between the right-hand sides of ( Σ θ0,ε ) and ( Σ θ0,ε ) is bounded by k ε for some constant k > 0; classical theorems on smooth dependence of solutions on "parameters" yield some constant c ′′ such that

d( x ε (λ), x ε (λ)) ≤ c ′′ ε for all λ ∈ [0, P(T )] . (61)
Then define

t = T(λ) = 1 ε λ 0 dℓ ω(θ 0 + ℓ ε , x ε (ℓ)) + ε g(θ 0 + ℓ ε , x ε (ℓ)) u ε (ℓ) (62)
and the controls t → u ε (t) by u ε (λ) = u ε (T(λ)); the solutions x ε (.) of ( 46) with these controls are given by x ε (λ) = x ε (T(λ)), and one therefore has

d(x ε (T•P(τ )), x 0 (τ )) < (c ′ + c ′′ )ε, τ ∈ [0, T ] . (63)
Now, on the one hand, (62) yields

T(P(τ )) = 1 ε P(τ ) 0 dℓ ω(θ 0 + ℓ ε , x ε (ℓ)) + ρ (64) with ρ = 1 ε P(τ ) 0 ω(ℓ) -ω(ℓ) ω(ℓ) + ε g(ℓ) u ε (ℓ) ω(ℓ) + ε g(ℓ) u ε (ℓ) dℓ ω(ℓ)
where ω(ℓ) stands for ω(θ 0 + ℓ ε , x ε (ℓ)), ω(ℓ) stands for ω(θ 0 + ℓ ε , x ε (ℓ)), and g(ℓ) stands for g(θ 0 + ℓ ε , x ε (ℓ)); using (61) and Lipschitz continuity of ω to bound the first term in the integral, this implies that |ρ| is bounded by a constant. On the other hand, one has, according to (56), τ = P(τ ) 0 dλ/ω( x 0 (λ)). Developing ω according to its definition (52), in which we add θ 0 + λ ε to θ without changing the integral due to periodicity, τ is also equal to 1 2π θ∈S 1 , 0≤λ≤P(τ ) dλ dθ/ω(θ 0 + λ ε + θ, x 0 (λ)). Finally, performing the change of variable λ = ℓ -εµ(θ), with µ(θ) defined in §2.4, yields

τ = P(τ )+εµ(θ) εµ(θ) 1 2π 2π 0 dθ ω(θ 0 + ℓ ε , x 0 (ℓ -εµ(θ)))
dℓ Using (59), the fact that |µ(θ)| < 2π and Lipschitz continuity of both x 0 and ω, we deduce from this and equation (64) that T•P(τ ) -τ ε ≤ |ρ| + k ′ for some constant k ′ and finally, using the fact that x ε (.) is Lipschitz continuous with constant 2 ε sup G /k ω , one has d(x ε (T • P(τ )), x ε ( τ ε )) < c ′′′ ε for some constant c ′′′ . This and equation (63) imply implies point 1 of the theorem, with c = c ′ +c ′′ +c ′′′ in (55).

For point 2, consider θ n (.), x n (.) : [0, T /ε n ] → S 1 × K a solution of system (46) with ε = ε n and some control u = u n (t). Following (49)-( 51) and setting λ = R n (t) (we write R n because R in (49) is constructed for system ( Σ θ0,εn ) and thus depends on n), one associates to these x and u a control λ → u n (λ) and a solution λ → x n (λ) of ( Σ θ0,εn ). The solutions λ → x n (λ) of ( Σ θ0,εn ) with same control and initial condition satisfy, for the same reasons as (61), d( x n (λ), x n (λ)) < c ′′ ε n for some constant c ′′ . By Theorem 3.7 (Point 2), the sequence ( x n ) is compact and subsequences converge to solutions λ → x 0 (λ) of (57), hence the same subsequences of ( x n ) converge as well, and, with τ = Q(λ)

∆ = λ 0 dℓ ω( x(ℓ)) , the maps τ → x n (Q -1 (τ )) = x n ((Q • R n ) -1 (τ ))
converge to a solution τ → x 0 (τ ) = x 0 (Q -1 (τ )) of the average system (53), with distance less than c ′ ε n for some constant c ′ . Using the same argument as in Point 1 for T•P(τ ), one gets a bound for

|(Q•R n ) -1 (τ ) -τ εn | and, for some constant c ′′′ , d( x n ((Q • R n ) -1 (τ )) , x n ( τ ε )) ≤ c ′′′ ε n . Point 2 is proved. 4.4.
Dimension of E(x). In §3.3, and in particular in Proposition 3.10, G can simply be replaced with G. It is however interesting to give a more intrinsic characterization of r(θ, x) and thus of dim E(x).

Proposition 4.4. If (46) and (48) represent the same control system, then

dim j∈N Range ∂ j G ∂θ j (θ, x) = -1 + rank {f 0 (θ, x)} ∪ ad j f0 f k (θ, x) , j ∈ N, 1 ≤ k ≤ m . (65)
Proof. Straightforward computation using the fact that f 0 = ∂/∂θ. Note that the right-hand side is r(θ, x). Proposition 3.10 applies, with this definition of r. In particular, the "full rank case" becomes: Proposition 4.5. If the vector fields f 0 and ad j f0 f k , 1 ≤ k ≤ m, j ∈ N span the whole tangent space of S 1 × M , then E(x) has a nonempty interior for all x. 4.5. The function H(x, p). Instead of ( 5), H has to be taken as follows, with ω defined in (52):

H(x, p) = ω(x) 1 2π 2π 0 p, G(θ, x) ω(θ, x) dθ . ( 66
)
The characterization of E(x) in Proposition 3.4 is unchanged. In the "full rank case", the results from §3.4 on the degree of differentiability apply without a change.

4.6. Application to the minimum time problem. As in §3.5, but for the Kepler system (46), let x 0 , x 1 be fixed, call T ε the minimum time such that, from some θ 0 , θ 1 , (θ 1 , x 1 ) can be reached from (θ 0 , x 0 ) in system (K ε ) (obviously T ε → +∞ as ε → 0) and T 0 the minimum time such that x 1 can be reached from x 0 in the average system (53). The equivalent of Theorem 3.18, with a similar proof, using Theorem 4.3, is: Theorem 4.6. In the full rank case, one has lim sup ε→0 εT ε ≤ T 0 (hence ε T ε is bounded as ε → 0). If, for all ε > 0 small enough, there is a minimizing solution (θ ε , x ε ) : [0, T ε ] → S 1 × M and they all remain in a common compact subset of M , then all accumulation points (as ε → 0) of the compact family (τ → x ε ( τ ε )) ε>0 in C 0 ([0, T 0 ], M ) are minimizing for the average system and lim ε→0 ε T ε = T 0 .

The Hamiltonian for minimum time for the average system is given by (66); one has to perform the time scaling described in §4.1 to have a result like Theorem 3.21 and the simple "commutation between averaging and writing Hamiltonian" noted in Remark 3.20. Let us translate in terms of (46) the sufficient condition for existence of a Hamiltonian flow given by Theorem 3.21:

Theorem 4.7. In the full rank case, assume that p, G(θ, x) and p, ∂G/∂θ(θ, x) do not vanish simultaneously outside {p = 0}, that θ → p, G(θ, x) vanishes at most once for each (x, p) ∈ T * M , p = 0, and that rank G(θ, x) = m for each (θ, x) ∈ S 1 ×M . Then (43), with H given by (66), has a unique solution for any initial condition.

The discussion that follows Theorem 3.21 also applies to the above; let us mention that, once it has been checked that, for each (x, p), p, G(θ, x) vanishes for at most one θ, the other conditions are guaranteed if (45) holds with G replaced by G or, in terms of the vector fields in (48), if, for all ξ = (θ, x), rank{f 0 (ξ), f 1 (ξ), . . . , f m (ξ), ad f0 f 1 (ξ), . . . , ad f0 f m (ξ)} = n + 1.

(67)

We prove in the next section that the above conditions are true for the planar control 2-body problem.

5. Application to the controlled 2-body system. In this section we study some properties of the planar control system and demonstrate that it satisfies the condition of Theorem 3.21 on the domain of non-degenerated elliptic orbits.

5.1. Planar control 2-body system. It is classically described by some first integrals of the free movement -here the semi-major axis a and the eccentricity vector (e x , e y )-and one angle L following the dynamics; we restrict to the set of non-degenerated elliptic orbits rotating in the direct sense, i.e. the state space is S 1 × M with M = {(a, e x , e y ) ∈ R 3 , a > 0 and e x 2 + e y 2 < 1}. The control u = (u t , u n ) is expressed in the tangential-normal frame and the system reads: The eccentricity e is the norm of the eccentricity vector, e = e 2

d dt     a e x e y L     = 1 a 3/2     0 0 0 w(e x , e y , L)     + √ a     2
x + e 2 y . Low thrust translates into u ≤ ε for a small ε.

Remark 5.1. This is indeed a "Kepler control system" of the type (46) except that ω = w/a 3/2 is, although strictly positive, not bounded from below by a positive constant on S 1 ×M . There is such a lowerbound if one replaces M by M c = {(a, e x , e y ) ∈ R 3 , a > 0 and e x 2 + e y 2 < c} with c < 1. Strictly speaking, the results of the paper have to be applied in M c, c < 1. However, Theorems 4.3 or 4.7, for instance, may be applied in M because each statement may ultimately be restricted to a compact subset of M , itself included in some M c, c < 1.

The Hamiltonian that both defines the average system according to [START_REF] Bombrun | Asymptotic behavior of time optimal orbital transfer for low thrust 2-body control system[END_REF] and yields the Hamiltonian system governing extremals for minimum time is given by (66). Since 2 + e y 2 < 1 and a > 0, and each (A, X, Y ) = (0, 0, 0), the vector (A X Y )G(e x , e y , L) vanishes for at most one angle L.

Proof. Removing denominators, the equations Aa a + Xa x + Y a y = 0 and Xb x + Y b y = 0 can be written:

2 e x A + 2(1 -e 2 )X cos L + 2 e y A +2(1 -e 2 )Y sin L = -(1 + e 2 )A -2 e x (1 -e 2 )X -2 e y (1 -e 2 )Y
-2 e x e y X + (e 2

x -e 2 y + 1)Y cos L+ (e 2

x -e 2 y -1)X + 2 e x e y Y sin L = 2 e y X -2 e x Y.

If ∆ = 2 e x A + 2(1 -e 2 )X 2 e y A + 2(1 -e 2 )Y -2 e x e y X + (e 2
x -e 2 y + 1)Y (e 2

x -e 2 y -1)X + 2 e x e y Y is nonzero, there is clearly at most one solution L. If ∆ = 0, there exists λ = 0 such that 2 e x A + 2(1 -e 2 )X = λ -2 e x e y X + (e 2

x -e 2 y + 1)Y , 2 e y A + 2(1 -e 2 )Y = λ (e 2

x -e 2 y -1)X + 2 e x e y Y , and there may be a solution to the system above only if

(1 + e 2 )A + 2 e x (1 -e 2 )X + 2 e y (1 -e 2 )Y = -2 λ(e y X -e x Y )
These three equations forms a linear system in (A, X, Y ), M (A, X, Y ) T = 0 with A brief computation gives det M = (1 -e) 3 (1 + e) 3 (λ 2 + 4), strictly positive when 0 ≤ e < 1. Hence M (A, X, Y ) T = 0 implies (A, X, Y ) = 0.

M =   2e x 2 (
Since the rank of G is obviously equal to 2 and the rank of {G, ∂G/∂L} equal to 3 for any (e x , e y , L), the hypotheses of Theorem 4.7 are satisfied by the planar control 2-body system, and it guarantees existence of a flow for the Hamiltonian system governing the extremals of minimum time for its average system.

Conclusion.

Attempting to formulate a control theory equivalent to the averaging theorems for ODEs naturally leads to, and justifies, the notion of average control system introduced in this paper. It has a conceptual importance as well as, for instance, applications to approximation of minimum time control.

Besides its definition and description, we gave results on its regularity and on the dimension of its velocity set ("number of inputs"). These are however mostly a starting point. The regularity of H has to be further explored when the conditions of Theorem 3.21 do not hold, see the last paragraph of §3.

It has already allowed us to give (with restrictions on the eccentricities, see Remark 5.1) a proof [START_REF] Bombrun | Asymptotic behavior of time optimal orbital transfer for low thrust 2-body control system[END_REF] that the minimum time between 2 ellipses grows like 1/ε for the planar 2-body problem. Here also, progress must be made. Explicit computation of the average system and its extremals for the 2-body problem has to be conducted.

Acknowledgements. The authors are indebted to Jana Němcová from Institute of Chemical Technology, Prague, for a careful proof-reading of the draft manuscript, and to two anonymous referees from this journal for extremely constructive reviews that make this paper considerably easier to read than the original submission.

Appendix. Proof of Proposition 3.16. Proof of Point 1. The integral in (37) is well defined (its integrand is bounded) and, by (36) and Lebesgue convergence theorem, it is continuous with respect to X and h. Let us prove that this dH is the derivative of H. Since V is smooth, one has

V(θ, X + h) -V(θ, X) - ∂V ∂X (θ, X).h ≤ k h 2 , ( 69
)
where ∂V ∂X (θ, X) is smooth with respect to (θ, X) and k is some local constant. Now, assuming V(θ, X) = 0, one has V(θ, X + h) -V(θ, X) = V(θ, X + h) -V(θ, X)

V(θ, X) V(θ, X) + a(θ, X, h) V(θ, X + h) -V(θ, X) 2 V(θ, X) + V(θ, X + h) with |a(θ, X, h)| ≤ 2. Hence, from (69) and (37), one has, for some local constant k ′ , H(X + h) -H(X) -dH(X).h h ≤ k ′ 2π 2π 0 h + V(θ, X + h) -V(θ, X) V(θ, X) + V(θ, X + h) dθ for h small enough. For fixed X and h → 0, the integrand in the right-hand side is bounded by 1 + h and converges to zero for θ outside the set {θ ∈ S 1 , V(θ, X) = 0}: by (36) and Lebesgue convergence theorem, the right-hand side tends to zero.

Let us now state two lemmas that are needed in the proof of Point 2. Lemma A.1. Assume that X ∈ Z and (38) is satisfied. There is a neighborhood U of X in O d and a smooth map χ : U → S 1 such that, for (θ, X) ∈ U , one has V(θ, X) = 0 only if θ = χ(X), and ∂V ∂θ ( χ(X), X) V( χ(X), X) = 0 , X ∈ U . Consider the map T : S 1 × O d → R given by T (θ, X) = ∂V ∂θ (θ, X) V(θ, X) . Let X be in Z; since V(χ( X), X) = 0, one has T (χ( X), X) = 0 and ∂T /∂θ(χ( X), X) = ∂V ∂θ (χ( X), X) 2 , nonzero from assumption (38.b): the implicit function theorem yields a unique map χ from a neighborhood U of X in O d to a neighborhood of χ( X) in S 1 such that θ = χ(X) solves T (θ, X) = 0; it must therefore coincide with χ in U ∩ Z and satisfies the lemma.

Lemma A.2. Assume that X ∈ Z and (38) is satisfied. There exist a neighborhood U of X in O d , local coordinates x 1 , . . . , x d defined on U , and smooth maps P : U → SO(m), α : U → R, and W : S 1 × U → R m such that, with

X I =     x 1 . . . x m-1     ,
V(θ, X) = P (X) X I α(X) (θ -χ(X))

+ (θ -χ(X)) 2 W (θ, X)

= P (X)

  X I 0   + (θ -χ(X)) W 1 (θ, X) (72)
with W 1 (θ, X)

=   0 m-1 α(X)   + (θ -χ(X)) W (θ, X) , (73)
in S 1 × U , where α is bounded from below: 0 < α 0 < α(X), X ∈ U . Furthermore, for a constant K 3 > 0, one has, for all (θ, X) ∈ S 1 × U ,

V(θ, X) ≥ K 3 X I 2 + α(X) 2 (θ -χ(X)) 2 , ( 74
)
and

X I = 0 ⇒ W 1 (θ, X) ≥ K 3 . (75)
Proof. The map X → ∂V ∂θ ( χ(X), X) is nonzero for X = X, hence it does not vanish on a sufficiently small neighborhood U of X, and one may write ∂V ∂θ ( χ(X), X) = P (X) 0 m-1 α(X) , α(X) > α 0 > 0 .

For i between 1 and m -1, ∂vi ∂θ ( χ(X), X) = 0 from (76), and v i ( χ( X), X) = 0 from Lemma A.1 and, using(38.a), the rank of the map X → (v 1 ( χ(X), X), . . . , v m-1 ( χ(X), X)) is m -1 at X = X: on a possibly smaller neighborhood U , there are local coordinates x 1 , . . . , x d such that v i (θ, X) = x i + (θ -χ(X))

2 W i (θ, X) for i ≤ m -1 and for some smooth W i ; substituting (76) and ( 77) in (70) implies v m ( χ(X), X) = 0, hence v m (θ, X) = α(X) (θ -χ(X)) + W m (θ, X) (θ -χ(X))

2 for a smooth W m ; (71) is proved.

Possibly restricting U to a subset with compact closure, W (θ, X) is bounded on S 1 × U ; if |θ -χ(X)| ≤ 1 2 α 0 / max W , then (74) holds with K 3 = 1 2 according to (71); on the set where |θ -χ(X)| ≥ 1 2 α 0 / max W , V does not vanish and hence ( X I 2 + α(X) 2 (θ -χ(X)) 2 ) 1/2 / V(θ, X) is bounded from below; (74) is proved, with K 3 smaller than this bound and than 1 2 . From (73), W 1 ( χ( X), X) = 0 because α does not vanish; from assumption (38.b) and (72) (where X I = 0 if X = X), W 1 (θ, X) = 0 if θ = χ( X), hence W 1 does not vanish on S 1 × { X}; it is therefore bounded from below on S 1 × U with U a small enough neighborhood of X: (75) holds with K 3 smaller than this bound.

Proof of Proposition 3.16 (Point 2). We use [-π, π] instead of [0, 2π] as an interval of integration. Let h ∈ R d , with h = 1. From (37), one has, for some constant K using bounds on the derivatives of the smooth V, |dH(X).h -dH(Y ).h| ≤ 1 2π

π -π ∂V ∂X (θ, X).h - ∂V ∂X (θ, Y ).h V(θ, X) V(θ, X) dθ + 1 2π π -π ∂V ∂X (θ, Y ).h V(θ, X) V(θ, X) - V(θ, Y ) V(θ, Y ) dθ ≤ K X -Y + K 2π π -π V(θ, X) V(θ, X) dθ - π -π V(θ, Y ) V(θ, Y ) dθ .
Finally, defining

V(ϕ, X) = V( χ(X) + ϕ, X) , W 1 (ϕ, X) = W 1 ( χ(X) + ϕ, X) , (78)
and making a different change of variables in the last two integrals, one has

dH(X).h -dH(Y ).h ≤ K X -Y + K 2π π -π V(ϕ, X) V(ϕ, X) - V(ϕ, Y ) V(ϕ, Y ) dϕ ≤ K X -Y + K π π -π V(ϕ, X) -V(ϕ, Y ) V(ϕ, X) dϕ (79)

Proposition 3 . 16 .

 316 Let d be a positive integer, O d an open subset of R d , V : S 1 × O d → R m a smooth map (C ∞ ), Z the subset of O d where V vanishes for some θ and H: O d → [0, +∞) the average of the norm of V:

  decreasing, Kamke uniqueness Theorem [15, chap. III, Th. 6.1] states that uniqueness of solutions holds if α 0 du ω(u) = +∞ for arbitrarily small α > 0. From Theorem 3.14, we are in this case with ω(u) = c u ln(1/u), and du ω(u) = -1

2 ( 1 -e 2 ) 3 / 2 ,

 21232 a a a (e x , e y , L) 0 2 a x (e x , e y , L) b x (e x , e y , L) 2 a y (e x , e y , L) b y (e x , e y , L) e x , e y , L) = (1 + e x cos L + e y sin L) a a (e x , e y , L) = 1 + e 2 + 2 e x cos L + 2 e y sin L √ 1 -e 2 , a x (e x , e y , L) = √ 1 -e 2 1 + e 2 + 2 e x cos L + 2 e y sin L (e x + cos L) , a y (e x , e y , L) = √ 1 -e 2 1 + e 2 + 2 e x cos L + 2 e y sin L (e x + cos L) , b x (e x , e y , L) = √ 1 -e 2 1 + e 2 + 2 e x cos L + 2 e y sin L × -2 e y + (e 2 x -e 2 y -1) sin L -2 e x e y cos L 1 + 2 e x cos L + 2 e y sin L , b y (e x , e y , L) = √ 1 -e 2 1 + e 2 + 2 e x cos L + 2 e y sin L × 2 e x + (e 2x -e 2 y + 1) cos L + 2 e x e y sin L 1 + 2 e x cos L + 2 e y sin L .

2π0. 5 . 2 .

 52 dL/w(e x , e y , L) = 2π, it can be expressed as H(a, e x , e y , p a , p ex , p ey ) = √ aH(e x , e y , ap a , p ex , p ey ) with H(e x , e y , A, X, Y ) = 1 2π 2π 0 (A X Y )G(e x , e y , L) , G(e x , e y , L) =   2 a a /w 0 2 a x /w b x /w 2 a y /w b y /w   Hamiltonian flow. Theorem 4.7 applies to this system. Indeed: Proposition 5.2. Fore each (e x , e y , a) with e x

( 70 )

 70 Proof. From (38.a), Z = {(θ, X)∈ S 1 × O d , V(θ, X) = 0} is a smooth submanifold of S 1 × O d and from (38.c), Z given by (35) a smooth submanifold of O d , both of dimension d + 1 -m, and the projection π : S 1 × O d → O d induces a diffeomorphism Z → Zwhose inverse is of the form x → (χ(x), x) with χ a smooth map Z → S 1 that satisfies, for all X ∈ Z: V(θ, X) = 0 if an only if θ = χ(x).

Define v 1

 1 , . . . , v m , smooth maps S 1 × U → R by    v 1 (θ, X) . . . v m (θ, X)    = P -1 (X) V(θ, X) .

  1 -e 2 + λe x e y ) -e 2 -λe x e y ) (1 + e 2 ) 2 e x (1 -e 2 ) + λe y 2 e y (1 -e 2 ) -λe x

	2e y	-λ(e 2 x -e 2 y -1)	-λ(e 2 x -e 2 y + 1) 2(1   .

Its relation to the limit case of (1) as ε → 0 is discussed in the next section.

Taking ε 0 small enough so that |ω(θ, x) + ε g(θ, x) u| > k ω /2 for x ∈ K, u ≤ 1, ε < ε 0 , one has dR/dt > ε k ω /2 hence R is strictly increasing and one-to-one, and

Then λ → x(λ) = x(R -1 (λ)) is a solution of the system

associated with the control λ → u(λ) = v(R -1 (λ))/ε. Except for the term εg u in the denominator, this is a fast oscillating system [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF] with G = G/ω. We now apply §3.

4.2. Average control system. The definition uses ω defined by

Definition 4.1 (Average control system of Kepler control systems). The average control system of the Kepler control system (46) is the differential inclusion

with E defined by (3) using

instead of (2). Solutions are defined as in Definition 3.2. Remark 4.2. This is almost Definition 3.2 applied to (51), which is equivalent to (46) via time changes, except :

(i) the term εg u in the denominator of (51) has been discarded, (ii) the right-hand side has been multiplied by ω(x). 1. Let x 0 (.) : [0, T ] → M be an arbitrary solution of (53) and θ 0 ∈ S 1 . There exist a family of measurable functions u ε (.) : [0, T ] → B m , indexed by ε > 0, and positive constants c, ε 0 , such that, if t → (θ ε (t), x ε (t)) is the solution of (46) with control u = u ε (t) and initial condition (θ ε (0), x ε (0)) = (θ 0 , x 0 (0)), it is defined on [0, T /ε] for ε smaller than ε 0 and

Let K be a compact subset of M , (ε n ) n∈N a decreasing sequence of positive real numbers converging to zero, and θ n (.), x n (.) : [0, T /ε n ] → S 1 × K a solution of system (46) for each n, with ε = ε n and some control u = u n (t),

Then the sequence τ → (x n (τ /ε n )) n∈N is compact for the topology of uniform convergence on [0, T ] and the limit of any converging sub-sequence is a solution x * (.) of the average differential inclusion (53).

where the last inequality uses the fact u

}, and also holds with V(ϕ, Y ) instead of V(ϕ, X) in the denominator. Now let us use Lemma A.2, let X = (x 1 , . . . , x d ) and Y = (y 1 , . . . , y d ) in these coordinates; from (72), one has, with W 1 defined by (78),

(81) Two cases are to be distinguished:

(i) If X I = Y I = 0, then ϕ factors out of V(ϕ, X) and V(ϕ, Y ) in (80) and the last term in (81) is zero: according to (75), the integrand in (79) is bounded by

and finally |dH(X).h -dH(Y ).h| ≤ K X -Y with a constant K that depends only on V, the open set U and the coordinates.

(ii) If X I = 0 (or Y I = 0, interchanging X and Y ), then (81), using (74), implies that the integrand in (79) is bounded by

but the same is also true replacing α(X) with α(Y ) and X I 2 with Y I 2 ; hence, since a -b 2 ≤ 4 max{ a 2 , b 2 }, the last term may be replaced by 2

K3

XI-YI 2 XI-YI 2 +4α0ϕ 2 , whose integral between -π and π is equal to

which is less than X I -Y I (k 1 +k 2 ln 1 XI-YI ) for some k 1 , k 2 when, say, XI-YI 2 √ α0 < 1. Finally, since X I -Y I is less than X -Y and u → u ln(1/u) is nondecreasing, less than X -Y (k 1 + k 2 ln 1 X-Y ). Cases (i) and (ii) do imply (39), possibly restricting U so that ln 1 X-Y ≥ 1.


Céline Blanc

Dominique Darmendrail

Benchmarking of Excavated Soil Reuse Management in Europe -Development of a French management framework

come

The consideration of the receptor site and soil qualities the basic principles of the national/regional approach (e.g. typology of acceptable / sustainable reuses, site specific or generic, cost-benefits analysis). The Netherlands with a generic approach (Values, Values and Values)

The United Kingdom with a specific approach (as the French legal framework for contaminated land & sites management)

vendredi 7 décembre 2012

D3E / DG > 6
The Netherlands : a generic approach Ensure materials traceability from the excavation site to the receiver site, Bring information on the status, location and quality of stocks, Connect soil producers / site owners and users of excavated soil, Enable the competent authorities to make their regulatory controls Produce useful indicators for monitoring the implementation.

> Further developments

Methodology for constrution materials reuse Regulatory obligations

>

  Linkage can be reused on a receiver sites:For road construction, As part of a development project (for which a building permit or an EIS is issued and under conditions)

Challenge with Excavated Soil Reuse in France > Demand of French Ministry of Environment to develop a national management framework with benefit from some experienced countries.

	Benchmarking of Excavated Soil Reuse Benchmarking of Excavated Soil Reuse
	Management in Europe Management in Europe
	> The Netherlands > EU Waste Framework Directive Application of soil
	Defines the legal Framework at EU level Benchmarking of Excavated Soil Implementation at a National Level: different options, different tools
	D3E / DG > vendredi 7 décembre 2012	Reuse Management in Europe & Development of a French management framework Blanc Céline, Darmendrail Dominique, BRGM, France	> 2

60 million tons/year primary (clean) sand 20 million tons/year reuse of lightly contaminated soil

	Benchmarking of Excavated Soil Reuse
	Management in Europe
	>
	soil
	5% treated before re-use
	20% used in construction material
	> France???
	vendredi 7 décembre 2012	> 3

Turnover 200 million euro/year > Flanders Belgium 20 million tons excavated soil per year from the excavated re-used soil: 75% re-used as

Work based on the ICCL questionnaire (including USA/Canada) > Assessment of the legal framework of European countries / regions Switzerland, Belgium/Flanders &Wallonia, the Netherlands, Finland, Spain /Basque Country, Germany, Italy and the United Kingdom Status of Excavated Soils

: waste, not waste, polluted, not polluted Specific legislation or not for excavated soils : soil legislation, waste legislation

some additional constraints: the French management policy on excavated soils needs to be in agreement :

	Benchmarking of Excavated Soil Reuse Benchmarking of Excavated Soil
	Management in Europe Reuse Management in Europe
	> With Local implementation
	Liability issues
	D3E / DG
	vendredi 7 décembre 2012 vendredi 7 décembre 2012	> 4 > 5

with the national policy on contaminated land management (site specific approach, risk based) with the Waste management policy (generic approach with guidelines values, hazardous substances contents source approach) > taking into consideration:

> Finally Focus on 2 countries:

> 7 Standstill principle: Applicable soil should be of better or the same quality as the soil it is applied upon Fit for use principle: Applicable soil should be of a quality competent for the function of the site it is applied upon Remediation criteria REMEDIATION Free use AGRICULTURE NATURE RESIDENCE INDUSTRY

	The Netherlands : traceability
	D3E / DG
	vendredi 7 décembre 2012 8	> 8

Soil application register Soil applier Competent authority Soil+ National authority Information on soil application Direct communication Access to the register

	The United Kingdom : a specific approach
	> Different methodologies (cf CoP v2):
	Use on the site of Origin
	Direct transfer
	Cluster/Hub site (if any pre-treatment/ treatment needed)

> Suitability for use Approach / reused soil quality Prescreening values (CLEA) Specific risk assesment if no guideline values or overtaking CLEA values Remedial Target Methodology :hydrogeological risk assessment for land contamination > Area approach / materials management plan (report)


Lei Yu

Pierre Maurel

Christian Barillot

Rémi Gribonval

Compressive Matched Filter for Cerebral Blood Flow Quantification with ASL: sampling diversity or repetition?

Keywords: ASL, CBF, Multi-TI, Diverse-TI, Compressive Matched Filter, Sparsity

Arterial Spin Labeling (ASL) is a noninvasive perfusion technique which allows the absolute quantification for Cerebral Blood Flow (CBF). The perfusion is obtained from the difference between images with and without magnetic spin labeling of the arterial blood and the captured signal is around 0.5-2% of the magnitude of the labeling images, so the noise is one of the main problems for further data analysis. Classical method, Mono-TI, for CBF quantification is averaging repetitions with only one Inversion Time (TI) -the time delay between labeling and acquisition to allow the labeled blood to arrive the imaging slice. It improves the robustness to noise, however, cannot compensate the variety of Arterial Arrival Time (AAT). In this paper, Diverse-TI is proposed to exploit different TI sampling instants (sampling diversity) to improve the robustness to variety of AAT and simultaneously average repetitions with each TI (sampling repetitions) to improve the robustness to noise. Generally, the sampling diversity is relatively small and can be considered as compressed measurements, thus the Compressive Matched Filter (CMF) enlightened from sparsity is exploited to directly reconstruct CBF and AAT directly from compressed measurements. Meanwhile, regarding the CBF quantification performance, the compromise between the sampling repetition and sampling diversity is discussed and the empirical protocol to determine the sampling diversity is proposed. Simulations are carried out to highlight our discussions.

Introduction

The Arterial Spin Labeling (ASL) is an MRI (Magnetic Resonance Imaging)based perfusion technique which uses the magnetically tagged water as a freely diffusible tracer to measure perfusion non-invasively. This blood water is first labeled with a radio-frequency pulse in the neck of the patient. After a delay, called Inversion Time (TI), which allows the labeled blood to arrive in the brain, a labeled image of the brain is acquired. A control image is also acquired without labeling and the CBF estimation is done on the difference between the control and labeled image. For a fixed TI t, the perfusion signal is generally described by a kinetic model introduced by Buxton et. al. in [START_REF] Buxton | A general kinetic model for quantitative perfusion imaging with arterial spin labeling. Magnetic resonance in medicine[END_REF]:

M ∆ (t) =      0, t ∈ [0, ∆) 2αf M 0b (t -∆)q p (t) exp(-t/T 1b ), t ∈ [∆, ∆ + τ ) 2αf M 0b τ q p (t) exp(-t/T 1b ), t ∈ [∆ + τ, ∞) (1)
where f is the CBF, α measures the labeling efficiency, M 0b is the equilibrium magnetization of arterial blood, ∆ is the Arterial Arrival Time (AAT) to the interesting slice, τ is the temporal width of the bolus, T 1b is the relaxation time in blood and q p (t) is considered to be approximately equal to 1. Specifically, in QUIPSS II [START_REF] W M Luh | QUIPSS II with thin-slice TI1 periodic saturation: a method for improving accuracy of quantitative perfusion imaging using pulsed arterial spin labeling[END_REF], a saturation pulse is given to the tagged bolus after a fixed period τ of the tagging time, thus the temporal width of the tagged bolus is known a priori.

The principle task of ASL analysis is quantifying the absolute value of CBF f . The traditional technique, which is called Mono-TI, only uses a single TI t after which the acquisition is gathered, where t is assumed to be bigger than ∆ + τ to guarantee a delay long enough to let the magnetic tagged blood arrive the interesting imaging slice. Then the quantification of CBF is a direct ratio between the magnetization difference M ∆ and the known terms, when t > ∆+τ .

However, two major problems prevent the efficiency of Mono-TI. At first, as the amplitude of the difference ASL signal is usually around 0.5-2% of the control image magnitude, its SNR is thus not sufficient for further analysis. Commonly, numbers of sampling repetitions are acquired for single TI (typically more than 30) and then averaged to improve the SNR. Secondly, the assumption t > ∆ + τ is not easy to be guaranteed in Mono-TI technique. Since ∆ is varying between different patients, ages and physical situations. Even excluding these factors, the AAT is still an uncertain value due to the presence of laminar and turbulence flow, complicated vessel networks and cardiac pulsations [START_REF] Wu | The effects of flow dispersion and cardiac pulsation in arterial spin labeling[END_REF]. This fact leads that the magnetic labeled blood may not have reached the imaging slice if the sampling time is too small and thus leads to underestimation of CBF. One possible way is giving the TI large enough, however, the magnetization difference M ∆ might be too small if TI is too big and thus leads to very low SNR.

In this paper, another measuring procedure with different sampling times is investigated for the CBF quantification, where the collected ASL data at different sampling times are captured during separate ASL RF pulse periods. It is different from the classical Multi-TI [START_REF] Reinoud | Arterial spin labeling perfusion MRI at multiple delay times: a correlative study with H(2)(15)O positron emission tomography in patients with symptomatic carotid artery occlusion[END_REF][START_REF] Macintosh | Multiple inflow pulsed arterial spinlabeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF][START_REF] Wang | Pediatric perfusion imaging using pulsed arterial spin labeling[END_REF] and thus called Diverse-TI to avoid ambiguities. Multi-TI is a more recent and confidential sequence than Mono-TI and is not currently available on most MR-Scanners, while Diverse-TI proposed in this paper can easily be produced since it only uses regular ASL sequences. Before further analysis on real ASL data, the first and also the essential step is giving a protocol to concretely design the measuring procedures. It is clear that some issues should be preferentially considered: (1) sampling repetition (for the same sampling time), which is intended to improve SNR; (2) sampling diversity (sampling at different times) which is intended to compensate for the inexact knowledge of parameters such as ∆, etc. For practical considerations, the most crucial criteria in real clinical studies is the total measuring time which is generally limited in a reasonable period. Thus the values of sampling repetition and sampling diversity can not be designed as large as possible to improve CBF quantification performance. A method to explicitly guide the design of these parameters is applaudable. However, to the best knowledge of the authors, the existing papers rarely focus on this point.

Consequently, as a preliminary study, the contribution of this paper is twofold. First, it formulates the considered problem as an instance of Compressive Matched Filter (CMF) [START_REF] Mark A Davenport | Signal Processing With Compressive Measurements[END_REF]. Then, it numerically investigates the tradeoffs between repetition and sampling diversity (for the same total measurement time) on CBF estimation with CMF. In a word, the main question investigated in this paper is: given a total measurement time, shall we favor repetitions or sampling diversity?

Diverse-TI Technique

Since parameters including τ, T 1b , α, M 0b and q p (t) are (or assumed to be) known a priori, the signal model ( 1) can be simplified as M ∆ (t) = f × g(∆; t) where the time related term g(∆; t) will be called "wave form" and can be defined by ( 1) with f = 1.

-Sampling Repetition

In practice, an important level of noise is affecting ASL measurements, hence the signal captured at time t can be expressed as follows:

y(t) = M ∆ (t) + (t) = f × g(∆; t) + (t)
where (t) is the noise term and is assumed to be Gaussian in this paper. In our setting, the measure at time t is repeated R times (sampling repetition), and averaging all measures divides the noise variance by R and thus it will improve the SNR according to the following equation:

Q R = 10 log 10 R + Q 1 (2)
where Q 1 is the SNR of one repetition and Q R is the SNR after R repetitions.

-Sampling Diversity Define g(∆) = {g(∆, iT s )} N i=1 where T s is a reference regular sampling interval and t 1 ...t M (possibly irregular) time instants located on this regular grid, where M N . With different TI (sampling diversity) t 1 ...t M , the ASL data are respectively collected y {y(t 1 ), ..., y(t M )} T ∈ R M . Then we can write the sampling model in the form of linear operation:

y = Φ(f × g(∆) + ) where Φ ∈ R M ×N is a sensing matrix, which verifies ΦΦ T = I M .
At last, the task is to find the CBF parameter f and the optimal ∆ that best match the observations y. To fulfill it, Compressive Matched Filter (CMF) [START_REF] Mark A Davenport | Signal Processing With Compressive Measurements[END_REF] is exploited in the following sections.

In the case considered in this paper, the CBF quantification problem can be formed as a multiple detection problem to distinguish the following hypothesis:

H i : y = Φ(f g(∆ i ) + ), with i ∈ {1, ..., d}
where f is the CBF value, ∆ 1 , ..., ∆ d are d possible AAT values that should be assumed before and ∼ N (0, σ 2 I N ) is the white noise term. In this case, the evidence under different hypothesis, for i ∈ {1, ...d}, can be written

p(y|H i ) = exp -1 2 (y -Φf g(∆ i )) T (σ 2 ΦΦ T ) -1 (y -Φf g(∆ i )) (2π) N/2 |σ 2 ΦΦ T | 1/2
Then the final detection is carried out by finding the hypothesis with the biggest conditional probability:

(f s , H s ) = arg max f,Hi,i∈{1,...,d} p(y|H i )
Thus by differentiating w.r.t. f , we can obtain its estimation

fi = y T (ΦΦ T ) -1 Φg(∆ i ) g(∆ i ) T Φ T (ΦΦ T ) -1 Φg(∆ i ) (3)
and by taking logarithm we obtain an equivalent test that simplifies to

H s = arg max Hi,i∈{1,...,d} y - 1 2 Φ fi g(∆ i ) T (ΦΦ T ) -1 Φ fi g(∆ i ) = arg max
Hi,i∈{1,...,d} 1 2

y T (ΦΦ T ) -1 Φg(∆ i ) 2 g(∆ i ) T Φ T (ΦΦ T ) -1 Φg(∆ i ) ΦΦ T =I

====== arg max

Hi,i∈{1,...,d} y,

Φg(∆ i ) Φg(∆ i ) 2 (4)
Above all, using (4) and (3), we can respectively obtain the estimation of arrival time ∆ = ∆ s and the estimation of CBF f = f s .

Analysis and Simulations

Without loss of generality, the time used for one single sampling repetition can be assumed to be T S , including the process of one time control imaging and tagging imaging. Then the total measuring time can be expressed as:

T = M • R • T S ( 5
)
where M indicates the level of sampling diversity i.e. the number of TI and R the number of sampling repetitions. For practical considerations, the most crucial criteria in real clinical studies is the total measuring time which is generally limited in a reasonable period. In one aspect, the sampling repetitions improve the signal SNR; while in the other aspect, the sampling diversity promote the robustness to the variety of AAT. Consequently, the compromise between the sampling diversity and sampling repetitions needs to be determined. In other words, given the maximum total time T , how to optimally choose R and M to reach better CBF quantification? To answer it, we first compare the extreme case when M = 1 to the other case when M > 1, then separately consider the case when M > 1.

Mono-TI technique v.s. Diverse-TI technique

In fact, M = 1 is exactly the Mono-TI technique. In order to well illustrate the performance of Mono-TI technique, a concrete simulation is carried on. In the following, the parameters of ASL model ( 1) are arbitrarily set as: M 0b = 100, α = 0.9, T 1b = 1.2s and τ = 1s, and Fig. 1(left) shows one example when AAT ∆ = 1 and CBF f = 10 mL/(100g)/min. Obviously, given constant noise variance σ, the highest SNR of single repetition Q 1 is reached when t i = ∆ + τ while it will decrease as t i increasing.

The assumption on TI t i > ∆ + τ is supposed to be verified. At first, t i can not be much smaller than ∆ + τ otherwise the tagged blood can not reached the interest slice on time, as shown by the red diamond point in Fig. 1(left) when t i = 1s; meanwhile, it can not be too big otherwise the captured signal will be too small which makes the SNR Q 1 very low, as shown by the green circle point in Fig. 1(left) when t i = 4.5s. Then the CBF f can be calculated directly by f = y(t i )/g(∆; t i ), however, the AAT ∆ can not be estimated in Mono-TI technique.

Then Mono-TI technique is carried out in the following procedure. Fix the CBF f = 10 mL/(100g)/min, give one AAT value ∆, sample at the time t i with noise corruption, and then one single repetition of the ASL signal is simulated. Without loss of generality, the time used for one single repetition is assumed T S = 5s (this value is set identically for all simulations). Then as an example, we can set the total time T = 30min, so according to [START_REF] Macintosh | Multiple inflow pulsed arterial spinlabeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF] the Mono-TI technique allows R = 360 repetitions. According to (2), the SNR can reach to Q R = 30.56dB with Q 1 = 5dB. However, this value is very sensible to the sampling location, or equivalently, the AAT value. This fact can be well illustrated by Fig. 1, where the result of Mono-TI is obtained by fixing t i = 1s, 3s or 4.5s, then ranging AAT ∆ from 0.1s to 3s (the performance shown in Fig. 1(right) is measured through the quantification SNR averaged over 1000 restarts). In Fig. 1, the result when t i = 1s shows the situation when TI is too short, t i = 4.5s shows the situation when TI is too long, even with proper TI t i = 3s, the performance of quantification might largely decrease when AAT and the TI are mismatch.

In the framework of Diverse-TI , following the CMF algorithm, we first assess a large interval where the value of AAT should locate. In this simulation, this interval is set to ∆ ∈ [0.1, 3]. We choose the number1 of sampling diversity to M = 10, and the sampling repetitions are set identically according to [START_REF] Macintosh | Multiple inflow pulsed arterial spinlabeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF] with the same total time as Mono-TI, then the CBF quantification performances through CMF algorithm with different AAT value are obtained, as shown in Fig. 1. Diverse-TI can give constant performance which is also better than Mono-TI with a majority of AAT. As a conclusion, in the case of M = 1, the CBF quantification quality might be better than M > 1, however, its performance is much sensitive to the sampling time, which is not easy to be designed practically. Moreover, as a supplement, Diverse-TI can also estimate the AAT.

Protocol to Sampling Diversity

In the last subsection, it has been shown that M = 1 is worse than M > 1, while in this subsection, we will propose a method to answer the question that how many sampling time locations are enough? Similarly, the total time T is the crucial parameter in clinical studies. Consequently, given a fixed total time T , the relationship between the CBF quantification performance and the number of sampling locations is essential to determine the required M .

In this simulation, the CBF f is fixed equal to 10mL/(100g)/min and the AAT ∆ is randomly choosing from the interval [0.1, 2], then the ASL signal are generated and corrupted by noises. Then to simulate the sampling procedure, the total time T is fixed respectively to 5, 10, 30, 60min and the number of sampling time locations M is varying from 2 to 30, thus the number of repetitions R is determined via [START_REF] Macintosh | Multiple inflow pulsed arterial spinlabeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF], then samples are collected. Since ∆ ∈ [0.1, 2] is assumed, the CMF can be used to carry out the CBF quantification and AAT estimation, and the SNR of CBF quantification and AAT estimation are respectively shown in Fig. 2 (bottom two figures). Meanwhile, the relationship between sampling repetition and sampling diversity can also be drawn according to [START_REF] Macintosh | Multiple inflow pulsed arterial spinlabeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF], as shown in Fig. 2 (upper figure). From the results, as expected, we can find that as M increases, the performance of Diverse-TI technique is improving before reaching an asymptotic value. From this simulation, we can find that the convergent point is around M = 8, as shown by the performance of CBF (or AAT) quantification in bottom figures of Fig. 2. Then according to [START_REF] Macintosh | Multiple inflow pulsed arterial spinlabeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF], the repetitions at one sample time R can be calculated. It is worth mentioning that this simulation is only an example of choosing the number of multiple sampling time locations, for concrete applications, the curve in Fig. 2 should be re-drawn with appropriate parameters, and then the best value should be re-selected.

Conclusion and Future Works

Through this preliminary study, it is shown that the new Diverse-TI technique is superior to the classical Mono-TI technique, since the first can give steady CBF performance and is capable to estimate the AAT. Furthermore, the sampling diversity of Diverse-TI technique has also been investigated where an experimental method is exploited and the protocol to design the sampling diversity, i.e. number of multiple sampling times, is proposed. Without loss of generality, from one exemplar simulation, we conclude that giving a total measuring time, the CBF quantification performance can be improved with the increase of the sampling diversity. Consequently, the sampling diversity should be designed large enough. On the other hand, when the sampling diversity is larger than a fixed value, the performance becomes constant. Thus, together with an upper bound resulted from the physical implementation constraints, the optimal region for sampling diversity can be well located through the proposed protocol.

The future works will consist in applying the parameter design protocol to guide the Diverse-TI technique in real ASL data acquisitions. Meanwhile, it is possible to extend CMF algorithm by considering additional priors to regularize the CBF estimation problem which might also improve the performance.

Fig. 1 .

 1 Fig. 1. (Left) Kinetic Model of ASL with ∆ = 1 (τ = 1s) and f = 10 mL/(100g)/min (dark solid curve) and its noisy version with average SNR 5dB (grey dashed curve). (Right) Comparison between Mono-TI and Diverse-TI (M = 10) methods for different AAT with total time T = 30.

Fig. 2 .

 2 Fig. 2. Performance Evolution with respect to different sampling diversities M via CMF using different total time: CBF quantification (bottom left), AAT (bottom right) and sampling repetitions (up).

This value is chosen according to empirical result, referring Fig.