To Thierry Dargent, for assistance with the many numerical emergencies and careful explainations -To Jean-Baptiste Caillau, for his suggestion which lead to including the bicircular model in Chapter 5 -To Jeremy Rouot, for all the help in Chapter 3 -To my parents, for

Introduction

In celestial mechanics, the 3-body problem is an ancient, classical problem of predicting the individual motions of celestial objects interacting with each other gravitationally. In space mechanics, where one of the bodies has a negligible mass and is controlled by a thrust force, the 3-body problem describes the motion of an artificial satellite in the Earth-Moon system. We will be interested in satellite transfers in this system, in the case of low thrust propulsion, which minimize transfer time. The movement of a satellite under the influence of two massive bodies can be divided, depending on satellite position relative to the other bodies, into two-body, perturbed two-body or threebody motion. There is also some interest in considering bicircular four-body motion, since the Sun's mass may play quite a significant role in satellite motion, to the same extent as the Moon's mass [START_REF] Roy | Luni-solar perturbations of an Earth satellite[END_REF]. This separation arises naturally for example from the separate planets' spheres of influence (SOI), or through the Hill's regions.

In fact, in the low thrust case, these different regions along a satellite trajectory permit very different techniques in both an analytic and a numerical sense. This is because when the satellite maximal thrust is a small parameter ǫ, the ratio between engine thrust and gravitational attraction of the three planets in the 3-body system may differ from one region to another. When this ratio is small, two time-scales are introduced into the dynamics which motivates the use of the technique of averaging in these regions. Also, each region has a different coordinate system, frame and model which best describe satellite motion within it. The existence of these separate regions lends a richness and complexity to the problem of low-thrust 3-and 4-body motion.

When Geffroy [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] framed the technique of averaging for low-thrust satellite transfers, the satellite motion was expressed in terms of two-body Keplerian motion. However, the veracity of the low-thrust assumption can be affected by the introduction of other primary bodies (the thrust becomes proportional to attraction to a primary, and so can no longer be considered as a 'low thrust') and even by other spatial effects such as the Earth's oblateness. A good question is, what effect does the wider space context, e.g. third-body gravity, have on the low-thrust satellite motion? The thesis directly addresses this aspect. We frame the low-thrust satellite motion within the wider space context, by studying specific aspects of not only satellite transfers within two-body motion but also of transfers where the motion is affected by naturally-occuring forces in the space environment -both in circumstances under which they can be considered perturbing forces and when they have a significant effect. Thus we extend the two-body problem naturally to the problem of three (and four) bodies by studying satellite transfers which initially require space perturbations to be considered, and finally which require us to consider a three and four-body model. We do this both by analytical studies of the perturbations (chapter 4) and by a numerical simulation (chapter 5) of an optimal low-thrust transfer from an Earth orbit to the L 1 Lagrange point, which takes place in a four-body context and under the influence of the low thrust.

Geffroy first introduced the concept of averaging within the Kepler problem 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]. This application to low-thrust satellite motion was exploited e.g. by Bonnard and Caillau and their collaborators ( [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF], [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF], [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF], [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF], [START_REF] Bonnard | Optimal Control with Applications in space and quantum dynamics[END_REF]) in a series of articles which discuss in detail aspects of the low-thrust energy-minimization problem of Kepler satellite transfers ('energy' here is a quadratic cost, the norm of the satellite thrust vector). The study of the energy-minimization optimal transfers is extensive, but the energy property is not a particularly physical property (although it has the advantage of being related to a Riemannian metric) and it is of interest to extend these techniques and directions of research considered here to other more physical optimization problems. Particularly, we choose to focus on time-optimal low-thrust transfers in the Kepler problem in our initial investigation on the effects of averaging in low-thrust two-body transfers.

When addressing non-Keplerian satellite motion, particularly perturbed two-body motion (where the perturbation is some other external force in the space environment besides the engine thrust) a difficulty encountered (e.g. in the description of lunar-perturbed satellite motion or in the ∆V studies carried out by Edelbaum in [START_REF] Edelbaum | Propulsion requirements for controllable satellites[END_REF]), is that simplifying assumptions are often made to obtain the cleanest expression of trajectories and perterbations which often have the additional affect of disallowing the inclusion of any further perturbations: for example, a quasi-circular assumption on satellite motion disincludes all perturbations on the satellite that are first-order in the eccentricity e. In this study we address this problem through the construction of low-thrust perturbed satellite models in which simplifications are done in such a way as to include the appropriate perturbations which makes the satellite motion increasingly accurate within this spatial region.

Previous studies have considered separation of the trajectories into various regions. However, not in any of the approaches we are aware of has the dominance of gravity over a low thrust along some arcs but not others been considered. For example, the approach given in [START_REF] Mingotti | Low-energy, low-thrust transfers to the moon[END_REF] consists of patching together two circular restricted 3-body problems, the 3D-Sun-Earth system and the Earth-Moon system, using invariant manifolds; firstly, of the Sun-Earth system, and then transferring to an invariant manifold of the Earth-Moon system until Moon capture. In [START_REF] Betts | Optimal low thrust trajectories to the moon[END_REF], the mission design for the Earth-Moon transfer includes decomposition of the trajectory into sub-sequences involving two-body transfer arcs. The transfer is divided into three phases, where the first phase consists of a two-body transfer around the Earth perturbed by the Moon, Sun and Earth's oblateness, the second phase is a free coast, and the third phase consists of two body transfer around the Moon, perturbed by the Earth and the Sun. A third thesis [START_REF] Daoud | Contribution au contrôle optimal du problème circulaire restreint des trois corps[END_REF] adresses optimal control of the circular restricted three-body problem (using no a-priori decomposition of the trajectory into distinct arcs) to compute optimal steering of the spacecraft on an Earth-Moon transfer. The minimization of time, energy, and fuel consumption is considered and the single shooting method is used. It is addressed via homotopy methods, which are used to continue solutions from 'easier' versions of the problem towards those for which the parameter values for the system are more in line with real values, but which are more difficult to solve.

Our study, in contrast with the previous given, involves both analytical and numerical analysis of the transition from the two-body to perturbed and multi-body models. Additionally, the use of averaging to simplify low-thrust transfers is emphasized in each kind of transfer, including the numerical study which takes place in a four-body context.

The organization of this thesis is as follows. In the first chapter, we provide the background to the two-body, three-body and bicircular four-body models, the perturbations considered, the averaging process and optimal control used in the thesis. The second chapter gives an overview of our contributions within the context of the background provided. In chapter 3, we consider the controlled Kepler two-body problem and its time-optimal solutions. Chapter 4 introduces the third-body and J 2 perturbations to the controlled Kepler problem and considers an energy-minimal problem (these two chapters 3 and 4 are reproductions of the publications [START_REF] Bonnard | Lunar and J 2 perturbations of the metric associated to the averaged orbital transfer[END_REF], [START_REF] Bonnard | Time versus energy in the averaged optimal coplanar Kepler transfer towards circular orbits[END_REF] which formed a part of this study). In the fifth chapter we discuss construction and numerical implementation of a transfer from low-Earth orbit to the L 1 Lagrange point.

Résumé en Français

En mécanique céleste, le problème ancien et classique des 3 corps est de prédire l'évolution des planètes sous l'action de la gravitation. En mécanique spatiale, le problème est de calculer la trajectoire d'un satellite artificiel, dont l'orbite osculatrice est caractérisée par cinq éléments orbitaux, sous l'attraction de deux planètes, la Terre et la Lune par exemple, le satellite étant contrôlé par des moteurs. Son mouvement est décrit par les équations de Gauss [START_REF] Battin | An introduction to the mathematics and methods of astrodynamics[END_REF]. L'objectif est de réaliser un transfert orbital ou un rendez-vous avec un autre engin spatial. Pratiquement le problème est souvent un problème de contrôle optimal où les critères sont le temps minimal ou la consommation de carburant. Le principe du maximum de Lev Pontryagin [START_REF] Pascoli | Eléments de mécanique céleste[END_REF] permet de sélectionner sous forme de conditions nécessaires, les candidats optimaux. Les contrôle optimaux sont calculés en maximisant un pseudo-Hamiltonien dépendant de l'état, du vecteur adjoint et du contrôle, pour calculer le Hamiltonien vrai dont les solutions produisent les trajectoires optimales. Avec la technologie à propulsion électrique le système est dit à poussée faible: le rapport entre la poussée maximale du satellite et le terme gravitationnel est petit. Dans le ces de deux corps, cette hypothèse permet de remplacer le vrai Hamiltonien par un Hamiltonien moyenné, voir par exemple [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF]. Les résultats classiques de convergence des trajectoires du moyenné vers le non moyenné pour les systèmes dynamiques sur des échelles de temps bien supérieures à nos temps de transfert permettent de justifier cette approximation. Notre travail permet une extension de ces résultats, notamment dans le cas de perturbations dues à une ou deux planètes (correspondant respectivement au problème de 2 ou 3 corps). Cet aspect est discuté dans les chapitres 1 et 2 de la thèse. Dans la suite de ce résumé on fait une synthèse de nos contributions présentées dans les chapitres 3, 4 et 5.

Comparaison de transferts coplanaires optimaux dans le cas de la minimisation du temps et de l'énergie L'application des techniques de moyennation pour le transfert orbital dans le cas des 2 corps a fait l'objet de nombreux travaux, citons ( [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF], [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF], [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF], [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF], [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF], [START_REF] Bonnard | Optimal Control with Applications in space and quantum dynamics[END_REF]) dans le cas de minimisation de la norme L 2 du contrôle (problème dit à énergie-minimale). L'objectif du chapitre 3 est de comparer ces résultats avec le cas du temps minimal. Dans ce cas le problème est très complexe, en particulier le Hamiltonien moyenné n'est plus explicite. On va donc considérer la situation la plus simple: le cas coplanaire et la cas d'un transfert vers un orbite circulaire, ce cas apparaissant souvent dans la pratique, pour un transfert vers un orbite géostationaire par exemple. Dans le cas du temps minimal le Hamiltonien moyenné est associé à une métrique de Finsler [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF], contrairement au cas de l'énergie où il est associé à une métrique Riemanienne. De plus le Hamiltonien moyenné n'est plus lisse sur le domaine physique, mais simplement C 1 , ce qui rend le calcul des trajec-7 toires très complexe. L'ensemble des singularités restant cependant assez petit pour pouvoir être évité dans de nombreuses missions, notament pour un transfert vers des orbites voisines. Notre principale contribution est de montrer que dans le cas du temps minimal le domaine physique est géodésiquement convexe, contrairement au cas de l'énergie où pourtant la métrique Riemanninne est plate (mais le domaine physique dans ces coordonnées n'est pas convexe). L'extension de ce résultat au cas coplanaire général et au cas non coplanaire est ouvert, mais le problème est complexe.

Transfert orbital perturbé par la Lune

Le problème des deux corps est seulement une approximation et un modèle plus complexe doit être utilisé pour certaines missions tenant compte de l'effet dit J 2 lié a l'aplatissement de la Terre et l'influence de la Lune par exemple. Dans le deux cas, les forces de perturbations sont conservatives. La méthode de Lagrange [START_REF] Battin | An introduction to the mathematics and methods of astrodynamics[END_REF] permet de modéliser ces perturbations. Notre contribution extraite de notre article [START_REF] Bonnard | Lunar and J 2 perturbations of the metric associated to the averaged orbital transfer[END_REF] est le suivante. On considerè avec une hypothèses de coplanarité l'influence de la perturbation lunaire moyenné sur le Hamiltonien moyenneé non perturbé, dans le cas du problème de l'énergie et qui est associée a une métrique Riemannienne (voir section précédente). En tenant compte de la perturbation, le Hamiltonien moyenné est associé a un problème de navigation de Zermelo. La problème est plus complexe que le cas Riemannien, mais il existe des études, notamment dans le contexte du contrôle optimal, voir par exemple [START_REF] Ronwell | Satellite lifetimes in nearly circular orbits for various Earth-Atmosphere models[END_REF], [START_REF] Roy | Luni-solar perturbations of an Earth satellite[END_REF]. Notre contribution est alors de présenter des résultats numériques concernant le calcul des lieux conjugués et de coupure obtenus en utilisant le code Hampath [START_REF] Chaplais | Averaging and deterministic optimal control[END_REF]. Ce travail permet de caractériser globalement l'optimalité des trajectoires. Ces études numériques restent dans un contexte académique mais peuvent aisément se généraliser en tenant compte de l'effet J 2 et d'un situation non coplanaire. Une question aussi intéressante est le problème de moyennation de toutes les perturbations et donc étudier l'effet des résonances. On l'étudie briévement dans notre travail dans le cas de la perturbation lunaire.

Transfert Terre-Lune et point de Lagrange L 1 Dans le cas d'un transfert Terre-Lune, le problème est associé au problème des trois corps et le calcul numérique de trajectoires optimales montre le passage par le point de Euler L 1 ou les 2 attractions se compensent. Cela conduit donc à observer que la trajectoire se décompose en trois phases: une première phase où la distance à la Terre est petite et où la perturbation lunaire est négligeable, une phase où la distance à la Lune est petite et où la perturbation terrestre est négligable et une phase intermédiaire au voisinage du point d'Euler L 1 . Pour les deux premières phases on peut remplacer le Hamiltonien par les Hamiltoniens moyennés respectifs qui permettent de réaliser une bonne approximation des trajectoires. Au voisinage du point d'Euler, on peut résoudre le problème de transfert en temps minimal, qui devient local et donc plus facile. On présente des résultats numériques décrivant la concatenation des trajectoires du moyenné avec la trajectoire temps minimum local à comparer avec une solution temps minimal calculée globalement par une méthode de tir. Modulo des conditions de recollement déduites du principe du maximum, ces solutions sont une bonne approximation des trajectoires temps minimals en restant bien plus simples à calculer numériquement, d'où l'intérêt de la méthode.

Enfin on discuté un contexte plus général où l'effet du Soleil est pris en compte et le modèle est un modèle 4 corps dit bi-circulaire. La problème est encore plus complexe numériquement pour appliquer une méthode de tir globale et les méthodes de moyennation sont donc très utiles sur certaines parties de la trajectoire, pour réduire selon le principe précédent la complexité des calculs numériques. The acceleration of a spacecraft (mass m) rotating around a celestial body of higher mass (mass M ) under the influence of gravity of the mass M , engine thrust and an external force F is

q = -µ q q 3 + u m + F m , u ≤ ǫ, q, u ∈ R 3 (1.1)
where q is the position of the satellite in an inertial frame with the origin at the primary celestial body (a vector joining the center of the satellite to the center of the Earth) and µ is the standard gravitational parameter µ = GM where G is the gravitational constant. The variable u is the force produced by the engines. It is the control, i.e. its direction and magnitude are to be decided to achieve some objective. The magnitude of this force is obviously limited by the capabilities of the engine; we denote this bound by ǫ > 0 . We shall see in section 1.1.5 to what extent it is 'small'. The force F is any external force on the satellite other than the control force and the gravity of mass M . The model in (1.1) assumes that the mass of the primary is concentrated at one point (or that it has spherical symmetry), hence the Newtonian gravity field.

Although dimensionless quantities are often used in constructing (1.1), the usual units are as follows: the entries of q are measured in meters (m), those of q in ms -2 , those of u and F in Newtons (N), the mass m in kilograms and of µ in m 3 s -2 .

Newtonian 2-body motion

We first recall the classical Keplerian motion, described by the equation

q = -µ q q 3 , (1.2)
that is, (1.1) with u = 0, F = 0. It is well-known that the following quantities (functions of (q, q)) are integrals of motion:

Energy K = 1 2 q 2 -µ q Angular momentum C = q × q Laplace vector L = q×C µ -q q .
(1.3)

Since K is a scalar and C and L are vectors of dimension 3, these are altogether 7 scalar quantities; they cannot be independent, and are indeed linked by the two relations

µ 2 ( L 2 -1) = 2 C 2 K L • C = 0. (1.4)
It is also well-known that the trajectory in (q, q) space projects on a conic section in R 3 with one locus at the origin for the position q that is -an ellipse if K < 0

-a parabola if K = 0 -a hyperbola if K > 0
that degenerates into a line segment if C = 0. These are the satellite orbits. We are mostly interested in the elliptic motions that take place in the elliptic domain

C = {(q, q) ∈ R 3 × R 3 , K < 0, C = 0}. (1.5)
It is interesting to make a change of coordinates on (q, q) such that the five first coordinates are independent first integrals of the motion (5 functions of the 7 scalars defined by (1.3)). If these coordinates are (x 1 , x 2 , ..., x 6 ), then (1.2) reads ẋ1 = 0, ..., ẋ5 = 0, ẋ6 = a function of (x 1 , x 2 , ..., x 6 ).

(

1.6)

There are many possible choices for these 5 first coordinates, often called 'orbital elements' because they define the orbit that would be described without control or perturbation. A possible choice is a = µ 2|K| the semi-major axis e = L the eccentricity (1.7)

and (Ω, i, ω) the Euler angles that send the bases {(1, 0, 0), (0, 1, 0), (0, 0, 1)} onto

L L , C C × L L , C C .
One may take as the last time-dependent coordinate x 6 the angle v that is between L and q in the plane orthogonal to C. Note that these do not make sense when L = 0, i.e. for e = 0, the circular orbits.

In terms of these orbital elements, orbits in the domain C are ellipses with radius

r = a(1 -e 2 )
1 + e cos v .

(1.8)

In the coordinates (a, e, i, ω, Ω, v), equation (1.2) reads

ȧ = 0, ė = 0, i = 0, ω = 0, Ω = 0, v = √ µ a 3/2 (1 -e 2 ) 3/2 (1 + e cos v) 2 .
(1.9)

Gauss equations

The Gauss equations are obtained by writing (1.1) in the coordinates (a, e, i, ω, Ω, v) (or another choice where the 5 first are another choice of orbital elements). This requires decomposing the acceleration u m or F m in a suitable frame. There are two common choices of frames used to describe the Gauss equations in the literature.

The radial-orthoradial frame is denoted by the radial, circumferential, and normal directions {i r , i θ , i c } where the radial direction is given by i r = q q , the orbit normal is given by the cross product i θ = q× q q× q and the third unit vector is perpendicular to both, i c = ir×ic ir×ic . In the tangential-normal frame {i t , i n , i c }, the tangential vector i t is in the plane of the osculating orbit and directed along the velocity vector, and i n is perpendicular to i c and i t . We suppose in the rest of this section that F is not present, and deal only with the control u. All external forces on the satellite are treated the same way, however, if perturbations are included, and one may read 'u + F ' instead of 'u' in the following equations if perturbing forces are taken into account.

We will make use of the following expressions of the Gaussian perturbation equations in this thesis. We remain in the elliptic domain, where these equations are defined. In the tangential-

normal frame {u t , u n , u c }, from [4, pp. 488-489],                          ȧ = 2a 2 √ 2a-r √ ar √ µ u t m ė = 2 √ ar(e+cos v) √ 2a-r √ µ u t m -r sin v a u n m i = r cos(v+ω) h u c m ω = √ ar ( 2e+ r cos v a ) e √ 2a-r √ µ u n m + 2 √ ar sin v e √ 2a-r √ µ u t m -r cot i sin(v+ω) h u c m Ω = r csc i sin(v+ω) h u c m v = h r 2 - √ ar ( 2e+ r cos v a ) e √ 2a-r √ µ u n m -2 √ ar sin v e √ 2a-r √ µ u t m , (1.10)
where h = aµ(1 -e 2 ), while in the radial-orthoradial frame {u r , u θ , u c }, they are

                 ȧ = 2a 2 e sin v h u r m + 2a 2 h p r u θ m ė = p sin v h u r m + (p+r) cos v+re h u θ m i = r cos θ h u c m ω = -p cos v he u r m + (p+r) sin v he u θ m + -r sin θ cos i h sin i u c m Ω = r sin θ h sin i u c m v = h r 2 + p cos v eh u r m -(p+r) sin v eh u θ m , (1.11)
where

θ = v + ω, p = a(1 -e 2 ).
If we assume that the orbits are circular, then e = 0, a singularity for the Gauss equations. To remove this singularity, we replace e and ω by the elements e x = e cos(ω + Ω), e y = e sin(ω + Ω). Furthermore, in the presence of an external force, for small e the orbital element a may change significantly during one period of the satellite's motion. Since in further sections (to apply averaging, as described in section 1.3.5) we will assume that the perturbed orbital elements are constant over an orbital period, it is desirable to use p = a(1 -e 2

x -e 2 y ) instead of a as one of the orbital elements because its change over one revolution is less than that of a for small eccentricity. The true longitude is the angular variable ℓ = v + ω + Ω which can be used to take the place of the variable v. The Gauss variational equations in these elements (p, e x , e y , i, Ω, ℓ) in the radialorthoradial frame ( [4, pp. 492-493]) are

                           ṗ = 2pu n √ p µ (1+e y cos ℓ+e x sin ℓ)m ėx = p µ u r m sin ℓ + u n (e x +sin ℓ(2+e y cos ℓ+e x sin ℓ)) m(1+e y cos ℓ+e x sin ℓ) u c m - e y √ p sin(ℓ-Ω) tan ( i 2 ) √ µ ėy = p µ -u r m cos ℓ + u n (e y +cos ℓ(2+e y cos ℓ+e x sin ℓ)) m(1+e y cos ℓ+e x sin ℓ) u c m + e x √ p sin(ℓ-Ω) tan ( i 2 ) √ µ i = √ p cos(ℓ-Ω) √ µ(1+e x cos ℓ+e y sin ℓ) u c m Ω = √ p sin(ℓ-Ω) csc i √ µ(1+e x cos ℓ+e y sin ℓ) u c m l = √ µ(1+e x cos ℓ+e y sin ℓ) 2 p √ p + √ p sin(ℓ-Ω) tan ( i 2 ) √ µ(1+e x cos ℓ+e y sin ℓ) u c m .
(1.12)

Occasionally it is also useful to use the variable n = µ/a 3 in the place of a or p. The Gauss equations (1.10), (1.11) or (1.12) may be written in a generic way as

ẋ = u 1 G 1 (x, ℓ) + u 2 G 2 (x, ℓ) + u 3 G 3 (x, ℓ) l = Q(x, ℓ) + u 3 g 3 (x, ℓ), (1.13)
where x = (a, e, ω, i, Ω), x = (p, e x , e y , i, Ω) or x = (n, e, ω, i, Ω), the variables (u 1 , u 2 , u 3 ) are the coordinates of the control vector in one of the frames , G 3 are three 5-dimensional vectors and Q, g 1 , g 2 , g 3 four scalars whose expression can be read in (1.10), (1.11) or (1.12). Using matrix notation, (1.13) may be written as

(i t , i n , i c ) or (i q , i θ , i c ) i.e. (u 1 , u 2 , u 3 ) = (u t , u n , u c ) or (u 1 , u 2 , u 3 ) = (u q , u θ , u c ) and G 1 , G
ẋ = G(x, ℓ)u, l = Q(x, ℓ) + Q 1 (x, ℓ)u, (1.14)
where G and Q 1 are 5 × 3 and 1 × 3 matrices, and u the control vector (the colums of G are G k ). Note that although G depends on m, since for our purposes m will always be constant (the derivative ṁ is proportional to u , so for u ≤ ǫ this assumption is feasible) we do not state m as a variable in (1.14).

These equations (1.10) can also be expressed using the so-called Lagrange equations when the force u takes the form of the gradient of a potential P. These equations are given in ( [4, p. 483])

1.1. MODEL in terms of the variables x = (a, e, ω, i, Ω), M , where M = nt + M 0 is the mean anomaly ( t is the time since epoch and M 0 is the mean anomaly at epoch). This variable increases uniformly from 0 to 2π radians during each orbit, but is not an angle; due to Kepler's second law it is proportional to the area swept by the focus-to-body line since the last periapsis.

                     ȧ = 2 na ∂P ∂M ė = 1-e 2 na 2 e ∂P ∂M - √ 1-e 2 na 2 e ∂P ∂ω i = cot i √ a √ 1-e 2 √ µ ∂P ∂ω Ω = cot i na 2 √ 1-e 2 ∂P ∂i ω = √ 1-e 2 na 2 e ∂P ∂e -cot i na 2 √ 1-e 2 ∂P ∂i Ṁ = n -2a na ∂P ∂a - ( -1+e 2 )
na 2 e ∂P ∂e .

(1.15)

1.1.3.

Planar case

If all transfers take place between orbits in the same plane, we may work in this plane only. Then q and q in (1.1) have dimension 2 instead of 3 -The orbital elements are (a, e, ω), (p, e x , e y ) or (n, e, ω) only.

-There is no out-of-plane control, i.e. one may use (1.10), (1.11) or (1.12), taking u c = 0 in the right-hand sides -They translate into (1.13) or (1.14) where u has dimension 2 only (i.e. u c = 0), the vectors

G 1 , G 2 , G 3 have dimension 3 and G is a 3 × 2 matrix.
In particular, (1.13) reads

ẋ = u 1 G 1 (x, ℓ) + u 2 G 2 (x, ℓ) l = Q(x, ℓ) (1.16)

Perturbations of the Newtonian 2-body problem

In equation (1.1), we add the term F different from the control because the Newtonian acceleration -µ q q 3 is only an approximation of the accelerating force that a spacecraft is subject to when the control is zero. Indeed it assumes that the distribution of the mass of the central body has spherical symmetry and that no other force acts on the spacecraft, however, the Earth (and most celestial bodies) does not have this symmetry; -µ q q 3 is only the dominant term of the gravitational potential of the so-called geocentric reference ellipsoid (described more fully in section 1.1.6) which more properly describes the gravitational potential by taking the actual shape of the Earth into account. Besides this, there is also the gravitational potential of other bodies present in the space environment which also affect the satellite motion, as well as features such as atmospheric drag, light pressure, Earth's magnetic field and the presence of cosmic particles which influence satellite motion. While Earth's gravity still dominates the satellite motion, these collective forces make up the perturbations of satellite motion, described by the force F .

The relative importance of these perturbations depends on the orbit. For example, the gravitational influence of third bodies acts under an inverse square law of the distance of the planet from the satellite, meaning that the further away the third body is, the less significant the third-body influence. Conversely, the oblateness of the Earth will affect the satellite to a greater degree the closer it is to the Earth's surface. Satellite orbits around the Earth are generally classified into the low-Earth orbits (LEO), medium-Earth orbits (MEO) and high-Earth orbits (HEO) according to their orbital altitude above the Earth's surface. A non-exhaustive table of perturbations which dominate the orbital motion in each altitude is given in table For our purposes, we define an external force on a satelite as a perturbation when the value of the ' thrust ratio' [START_REF] Dargent | Initial and final boundaries transformation when solving optimal control problem with averaging techniques and application to low thrust orbital transfer[END_REF], the ratio of the acceleration F m of the perturbation and the Earth gravitational acceleration µ r 2 , satisfies the inequality

F r 2 µ m ≤ ǫ ratio , (1.17)
where ǫ ratio is a suitably small number (we set ǫ ratio = 0.0005 for our purposes; [START_REF] Dargent | Initial and final boundaries transformation when solving optimal control problem with averaging techniques and application to low thrust orbital transfer[END_REF] uses 0.001).

We note that at some altitude within the HEO orbits, the n-body perturbations cease to be perturbations as the satellite enters the zone where full three-body motion is appropriate due to the increase of the value of the thrust ratio (through the increase of r in (1.17)). Of these perturbations, in this thesis we consider only LEO and HEO perturbations (the MEO case can be thought of as a mix between the two) and concentrate on the oblateness perturbation of LEO (specifically, the so-called 'J 2 perturbation', which we discuss in section 1.1.6) and the thirdand forth-body perturbations of HEO. These are generally considered the dominant perturbations in these zones. Further, the third and forth body perturbation can be considered an approximation (appropriate while the thrust ratio (1.17) remains low) of full four-body motion, which gives a way to bridge the gap between full three and four body motion and the two-body problem. Finally, we note that orbits such as the highly-elliptic Molniya orbit, which has an orbit with an inclination of 63.4 degrees, a perigee altitude of 1000km -qualifying as LEO -and an apogee altitude of 40000km, within the range of HEO (named after the first successful satellite to use an orbit with these specifications, Molniya 1-01, 1965) can experience both LEO and HEO perturbations at once, so a superposition of these perturbations also has some interest.

Electric propulsion and low thrust

In equation (1.1), we made use of the condition u < ǫ without justification of the fact that ǫ is small. In this thesis, we consider the case of electrically powered spacecraft. Electric propulsion engines electrically accelerate the propellant to higher energies and expel it at high speed. These are the industry's current state-of-the-art and have been used in many recent missions such as the SMART-1 mission (2006), FalconSat-3 (2007), Boeing-702SP (2015) and Microscope (2016).

The specific impulse I sp is a measure of the impulse produced per unit of propellant expended. Electric propulsion engines tend to have a high specific impulse (related to the high exhaust veolocity) and a low fuel consumption rate, which makes satellites equipped with these engines much more fuel-effective.

In our computations, we will define a 'weak thrust' by replacing F /m with u /m in the inequality (1.17). This gives the restriction on thrust magnitude

u ≤ ǫ ratio µ m r 2 . ( 1.18)
Again, we will set ǫ ratio = 0.0005. The inequality in (1.18) indicates that the thrust is a perturbation the way we have defined it in section 1.1.4. We will use the notation for the thrust vector

u = u max v, v ≤ 1
where u max = ǫ, a small parameter.

The J 2 perturbation

Although in equation (1.2), we assumed that the gravitational force acts as if the entire mass of the Earth is concentrated at the center, i.e. as if the Earth is spherical, however, in reality the Earth has an oblate ellipsoid shape. In the analysis of motion particularly of LEO satellites, the assumption of a spherical Earth is no longer valid, and it is necessary to consider the Earth's gravitational potential as that of a body with internal structure rather than a point mass. Here we give a brief description (based mainly on [START_REF] El'iasberg | Introduction to the Theory of Flight of Artificial Earth Satellites (Vvedenie Y Teoriu Poleta Iskusstvennykh Sputnikov Zemli)[END_REF] and also on [START_REF] Kozai | The Earth gravitational potential derived from satellite motion[END_REF]) of the construction of a geopotential model of the Earth as an oblate ellipsoid, sufficient to explain the 'J 2 -potential' which we will use in later chapters.

The geocentric reference ellipsoid is an oblate ellipsoid of revolution with an axis along the Earth's spin axis, so that the gravity field of the reference ellipsoid is symmetric about the Earth's axis. Since the mass distribution is the same with respect to the axis of rotation, the potential does not depend on the geocentric longitude, and in terms of the spherical harmonic functions, the normal potential takes the form

U (r, θ) = - µ r 1 + Nz n=2 J n P 0 n (sin θ) ( r R ) n , (1.19)
where J n are dimensionless coefficients and P 0 n (x) are the Legendre polynomials. The functions

P 0 n (sin θ) r n+1 n = 0, 1, 2, (1.20)
are known as the zonal harmonics.

To make the normal field represent the gravitational field of the Earth, the coefficients J n are measured by comparing true satellite motion with that predicted by a spherical potential, using the tracking of satellite arrays. This was first done by O'Keefe, Eckels and Squires [START_REF] O'keefe | The gravitational field of the earth[END_REF] who discovered the value of the J 2 -coefficient by comparing the true motion of the satellite Vanguard 1 to that predicted by classical geodesy with a spherical model. Jefferys set this dimensionless coefficient at 1.632 × 10 -3 in 1959, whereas we now know the J 2 value to be at around 1.0833 × 10 -3 [START_REF] Jeffreys | The earth. Its origin, history and physical constitution[END_REF], [START_REF] Kozai | The Earth gravitational potential derived from satellite motion[END_REF]. The subsequent values J 3 , J 4 ,... for the Earth are of the order of 10 -6 or less.

This potential (1.19) alone is generally sufficient to describe the Earth's deviation from a sphere. However, the geocentric reference ellipsoid is still an approximation to the Earth's figure. Geoid undulations, i.e. the deviations of the geoid from the reference surface, are introduced for the exact descriptions of the geoid figure. These result in additional terms to represent this deviation, known as the gravity anomalies, described by the tesseral harmonics. However, we are interested in the geopotential of an oblate Earth; thus we consider the model (1.19).

In fact, for our work it is sufficient to set all other terms besides J 2 in (1. [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF] to zero. This gives the potential

U (r, θ) = µ r - 3 2 R 2 r 3 J 2 (sin 2 θ - 1 3 ) R 1 (1.21)
where R 1 is the perturbing potential due to the Earth's oblateness in terms of the Kepler orbital elements. Using the fact that for the geocentric latitude, sin θ = sin(ω + v) sin i in the Kepler orbital elements,

R 1 = R 2 (1 + e cos v) 3 2a 3 (-1 + e 2 ) 3 J 2 -1 + 3 sin 2 i sin 2 (v + ω) . (1.22)
The averaged perturbation can be computed by

R1 = 1 2π 2π 0 R 1 dv, (1.23)
which gives

R1 = 3 2 R 2 J 2 a 2 (1 -e 2 ) 3/2 1 3 - 1 2 sin 2 i . (1.24)

Circular restricted 3-body problem

The circular restricted 3-body problem ( [START_REF] Battin | An introduction to the mathematics and methods of astrodynamics[END_REF], [START_REF] Szebehely | Theory of orbits in the restricted problem of three bodies[END_REF]) is three-body motion (the two primaries are the Earth and Moon) simplified by the assumptions -The two primaries move in circular orbits around their common barycenter -The satellite (third body) moves in the same plane as the two primaries Note that the previous assumptions mean that the orbit is planar, and we will further restrict to a two-dimensional transfer q ∈ R 2 . The circular restricted 3-body problem is generally expressed in coordinates which make use of a standard normalization fixing the Earth-Moon distance to 1. However, for the sake of consistency with, for example, perturbed 2-body problems, we will use non-normalized coordinates, where µ E is the standard gravitational parameter of the Earth µ M is the standard gravitational parameter of the Moon c 1 is the distance of the Earth from the Earth-Moon barycenter c 2 is the distance of the Moon from the Earth-Moon barycenter

-d 1 is the Earth-Moon distance (d 1 = c 1 + c 2 ) -ω b = µ E +µ M d 3 1
is the angular speed of the Moon Two frames are relevant to the circular retricted 3-body problem.

'Inertial' frame: It is the Cartesian frame centered at the Earth with fixed x and y directions (this frame is not really inertial for it has a translation movement with respect to the true inertial frame centered at the Earth-Moon barycenter). We denote by ( X, Y ) the coordinates in this frame Synodic frame: It is centered at the Earth-Moon barycenter and rotates with the angular speed of the Moon, so that coordinates of the Moon are (c 2 , 0) at all times. We denote by (X, Y ) the coordinates in this frame.

The coordinates ( X, Y ) and (X, Y ) are related by

X Y = cos(ω b t) -sin(ω b t) sin(ω b t) cos(ω b t) X Y - c 1 0 . (1.25)
In these coordinates, the satellite acceleration in the circular restricted 3-body problem is described by

Ẍ = 2ω b Ẏ + ω 2 b X -µ E X + X M r 3 1 -µ M X -X E r 3 2 (1.26) Ÿ = -2ω b Ẋ + ω 2 b Y -µ E Y r 3 1 -µ M Y r 3 2
.

(1.27)

where

r 1 = (X + c 1 ) 2 + Y 2 (1.28) r 2 = (X -c 2 ) 2 + Y 2 , (1.29)
and the position of the Earth and Moon in synodic coordinates are

(X E , Y E ) = (-c 1 , 0) (1.30) (X M , Y M ) = (c 2 , 0). (1.31) (1.32)
The Lagrange points are positions which occur in the problem of three bodies where the satellite can maintain an equilibrium position (in synodic coordinates) relative to the two primaries (e.g. [START_REF] Szebehely | Theory of orbits in the restricted problem of three bodies[END_REF]). This is because the Lagrange points mark positions where the combined gravitational pull of the two primaries are precisely equal. There are five such points, labeled L 1 to L 5 . The first three are on the line connecting the two primaries (known as collinear points) and the last two, L 4 and L 5 , form an equilateral triangle with the two primaries (known as triangular points). The configuration of these five points is illustrated in figure 1.1 (adapted from [START_REF] Szebehely | Theory of orbits in the restricted problem of three bodies[END_REF]). In our study, we will focus on the collinear points, particularly on the L 1 point, because it lies between the Earth and Moon and so is of interest on any Earth-Moon transfer. Determining the position of collinear Lagrange point in synodic coordinates involves solving for X in (1.26), where Ẍ = Ÿ = Ẋ = Ẏ = 0, i.e. the fifth-degree polynomial

-(1 -µ)(1 -µ + X) 2 + µ(µ + X) 2 + X(1 -µ + X) 2 (µ + X) 2 , (1.33)
where µ = µ M µ M +µ E ; the position of L 1 is the only root in the interval [0, d 1 ]. In order to be sure to use consistent constants without solving (1.33), we set the L 1 -point location to X = d 1 β, where β = 0.8491 (chosen from the values for the L 1 position determined 1.1. MODEL using the roots of the characteristic equation (1.33) given in [START_REF] Szebehely | Theory of orbits in the restricted problem of three bodies[END_REF]) and use (1.33) to solve for the standard gravitational parameter µ M of the Moon in terms of the values β and µ E . Note that µ in (1.33) is a distance measure in the standard normalized synodic coordinates (and in fact gives the position of the Earth in these coordinates) and the unknown distance X from the Earth to the Lagrange point L 1 is then given by X = µ -β in these coordinates, hence substituting X = µ -β into equation (1.33) gives

β 2 (1 -β 2 )(µ -β) + (1 -µ)(1 -β) 2 -µβ 2 = 0, (1.34) hence µ = µ M µ M + µ E = (1 -β) 2 (1 -β 3 ) 1 -2β + β 2 + 2β 3 -β 4 .
(1.35)

Thus, for β = 0.8491 and µ E the standard gravitational parameter µ E = 3.9860044189×10 14 m 3 /s 2 , we may solve for µ M from (1.35). We in turn determine the corresponding value of the Moon's mass M M by M M = µ M /G (where G is the gravitational constant) so that all choices of variables correlate throughout all of the numerical computations in chapter 5.

Bicircular four-body motion

The bicircular problem ( [START_REF] Gómez | II Fundamentals: the case of triangular libration points[END_REF]) is full four-body motion (the three primaries are the Sun, Moon and Earth) simplified using the hypotheses (a) Two primaries (Earth and Moon) move in circular orbits around their mutual center of mass and are coplanar.

(b) The third primary (Sun) is in a circular orbit around the center of mass of the system formed by the first two primaries, and its orbit is coplanar with the orbits of those primaries.

The configuration of the three primaries based on (a) and (b) is shown in figure 1.2 (adapted from [START_REF] Gómez | II Fundamentals: the case of triangular libration points[END_REF]). We will further impose that the satellite moves in the same plane as the Earth, Moon and Sun, as in the restricted three-body problem. Thus, although the equations for the bicircular motion in [START_REF] Gómez | II Fundamentals: the case of triangular libration points[END_REF] are given in dimension three, we will set the third component to zero and work in the variables (X, Y, Ẋ, Ẏ ).

In [START_REF] Gómez | II Fundamentals: the case of triangular libration points[END_REF], the standard normalization fixing the Earth-Moon distance to 1 is also used to express the equations of motion of the bicircular problem in the synodic coordinates, but we again choose to re-express these equations in non-normalized coordinates, where µ E is the standard gravitational parameter of the Earth µ M is the standard gravitational parameter of the Moon µ S is the standard gravitational parameter of the Sun µ S R S is the distance between the Sun and the Earth-Moon barycenter The arrangement of the four bodies in the bicircular problem (adapted from [START_REF] Gómez | II Fundamentals: the case of triangular libration points[END_REF]).

In terms of these non-normalized coordinates the equations of motion in the synodic frame are

Ẍ = 2ω b Ẏ + ω 2 b X + µ S R 3 S -µ E X + X M r 3 1 -µ M X -X E r 3 2 -µ S X -X S r 3 3 (1.36) Ÿ = -2ω b Ẋ + ω 2 b Y - µ S R 3 S -µ E Y r 3 1 -µ M Y r 3 2 -µ S Y -Y S r 3 3 , (1.37)
where

(X E , Y E ) = (-c 1 , 0) (1.38) (X M , Y M ) = (c 2 , 0) (1.39) (X S , Y S ) = (A S cos((ω S -ω b )t -α 0 S ), A S sin((ω S -ω b )t -α 0 S )) (1.40)
and the variables

-ω S = µ S +µ M +µ E A 3 s
is the angular speed of the Sun α 0 S is the Sun position at epoch and the three distances from the satellite to the Earth, Moon and Sun, respectively, are given by

r 1 = (X + c 1 ) 2 + Y 2 (1.41) r 2 = (X -c 2 ) 2 + Y 2 (1.42) r 3 = (X -X S ) 2 + (Y -Y S ) 2 . (1.43)
Notice that if we set µ S = 0, the equations (1.36)-(1.37) describe the circular restricted three-body problem in the synodic frame (equations (1.26) -(1.27)). We will use these equations to describe the satellite motion in the numerical simulations described in chapter 5.

Control objectives

In the most classical case, we assume Newtonian force around the Earth (Earth is a point mass) and zero control: the satellite then describes an orbit, characterized by the five orbital elements (section 1.1.1). However, for a controlled spacecraft u = 0, we may design a control u to alter the spacecraft's position by one or more discrete changes in velocity for the purpose of fulfilling certain mission objectives. To perform a satellite transfer or rendezvous (the difference between transfer and rendezvous is described for our purposes in section 1.3.2), we require to design a control which moves the satellite from its initial point to the chosen final point. These objectives may also include a minimization criterion to be satisfied. In this case, designing the control for the transfer or rendezvous problem becomes a problem of optimal control.

Optimal control

Optimal control is the use of the thrust-force u to displace a satellite from a fixed initial point q(0) to a fixed final point or condition q(t f ) ∈ M 1 in such a way as to minimize some possible cost L. These costs include -Transfer time t f -t 0 -Fuel consumption

t f 0 |u|dt -Quadratic criterion t f 0 u 2 dt.
(1. [START_REF] Kechichian | Reformulation of Edelbaum's low-thrust transfer problem using optimal control theory[END_REF] In the case of fuel consumption or the quadratic criterion, final time t f has to be fixed. While the quadratic criterion is sometimes called the 'energy' of the transfer (e.g. in [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF], [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF]), it is not a particularly physical criterion in that it does not directly represent any property of the satellite. However, it is convenient due to the fact that minimizing such a criterion gives rise to a Hamiltonian associated to a Riemannian metric, ( [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF], [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF] ).

Control methods

The Pontryagin Maximum Principle

We refer to [START_REF] Pontryagin | The mathematical theory of optimal processes Translated from the Russian by[END_REF] and [START_REF] Bonnard | Optimal Control with Applications in space and quantum dynamics[END_REF]. Consider the generic control system

ẋ = f (x, u, t), x ∈ X
where X is an n-dimensional smooth manifold, t ∈ R and u ∈ B ⊆ R k , the control set (a set of bounded, measureable mappings valued in the control domain B). An optimal control problem on X associated with the control system ẋ = f (x, u, t) is the problem of finding relative to the given points x 0 , x t f the trajectory x(•) and control u(•) such that

ẋ = f (x, u, t), x ∈ X , (u 1 , u 2 , ..., u ℓ ) ∈ R ℓ x(0) ∈ M 0 x(t f ) ∈ M 1 J = t f 0 L(u(t))dt → min.
(1.45)

For example, (1.45) describes the minimum time problem when L(u) = 1 and t f is free, and the quadratic criterion is the case in which L(u) = u 2 and t f is fixed.

CHAPTER 1. PRELIMINARIES

The Hamiltonian of the optimal control problem (1.45) is the function

H(x, λ, u, t, λ 0 ) = λf (x, u, t) + λ 0 L(u), (1.46)
where λ is a vector of costate variables (the adjoint vector) of the same dimension as the state variables x, and λ 0 is either 0 or -1. The Pontryagin's maximum principle can be written as Theorem 1.1 ( [START_REF] Bonnard | Optimal Control with Applications in space and quantum dynamics[END_REF]). For the pair (x, u) : [0, t f ] → X × B to be a solution of the problem (1.45), it is necessary that there exists an adjoint vector function λ : [0, t f ] → R n such that the following equations are satisfied for almost every t by the triple (x(t), λ(t), u(t)):

1. ẋ = ∂H ∂λ (x(t), λ(t), t, u(t), λ 0 ), λ = -∂H ∂x (x(t), λ(t), t, u(t), λ 0 ) 2. λ 0 (t) is constant: λ 0 (t) ≡ 0 or λ 0 (t) ≡ -1 3. H(x(t), λ(t), t, u(t), λ 0 ) = H(x(t), λ(t), t, λ 0 ) (1.47)
where H(x, λ, t, λ 0 ) = max v∈B H(x, λ, v, t, λ 0 ).

(1.48)

The following initial and final conditions are satisfied:

x(0) ∈ M 0 , x(t f ) ∈ M 1 , λ(0) ⊥ T x(0) M 0 , λ(t f ) ⊥ T x(t f ) M 1 .
Moreover, d dt H(x(t), λ(t), t, λ 0 ) = ∂H ∂t (x(t), λ(t), t, λ 0 ).

(1.49)

In the case that the function f in (1.45) does not depend explicitly on time, the condition (1.49) has the simpler form,

H(x(t), λ(t), t, λ 0 ) on [0, t f ] is constant (1.50)
If the final time t f is free, we have the additional condition

H(x(t f ), λ(t f ), t f , λ 0 ) = 0. (1.51)
The Hamiltonian vector field with Hamiltonian H takes the form

- → H λ 0 (x, λ, t) = ∂H ∂λ , - ∂H ∂x (x, λ, λ 0 ).
The Hamilton's equations are the two components of this field

ẋ = ∂H ∂λ (x, λ, t, λ 0 ), λ = - ∂H ∂x (x, λ, t, λ 0 ). (1.52)
The curve (x(t), λ(t)) is known as an extremal pair. The extremal pairs corresponding to λ 0 ≡ -1 are normal extremal pairs, while those corresponding to λ 0 ≡ 0 are abnormal extremal pairs. The projection x(•) of the extremal pair is a geodesic or extremal trajectory.

Applying the maximum principle to the Gauss equations

Let us now apply the maximum principle to the equations (1.13). We will replace the term x in theorem 1.1 by (x, ℓ), where x is the five-dimensional variable x = (a, e, i, ω, Ω), x = (p, e x , e y , i, Ω) or x = (n, e, ω, i, Ω) and ℓ is the true longitude. Correspondingly, we replace λ in this theorem with (λ, λ ℓ ). In this section we will also focus exclusively on the case that F does not depend on time t; i.e. the Hamiltonian H is not time-dependent. Since F does not play a role in the definition of the optimal controls, and since t is included in our definition the maximum principle in theorem 1.1, no confusion results by neglecting the dependence of F on t in our notation in this section.

The criterion to minimize is

t f 0 L(u)dt.
We will concentrate on the costs L(u) = u 2 and L(u) = 1 (the quadratic 'energy' criterion and time).

The following result has been determined for extremals of the energy-minimum problem:

Lemma 1.2. [START_REF] Bonnard | Optimal Control with Applications in space and quantum dynamics[END_REF] Abnormal extremals do not exist for the Hamiltonian systems associated to minimum energy Kepler problem in two-input case.

Thus in our investigation of extremals for the time and energy-minimal problems with one and two inputs, we will similarly restrict to studying the normal extremals in the time-minimal case.

The Hamiltonians associated to time-and energy-minimal Hamiltonians then have the form

H(x, ℓ, λ, λ ℓ , t) = -L(u)+ k u k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ))+ k F k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ))+λ ℓ Q(x, ℓ), (1.53)
where F k are the components of the perturbing force F in the same frame as the control, and

H k (x, ℓ, λ) = λ, G k (x, ℓ) , (1.54)
where G is defined by equation (1.14).

Remark 1.3. The number of controls cannot be more than the dimension of q, but it may sometimes be less, for instance if technology forces the thrust to be only tangential, for example in sections 2.1.2 and 3.4. Remark 1.4. In the minimizing case, one may use the constraint u ≤ 1, i.e. k u 2 k ≤ 1. In the 'energy' case, it is simpler to put no constraint on u (this gives a smooth dependence of the minimizing u * on x, ℓ, λ, λ ℓ ) and fix a 'sufficiently large' t f , so that the small parameter is not ǫ but 1 t f . Let us detail more what happens when L(u) = u 2 and L(u) = 1 (section 1.2.1). The maximization (1.48) is explicit; consequently, the Hamilton's equations (1.52) become

ẋ = ∂H ∂λ (x, ℓ, λ, λ ℓ , u) l = ∂H ∂λ ℓ (x, ℓ, λ, λ ℓ , u) λ = -∂H ∂x (x, ℓ, λ, λ ℓ , u) λℓ = -∂H ∂ℓ (x, ℓ, λ, λ ℓ , u).
(1.55)

If L(u) = u 2 , and there is no bound on u, then the maximum is reached for

u k = u * k (x, ℓ, λ, λ ℓ ) = 2 (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) , (1.56)
and

H(x, ℓ, λ, λ ℓ ) = (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 + F k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)). (1.57)
If L(u) = 1 (the minimum time case), with the bound

u 2 k = ǫ, then the maximum is reached, if (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 = 0 (1.58)
for

u k = u * k (x, ℓ, λ, λ ℓ ) = ǫ (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 (1.59)
and

H(x, ℓ, λ, λ ℓ ) = -1 + ǫ k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 + k F k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)).
(1.60) We discuss the surface on which the denominator

k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2
in the equation (1.59) goes to zero in the result (2.17) in section (2.1.2). The existence of singular extremals in the time-minimal case is covered in [START_REF] Bonnard | Geometric optimal control of elliptic Keplerian orbits[END_REF]. There do exist singular extremals in the single-input time-minimal case; however, since they do not correspond to optimal curves, we choose not to pursue a study of such extremals. The boundary and transversality conditions we use will be that x(0) and ℓ(0) are given (initial conditions). For final conditions, we will consider either a rendezvous or a transfer: for a rendezvous, the final values of both x(t f ) and ℓ(t f ) are given, and for a transfer, x(t f ) is given and ℓ(t f ) is free, which implies the transversality condition λ ℓ = 0. Note that the final condition may be looser: F(x(t f ), λ(t f )) = 0, F : R 6 → R r , r < 6 (4 instead of 6 in the planar case). This implies the transversality condition

(λ(t f ), λ ℓ (t f )) • dF(x(t f ), ℓ(t f )) = 0.
If, for instance, we are doing a transfer (not a rendezvous) in minimum time, then x(0), x(t f ) are fixed and ℓ(t f ), as well as t f , are free. Hence one has to add the transversality conditions H(x, ℓ, λ, λ ℓ ) = 0, λ ℓ (t f ) = 0 and also λ ℓ (0) = 0 if ℓ(0) is also left free.

Conjugate points

For the next two sections we will return to the notations (x, λ) of section 1.3.1 rather than (x, ℓ, λ, λ ℓ ).

Definition 1.5. Let z = (x, λ) be a reference extremal solution of

- → H on [0, t f ]. The variational equation δz(t) = d - → H(z(t))δz(t)
is called the Jacobi equation. A Jacobi field is a non trivial solution δz = (δx, δλ) of the Jacobi equation and it is said to be vertical at time t if δx(t) = 0. Definition 1.6. We define the exponential mapping exp(t H)(x 0 , λ 0 ) = (x(t, x 0 , λ 0 ), λ(t, x 0 , λ 0 )) and its projection exp x 0 ,t (λ 0 ) = Π(x(t, x 0 , λ 0 ), λ(t, x 0 , λ 0 ))

where Π is the projection on the first element, and λ 0 can be restricted to the sphere |λ 0 | = 1.

If z = (x, λ) is the reference extremal, a time t c > 0 is said to be conjugate to 0 if the mapping λ 0 → exp x 0 ,t (λ 0 ) is not of rank n -1 at t = t c (with n = dimQ) and the associated point x(t c ) is said to be conjugate to x 0 . We denote by t 1c the first conjugate time and C(x 0 ) is the conjugate locus formed by the set of first conjugate points.

Testing conjugacy An algorithm can be deduced which is implemented in the Hampath Code [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF]. Let z(t) = (x(t), λ(t)) be the reference extremal and consider the vector space of dimension n -1 generated by the Jacobi fields δz i = (δx i , δλ i ), i = 1, ..., n -1 vertical at t = 0 and such that δλ i (0) is orthogonal to λ 0 . At a conjugate time t c , one has rank[δx 1 (t c ), ..., δx n-1 (t c )] < n -1 or equivalently, det[δx 1 (t c ), ..., δx n-1 (t c ), ẋ(t c )] = 0.

The shooting method

The Pontryagin maximum principle (theorem 1.1) gives a Hamiltonian o.d.e. (1.52) in the state and adjoint variable, and conditions both at initial time t = 0 and final time t = t f , for the solution to be relevant to our optimal control problem with initial and final conditions. This is a boundaryvalue problem. It is well-known that the two-point boundary value problem is more complex than a Cauchy problem, i.e. an o.d.e. with prescribed initial conditions, which is known to have a unique solution (Cauchy-Lipschitz theorem) which can be computed numerically using more or less sophisticated integration schemes, for example Runge-Kutta 4,5. Often in optimal control the right hand side is not Lipschitz and so the Cauchy-Lipschitz theorem does not apply; special care in the numerics to account for this is sometimes necessary. Shooting methods are the most popular way to solve the boundary-value problem numerically.

One looks for solutions of (1.52), (x(t), λ(t)) with initial and final conditions, for instance

x(0) = x 0 F(x(t f ), λ(t f )) = 0, (1.61)
where F has a different form depending on what is specified. Obviously, (x(t f ), λ(t f )) is a function of x(0), λ(0):

(x(t f ), λ(t f )) = exp(t f H)(x(0), λ(0))
through a Cauchy problem. The equation (1.61) can be written as an equation with respect to the (unknown) final time t f and λ(0) = λ 0 , i.e. one has to solve S(t f , λ 0 ) = 0 where S : (t f , λ 0 ) → F(exp(t f H))(x 0 , λ 0 ).

(1.62) This is the so-called 'shooting equation' associated to the optimal control problem (1.45) where the manifolds M 0 and M 1 are replaced by the conditions (1.61). For an 'initial guess' (t f , λ 0 ), S(t f , λ 0 ) can be computed for instance using a Runge-Kutta scheme. For solving S(t f , λ 0 ) = 0 numerically in this study, we use the matlab function tdfsolve developed by Thierry Dargent, which is an adaptation of the matlab fsolve function to include a constant parameter in the input function, which makes use of the Newton's method with either Cauchy steps (a step along the steepest descent direction) or Gauss-Newton steps (or a convex combination of these) to determine the roots of S starting from the given initial guess.

Averaging

Averaging is a classical technique to treat perturbations of conservative systems (e.g. integrable Hamiltonian systems) that admit periodic or quasi-periodic trajectories [START_REF] Arnold | Supplementary chapters to the theory of ordinary differential equations[END_REF]. The general idea is that it is possible to select coordinates such that the dynamics of a perturbed system with periodic or quasi-periodic trajectories have a 'fast' and 'slow' component in the sense that the two dynamics develop in different time scales. These dynamics are often referred to as 'fast' and 'slow' dynamics, respectively. This existence of two time scales makes solving such a system more complex both analytically and numerically; analytically, the equations must be solved with respect to both time-scales, while numerically, the fast time scale introduces 'rapid oscillations' into the solution, which can cause problems for numerical solvers by introducing plural local solutions which pose a problem when seeking to identify an optimal trajectory (we discuss this further in section 2.3).

The underlying principle of averaging is to eliminate the fast time to obtain a simpler system. Different specific cases of averaging have been extensively studied ( [START_REF] Bogoliubov | Les methodes asymptotiques en theorie des oscillators nonlineares[END_REF], [START_REF] Volosov | Some types of computations in the theory of nonlinear oscillators which involve averaging[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]); here we give a summary of one general and a more specific case which can be applied to Keplerian systems, mainly following [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].

Generalities

Consider any dynamical system which may be written in the form

dx dt = ǫf (ǫt, t, x, ǫ) = ∞ i=1 ǫ i f i (ǫt, t, x) (1.63)
where x is a vector of dimension n and ǫ is a small parameter compared to the size of the others. This system is by definition non-autonomous (since the second member of the differential equation depends explicitly on time) and more precisely it depends on a 'slow time' ǫt and a 'fast time' t. Note that the infinite series is generally limited to some number of terms, and the others are supposed zero. We will consider the case where the function f is periodic, because such a system is similar to the one which occurs in the perturbed two-body problem.

We search for a way to establish a system of the form

dx dt = ǫF (ǫt, x, ǫ) = ∞ i=1 ǫ i F i (ǫt, x) (1.64)
where x is related to x by the asymptotic expansion

x = x + ǫU (ǫt, t, x, ǫ) = x + ∞ i=1 ǫ i U i (ǫt, t, x).
(1.65)

The development (1.65) is interpreted as the decomposition of the 'real' motion described by x to an 'averaged' motion described by x, and weak perturbations (corresponding to the oscillations in this periodic case) of this averaged motion described by U . We will consider only the first term in the infinite series (1.64) and (1.65), an asymptotic approximation of order 1 of the system (1.63). We may then determine the terms F 1 and U 1 terms in the 'classical' way (e.g. [START_REF] Bogoliubov | Les methodes asymptotiques en theorie des oscillators nonlineares[END_REF], [START_REF] Volosov | Some types of computations in the theory of nonlinear oscillators which involve averaging[END_REF]), by substituting x in the expression (1.65) into (1.63) where the successive powers of ǫ are developed in (1.64) until the order k; we then equate the terms in ǫ 1 between the two equations. This first asymptotic approximation of the initial system, is the average we seek:

dx dt = ǫf 1 (ǫt, t, x), (1.66)
In the case that f is 2π-periodic with respect to t, this term turns out to be equivalent to the classical average of an 2π-periodic function φ(τ, x),

φ(τ, x) = 1 2π 2π 0 φ(τ, t, x)dt. (1.67)
One may use averaging on these o.d.e.'s of a more specific form

dx dt = ǫf (x, y, ǫ) = ∞ i=1 ǫ i f i (x, y) dy dt = g(x, y, ǫ) = g 0 (x, y) + ∞ i=1 ǫ i g i (x, y) (1.68)
where the zero-order term in ǫ in the equation of the 'fast movement' y does not depend exclusively on the movement x. In this case, the 'fast' variable is y, and the 'slow' variable is x. We will further assume that the fast movement is a scalar; the method of averaging consists of replacing the independent variable t by y in (1.63) and dividing the equation associated with x by that associated with y (under the qualification that g(x, y, ǫ) is nonzero) and considering τ = ǫt as a new (slow) state variable. We obtain in this way the standard system

dx dy = ǫ f (x,y,ǫ) g(x,y,ǫ) = ∞ i=1 ǫ i r i (x, y) dτ dy = ǫ 1 g(x,y,ǫ) = g 0 (x, y) + ∞ i=1 ǫ i q i (x, y) (1.69)
for which the functions r 1 (x, y) f 1 (x,y)

g 0 (x,y) r 2 (x, y) f 2 (x,y) g 0 (x,y) -f 1 (x,y)g 1 (x,y) g 0 (x,y) 2 . . . q 1 (x, y) 1 g 0 (x,y) q 2 (x, y) -g 1 (x,y) g 0 (x,y) 2 . . . (1.70)
and the averaged system has the form [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]

dx dȳ = ǫf 1 (x)/g 0 (x) (1.71) dτ dȳ = ǫ1/g 0 (x). (1.72)
We then obtain the derivative dx/dt as

dx dt = f 1 (x)/g 0 (x) 1/g 0 (x) . (1.73)
Averaging in optimal control of the Kepler problem Conforming the system (1.55) to the form of equation (1.68) is not straightforward. Obviously, in (1.55), x is slow and ℓ is 'fast', but it is not so easy to sort out which components of the adjoint vector (λ, λ ℓ ) are fast or slow. Hence the y-variable in (1.68) must contain at least ℓ and the x variable in (1.68) contains x from (1.55) and part of λ and λ ℓ . We do not give any details of the proof but, to compute the average system that accounts for the movement of (x, λ) one may assume that x and λ are slow variables and take λ ℓ equal to zero; this amounts to saying that the average system of (1.55) for transfer (unprescribed ℓ(0) and ℓ(t f )) is the Hamiltonian equation associated with the Hamiltonian

H(x, λ) = 1 2π 2π 0 ̟(x, ℓ)H(x, ℓ, λ, 0)dℓ, (1.74)
where

̟(x, ℓ) = 1 Q(x, ℓ) 1 2π 2π 0 dℓ Q(x, ℓ)
.

(1.75) Note: if Q does not depend on ℓ, one would have ̟(x, ℓ) = 1. For example, we may replace ℓ with the mean eccentric anomaly E -e sin E, for which E -e sin E = nt; then replacing y = ℓ by y = E-e sin E n , we have ẏ = 1 which is independent of E and so in this case ̟(x, y) = 1.

Chapter 2

Contributions

The Kepler motion described in (1.2) is a very particular constrained case of satellite motion in the space environment. Within this study we investigate the satellite motion in a wider space context, by studying specific aspects of not only satellite transfers with dynamics (1.2) but also when these dynamics are affected by naturally-occuring forces F in the space environment; both when these forces are perturbations (described in sections 1.1.4 -1.1.6) and when the external forces have a more significant effect (as in the models in sections 1.1.7 and 1.1.8). Thus we extend the two-body problem (1.2) naturally to the problem of three (and four) bodies by studying satellite transfers which initially include perturbations, and finally which require us to consider a three and four-body model. In this section we will detail the contributions we made to each kind of lowthrust, time-optimal transfer (with two-body dynamics, perturbed two-body dynamics and three and four body dynamics) within the course of this study.

Planar two-body Kepler transfers

Here we investigate controlled Kepler motion, the case where no perturbation is present in the model (1.1) (F = 0) for planar transfers (q is 3-dimensional in (1.1)). We are concerned with minimizing time or the quadratic 'energy' criterion (L(u) = 1 or L(u) = 1 2 u 2 ) where the dynamics are given by (1.13) for k = 2 (planar case) and x is (n, e, ω) . As we have stated in section 1.3.2, for the sake of comparision between these two minimization problems, we will concern ourselves with only normal extremals; thus, according to equation (1.46), the Hamiltonian is then

H(x, ℓ, λ, λ ℓ , λ 0 , u) = -L(u) + k u k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ, λ)) + λ ℓ Q(x, ℓ) (2.1)
where

H i = λ, G i for G i the fields in (1.16).
We are interested in transfer rather than rendezvous, hence x(0) and incidentally ℓ(0) are fixed,

x(t f ) is fixed but ℓ(t f ) is free, hence λ ℓ (t f ) = 0 in the Pontryagin maximum principle.
In the minimum time case, the minimizing control (refer to theorem (1.1)) is given by (1.59) and the maximized Hamiltonian (which we denote by H time ), is given by (1.60) where F = 0:

H time (x, ℓ, λ, λ ℓ ) = -1 + ǫ k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 + λ ℓ Q(x, ℓ). (2.2)
Note that if the number of controls is only 1, the square root is simplified to

|H 1 (x, ℓ, λ) + λ ℓ g 1 (x, ℓ) + λ ℓ Q(x, ℓ)|
(the control dimension may be lower than the manifold dimension; see remark 1.3) . In the energy case, the (unbounded) minimizing control is given by (1.56), and the minimized Hamiltonian (which we denote by H energy ) has the form of (1.57) where F = 0:

H energy (x, ℓ, λ, λ ℓ ) = (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 + λ ℓ Q(x, ℓ). (2.3)
Let us call M the manifold where x evolves (in coordinates, it is R 3 in the planar case, and R 5 in the 3-D case). We restrain ourselves to the elliptic domain C (1.5), hence ℓ evolves on the bounded set S 1 (in non-negative energy, some values of ℓ send q, in (1.1), to infinity). These two Hamiltonians (2.2) and ( 2.3) give rise to Hamiltonian flows on T * (M × S 1 ), that can be expressed in the coordinates (x, ℓ, λ, λ ℓ ). In these equations

• ℓ is a fast variable, x is slow. We can also show that with some normalization, λ ℓ is small: when the final time t f is unknown, then H time (x, ℓ, λ, λ ℓ ) = 0 and by homogeneity we may express equation (2.2) as

0 = -ǫ + ǫ k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 + λ ℓ Q(x, ℓ). (2.4)
Thus, solving for λ ℓ gives

λ ℓ = ǫ (1 - k (H k (x, ℓ, λ) + λ ℓ g k (x, ℓ)) 2 ) Q(x, ℓ) (2.5)
and we can see that λ ℓ is small in this case. A similar argument can be carried out in the energy case. Since λ ℓ is zero at final time, we may assume that since λ ℓ is small, then it can be set to zero on trajectories of interest.

Finally, because λ ℓ is small, then λ is slow.

This remark is only heuristic. We obtain the averaged Hamiltonians as described in section 1.3.5, using equation (1.74), to give

Henergy (x, ℓ) = (1 -e 2 ) 3/2 2π 2π 0 k=1,2 H k (x, λ, ℓ) 2 dℓ (1 + e cos(ℓ -ω)) 2 (2.6) Htime (x, ℓ) = (1 -e 2 ) 3/2 2π 2π 0 k=1,2 H k (x, λ, ℓ) 2 dℓ (1 + e cos(ℓ -ω)) 2 .
(2.7)

Now we are concerned with studying the Hamiltonian flow associated with the Hamiltonians Henergy and Htime .

The quadratic case

The study [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF] by Bonnard, Caillau and Dujol discussed the properties of the Hamiltonian (2.3) and the flow of the associated Hamiltonian system. The first noteable aspect is that Henergy is smooth (real analytic). This is clear from the fact that it is a quadratic function. In the coordinates x = (n, e, ω), the Hamiltonian Henergy has the form

Henergy (x, λ) = 18n 2 λ 2 n + 5(1 -e 2 )λ 2 e + (5 -4e 2 ) λ 2 ω e 2 .
(2.8)

Secondly, they determine that it is associated to a Riemannian metric on M ,

g = 1 9n 1/3 dn 2 + 2n 5/3 5(1 -e 2 ) de 2 + 2n 5/3 5 -4e 2 dω 2 . (2.9)
Because ω is a cyclic variable (the Hamiltonian, and hence the right hand side of the Hamiltonians equations, does not depend on ω, and λ ω = 0 ⇒ ω = 0), one may consider a 'subproblem' in the variables (n, e) only. Geometrically, the condition λ ω = 0 is the transversality condition for a transfer towards a 'circular orbit', or more correctly any orbit for which the angle of the pericenter is unprescribed. Setting λ ω = 0 in (2.3), the associated metric becomes

g = 1 9n 1/3 dn 2 + 2n 5/3 5(1 -e 2 ) de 2 , (2.10)
which the authors discover is a Liouville metric with a linear first integral, and so the geodesic flow can be integrated using elementary functions. This integrability is a rare feature which immediately makes the flow totally understood. By choosing appropriate variables, indeed, this metric is a flat metric, and the geodesics are straight lines in the plane with polar coordinates (n, ψ) where 3nλ n = ρ cos ψ, cos ϕλ e = ρ sin ψ and sin ψ = e (the variable ρ, due to the fact that Henergy is homogeneous of degree 1, plays no role in the evolution of the optimal trajectories and so it is neglected). This gives a very clean and simple description of the flow. This flow however has a surprising feature; the elliptic domain is not geodesically convex for the Riemannian metric in question. Namely there are points (n int , e int ), (n final , e final ) such that no geodesic (or extremal curve in the optimal control language) contained within the elliptic domain joins them.

Structural properties of the flow: conjugate points

From the summary of the energy-minimizing flow in the circular restricted case given in (2.1.1) it is obvious that the cut and conjugate loci are empty, since the metric (2.10) in the right coordinates is the flat Euclidean one. When the transfers are not restricted to the circular case, the result is less clear, but the conjugate points in this case are discussed in [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF]. Here, they first use a change of variables n = (5ρ/2) 6/5 , e = sin r to transform the metric in (2.9) to a metric of the form

ḡ = dρ 2 + (ρ 2 /c 2 )g, (2.11)
where g has the form g = dr 2 + m 2 (r)dθ 2 .

(2.12)

They note that by homogeneity, it is possible to restrict the optimality analysis to the metric g with r ∈ [0, π/2]. This metric can be extended to an analytic metric on a two-sphere of revolution, where (r, θ) are spherical coordinates. The metric g turns out to be a specific case (q = 4) of the more general metric they wish to study in [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF]:

g q = dr 2 + m 2 q (r)dθ 2 , (2.13)
with m q (r) = √ q + 1 sin r/ 1 + q cos 2 r where q ≥ 0. The space M q is the two-sphere S 2 , and (g q , M q ) is a Riemannian metric. For this more general metric, the result on the cut and conjugate locus is proved:

Theorem 2.1. [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF] If q > 0, then for each point m of M q distinct from a pole, the cut locus of m is a subarc of the antipodal parallel to m and the first conjugate locus of m has exactly four cusps.

Since the metric g in (2.12) is the restriction of the metric (2.9) to the two-sphere of revolution in the case q = 4, then this theorem serves to describe the cut and conjugate loci of the optimal solutions for the unrestricted two-input case.

Time-minimal case

Studying the Hamiltonian flow associated to Htime is somewhat more difficult than Henergy for many reasons. Firstly, due to the presence of the square root in the Hamiltonian Htime , it is not smooth -this loss of analyticity occurs at points (x, ℓ) such that all functions H k in (2.2) vanish for at least one value of ℓ. This loss of smoothness is a key difference between the two Hamiltonians Htime and Henergy . Another difficulty, which requires fundamental difference of approach between the two problems, is the fact that the integral cannot be explicitly computed for the case of Htime as in the quadratic case.

Finsler geometry

There is a common feature between the Hamiltonian flows of Htime and Henergy : the variational problem associated with the average Hamiltonian comes from a rather classical geometry, namely Finsler geometry. Finsler geometry is the geometry of the Finsler metric, which is derived in much the same way as the Riemannian metric, as the data in each tangent space of a certain norm. In the Riemannin case, this norm is a positive quadratic form, whereas for a Finsler metric [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF], it is the data in the same tangent spaces of a norm that does not necessarily derive from an inner product (unit ball in the tangent space are not ellipsoids).

In our case, the Finsler norm itself is implicitly given by

v x = max H(x,ℓ)≤1 λ, v , (2.14)
and cannot be explicitly used as in [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF].

Non-smoothness

The points (n, e, ω, ℓ, λ n , λ e , λ ω ) where Htime is not smooth is obtanied by eliminating ℓ in

H 1 (n, e, ω, ℓ, λ n , λ e , λ ω ) = 0 (2.15) H 2 (n, e, ω, ℓ, λ n , λ e , λ ω ) = 0 (2.16)
in the case of two controls.

In the case of the transfer to circular orbits (λ ω = 0) (see section 2.1.1), this yields a hypersurface in T * M described by (3.34) where sin ϕ = e and ψ is the polar angle of (λ n , λ e ):

S = (ψ, ϕ), tan ψ = 1 + sin ϕ 2 cos ϕ ∪ (ψ, ϕ) , tan ψ = -1 + sin ϕ 2 cos ϕ .
(2.17)

In the case of one control, the situation is a bit different. The integrand k=1,2 H 2 k becomes k=1,2 H 1 , and it turns out that the equation

H 1 (n, e, λ n , λ e , ℓ) = 0 (2.18)
has, for fixed (n, e, λ n , λ e ) either no solution (this happens on an open set) or two distinct solutions (this happens on another open set) or one double solution (which happens on a closed set that is the common boundary of these two open sets). Htime is smooth on the two above-mentioned open sets and nonsmooth (although C 1 ) on the above-mentioned closed set. When reducing to the case of transfer to circular orbits in the variables (n, e), that closed subset is the same smooth hypersurface S given by (3.34).

Convexity

One point that we established is the convexity of the elliptic domain. The main result in the paper [START_REF] Bonnard | Time versus energy in the averaged optimal coplanar Kepler transfer towards circular orbits[END_REF] reproduced in chapter 3 is that, contrary to the quadratic case, the elliptic domain is geodesically convex for our problem. Namely, given initial and final values of (n, e), there is one extremal curve (projection of a solution of the Hamiltonian equation) that remains in the elliptic domain and joins them. This is true for one control and two controls.

Recall that averaging as defined in (1.74) makes sense only in the elliptic domain (this can be seen from the fact that (1.74) is defined for the case of periodic quasi-periodic solutions of the orbital motion, which we have already stated occur only in the elliptic domain ( orbits are ellipses, as opposed to hyperbolas and parabolas) and also since the function 1/Q(x, ℓ) in the denominator in (1.74) tends to zero as the integral leaves the elliptic domain). Thus convexity of the elliptic domain of the flow of Htime is truly convexity over the whole region where Htime (and so its flow) is defined.

Theorem 2.2. For any (n int , e int ) and (n final , e final ) in the elliptic domain X = {(n, e), 0 < n < +∞, -1 < e < 1}, there exist a time t f ≥ 0 and a solution t → (n(t), e(t), p n (t), p e (t)) of the associated Hamiltonian system

ṅ = ∂ Htime ∂p n , ė = ∂ Htime ∂p e , λn = - ∂ Htime ∂n , λe = - ∂ Htime ∂e defined from [0, t f ] to X , such that (n(0), e(0)) = (n int , e int ) and (n(t f ), e(t f )) = (n final , e final ).
Another interesting point in [START_REF] Bonnard | Time versus energy in the averaged optimal coplanar Kepler transfer towards circular orbits[END_REF] is that we were able to look at the quadratic problem in the same coordinates as the time-minimal problem, and understand why our arguments for convexity in the minimum-time case fail when applied to the minimum-energy case. Making the symplectic change of variables from (e, ω, λ e , λ ω ) to (n, ϕ, λ n , λ ϕ ), where e = sin ϕ (forπ 2 < ϕ < π 2 ) and λ ϕ = (1 -e 2 )λ e , in the Hamiltonian Henergy one can see from the phase portrait of these solutions shown in figure 3.3 that due to a continuum of equilibria that form a 'barrier', the maximum possible variation of the variable ϕ for some solution is 2 5 π. Thus, if |ϕ 0 | > ( 2 5 -1 2 )π, there are some values of ϕ that cannot be reached by any solution starting from the line {ϕ = ϕ 0 }. This obstruction disappears in the minimum time case, where the phase portrait is not as degenerate, with only two equilibrium points (figures 3.1, 3.2).

The convexity result demonstrates that, in contrast to the energy-minimization problem under the transfer to circular orbits, any final state may be reached time-optimally from any initial state both under two-input and one-input control in only the tangential direction.

An interesting further direction of study would be to remove the simplification involved by considering only the case of circular transfers, i.e. to consider a non-zero λ ω . This increases the dimension of the solution space, but would provide a more physically-realistic solution. A further generalization is in the non-coplanar case, i.e. the case where three control directions are applied.

Planar and non-planar perturbed two-body transfers

So far, we have only taken into account the Earth's gravity, i.e. F = 0 in the equation (1.1). As explained in sections 1.1.4 and those following, there are perturbations due to the Earth's asymmetry (J 2 perturbations) and other celestial bodies (as well as smaller perturbations due to light pressure, cosmic particles, etc. which we do not consider) which influence the satellite's motion. Let us investigate the Moon (third body) and J 2 perturbations.

Lunar perturbation of the metric associated to the averaged orbital transfer

We will consider the satellite dynamics described by

                                       ȧ = 2a 2 √ 2a-r √ ar √ µ ut m + 2 √ a √ µ ∂R ∂M ė = 2 √ ar(e+cos v) √ 2a-r √ µ ut m -r sin v a un m + 1-e 2 √ ae √ µ ∂R ∂M - √ 1-e 2 √ ae √ µ ∂R ∂ω i = r cos(v+ω) h uc m + cot i √ a √ 1-e 2 √ µ ∂R ∂ω Ω = r csc i sin(v+ω) h uc m + cot i √ a √ 1-e 2 √ µ ∂R ∂i ω = √ ar(2e+ r cos v a ) e √ 2a-r √ µ un m + 2 √ ar sin v e √ 2a-r √ µ ut m -r cot i sin(v+ω) h uc m + √ 1-e 2 √ ae √ µ ∂R ∂e - cot i √ a √ 1-e 2 √ µ ∂R ∂i Ṁ = n + √ 1-e 2 (-1+e 2 ) cos v e(1+e cos v) µ(1+e 2 +2e cos v) a-ae 2 un m - 2 √ 1-e 2 (a(e 2 -e 4 )+p+ep cos v) sin v ep(1+e cos v) µ(1+e 2 +2e cos v) a-ae 2 ut m -2a 2 a 3/2 √ µ ∂R ∂a + a(-1+e 2 ) a 3/2 e √ µ ∂R ∂e , (2.19
) which can be put in the form similar to (1.13),

ẋ = G(x, M )u + 3 k=1 G k (x, M )F k Ṁ = g 0 (x, M ) + 3 k=1 u k g k (x, M ) + 3 k=1 F k g k (x, M ) (2.20)
in which the colums of G are given by (2.19), where R is the lunar perturbing potential [START_REF] Pascoli | Eléments de mécanique céleste[END_REF] and F = ∇R, decomposed in the same frame as u (the tangential-normal frame). This lunar perturbing potential has the form

R(q, q ′ ) = µ M M E 1 |q -q ′ | - q -q ′ (r ′ ) 3 , (2.21)
where q is the satellite position q ′ is the Moon position µ M is the standard gravitational parameter of the Moon -M E is the mass of the Earth r ′ is the distance of the Moon from the Earth center, r ′ = q ′ .

In order to obtain an expression of R in terms of the orbital elements x, x ′ , M, t, (for x = (n, e, ω, i, Ω), the satellite orbital elements, and x ′ = (n ′ , e ′ , ω ′ , i ′ , Ω ′ ), the Moon orbital elements) we make the development given in chapter 4 in the case where the spacecraft is not forced to be in the same plane as the orbit of the Moon around the Earth. This development takes the form (4.5).

If we assume that the Lunar orbit is circular around the Earth, we obtain a significant simplification of (4.5). In the case of a circular orbit, r ′ and n ′ (the Moon mean motion) are fixed, and R(q, q ′ ) is instead R(n, n ′ , e, ω, M, t), where the Moon's rotation around the Earth is periodic with respect to time t with period 2π n ′ . We note that there are two angular motions involved in the Lunar perturbation: the angular variable of the satellite rotation around the Earth, M , and the angular motion of the Moon, given by t. Thus if we wish to simplify the equations of motion (2.19) to remove the angular variables, it becomes necessary to integrate with respect to both these variables, giving rise to a so-called 'double-averaging'

R = n ′ 2π 1 2π 2π n ′ t=0 2π M =0 R(n, n ′ , e, ω, M, t) dM dt. (2.22)
Note that in the paper [START_REF] Bonnard | Lunar and J 2 perturbations of the metric associated to the averaged orbital transfer[END_REF] (in chapter 4), the variable M ′ is used instead of t, and expansions are valid even if the Moon eccentricity ρ ′ is nonzero. Carrying out this computation gives rise to the doubly-averaged Lunar potential,

R = n ′ 4n 4/3 1 + 3 2 e 2 .
(2.23)

Satellite motion under this third-body perturbation is obtained from the Lagrange equations (1.15) with R used in the place of R.

The solution (n(t), e(t), ω(t)) of the Hamiltonian vector field of the free system (deduced from the Lagrange equations applied to the potential (2.23), which is the free-system Hamiltonian) is determined numerically and plotted in figure 4.1-4.2. Both the non-averaged perturbation and the double-averaged perturbation (2.23) (associated to the Hamilton's equations of the free system with and without averaging) are used, to compare the averaged and the non-averaged case. Note that the averaged solution is not fully 'centered' within the averaged solution: this derives from the fact that, when the same initial conditions are used for both the averaged and non-averaged solution (the classical assumption), we solve a problem 'close to' the averaged problem, but which is not actually the averaged problem associated to the non-averaged problem. This dilemma is discussed in detail in the paper [START_REF] Dargent | Initial and final boundaries transformation when solving optimal control problem with averaging techniques and application to low thrust orbital transfer[END_REF], and a way of constructing 'correct' initial conditions such that the averaged problem solved is actually the average of the associated non-averaged problem.

Perturbed energy-minimal transfer

We consider the control problem

ẋ = G(x, M )u + G ′ (x, M )F (x, M, t), x ∈ C, (u 1 , u 2 , ..., u ℓ ) ∈ R ℓ Ṁ = g 0 (x, M ) + 3 k=1 u k g k (x, M ) + 3 k=1 F k (x, M )g k (x, M ) (x(0), M 0 ) = x 0 , x(t f ) = x f J = t f 0 |u| 2 dt → min (2.24)
where C is the elliptic domain (1.5), G, g k and F (x, M, t) are given by (2.20). The Hamiltonian associated to the optimal control problem (2.24) would be given by equation

H 2 k + H k F k (x, M, t).
(2.25)

In the case that F = 0, the Hamiltonian in the coordinates x = (n, e, ω) is given by (2.8), which we will denote by Henergy .

Averaging would give, according to (1.74), a non-homogeneous Hamiltonian

Henergy + H k F k . (2.26)
We prefer the homogeneous

d Henergy + H k F k , (2.27)
where d is a scaling parameter associated to the maximal control magnitude, firstly because it is homogeneous, and secondly, because it derives from a Zermelo problem.

A Zermelo problem essentially involves a perturbation in one direction of a control system. It derives from the classical 'Zermelo's navigation problem', which considers a boat navigating on a body of water affected by a current, originating from a point O and transferred to a destination point D. The boat is capable of a certain maximum speed, and the problem is a time-minimal transfer from the origin O to the point D. This problem extends to n-dimensional Riemannian manifolds; formally, Definition 2.3. A Zermelo navigation problem on a n-dimensional Riemannian manifold (X, g) is a time minimal problem associated to the system

dx dt = F 0 (x) + n i=1 u i F i (x) (2.28)
where F i forms an orthonormal frame for the metric g, and |u| ≤ 1. The field F 0 represents the current of magnitude |F 0 | g .

Applying the maximum principle to the Zermelo problem defines a Hamiltonian which is homogeneous in λ and of the form

H = H 0 + d H 1 (2.29)
where H 0 is linear in λ, and H 1 is quadratic with respect to λ. Conversely, one can associate to a Hamiltonian of this form, like (2.27), a Zermelo navigation problem.

Computations of cut and conjugate loci of the controlled coplanar system are provided (in figures 4.3-4.7) using the software Hampath (developed in [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF]). The solutions are expressed in (n, ψ, ω) coordinates, where e = cos ψ. In these numerical results, various initial points (ψ 0 , ω 0 ) are used and the extremal trajectories are computed for different values of d in equation (2.27). These extremals are shown in figures 4.3-4.7, where the extremal trajectory in the perturbed case is represented in dash-dot line, and the unperturbed case with a solid line. In figure 4.8, the time evolution of the determinant of the matrix (δx 1 (t), δx 2 (t), ẋ(t)) is shown. This is not as conclusive as the precise results in [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF], it is in a sense preliminary.

The Zermelo geometry is more complex than Riemannian, but there are some systematic studies: see for example [START_REF] Serres | On Zermelo-like problems: Gauss-Bonnet inequality and E. Hopf theorem[END_REF], [START_REF] Serres | On the curvature of two-dimensional optimal control systems and Zermelo's navigation problem[END_REF].

Three and four-body transfers

In the case of orbital transfers constrained to the region of space where perturbing forces in the space environment are low (i.e. the chosen value of the thrust ratio in equation (1.17) is satisfied), the perturbed two-body motion (1.2) is a good first approximation. However, if a transfer is such that the external third and fourth-body perturbations makes significant changes over any arc in the orbit, it is necessary to consider (at least for a section of the trajectory) a model such as a full three or four-body problem.

Spheres of influence (SOI) are spheroid-shaped regions around a celestial body where the primary gravitational influence on an orbiting object is that body. In space missions that visit the SOIs of different celestial bodies, it is interesting to use different techniques and coordinate in each of these SOIs, and also in regions outside the SOIs. A transfer from an Earth orbit to a Moon orbit, for instance, has clearly three phases: one where the Earth gravity is preponderant, one where they have the same order of magnitude, and one where the gravity of the Moon is preponderant.

Such transfers have already been realised (e.g. SMART-1 mission (2006), Boeing-702SP (2015) ) and studied, for example in [START_REF] Daoud | Contribution au contrôle optimal du problème circulaire restreint des trois corps[END_REF], where an minimizing Earth-Moon transfer is constructed numerically within the restricted three-body problem in synodic coordinates. The goal here is to try to use averaging in the regions that are close to a celestial body, and join these sections with others where averaging is not suitable.

If we do not use averaging in the SOI of a primary, the existence of two temporal time-scales (the so-called "fast" and 'slow' time discussed in section 1.3.5) requires a very fine discretization grid in the numerical integration scheme.

-The presence of strong oscillations (or weak perturbations by the way of formulating the problem) can also cause the existence of plural local solutions (in the sense of optimization) [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] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation : resolution par homotopie différentielle[END_REF]. The presence of these local solutions, more closely spaced when the oscillations are strong, create a strong numerical instability.

-Conventional methods are designed simply to seek a local solution of the problem, and so you are never sure if the solution obtained corresponds to the overall solution (i.e. the lowest value to minimize the condition). Thus averaging in the numerical solution over at least one arc of the trajectory increases our chances of avoiding these local solutions and also makes for a much faster convergence (since the small oscillations are removed, making the computation more numerically stable).

We use this method of averaging along suitable arcs on a more modest goal than an Earth-to-Moon transfer; in fact, we target the L 1 -Lagrange point from an Earth-centered orbit. The method and numerical experiments are detailed in chapter 5. The underlying concept is that of a multiple shooting method. The idea is to (1.) First solve a minimum-time problem from the Earth orbit with partial final conditions in terms of distance to the Earth and/or mechanical energy that uses averaging.

(2.) Then solve another minimum-time problem from (part of) these initial conditions to the Lagrange point. These two steps yield two trajectories (together with two sets of initial adjoint vectors) that are minimum-time individually (only in the average sense for the first one), but the concatenation of the two has a discontinuity in the state and adjoint space and carries no extremality property.

(3.) We then solve a third combined problem, taking as initial guess for the adjoint vectors the result of the first two problems, that is a combination of the two above whose shooting function is zero if and only if the two pieces of extremal join continuously at the point where we switch from one to the other.

The numerical experiments are conclusive in that the first two steps make the convergence of the third one easier and the comparison with a minimum-time solution without averaging shows very little mismatch (see figure 5.7).

We will use the dynamics of the bicircular problem (section 1.1.8) instead of that of the circular restricted three-body problem (section 1.1.7) in the 'non-averaged' arc (2.). We do this because, although the Sun perturbation is relatively insignificant at L 1 compared to the triangular Lagrange points, however, due to the low thrust of the satellite, this effect may firstly play a role comparable to the engine thrust as the satellite moves out of the Earth SOI. Secondly, when we consider the perturbations of the Kepler orbits (from the approximations given in [START_REF] Roy | Luni-solar perturbations of an Earth satellite[END_REF]),

F Sun ∼ 2 m S m E r r S 3 1 - 3 2 r R S , (2.30)
F Moon ∼ 2 m M m E r r M 3 1 - 3 2 r R M (2.31)
(where for the satellite position q in the plane, r = q 2 1 + q 2 2 , r s = (q S 1 -q 1 ) 2 + (q S 2 -q 2 ) 2 and r M = (q M 1 -q 1 ) 2 + (q M 2 -q 2 ) 2 ) and the plots in figure 2.3 of the the behaviour of these perturbations as r changes, we note that the Lunar and Solar perturbations increase comparably and have the same order of magnitude from approximately 8 Earth radii from the Earth's center until 21 Earth radii, where the two values are equal. Since the L 1 point is more than 50 Earth radii from Earth's center, there are sections of the trajectory where they each dominate the other and so should both be considered to give a totally accurate model in the low-thrust case. Time versus energy in the averaged optimal coplanar Kepler transfer towards circular orbits

This chapter is a reproduction of the paper 'Time versus energy in the averaged optimal coplanar Kepler transfer towards circular orbits' which is authored by Bernard Bonnard, Helen Henninger, Jean-Baptiste Pomet and Jana Nemčova. It appeared in Acta Appl. Math. [START_REF] Bonnard | Time versus energy in the averaged optimal coplanar Kepler transfer towards circular orbits[END_REF] .

Introduction

We consider the controlled Kepler equation describing orbital transfers with low thrust engines, that we normalize as q = -q q 3 + u;

(3.1)
the control is constrained by u ≤ ε, where ε is a small parameter. The phase space, or state space, is the one with coordinates (q, q). Let K = 1 2 q 2 -1/ q be the mechanical energy of the uncontrolled system and X be the elliptic domain:

X = {K < 0, q ∧ q = 0} .
For the free motion (u = 0), the solutions that lie in X are ellipses -or more precisely closed curves that project on the q component as ellipses-and they form a foliation of X.

In this domain, we may chose coordinates (x, l) where x is made of independent first integrals of the uncontrolled motion (so that x describes the geometry of the ellipses) and the "longitude" l defines the position of the spacecraft on this ellipse; (q, q) can be expressed in terms of (x, l) and vice versa. Restricting to the coplanar case, where q and q have dimension 2 and x has dimension 3, the system can be written as

ẋ = i=1,2 u i F i (x, l) , l = Ω(x, l) ,
where the control u = (u 1 , u 2 ) is the coordinates of the original acceleration u in some frame F 1 , F 2 , e.g., the tangential/normal frame (the vector fields F 1 , F 2 are another basis of the distribu-47 tion spanned by ∂/∂ q1 , ∂/∂ q2 in the original cartesian coordinates). In these coordinates, the free motion is ẋ = 0, l = Ω(x, l); there may be a control term in l too but we neglect it for clarity.

The energy minimization problem is the one of minimizing a quadratic criterion u 2 dt for fixed initial and final value of x, and free l; it was analyzed from the averaging point of view in a series of articles [START_REF] Edelbaum | Optimum low-thrust rendezvous and station keeping[END_REF][START_REF] Edelbaum | Optimum power-limited orbit transfer in strong gravity fields[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 | Geodesic flow of the averaged controlled Kepler equation[END_REF]. The Pontryagin maximum principle yields (for any type of cost: energy, final time or others) an Hamiltonian on the cotangent bundle of the state space with the property that a minimizing trajectory must be the projection of an integral curve of the Hamiltonian vector field. For energy minimization, this Hamiltonian is

H(x, p, l) = 1 2 (H 1 (x, p, l) 2 + H 2 (x, p, l) 2
) where H i (x, p, l) = p, F i (x, l) are the Hamiltonian lifts of the vector fields F i and p is the vector of costate variables of the same dimension as the state vector.

As the bound ε tends to zero, the time needed to reach a given orbit tends to infinity. During this very long time, the variable x move slowly because the control is small while variables like l move fast thanks to the term Ω; this yields ill conditioned integration if numeric methods are used. It may be shown that there is an average Hamiltonian

H(x, p) = 1 2π 2π 0 ̟(x, l)H(x, p, l)dl ,
with ̟ some weight function to be determined, that eliminates the fast variable l and whose Hamiltonian flow gives a remarkably good approximation of the movement of x in the original system if ε is indeed small. It sometimes leads to explicit formulas, and is anyway much better conditioned numerically because the fast variable has been eliminated. We shall recall briefly these facts but are more interested in studying qualitatively this new Hamiltonian. We refer the reader to [1, §52] (although no control is considered there) for details on this approximation and its validity. It turns out that it is quadratic definite positive with respect to p and hence derives from a Riemannian metric on X; furthermore, the coefficients of this metric can be explicitly computed. In the coplanar case the geodesic flow is Liouville integrable and the metric associated to a subproblem related to transfer from an arbitrary orbit (in X) to a circular one is even flat: in suitable coordinates the minimizing solutions are straight lines [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF]. Moreover this result is still true if the thrust is oriented only in the tangential direction [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF].

The same averaging technique can be applied in the minimum time case. The non averaged Hamiltonian reads H 2 1 (x, p, l) + H 2 (x, p, l) 2 and again an averaged Hamiltonian may be constructed:

H(x, p) = 1 2π 2π 0 ̟(x, l) H 2 1 (x, p, l) + H 2 (x, p, l) 2 dl .
Like in the energy case, this Hamiltonian derives from a metric on X, i.e. the data of a norm on each tangent space to X; however, unlike in the energy case and as observed in the article [START_REF] Bombrun | The averaged control system of fast oscillating control systems[END_REF], these norms are not associated with inner products on these tangent spaces -this defines a Finsler metric [START_REF] Bao | An introduction to Riemann-Finsler geometry[END_REF], not necessarily Riemannian-and are not everywhere smooth. Technical problems involved in going from Riemannian to non smooth Finsler geometry make the computations of time minimal transfer towards circular orbits a complicated problem.

The objective of this article is to make a preliminary qualitative description of the time minimum transfers and to compare them with the energy minimum ones: section 3.2 recalls the equations and the computation of the average Hamiltonians; section 3.3 recalls the results from [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF][START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF] on the minimum energy problem; section 3.4 provides a new analysis of the minimum time problem, for transfers to circular orbits, and in particular proves that the elliptic domain is geodesically convex in this case; section 3.5 explains why that proof fails in the minimum energy problem, which is consistent with the non-convexity mentioned in [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF].

Preliminaries

Hamiltonian formalism, Pontryagin maximum principle

The goal of this paper is to study some Hamiltonian systems associated to optimal control problems. For the sake of self containedness, let us sketch the relation to the optimal control problems.

Consider the smooth control system ẋ = f (x, u, t) for x ∈ X, an n-dimensional manifold, t ∈ R and u ∈ B ⊂ R ℓ .

An optimal control problem on X associated with the control system ẋ = f (x, u, t) is, for instance, the problem of finding relative to the given points x 0 , x T the trajectory x(•) and control u(•), and possibly the final time T if it is not specified, such that

ẋ = f (x, u, t), x ∈ X, (u 1 , u 2 , ..., u ℓ ) ∈ B ⊂ R ℓ x(0) = x 0 , x(T ) = x T J = T 0 L(x(t), u(t))dt → Min. (3.2)
We call "minimum time" the problem where L(x, u) = 1 and T is free, and "minimum energy" the one where T is fixed and L(x, u) = u 2 . The Hamiltonian of the optimal control problem (3.2) is the function

H(x, p, u, p 0 , t) = p 0 L(x, u) + p, f (x, u, t)
where p is a vector of costate variables (the adjoint vector) of the same dimension as the state variables x(t), and p 0 is either 0 or -1. The Pontryagin maximum principle [START_REF] Pontryagin | The mathematical theory of optimal processes Translated from the Russian by[END_REF] (see also [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF]Chap. 6] for applications to the problems we consider here) is a powerful necessary condition for optimality, that states the following: if (x(•), u(•)) is an optimal trajectory-control pair of the above optimal control problem on a time interval [0, T ], then it can be lifted to a parameterized curve t → (x(t), p(t)) on the cotangent bundle T ⋆ X (p is the adjoint vector, or the vector of costate variables) that satisfies, for almost all time and either for p 0 = 0 or for p 0 = -1,

ẋ(t) = ∂H ∂p (x(t), p(t), u(t), p 0 , t) = f (x(t), u(t), t) ṗ(t) = - ∂H ∂p (x(t), p(t), u(t), p 0 , t) (3.3)
and, for almost all t, H(x(t), p(t), u(t), p 0 , t) is the maximum of H(x(t), p(t), u, p 0 , t) with respect to u ∈ B. The solutions where p 0 = 0 are called abnormal. Let us assume p 0 = -1.

In the problems we consider here, we are in the nice situation where for all (x, p, t), or almost all (x, p, t), there is a unique u ⋆ (x, p, t) such that

H(x, p, t) = H(x, p, u ⋆ (x, p, t), -1, t) = max u∈B H(x, p, u, -1, t)
(the second equality is a property of u ⋆ (x, p, t); the first equality is the definition of H from H and u ⋆ ). In that case, one may sum up the above in the following way: if (x(•), u(•)) is an optimal trajectory, then x(.) may be lifted to a solution of the Hamiltonian vector field associated to H on T ⋆ X:

ẋ = ∂H ∂p (x, p, t), ṗ = - ∂H ∂x (x, p, t). (3.4)
The situation is even nicer if u ⋆ is a smooth function of x, p, t; if not, one must be careful about existence and uniqueness of solutions of solutions to this differential equation. We kept the above time-varying system because we will encounter time-periodic Hamiltonians that we average with respect to time, or with respect to a variable that we may view as a new time.

Coordinates

First of all, we recall the equations describing the planar controlled Kepler problem in the elliptic case (mechanical energy K is negative).

If we chose as coordinates (n, e x , e y , l) where n is the mean movement (n = 1/a 3 = (-2K) 3/2 ; a is the semi-major axis), (e x , e y ) are the coordinates of the eccentricity vector in a fixed frame and l is the "longitude", or the polar angle with respect to a fixed direction, then the elliptic domain is given by {n > 0, e x 2 + e y 2 < 1}. The control system is described by the Gauss equations, where u t , u n are the coordinates of the control in the tangential-normal frame: ṅ = -3n 2/3 1 + 2 (e x cos l + e y sin l) + e x 2 + e y 2

1 -e x 2 -e y 2 u t (3.5a)

ėx = n -1/3 1 -e x 2 -e y 2
1 + 2 (e x cos l + e y sin l) + e x 2 + e y 2

× 2 (cos l + e x ) u t -sin l + 2 e y + 2e x e y cos l -(e x 2 -e y 2 ) sin l

1 -e x 2 -e y 2 u n (3.5b)

ėy = n -1/3 1 -e x 2 -e y 2
1 + 2 (e x cos l + e y sin l) + e x 2 + e y 2

× 2 (sin l + e y ) u t -cos l + 2 e x + (e x 2 -e y 2 ) cos l + 2e x e y sin l

1 -e x 2 -e y 2 u n (3.5c) l = n (1 + e x cos l + e y sin l) 2 (1 -e 2 ) 3/2 . (3.5d)
Instead of e x , e y , it will be more convenient to use the eccentricity e and the argument of the pericenter ω (not defined if e = 0), defined by e x = e cos ω, e y = e sin ω .

(3.6)

The equations become:

ṅ = - 3n 2/3 √ 1 -e 2 √ 1 + 2e cos v + e 2 u t (3.7a) ė = √ 1 -e 2 3 √ n 1 √ 1 + 2e cos v + e 2 2(e + cos v) u t -sin v 1 -e 2 1 + e cos v u n (3.7b) ω = √ 1 -e 2 e 3 √ n 1 √ 1 + 2e cos v + e 2 2 sin v u t + 2e + cos v + e 2 cos v 1 + e cos v u n (3.7c) l = n (1 + e cos v) 2 (1 -e 2 ) 3/2 . (3.7d)
The angle v is the true anomaly

v = l -ω. (3.8)
In these coordinates, the elliptic domain is

X = {x = (n, e, ω) , n > 0, 0 ≤ e < 1, ω ∈ S 1 } . (3.9)
Remark 3.1 (Transfer towards a circular orbit). In the transfer "towards a circular orbit" (or merely if we do not take into account the direction of the semi-major axis during the transfer), we may use these coordinates although they are singular at e = 0, because the variable ω may simply be ignored; this is possible because it is a cyclic variable, i.e. it does not influence the evolution of the other variables (n, e, v). In the variables (n, e), the elliptic domain is:

X = {(n, e), 0 < n < +∞, -1 < e < 1} . (3.10)
The fact that negative values of e are allowed comes from identifying (-e, ω) with (e, ω + π), or, equivalently, considering that (e x , e y ) (see (3.6)) lies on a line of fixed arbitrary direction instead of a half-line. This line may for instance be {e y = 0}, and X is then identified with {(n, e x , e y ), n > 0, -1 < e < 1, e = e x , e y = 0}.

Equations (3.7a)-(3.7d) read: ẋ = 1≤i≤2 u i F i (x, l), l = Ω(x, l) (3.11)
where u 1 , u 2 stand for u n , u t , x = (n, e, ω), the vectors F 1 , F 2 are readily obtained from (3.7a)-(3.7c), and

Ω(x, l) = n (1 + e cos(l -ω)) 2 (1 -e 2 ) 3/2 . (3.12)
One way to introduce averaging is to use the so-called "mean eccentric anomaly". The eccentric anomaly is E, related to e and v by

tan v 2 = 1 + e 1 -e tan E 2 (3.13)
and the mean eccentric anomaly is E -e sin E; the Kepler equation (third Kepler law) implies that, when the control is zero, E -e sin E = n t , t = 0 being the time at the pericenter. Introducing (see for instance [14, sec. 3.6.3])

x 0 = (E -e sin E)/n , one has ẋ0 = 1 if u = 0, i.e. the variable x 0 behaves like time modulo an additive constant; this is an implementation of the flow-box theorem. In the coordinates (x, x 0 ), the system becomes

ẋ = i=1,2 u i F i (x, x 0 ), ẋ0 = 1 + i=1,2 u i G i (x, x 0 ).
Due to the implicit relation between E and x 0 , the practical derivation of such equations is complicated, but they will be useful in formally identifying averaging with respect to l ∈ [0, 2π] and averaging with respect to t ∈ [0, 2π/n].

We define the Hamiltonian lifts (i = 1, 2):

H i (x, p, l) = p, F i (x, p, l) , H i (x, p, x 0 ) = p, F i (x, p, x 0 ) . (3.14)

Averaging

Using the previous equations and rescaling the control with u = εv to introduce the small parameter, the trajectories parameterized by x 0 are solutions of

dx dx 0 = ε i=1,2 v i F i (x, x 0 ) 1 + ε i=1,2 v i G i (x, x 0 )
, which is approximated for small ε by

dx dx 0 = ε i=1,2 v i F i (x, x 0 ).
For this system, we consider the following minimization problems:

• Energy : min v ε 2 x 0 0 i=1,2 v 2 i dt • Time : min v x 0 , v ≤ 1.
Applying the Pontryagin maximum principle leads to the following respective Hamiltonians (normal case in the energy minimization problem),

H e (x, p, x 0 ) = i=1,2 H i (x, p, x 0 ) 2 , H t (x, p, x 0 ) = i=1,2 H i (x, p, x 0 ) 2 , (3.15)
where the lifts H i , defined by (3.14), are periodic with respect to x 0 with period 2π/n.

H i (x, p, l) 2 dl (1 + e cos(l -ω)) 2 (3.18) H t (x, p) = (1 -e 2 ) 3/2 2π 2π 0 i H i (x, p, l) 2 dl (1 + e cos(l -ω)) 2 . ( 3.19)
Remark 3.3. In the original system, the control is "small" (parameter ε). The average system that we study in the next sections can be seen as a limit as ε → 0.

The smaller ε is, the better the average system approximates the real system, but neither the results of this paper not any analysis or simulation in the next sections depend on the size of ε, that is on the magnitude of the thrust.

Singularities.

Let us explain how the non smoothness is a result of the averaging of singularities of a control system. Consider the time minimal control problem for a generic smooth system of the form

ẋ = F 0 (x) + i=1,m u i F i (x), u ≦ 1.
Moreover assume for simplicity that the control distribution D = span{F 1 , . . ., F m } is involutive. From the maximum principle in this case, the extremal control is defined by

u i = H i (x,p) √ i H 2 i (x,p)
where H i (x, p) are the Hamiltonian lifts of F i (x). More complicated extremals are related to the switching surface Σ : H i = 0. Observe that in the single-input case the control is given by u 1 = sign H 1 (x, p) and meeting the surface Σ transversally corresponds to a regular switching. This can be generalized to the multi-input case. More complicated singularities can occur in the non transversal case, for instance in relation with singular trajectories of the system (contained by definition in the surface Σ) [START_REF] Bonnard | Optimal Control with Applications in space and quantum dynamics[END_REF].

The analysis of the averaged systems for minimum energy

First of all we recall the results from the energy case [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF] . The energy minimization problem is expressed as

l f 0 u 2 1 (t) + u 2 2 (t) dt → Min,
where we fix the final cumulated longitude l f (this is slightly different from fixing the transfer time).

The coplanar energy case

In this case the averaged system can be computed explicitly by quadrature, and we have the following proposition.

Proposition 3.4. In the coordinates (n, e, ω) the averaged Hamiltonian (up to a positive scalar) is given by

H e = 1 n 5/3 [18n 2 p 2 n + 5(1 -e 2 )p 2 e + 5 -4e 2 e 2 p 2 ω ] (3.20)
where the singularity e = 0 corresponds to circular orbits. In particular (n, e, ω) are orthogonal coordinates for the Riemannian metric associated to H, namely

g = 1 9n 1/3 dn 2 + 2n 5/3 5(1 -e 2 ) de 2 + 2n 5/3 5 -4e 2 dω 2 .
Further normalizations are necessary to capture the main properties of the averaged orbital transfer.

Proposition 3.5. In the elliptic domain we set r = 2 5 n 5/6 , ϕ = arcsin e and the metric is isometric to

g = dr 2 + r 2 c 2 (dϕ 2 + G(ϕ)dω 2 )
where c = 2/5 and G(ϕ) = 5 sin 2 ϕ 1+4 cos 2 ϕ .

Transfer towards circular orbits

As noticed in Remark 3.1, for such transfers we may ignore the cyclic variable ω and allow negative e. In this case, the elliptic domain is the X given by (3.10). The metric above then reduces to g = dr 2 + r 2 dψ 2 , with ψ = ϕ/c defined on the domain {(r, ψ), 0 < r < +∞, -π 2c < ψ < π 2c }; it is a polar metric isometric to the flat metric dx 2 + dz 2 if we set x = r sin ψ and z = r cos ψ. Flatness in the original coordinates can be checked by computing the Gauss curvature. We deduce the following theorem: with c = 2/5. Since c < 1, the domain is not (geodesically) convex and the metric is not complete.

Remark 3.7 (Tangential thrust). The properties of theorem 3.6 are still true when the thrust is only in the tangential direction except that the metric has a singularity at e = 1. The formula is

g = 1 9n 1/3 dn 2 + (1 + √ 1 -e 2 )n 5/3 4(1 -e 2 ) 1 √
1 -e 2 de 2 + e 2 dω 2 .

We may slightly twist the previous coordinates using e = sin ϕ 1 + cos 2 ϕ to get the normal form

dr 2 + (r 2 /c t )( dϕ 2 + G t (ϕ) dω 2 ), c t = c 2 = 2/5, G t (ϕ) = sin 2 ϕ( 1-(1/2) sin 2 ϕ 1-sin 2 ϕ ) 2 .

3.4

The analysis of the averaged systems for minimum time

The Hamiltonian

We compute H t according to (3.19). The functions H i , i = 1, 2 depend on n, e, ω, p n , p e , p ω , l.

Since we only consider transfer towards a circular orbit, we set p ω = 0 and define h 1 , h 2 by h i (n, e, p n , p e , v) = H i (n, e, ω, p n , p e , 0, ω + v) .

(3.21)

The right-hand side does not depend on the cyclic variable ω, see Remark 3.1. From here on we will use the subscripts 1 and 2 to denote respectively the tangential and normal directions, rather than t and n, for ease of notation. From (3.7), we get

h 1 = n -1/3 -3n p n √ 1 + 2e cos v + e 2 √ 1 -e 2 + 2p e (e + cos v) √ 1 -e 2 √ 1 + 2e cos v + e 2 (3.22a) h 2 = -n -1/3 p e sin v (1 -e 2 ) 3/2 (1 + e cos v) √ 1 + 2e cos v + e 2 (3.22b)
Note that ω does not vary in the integral; the integrand has period 2π with respect to either l or v. This allows us to make the change of variable l = ω + v in the integral in (3.19). In the full control case (both tangential and normal control), the sum in (3. In order to highlight some properties of these Hamiltonians, we perform a canonical change of coordinates (n, e, p n , p e ) → (λ, ϕ, p λ , p ϕ ): n = e 3λ , e = sin ϕ , p n = p λ 3n , p e = p ϕ cos ϕ followed by taking (ρ, ψ) as polar coordinated for the adjoint vector (p λ , p ϕ ); we shall never use again the notations λ, p λ , p ϕ and directly write the change as

3 n p n = ρ cos ψ , √ 1
H t (n, sin ϕ, ρ cos ψ 3n , ρ sin ψ cos ϕ ) = ρ n -1/3 L(ψ, ϕ) (3.26)
H 1 t (n, sin ϕ, ρ cos ψ 3n , ρ sin ψ cos ϕ ) = ρ n -1/3 M (ψ, ϕ) (3.27)
with L and M some functions C → R, where C is the cylinder

C = {(ψ, ϕ), ψ ∈ (R/2πZ), ϕ ∈ R, - π 2 < ϕ < π 2 } = R/2πZ × (- π 2 , π 2
) .

(3.28)

The expressions of L and M are, taking the eccentric anomaly E as the variable of integration instead of v (see (3.13); in particular, dv/(1 + e cos v) = dE/ √ 1 -e 2 ) and restricting the interval of integration from [0, 2π] to [0, π] because the integrand depends on cos E only:

L(ψ, ϕ) = 1 π π 0 I(ψ, ϕ, E) dE , (3.29)
I(ψ, ϕ, E) = α 1,1 (ϕ, cos E) cos 2 ψ + 2 α 1,2 (ϕ, cos E) cos ψ sin ψ + α 2,2 (ϕ, cos E) sin 2 ψ , (3.30
)

α 1,1 = 1 -sin 2 ϕ cos 2 E , α 1,2 = -2 cos ϕ (1-sin ϕ cos E) cos E , α 2,2 = (1-sin ϕ cos E) 1 -3 sin ϕ cos E + 3 cos 2 E -sin ϕ cos 3 E (3.31)
and

M (ψ, ϕ) = 1 π π 0 J(ψ, ϕ, E) dE , (3.32)
J(ψ, ϕ, E) = 1 -sin ϕ cos E 1 + sin ϕ cos E (2 cos ϕ sin ψ -sin ϕ cos ψ) cos E -cos ψ . (3.33)
In the sequel we take advantage of the double homogeneity with respect to ρ and n displayed in (3.26) and (3.27).

Singularities of the Hamiltonian in the single-input and two-input cases

According to (3.26) and (3.27), the Hamiltonians H t and H 1 t , have the same degree of smoothness as, respectively the maps L and M . They are both continuously differentiable on C, but their differentials are not locally Lipschitzcontinuous on the set S; we have the following moduli of continuity of the differentials: in a neighborhood of a point ξ = (ψ, ϕ) ∈ S, in some corrdinates and for a "small" δ,

dL(ξ + δ) -dL(ξ) ≤ k δ ln(1/ δ ) (3.35) dM (ξ + δ) -dM (ξ) ≤ k δ 1/2 (3.36)
Proof. The set S is the set of points (ψ, ϕ) such that I(ψ, ϕ, E) vanishes for some value of E; hence the integrand in (3.29) is real analytic on C \ S and so is L. The degree of regularity (3.35) for L at points in S is given in [START_REF] Bombrun | The averaged control system of fast oscillating control systems[END_REF].

Let us now treat M . It turns out that S is also the border between the region

R 1 = {(ψ, ϕ) ∈ C : -1 + sin ϕ 2 cos ϕ < tan ψ < 1 + sin ϕ 2 cos ϕ } (3.37)
where the sign of J(ψ, ϕ, E) does not depend on E and the region

R 2 = {(ψ, ϕ) ∈ C : tan ψ < -1 + sin ϕ 2 cos ϕ or 1 + sin ϕ 2 cos ϕ < tan ψ} (3.
M (ψ, ϕ) =    -sign cos ψ π π 0 J(ψ, ϕ, E) dE on R 1 , sign P (ψ,ϕ) π arccos R(ψ,ϕ) 0 J(ψ, ϕ, E) dE - π arccos R(ψ,ϕ) J(ψ, ϕ, E) dE on R 2 . (3.40) It is therefore clear that M is real analytic on C \ S = R 1 ∪ R 2 .
The singularity of M on S is not of the type treated in [START_REF] Bombrun | The averaged control system of fast oscillating control systems[END_REF], but it is clear above that the restriction of M to R 1 has a real analytic continuation through S while its restriction to R 2 , on the contrary, behaves like a square root in a neighborhood of S, whence (3.36).

The properties of the differential of the Hamiltonian are important because it is the right-hand side of the Hamiltonian equation. Studying these singularities more precisely is an interesting program that is not yet carried out.

The Hamiltonian flow

Let us now study the solutions of the Hamiltonian equation associated with the minimum time problem in the full control or single control (tangential thrust) cases, namely: ). Specifically, we establish geodesic convexity of the elliptic domain X (see (3.10)), i.e. any two points in X can be joined by a extremal curve. This is contained in the following result: Theorem 3.9 (geodesic convexity). For any (n 0 , e 0 ) and (n 1 , e 1 ) in X , there exist a time T ≥ 0 and a solution [0, T ] → X , t → (n(t), e(t), p n (t), p e (t)) of (3.41) (resp. of (3.42)) such that (n(0), e(0)) = (n 0 , e 0 ) and (n(T ), e(T )) = (n 1 , e 1 ).

ṅ = ∂H t ∂p n , ė =
In order to ease the proof, let us write (3.41) and (3.42) in other coordinates. where a, b, c are given by 1 :

a(ψ, ϕ) = -L(ψ, ϕ) sin ψ -L ϕ (ψ, ϕ) cos ψ , b(ψ, ϕ) = L(ψ, ϕ) sin ψ + L ψ (ψ, ϕ) cos ψ , (3.45)
c(ψ, ϕ) = L(ψ, ϕ) cos ψ -L ψ (ψ, ϕ) sin ψ (resp. given by:

a(ψ, ϕ) = -M (ψ, ϕ) sin ψ -M ϕ (ψ, ϕ) cos ψ , b(ψ, ϕ) = M (ψ, ϕ) sin ψ + M ψ (ψ, ϕ) cos ψ , (3.46) c(ψ, ϕ) = M (ψ, ϕ) cos ψ -M ψ (ψ, ϕ) sin ψ )
and the evolution of ρ is given by:

ρ(τ ) = ρ(0) n(0) n(τ ) -1/3 L(ψ(0), ϕ(0)) L(ψ(τ ), ϕ(τ )) (3.47) (resp. ρ(τ ) = ρ(0) n(0) n(τ ) -1/3 M (ψ(0), ϕ(0)) M (ψ(τ ), ϕ(τ ))
).

(3.48)

The "time" τ is related to the real time t by ∂pe ), we obtain the latter as linear combinations of L(ψ, ϕ), L ϕ (ψ, ϕ), L ψ (ψ, ϕ) (resp. of M (ψ, ϕ), M ϕ (ψ, ϕ), M ψ (ψ, ϕ)) with coefficients depending on n, ϕ, ρ, ψ; substituting these expressions into (3.51) gives

t = n(τ ) 1/3 cos ψ(τ ) L(ϕ(τ ), ψ(τ )) - n(0) 1/3 cos ψ(0) L(ϕ(0), ψ(0)) . (3.49) (resp. t = n(τ ) 1/3 cos ψ(τ ) M (ψ(τ ), ϕ(τ )) - n(0) 1/3 cos ψ(0) M (ψ(0), ϕ(0)) ). ( 3
ψ = n -1/3 (-L sin ψ -L ϕ cos ψ) , φ = n -1/3 (L sin ψ + L ψ cos ψ) , ṅ = -3 n 2/3 (L cos ψ -L ψ sin ψ) (resp. ψ = n -1/3 (-M sin ψ -M ϕ cos ψ) , φ = n -1/3 (M sin ψ + M ψ cos ψ) , ṅ = -3 n 2/3 (M cos ψ -M ψ sin ψ) ).
With the new time τ given by (3. The first two equations in (3.44) form an autonomous system of equations in the two variables (ψ, ϕ) ∈ C that will be the core of our analysis; the third one may be integrated and yields n(τ ):

n(τ ) = n(0) exp -3 τ 0 c(ψ(σ), ϕ(σ))dσ .
(3.52)

The variable ρ (the magnitude of the adjoint vector) plays no role in the evolution of the other variables, in particular the state (n, e) (e = sin ϕ); this is a well-known consequence of the Hamiltonian being homogeneous of degree 1 with respect to the adjoint vector and is anyway obvious from (3.44).

Let us now gather some properties of the maps a, b, c, i.e. of the differential equation (3.44), that are valid both for a, b, c given by (3.45) and for a, b, c given by (3.46); they contain all the information to prove Theorem 3.9.

Proposition 3.11. The maps a, b, c given by (3.45) satisfy the following properties with σ = 0. The maps a, b, c given by (3.46) satisfy the same properties with σ = arctan 1 2 .

1. Symmetries. For all (ψ, ϕ) in C,

a(ψ + π, ϕ) = -a(ψ, ϕ) , b(ψ + π, ϕ) = -b(ψ, ϕ) , c(ψ + π, ϕ) = -c(ψ, ϕ) , a(-ψ, -ϕ) = -a(ψ, ϕ) , b(-ψ, -ϕ) = -b(ψ, ϕ) , c(-ψ, -ϕ) = c(ψ, ϕ) .
(3.53)

2. Uniqueness of solutions. The following differential equation on C:

ψ = a(ψ, ϕ) , φ = b(ψ, ϕ) (3.54)
has, for any (ψ o , ϕ o ) ∈ C, a unique solution t → (ψ(t), ϕ(t)) such that (ψ(0), ϕ(0)) = (ψ o , ϕ o ), defined on a maximum open interval of definition (τ -, τ + ). In this interval, τ -< 0 < τ + where τ -is such that either τ -= -∞ or ϕ(τ -) = ± π 2 , and τ + is such that either τ + = +∞ or ϕ(τ + ) = ± π 2 . This defines a flow Φ from an open subset of C × R to C such that the above unique solution is

t → Φ(ψ o , ϕ o , t) .
(3.55)

Sign and zeroes of b.

There exists a continuous map

Z b : [0, π 2 ] → (- π 2 , 0] (3.56)
continuously differentiable on the open interval (0, π 2 ), such that

Z b (0) = -σ (3.57) and b(ψ, ϕ) = 0, ϕ ≥ 0 ⇔      either ψ = Z b (ϕ), or ψ = π + Z b (ϕ), or ϕ = 0 and ψ ∈ [-σ, σ] ∪ [π -σ, π + σ] . (3.58) Furthermore, b(ψ, ϕ) > 0, ϕ ≥ 0 ⇔ either ϕ > 0 and Z b (ϕ) < ψ < π + Z b (ϕ) , or ϕ = 0 and σ < ψ < π + σ , (3.59)
4. Sign and zeroes of a. One has

0 < ϕ < π 2 ⇒ a(Z b (ϕ), ϕ) > 0 , a(0, 0) = 0 ,
and, if σ > 0 , -σ ≤ ψ < 0 ⇒ a(ψ, 0) > 0 .

(3.60)

5.

Hyperbolic saddle point at (0, 0). The maps a and b are smooth in a neighborhood of (0, 0) and a(0, 0) = b(0, 0) = 0 , ∂a ∂ψ (0, 0) ∂b ∂ϕ (0, 0) -∂a ∂ϕ (0, 0) ∂b ∂ψ (0, 0) < 0 .

(3.61)

6.

Values of c at equilibria. c(0, 0) = 1 , c(π, 0) = -1 .

(3.62)

7.

Stable and unstable manifolds of (0, 0). There exists continuous maps

S : [0, π 2 ] → [-π, 0] , U : [0, π 2 ] → [0, π] , (3.63)
continuously differentiable on the open interval (0, π 2 ), and a number σ with

U (0) = 0, S(0) = -σ, σ ≥ 0 , (3.64)
such that the stable and unstable manifolds of (0, 0) are described by

S 0 = {(S(ϕ), ϕ), 0 ≤ ϕ < π 2 } ∪ [-σ, σ]×{0} ∪ {(-S(-ϕ), ϕ), -π 2 < ϕ ≤ 0} U 0 = {(U (ϕ), ϕ), 0 ≤ ϕ < π 2 } ∪ {(-U (-ϕ), ϕ), -π 2 < ϕ ≤ 0} (3
.65) Furthermore, the zeroes of b are positioned with respect to the stable and unstable manifolds so that the maps S, U, Z b satisfy:

0 < ϕ < π 2 ⇒ S(ϕ) < Z b (ϕ) < 0 < U (ϕ) . (3.66)
Proof. See section 3.7.

The following theorem is almost independent of the rest of the paper: it states that for any a, b, c that satisfy the seven conditions established in Proposition 3.11, the differential equation (3.54) has some properties (that will lead to geodesic convexity); the conditions are of course much more general than the two cases considered in Proposition 3.11. Theorem 3.9 will be easily deduced from Theorem 3.12. Theorem 3.12. If a, b, c satisfy the properties of Proposition 3.11, i.e. (3.53) through (3.62), then, for any ϕ 0 and ϕ 1 in the interval (-π/2, π/2) and any λ ∈ R, there exists τ fin ≥ 0 and a solution (ψ(.), ϕ(.)) : [0, τ fin ] → C of (3.54) such that

ϕ(0) = ϕ 0 , ϕ(τ fin ) = ϕ 1 , τ fin 0 c(ψ(τ ), ϕ(τ ))dτ = λ . (3.67)
Proof of Theorem 3.12. See section 3.8. Proof of Theorem 3.9. Pick n 0 , e 0 , n 1 , e 1 ; according to Proposition 3.11, Theorem 3.12 applies to a, b, c defined either by (3.45) or by (3.46). Take ϕ 0 = arcsin e 0 , ϕ 1 = arcsin e 1 , λ = -1 3 ln n 1 n 0 and apply this theorem. Use (3.52) to get n(τ ) and (3.47) or (3.48) to get ρ(τ ) (with some arbitrary ρ(0), for instance ρ(0) = 1) and finally (3.25) to get e(τ ), p n (τ ), p e (τ ) from ψ(τ ), ϕ(τ ), n(τ ), ρ(τ ). Apply the time reparametrization (τ t) given by (3.49) or (3.50), T being deduced from τ fin in the same way. According to Proposition 3.10, the obtained t → (n(t), e(t), p n (t), p e (t)) satisfies the conclusions of Theorem 3.9. A numerical simulation of the phase portrait of the differential equation (3.54) (or the first two equations in (3.44)) is displayed in Figure 3.1 in the "full control case" where a and b are given by (3.45) and in Figure 3.2 in the "tangential thrust case" where a and b are given by (3.46). This is supposed to be a phase portrait on the cylinder C (for instance, identify {ψ = π 2 } with {ψ = 3π 2 }). The thick trajectories are the stable and unstable manifolds of (0, 0) and (π, 0); the other thick curve is the set of zeroes of b(ψ, ϕ) (i.e. the isocline { φ = 0}). One may check visually the properties established in Proposition 3.11; in particular the unstable manifold of (0, 0) is, in both cases, a graph ϕ → ψ while the stable manifold is also such a graph in the full control case (Figure 3.1) but not in the tangential thrust case (Figure 3.2) where it comprises a segment of the ψ-axis.

Simulations

It can be seen that in both cases, the cylinder C is divided into six regions by these invariant manifolds: one region (called F in section 3.8) where all trajectories go "up" (ϕ is monotone increasing), one (called F + in section 3.8) where all trajectories go "down", and four other regions (called E, E ♯ , E + and E +♯ in section 3.8) where all trajectories cross once the isocline { φ = 0} so that they go up and then down or down and then up. This is exploited in the proof of Theorem 3.12. The generic figure 3.7 is a drawing used to support that proof, that figures in an illustrative manner the features contained in the assumptions of Theorem 3.12, and that can also be observed in the numerical simulations of the two cases that we are really interested in (Theorem 3.9).

Comparison between the minimum-energy and minimumtime cases from the convexity point of view

In section 3.3 we recalled some results from [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF] (and previous work by the same authors); in particular, Theorem 3.6 states that the elliptic domain is not geodesically convex for the energy minimization problem, i.e. some pairs of points in E cannot be joined by a geodesic. In that case, in suitable coordinates ((n 5/6 , 5/2ϕ) as polar coordinates), geodesics are straight lines hence geodesic convexity reduces to usual (affine) convexity, thus the simplest way to see this non convexity is to determine the shape on the elliptic domain in these polar coordinates.

Here we try to explain why convexity holds in the minimum-time case and not in the minimumenergy case. Using the coordinates from Theorem 3.6 for the time-minimizing problem does not seem to shed any light. Rather, we explain how the proof of convexity that we made in the minimum-time case fails when applied to the minimum-energy case.

When p ω = 0, the Hamiltonian in the minimum-energy case is given by (3.20) and can be written as follows is a first integral (it is smooth at ψ = π 2 ). These solutions are drawn on Figure 3.3. It is clear that, on a solution, the maximum possible variation of the variable ϕ is 2/5 π; this implies that, if

H = n -5/3 2(3n p n ) 2 + 5p 2 ϕ ( 3
|ϕ 0 | > ( 2/5 - 1
2 )π, there are some values of ϕ that cannot be reached by any solution starting from the line {ϕ = ϕ 0 }.

Conclusion and open problems

We have studied the average minimum time problem as described in section 3.4.1. This is a reduced subproblem of the planar transfer problem: the state has dimension 2, whereas it would have dimension 3 in the real planar problem (we have set p ω = 0; this imposes that the cyclic variable is constant along transfers) and dimension 5 in the full problem where the plane containing the orbits is not fixed.

In [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF][START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF], the energy problem in full dimension is treated; the planar case is integrable (but only the reduced planar case is flat); the full problem is not integrable but extremals may still be computed explicitly. Studying minimum time in higher dimension is an interesting program.

Concerning the reduced problem considered here, the main contribution of the paper is to prove geodesic convexity of the elliptic domain (any two points in the domain may be joined by an extremal trajectory). On the one hand, it is not clear that this result holds true in higher dimension, and on the other hand, in the present small dimension, optimality and/or uniqueness of the extremal trajectories has not been studied.

Finally the singularities of the Hamiltonian have been investigated roughly, mostly to ensure existence of a Hamiltonian flow. It would be interesting to better understand their nature and their role, in particular the singularities they cause on the balls of small radius for the metric.

Proof of Proposition 3.11

Let us prove that the seven points in Proposition 3.11 are satisfied by a, b, c given by (3.45) (full control case) and also by a, b, c given by (3.46) (tangential thrust case).

1. Symmetries. Equations (3.30) and (3.31) imply

I(π + ψ, ϕ, E) = -I(ψ, ϕ, E), I(-ψ, -ϕ, π -E) = I(ψ, ϕ, E) while (3.33) implies J(ψ + π, ϕ, E) = -J(ψ, ϕ, E), J(-ψ, -ϕ, π -E) = J(ψ, ϕ, E).
Substituting in (3.29) and (3.32) yields, using the change of variable E → π -E in the integrals L(-ψ, -ϕ) and M (-ψ, -ϕ),

L(π + ψ, ϕ) = L(-ψ, -ϕ) = L(ψ, ϕ), M (π + ψ, ϕ) = M (-ψ, -ϕ) = M (ψ, ϕ).
This yields identities (3.53) with a, b, c given either by (3.45) or by (3.46).

Uniqueness of solutions.

This follows from the classical Cauchy-Lipschitz theorem away from S (see Proposition 3.8). On S, -in the full control case (a, b given by (3.45)), as seen in [START_REF] Bombrun | The averaged control system of fast oscillating control systems[END_REF], the regularity properties (3.35) of the right hand side of (3.54) guarantee the existence and uniqueness of solutions to the Cauchy problem (Kamke uniqueness Theorem [40, chap. III, Th. 6.1]), -in the tangential thrust case (a, b given by (3.46)), the same argument does not apply but one may check that the derivative of tan ψ -1 2 (±1 + sin ϕ)/ cos ϕ along ψ = a(ψ, ϕ), φ = b(ψ, ϕ) is nonzero along the curve tan ψ = 1 2 (±1 + sin ϕ)/ cos ϕ, hence the vector field is transverse to S and this implies uniqueness of solutions starting from a point in S (see e.g. [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF]). On the other hand, the derivative of b(ψ, ϕ) with respect to ψ is given by

∂b ∂ψ (ψ, ϕ) = cos ψ π π 0 (1 -sin 2 ϕ cos 2 E) 4 sin 2 E I(ψ, ϕ, E) 3/2 dE . (3.72)
The integrals in both equations are positive. Hence for any fixed ϕ ≥ 0, b(ψ, ϕ) is increasing with respect to ψ on (-π 2 , π 2 ) and decreasing on ( π 2 , 3π 2 ); according to (3.71), it is positive in (0, π 2 ] and negative in (π, 3π 2 ] (identified with (-π, -π 2 ]), hence it must vanish for a unique value of ψ betweenπ 2 and 0, that we call Z b (ψ), thus defining Z b : [0, π 2 ) → (-π 2 , 0]. It also vanishes for a unique value of ψ between π 2 and π that must be equal to π + Z b (ψ) according to (3.53). According to (3.71), Z b (0) = 0 and Z b (0) < 0 if ϕ > 0. Proposition 3.8 says that L, and hence b, are smooth away from S. The part of S that is contained in the square (-π 2 , 0] × [0, π 2 ) is the curve {tan ψ = -1+sin ϕ 2 cos ϕ , 0 ≤ ϕ < π 2 }. We claim that γ(ϕ) = b arctan -1+sin ϕ 2 cos ϕ , ϕ does not vanish between 0 and π 2 ; this is numerically checked by plotting, on Figure 3 The derivative of b(ψ, ϕ) with respect to ψ is zero in R 1 because the above does not depend on ψ, and in R 2 it is given by

b ψ (ψ, ϕ) = sign P (ψ, ϕ) 8 cos 2 ϕ R(ψ, ϕ) (1 -R(ψ, ϕ) sin ϕ) 1 -R(ψ, ϕ) 2 1 -R(ψ, ϕ) 2 sin 2 ϕ (2 sin ϕ cos ψ -cos ϕ sin ψ) 2 . (3.75)
Since |R| < 1 on R 2 , all factors are positive except R(ψ, ϕ), hence b ψ vanishes in R 2 at points where R vanishes, and this is exactly, according to (3.39), on the lines {ψ = ± π 2 }, so that b ψ (ψ, ϕ) has the sign of P (ψ, ϕ)R(ψ, ϕ), i.e. of cos ψ. Hence, for fixed ϕ ≥ 0, b(ψ, ϕ) is -minimum and negative for ψ = -π/2, -increasing for ψ in (-π 2 , arctan -1+sin ϕ

2 cos ϕ ),
-constant, positive if ϕ > 0 and zero if ϕ = 0, for ψ in [arctan -1+sin ϕ 2 cos ϕ , arctan 1+sin ϕ 2 cos ϕ ], -increasing, hence positive, for ψ in (arctan 1+sin ϕ 2 cos ϕ , π

2 ), -maximum and positive for ψ = π 2 . Hence, for any ϕ > 0, there is a unique ψ, -π 2 < ψ < arctan -1+sin ϕ

2 cos ϕ such that b(ψ, ϕ) = 0; we call it Z b (ϕ), thus defining Z b : (0, π 2 ) → (-π 2 , 0), satisfying -π 2 < Z b (ϕ) < arctan -1+sin ϕ 2 cos ϕ satisfying (3.58)-(3.59) (situation on [ π 2 , 3π 2
] by symmetry, see (3.53)). Since (3.74) is valid also on S by continuity, b does not vanish on S except at ϕ = 0, hence b is smooth when it vanishes, away from ϕ = 0; since we also proved that b ψ is nonzero at these points, the inverse function theorem implies that Z b is smooth on the open interval (0, π 2 ); also the monotonicity argument shows that

lim ϕ→0 Z b (ϕ) = -arctan 1 2 , hence Z b is defined [0, π 2 ) → (-π 2 , 0) with Z b (0) = -arctan 1 2 . Since our considerations above for ϕ ≥ 0 imply that b(0, ϕ) is zero if and only if ψ ∈ [-arctan 1 2 , arctan 1 2 ] ∪ [π -arctan 1 2 , π + arctan 1 2 ],
we have proved point 3 of the proposition with σ = arctan 1 2 . 4. Sign and zeroes of a. In (3.60), the part saying that a(ψ, 0) > 0 for -σ ≤ ψ < 0 needs no proof in the full control case because σ = 0 and is easy in the tangential thrust case because, from (3.40) and (3.46), a(ψ, 0) = (1 + cos ψ) sin ψ in R 2 .

We give numerical evidence that a(Z b (ϕ), ϕ) is positive if 0 < ϕ < π 2 . Note that the map Z b can only be determined numerically as the zero ψ = Z b (ϕ) of b(ψ, ϕ) = 0 betweenπ 2 and 0 for fixed ϕ, but the determination is very reliable for b is monotonous with respect to ψ = Z b (ϕ) in the considered region; see point 3. is the curve where R(ψ, ϕ) = -1, the border between R 1 and R 2 ).

5. Hyperbolic saddle. Smoothness around the origin follows from Proposition 3.8. It is clear in both cases that a(0, 0) = 0, b(0, 0) = 0. The computation of the Jacobians is easy (in the tangential thrust case it takes place in the region R 1 with cos ψ ≥ 0, see (3.74)).

In the full control case,

∂a ∂ψ (0, 0) = -2 , ∂a ∂ϕ (0, 0) = 1 2 , ∂b ∂ψ (0, 0) = 1 2 , ∂b ∂ϕ (0, 0) = 1 . (3.76)
The eigenvalues of the Jacobian -2 1/2 In the tangential thrust case,

∂a ∂ψ (0, 0) = -2, ∂a ∂ϕ (0, 0) = 1 2 , ∂b ∂ψ (0, 0) = 0, ∂b ∂ϕ (0, 0) = 1. (3.77)
The eigenvalues of the Jacobian -2 1/2 0 1 at (0, 0) are 1 (unstable) and -2 (stable), associated to the eigenvectors (1/6, 1) and (1, 0), respectively.

6.

Values of c at equilibria. One deduces L(0, 0) = L(π, 0) = 1 from (3.29), (3.30), (3.31) and M (0, 0) = M (π, 0) = 1 from (3.32), (3.33). According to (3.45) and (3.46), this implies c(0, 0) = 1, c(0, π) = -1 in both cases.

7. Stable and unstable manifolds of (0, 0). 7.1 full control case. The unstable manifold is the union of the equilibrium (0, 0) and two solutions that tend to (0, 0) as time τ tends to -∞ and are, according to (3.76), both tangent to the line {ψ = ( √ 10 -3)ϕ} at (0, 0). One of the solutions approaches with positive ψ and ϕ and the other with negative ψ and ϕ. We consider only the first one and call it τ → ( ψ(τ ), φ(τ )), defined on the time interval (-∞, τ + ), τ + ≤ +∞; the result for the other one follows by symmetry. Let D 1 be the rectangle

D 1 = {(ψ, ϕ) ∈ C, 0 < ψ < π, 0 < ϕ} .
On the one hand, we have lim τ →-∞ ψ(τ ) = lim τ →-∞ φ(τ ) = 0, lim τ →-∞ ( ψ(τ )/ φ(τ )) = √ 10 -3, hence ( ψ(τ ), φ(τ )) ∈ D 1 for τ close enough to -∞. On the other hand, the border of D 1 is made of the two equilibria and three segments

{ψ = 0, 0 < ϕ < π 2 } , {ψ = π, 0 < ϕ < π 2 } , {0 < ψ < π, ϕ = 0} .
A short computation shows that

a(0, ϕ) = cos ϕ sin ϕ π π 0 cos 2 E dE 1 -sin 2 ϕ cos 2 E , a(π, ϕ) = -a(0, ϕ) , (3.78)
hence a is positive on the first segment and negative on the second one; according to the proof of Point 3 above, b is positive on the last one; hence solutions starting on these segments all enter D 1 . This proves positive invariance of D 1 (solutions may "exit" through the segment {ϕ = π 2 }, but they are no longer defined). Hence the solution remains in D 1 for all time in the open interval (-∞, τ + ). According to Point 3, b(ψ, ϕ) is positive on D 1 . This solution cannot remain in a compact subset of C for all time because then it would have a non-empty ω-limit set that would have to be a union of equilibria and periodic solutions by Poincaré-Bendixon Theorem, but the fact that φ > 0 in D 1 prevents periodic solutions from existing and the only equilibria are (0, 0) and (π, 0), that cannot be approached because φ(t) cannot become small for φ(τ ) is increasing. Hence necessarily, lim t→τ + φ(t) = π 2 . We have established that the parametrized curve t → ( ψ(t), φ(t)), -∞ < t < τ + defines the graph of a function ψ = U (ϕ), (0, π 2 ) → (0, π); it is continuously differentiable from the implicit function theorem: since the right-hand side of the differential equation is continuous, the parameterized curve is continuously differentiable, and we saw that the derivative of ϕ with respect to the parameter (time) remains positive (again because b > 0 in D 1 ).

Let us turn to the stable manifold. It is the union of the equilibrium (0, 0) and two solutions that tend to (0, 0) as time τ tends to +∞. According to the proof of point 5, both solutions are tangent to the line {ψ = -( √ 10 + 3)ϕ} at the origin. We consider the solution that approaches (0, 0) with negative ψ and positive ϕ, and call it τ → ( ψ(τ ), φ(τ )); the result for the other one follows by symmetry. The proof is now very similar to the one for the unstable manifold, reversing time and replacing D 1 by the domain

D 2 = {(ψ, ϕ) ∈ C , -π < ψ < Z b (ϕ) , ϕ > 0} .
Firstly, b(ψ, ϕ) is negative in this domain. Secondly, the solution is in this domain for τ large enough (obviously ψ is negative and ϕ is positive, and it is on the right side of ψ = Z b (ϕ) because ϕ tends to zero so the solution must spend some time in the region where b < 0). Thirdly, the domain D 2 is negatively invariant: its border is made of the equilibria, the segments {ψ = -π, 0 < ϕ < π 2 }, {-π < ψ < 0, ϕ = 0} and the curve {ψ = Z b (ϕ), 0 < ϕ < π 2 }. Solutions which start on the segments leave D 2 because of (3.78) and the fact that that b(ψ, 0) < 0 if -π < ϕ < 0 (see Point 3 above). Solutions which start on the curve leave the domain D 2 because of Point 4 above (at these points,

Z b is differentiable, Z b ( ψ(τ ), φ(τ )) = 0, a( ψ(τ ), φ(τ )) > 0, b( ψ(τ ), φ(τ )) = 0).
This with the second point implies that the solution is in the domain for all time. The end of the proof, i.e. definition of the continuously differentiable S is exactly the same as the previous proof, only with b < 0 instead of b > 0. Moreover, we get for free that -π < S(ϕ) < Z b (ϕ) from the definition of D 2 ; this and the above implies (3.66).

a (0, ϕ) = cos ϕ sin ϕ π π 0 cos 2 E dE 1 -sin 2 ϕ cos 2 E , b(ψ, 0) =    sign(sin ψ) 2 π 4 -cot 2 ψ if tan ψ ≥ 1 2 (i.e. in R 2 ) , 0 if tan ψ ≤ 1 2 (i.e. in R 1 ) .
(3.80)

The unstable manifold comprises the equilibrium point (0, 0) and two solutions that tend towards it as time τ tends to -∞, tangent, according to (3.77), to the line {ψ = ϕ/6} at (0, 0). Thus either both ϕ and ψ are positive as they approach (0, 0) or they are both negative. Of the two solutions, we will only consider the one where ψ, ϕ are both positive, and call it τ → ( ψ(τ ), φ(τ )). Obviously, ( ψ(τ ), φ(τ )) is in the rectangle

D 1 = {(ψ, ϕ) ∈ C : 0 < ψ < π/2, ϕ > 0}
for τ negative large enough. From Point 3.2 above, b(ψ, ϕ) > 0 in the whole of D 1 . From (3.79) and (3.80), a(0, ϕ) > 0, a(π/2, ϕ) < 0, and b(ψ, 0) ≥ 0 for 0 < ψ < π, ϕ > 0, thus D 1 is positively invariant. The rest of the proof follows exactly the same argument as for the case of the unstable manifold in Point 7.1 above.

We now consider the stable manifold. It comprises the equilibrium point (0, 0) and two solutions that tend towards it as time τ tends to +∞. They are both tangent at (0, 0) to the stable eigenvector, i.e. (see (3.77) and the sequel) to the line {ϕ = 0}. Since φ, i.e. b(ψ, ϕ), is zero on the segment {ϕ = 0, -arctan 1 2 ≤ ψ ≤ arctan 1 2 }, these solutions follow this segment. We examine the one that approaches (0, 0) with negative ψ, the other one follows by symmetry. Call τ → ( ψ(τ ), φ(τ )) the solution such that ψ(0) = -arctan 1 2 , φ(0) = 0. One has φ(τ ) = 0 for all positive τ and ψ(τ ) is increasing for positive τ and tends to zero as τ → +∞. Define the domain

D 2 = {(ψ, ϕ) ∈ C : - π 2 < ψ < Z b (0), ϕ > 0} .
The solution is outside D 2 for positive τ but on its border at τ = 0. From (3.79), a(-π 2 , ϕ) > 0 for all ϕ between 0 and π 2 ; from (3.79), b(ψ, 0)

< 0 for ψ ∈ (-π 2 , -arctan 1 2 ) (-arctan 1 2 is Z b (0))); from point 4 above, a(Z b (ϕ), ϕ) > 0 if ψ ∈ (0, π 2
). Hence D 2 , as well as its topological closure, are negatively invariant. Since ( ψ(0), φ(0)) is on the boundary of D 2 , one has ( ψ(τ ), φ(τ )) ∈ D 2 for all τ ∈ (τ -, 0) where (τ -, +∞) is the maximal interval of definition of the solution we consider. From (3.59) and (3.58), b(ψ, ϕ) < 0 for all (ψ, ϕ) in D 2 . Then, following the same argument as in the proof concerning the stable manifold in Point 7.1 above, we obtain that the restriction to negative times of the solution is the graph ϕ = S(ψ) where S : [0, π 2 ) → (-π 2 , 0] is continuously differentiable on (0, π 2 ) and S(0) = -arctan 1 2 . We already noticed that the other part of the solution covers the segment {ϕ = 0, -arctan 1 2 ≤ ψ < 0}. This ends the proof of point Point ((3.63) to (3.66)) in the tangential thrust case.

3.8 Proof of Theorem 3.12 Theorem 3.12 and this section are independent of the rest of the paper: here we only refer to the seven conditions ranging from equation (3.53) to equation (3.66). Lemma 3.13. Assume that a, b, c satisfy (3.53) (i.e. point 1).

If τ → (ψ(τ ), ϕ(τ )) ∈ C is a solution of (3.54) defined on the time interval [0, τ fin ], then τ → (ψ ♯ (τ ), ϕ ♯ (τ )) and τ → (ψ + (τ ), ϕ + (τ )) with ψ ♯ (τ ) = -ψ(τ ) , ψ + (τ ) = ψ(τ fin -τ ) + π , ϕ ♯ (τ ) = -ϕ(τ ) , ϕ + (τ ) = ϕ(τ fin -τ ) , (3.81)
are also solutions of (3.54) on the same time interval [0, τ fin ] and they satisfy Proof. This is straightforward.

We also use the "♯" and "+" notation to denote the transformations in C:

(ψ, ϕ) ♯ = (-ψ, -ϕ) , (ψ, ϕ) + = (π + ψ, ϕ) . (3.83)
Let us make further constructions and remarks on the conditions (3.53)-(3.66) before proceeding with the proof per se.

Stable and unstable manifolds. Equation (3.61) implies that the Jacobian of the vector field at the equilibrium (0, 0) has two real eigenvalues, of which one is positive and the other negative; i.e. (0, 0) is an hyperbolic saddle (see e.g. [41, section 8.3]). Thus it has a stable manifold S and an unstable manifold U 0 ; these are curves passing through (0, 0) tangent to the corresponding eigenvectors. Their existence is a consequence of (3.61) but point 7 assumes a more specific description.

The number σ. Everything may be stated in a much simpler if σ = 0: in particular in Points 3 and 7, S and Z b may be continued into continuous even maps (-π 2 , π 2 ) → R, and, for instance, the equations of S 0 and U 0 , instead of (3.65), as {ψ = S(ϕ), -π 2 < ϕ < π 2 } and {ψ = U (ϕ), -π 2 < ϕ < π 2 }. We would have preferred this simpler formulation but we do not assume σ = 0 because the proof of Theorem 3.9 in the tangential case uses Theorem 3.12 with a nonzero σ (σ = arctan 1 2 ). However, in order to avoid considering positive and negative ϕ's as different cases in (3.65), (3.58) and (3.59), we define the functions

U 0 : [-π 2 , π 2 ] → [-π, π], S 0 : [-π 2 , π 2 ] \ {0} → [-π, π], Z 0 b : [-π 2 , π 2 ] \ {0} → (-π 2 , π 2
) after S, U , Z b ; these functions are odd and coincide with the former on (0, π 2 ]:

S 0 (-ϕ) = -S 0 (ϕ), U 0 (-ϕ) = -U 0 (ϕ), Z 0 b (-ϕ) = -Z 0 b (ϕ) , (3.84) 0 < ϕ ≤ π 2 ⇒ S 0 (ϕ) = S(ϕ), U 0 (ϕ) = U (ϕ), Z 0 b (ϕ) = Z b (ϕ) . (3.85)
Then, the description (3.65) of the stable and unstable manifolds of (0, 0) may be replaced by

S 0 = {(ψ, ϕ) ∈ C, ϕ = 0 and ψ = S 0 (ϕ)} ∪ [-σ, σ]×{0} , U 0 = {(ψ, ϕ) ∈ C, ψ = U 0 (ϕ)} . ( 3
(ψ, ϕ) = 0 ⇔ ϕ = 0 and ψ ∈ [-σ, σ] ∪ [π -σ, π + σ] or ϕ = 0 and ψ ∈ {Z 0 b (ϕ)} ∪ {Z π b (ϕ)} .
The equilibrium point (π, 0). The "+" symmetry (see (3.81)) obviously maps (0, 0) to (π, 0), the stable manifold of (0, 0) to the unstable manifold of (π, 0), the unstable one to the stable one and the set of zeroes of b to itself. Define S π , U π , Z π b after S 0 , U 0 , Z 0 b by

S π (ϕ) = π + U 0 (ϕ), U π (ϕ) = π + S 0 (ϕ), Z π b (ϕ) = π + Z 0 b (ϕ) . ( 3
0 < ϕ < π 2 ⇒ S 0 (ϕ) < Z 0 b (ϕ) < 0 < U 0 (ϕ) , U π (ϕ) < Z π b (ϕ) < π < S π (ϕ) , -π 2 < ϕ < 0 ⇒ U 0 (ϕ) < 0 < Z 0 b (ϕ) < S 0 (ϕ) , S π (ϕ) < π < Z π b (ϕ) < U π (ϕ) .
(3.88)

The stable and unstable manifolds of (π, 0) are: Invariant regions of C. The stable and unstable manifolds S 0 , U 0 , S π , U π , that intersect at the equilibria (0, 0) and (π, 0) are invariant sets that divide the cylinder C into six open regions:

S π = {(ψ, ϕ) ∈ C, ψ = S π (ϕ)} , U π = {(ψ, ϕ) ∈ C, ϕ = 0 and ψ = U π (ϕ)} ∪ [-σ, σ]×{0} . ( 3
F = {(ψ, ϕ), 0 < ϕ < π/2 and S 0 (ϕ) < ψ < U 0 (ϕ)} , (3.92)
F + = {(ψ, ϕ), 0 < ϕ < π/2 and U π (ϕ) < ψ < S π (ϕ)} , (3.93)
F ♯ = {(ψ, ϕ), -π/2 < ϕ < 0 and U 0 (ϕ) < ψ < S 0 (ϕ)} , (3.94)
F ♯+ = {(ψ, ϕ), -π/2 < ϕ < 0 and S π (ϕ) < ψ < U π (ϕ)} , (3.95) E = {(ψ, ϕ), S 0 (ϕ) < ψ < S π (ϕ) if ϕ ≤ 0 , U 0 (ϕ) < ψ < U π (ϕ) if ϕ ≥ 0 } , (3.96)
E + = {(ψ, ϕ), U π (ϕ) < ψ < U 0 (ϕ) if ϕ ≤ 0 , S π (ϕ) < ψ < S 0 (ϕ) if ϕ ≥ 0 } . (3.97)
E is self-symmetric under the "♯" symmetry and E + is its own image under the "+" symmetry; F + , F ♯ , F ♯+ are the images of F by the "♯" and "+" symmetries; see Lemma 3.13. These regions are represented in Figure 3.7.

We now state and prove two preliminary lemmas and give the proper proof of Theorem 3.12.

Lemma 3.14. Assume that a, b satisfy assumptions (3.53) to (3.60), and consider a solution t → (ψ(t), ϕ(t)) of (3.54) defined on [0, τ fin ], τ fin > 0.

1. If it starts in E or in the upper part of the unstable manifold U 0 , ϕ(.) is monotonic increasing.

2. If it starts in F , then -if it starts in {(ψ, ϕ), Z 0 b (ϕ) ≤ ψ < U 0 (ϕ)}, it remains in this part of E and ϕ(.) is monotonic increasing, -if it starts in {(ψ, ϕ), S 0 (ϕ) < ψ < Z 0 b (ϕ)}, either it remains in this part of E and ϕ(.) is monotonic decreasing, or there is some τ , 0 < τ < τ fin such that t → ϕ(t) is monotonic decreasing for t between 0 and τ , minimum for t = τ and monotonic increasing for t between τ and τ fin , with ψ(τ ) -Z 0 b (ϕ(τ )) = b(ψ(τ ), ϕ(τ )) = 0. 3. If it starts in the upper part of the stable manifold S 0 , ϕ(.) is monotonic non-increasing.

The behavior in the regions E + , F + , F ♯ , F ♯+ and on the other pieces of stable or unstable curves are obtained by symmetry; see Lemma 3.13. For any number f , 0 < f < π 2 , let For any number f , 0 < f < π 2 , there are two neighborhoods Ω 0 and Ω π of (0, 0) and (π, 0) respectively and, a number

C f = {(ψ, ϕ) ∈ C , |ϕ| < f }. ( 3
T (f ) ≥ 0 such that (ψ, ϕ) ∈ Ω 0 ⇒ c(ψ, ϕ) > 1 2 , (ψ, ϕ) ∈ Ω π ⇒ c(ψ, ϕ) < -1 2 , ( 3
meas {t ∈ [0, τ fin ] , (ψ(t), ϕ(t)) / ∈ Ω 0 ∪ Ω π } ≤ T (f ) . (3.101)
The left-hand side (meas stands for the Lebesgue measure of a subset of R) is simply the time spent by the solution outside the neighborhoods Ω 0 and Ω π of (0, 0) and (π, 0); this bound depends on f < π 2 because a and b could tend to zero as ϕ tends to ± π 2 .

Proof of Lemma 3.15. To deal with the case where the stable or unstable manifolds contain a segment of the ψ-axis, we must account for both the cases where σ = 0 and where σ = 0. We define the set K ε (ε > 0) as follows:

K ε = ∅ if σ = 0 (3.102)
and, if σ = 0,

K ε = [ε, σ + ε]×[-kε, kε] ∪ [-σ -ε, -ε]×[-kε, kε] ∪ [π + ε, π + σ + ε]×[-kε, kε] ∪ [π -σ -ε, π -ε]×[-kε, kε] (3.103) with 0 < k ≤ min{ 1, 1 2 ∂a ∂ψ (0, 0) ∂a ∂ϕ (0, 0) }. (3.104)
Since ∂b ∂ψ (0, 0) is zero, (3.61), (3.59) and (3.60) imply ∂b ∂ϕ (0, 0) > 0, ∂a ∂ψ (0, 0) < 0, and (3.66) implies ∂a ∂ϕ (0, 0) ≥ 0. The slope of the tangent to the curve a = 0 is -∂a ∂ψ (0, 0) ∂a ∂ϕ (0, 0), and the slope of the unstable manifold {ψ = U 0 (ϕ)} (i.e. the slope of the eigenvector corresponding to the negative eigenvalue ∂a ∂ψ (0, 0)) is ∂b/∂ϕ(0, 0) -∂a ∂ψ (0, 0) ∂a ∂ϕ (0, 0), larger than the previous slope. Hence, for some open ball B around the origin,

|ϕ| ≤ k|ψ|, ψ = 0, (ψ, ϕ) ∈ B ⇒ |a(ψ, ϕ)| = 0 and |ψ -U 0 (ϕ)| = 0 . (3.105)
This implies that a does not vanish on B ∩ ([ε, σ + ε] × [-kε, kε]). On the compact segment {(ψ, 0), 0 ≤ ψ ≤ σ, (ψ, 0) / ∈ B}, |a(ψ, ϕ)| has a positive lower bound a; hence, for ε small enough, a(ψ, ϕ) is larger that 1 2 a/2 on the part of the compact rectangle [ε, σ + ε] × [-kε, kε] that is outside B. Gluing the piece inside B and the piece outside B together, we get that, for ε small enough, a does not vanish on the compact rectangle [ε, σ + ε] × [-kε, kε]. By symmetry, i.e. from (3.53) and (3.103), we get that, for ε small enough, a does not vanish on K ε .

(3.106) By a similar argument, since ψ -U 0 (ϕ) does not vanish on the compact segment {(ψ, 0), 0 ≤ ψ ≤ σ, (ψ, 0) / ∈ B}, (3.105) implies that it also does not vanish on K ε for ε small enough either, hence

U 0 ∩ K ε = U π ∩ K ε = ∅ (3.107)
for ε small enough. Define the neighborhoods Ω ε 0 and Ω ε π of (0, 0) and (π, 0) as:

Ω ε 0 = {(ψ, ϕ), |ψ| < ε , |ϕ| < ε} , Ω ε π = (Ω ε 0 ) + = {(ψ, ϕ), |ψ -π| < ε , |ϕ| < ε} . (3.108)
Consider the two distinct solutions going through ( π 2 , 0) and (-π 2 , 0); for ε small enough they cross neither Ω ε π nor Ω 0 π , and hence they separate C into two regions, one containing Ω ε π and the other one Ω ε π . Hence, for ε small enough, no solution can cross both Ω ε π and Ω ε 0 . From (3.62), it is also clear that, for ε small enough,

min (ψ,ϕ)∈Ω ε 0 c(ψ, ϕ) > 1 2 and max (ψ,ϕ)∈Ω ε π c(ψ, ϕ) < - 1 2 .
Let us now fix some ε small enough that this is true, (3.106) and (3.107) hold and no solution can cross both Ω ε π and Ω ε 0 . Take

Ω 0 = Ω ε 0 , Ω π = Ω ε π , K = K ε (3.109)
for this fixed value of ε. With this choice, one has

U 0 ∩ K = U π ∩ K = ∅ , (3.110)
and there is some a > 0 such that

(ψ, ϕ) ∈ K ⇒ |a(ψ, ϕ)| > a (
Σ 0 = {Φ(Z 0 b (ϕ), ϕ, t), -π 2 < ϕ < π 2 , ϕ = 0, -1 2 < t < 1 2 } ∩ C f , Σ π = (Σ 0 ) + = {Φ(Z π b (ϕ), ϕ, t), -π 2 < ϕ < π 2 , ϕ = 0, -1 2 < t < 1 2 } ∩ C f .
(3.112)

Note: if, for the initial condition (ψ o , ϕ o ) = (Z π b (ϕ), ϕ), either τ ->-1 2 or τ + < 1 2 (see (3.54)-(3.61)), then Φ(Z π b (ϕ), ϕ, t)
is not defined up to -1 2 or 1 2 ; we however kept, for the sake of simplicity, "-

1 2 < t < 1 2 " instead of "max{-1 2 , τ -} < t < min{ 1 2 , τ + }". The topological closure of C f \ Ω 0 ∪ Ω π ∪ K ∪ Σ 0 ∪ Σ π does
not contain any zero of b, and is obviously compact. Hence b has a positive lower bound m(f ) on that compact set:

(ψ, ϕ) ∈ C f \ Ω 0 ∪ Ω π ∪ K ∪ Σ 0 ∪ Σ π ⇒ |b(ψ, ϕ)| > m(f ) > 0. (3.113)
Now consider a solution [0, τ fin ] → C f , and partition [0, τ fin ] as follows:

[0, τ fin ] = I 0 ∪ I 1 ∪ I 2 ∪ I 3 with I 0 = {t ∈ [0, τ fin ] , (ψ(t), ϕ(t)) ∈ Ω 0 ∪ Ω π } , I 1 = {t ∈ [0, τ fin ] , (ψ(t), ϕ(t)) ∈ K } , I 2 = {t ∈ [0, τ fin ] , (ψ(t), ϕ(t)) ∈ Σ 0 ∪ Σ π and t / ∈ I 0 ∪ I 1 } , I 3 = {t ∈ [0, τ fin ] , (ψ(t), ϕ(t)) / ∈ Ω 0 ∪ Ω π ∪ K ∪ Σ 0 ∪ Σ π } . (3.114) Obviously, meas {t ∈ [0, τ fin ] , (ψ(t), ϕ(t)) / ∈ Ω 0 ∪ Ω π } = meas I 1 + meas I 2 + meas I 3 . (3.115)
Either the solution is one of the two equilibria or it stays in one of the stable or unstable manifolds or in one of the six regions E, E + , F , F + , F ♯ , F ♯+ . According to Lemma 3.13, and seen that the neighborhoods are invariant by the ♯ symmetry and exchanged by the + symmetry, it is enough to prove the property for solutions in the regions E and F , the equilibrium (0, 0) and the upper parts of its stable and unstable manifolds S 0 and U 0 .

In order to bound meas I 1 if σ = 0 (it is zero if σ = 0), let us prove that either the solution does not cross K or there are times t 1 , t 2 , 0 ≤ t 1 < t 2 ≤ τ fin such that the solution is in K on the time interval [t 1 , t 2 ] and outside K on [0,

τ fin ] \ [t 1 , t 2 ]. (3.116)
• This is obvious if the solution is an equilibrium or is on the unstable manifold U 0 , that do not cross K (see (3.110)).

• In E: among the rectangles in (3.103), only

[ε, σ + ε] × [-kε, kε] and [π -σ -ε, π - ε]×[-kε, kε
] intersect E; a solution that lies in E cannot cross both rectangles because the solution passing through (0, π 2 ) separates them. Consider a solution that crosses one of them, say [ε, σ + ε]×[-kε, kε] (the situation around [π -σ -ε, π -ε]×[-kε, kε] is similar). Since ψ is negative in the rectangle (a does not change sign according to (3.106) and a(0, σ) < 0 according to (3.60)) and φ is positive in E, a solution may only exit through the top or lefthand edge; if it exits through the top edge {ϕ = kε}, it will not enter again because ϕ will remain larger than kε. If it exits through the left-hand edge, the fact that φ > 0 in E only allows it to enter again through the same edge, but this is impossible because ψ > 0 on this edge; this proves (3.116).

• In F : from (3.110), [-σ -ε, -ε]×[-kε,
kε] is the only rectangle in (3.103) that intersects F ; hence the solution may only cross this rectangle. Since, according to (3.106), a does not change sign in the rectangle and, according to (3.60), a(0, σ) > 0, a is positive on the rectangle, then the vector field points inwards on the left-hand edge {ψ = -σ -ε} and outwards on the right-hand edge {ψ = -ε}. The bottom edge is not in F . A solution may only exit through the top or right-hand edge; if it exits through the top edge, it means that ϕ(t) is increasing at the exit time, and so, according to Lemma 3.14, it will continue increasing and cannot go back to the rectangle. Also, if it exits through the right-hand edge, re-entering the rectangle through the top or left-hand edge would require ϕ to increase and reach at least kε, making it impossible to reach the rectangle again because ϕ will continue increasing. This proves (3.116) for solutions that remain in F .

• A solution in the upper part of the stable manifold S 0 also satisfies (3.116) because it enters the rectangle through the left-hand edge or the top edge and exits it through the right-hand edge and then goes asymptotically to (0, 0).

We have proved (3.116) for any solution. Either I 1 = ∅ or I 1 = [t 1 , t 2 ] and connectedness implies that the solution stays in a single rectangle of K. Using (3.111), ψ(t) varies monotonically in [t 1 , t 2 ], and so its variation is at most σ: a(t 2 -t 1 ) < σ. This yields

meas I 1 ≤ σ/a . (3.117)
Solutions in E or in the stable and unstable manifolds cross neither Σ 0 nor Σ π ; hence I 2 = ∅ for these solutions. Solutions in F may cross Σ 0 but not Σ π , and they cannot enter Σ 0 again after leaving it because the region between Σ 0 and the unstable manifold U 0 is invariant in positive time. They stay in Σ 0 on a time-interval of length at most 1 from the definition (3.112). Hence, for any solution, meas I 2 ≤ 1 .

(3.118) Using Lemma 3.14, for any solutions in E or F or the upper part of S 0 or U 0 , the total variation of ϕ on the interval [0, τ fin ] is at most π -2ℓ. The inequality (3.113) then implies meas I 3 ≤ (π -2ℓ)/m(f ). Proof of Theorem 3.12. From Lemma 3.13, the "+" symmetry allows one to interchange ϕ 1 and ϕ 0 while the "♯" symmetry changes their sign: we may assume

- π 2 < ϕ 0 ≤ ϕ 1 < π 2 and ϕ 1 ≥ 0 (3.120)
in the proof without loss of generality. We distinguish four cases.

Case a:

-π 2 < ϕ 0 ≤ 0 < ϕ 1 < π 2 .
In this paragraph, by convention, if ϕ 0 = 0, then S 0 (ϕ 0 ) stands for σ and S π (ϕ 0 ) stands for π -σ.

(3.121)

The solutions such that ϕ(0) = ϕ 0 , ϕ(τ fin ) = ϕ 1 , must be in the region E (see Lemma 3.14 and Figure 3.7), and satisfy S 0 (ϕ 0 ) < ψ(0) < S π (ϕ 0 ). For any χ in the open interval (S 0 (ϕ 0 ), S π (ϕ 0 )), let (ψ χ (.), ϕ χ (.)) be the unique solution to the Cauchy problem (3.54) with initial condition

ψ χ (0) = χ , ϕ χ (0) = ϕ 0 . (3.122)
It is -see (3.55)-continuous with respect to χ and continuously differentiable with respect to τ . Since this solution is in E, ϕ χ is an increasing function of time, and so there is a unique time

τ χ fin > 0 such that ϕ χ (τ χ fin ) = ϕ 1 . (3.123) Since b( ϕ χ (τ χ fin ) , ψ χ (τ χ fin ) ) > 0, there is a constant k > 0 such that |ϕ χ (t ′ ) -ϕ χ (t ′′ )| > k|t ′ -t ′′ | for t ′ , t ′′ in a neighborhood of τ χ
fin . This implies that τ χ fin depends continuously on χ. This allows us to define a continuous map Λ : (S 0 (ϕ 0 ), S π (ϕ 0 )) → R by

Λ(χ) = τ χ fin 0 c(ψ χ (s), ϕ χ (s))ds. (3.124)
All we need to prove is that, for any λ ∈ R, there is at least one χ in (S 0 (ϕ 0 ), S π (ϕ 0 )) such that Λ(χ) = λ, i.e. that Λ is onto. Since Λ is continuous, it is sufficient to prove that

lim χ→S 0 (ϕ 0 ) χ>S 0 (ϕ 0 ) Λ(χ) = +∞ , lim χ→S π (ϕ 0 ) χ<S π (ϕ 0 ) Λ(χ) = -∞ . (3.125)
The solution t → (ψ S 0 (ϕ 0 ) (t), ϕ S 0 (ϕ 0 ) (t)) of (3.54) with initial condition (S 0 (ϕ 0 ), ϕ 0 ) is on the stable manifold of (0, 0); it is defined on [0, +∞); ϕ S 0 (ϕ 0 ) (t) is negative for all time and tends to zero as t → +∞. By continuity with respect to initial conditions, the solutions (ψ χ (.), ϕ χ (.)) starting from (χ, ϕ 0 ) with χ close enough to S 0 (ϕ 0 ) are also defined on [0, τ ] for arbitrarily large fixed τ > 0, and converge uniformly to (ψ S 0 (ϕ 0 ) (.), ϕ S 0 (ϕ 0 ) (.)) on the compact interval [0, τ ] as χ → S 0 (ϕ 0 ). This proves that, for χ close enough to S 0 (ϕ 0 ), ϕ χ (τ ) < 0 and hence τ fin > τ . The situation near S π (ϕ 0 ) being similar, we have proved that

lim χ→S 0 (ϕ 0 ) χ>S 0 (ϕ 0 ) τ χ fin = lim χ→S π (ϕ 0 ) χ<S π (ϕ 0 ) τ χ fin = +∞ . (3.126)
Now define Ω 0 and Ω π according to Lemma 3.15, and τ large enough that (ψ S 0 (ϕ 0 ) (τ ), ϕ S 0 (ϕ 0 ) (τ )) is in Ω 0 . For χ close enough to S 0 (ϕ 0 ), we also have (ψ χ (τ ), ϕ χ (τ )) ∈ S 0 (ϕ 0 ); hence, according to (3.100), solutions (ψ χ (.), ϕ χ (.)) with χ close enough to S 0 (ϕ 0 ) never cross Ω π . The interval [0, τ χ fin ] is partitioned into times t such that (ψ χ (t), ϕ χ (t)) ∈ Ω 0 and also (ψ

χ (t), ϕ χ (t)) / ∈ Ω 0 ∪ Ω π . Since, according to (3.101) ( with f = max{|ϕ 0 |, |ϕ 1 |}) meas{t ∈ [0, τ χ fin ] , (ψ χ (t), ϕ χ (t)) / ∈ Ω 0 ∪ Ω π } ≤ T (f ) and meas {t ∈ [0, τ χ fin ] , (ψ χ (t), ϕ χ (t)) ∈ Ω 0 } ≥ τ χ fin -T (f ), then (3.124) and (3.99) imply Λ(χ) ≥ 1 2 (τ χ fin -T (f )) -c f T (f ) with c f = max (ψ,ϕ)∈C, |ϕ|≤f |c(ψ, ϕ)| , (3.127)
and this does imply, using (3.126), the first limit in (3.125). Similarly, for χ close enough to

S π (ϕ 0 ), one has meas{t ∈ [0, τ χ fin ] , (ψ χ (t), ϕ χ (t)) ∈ Ω π } ≥ τ χ fin -T (f ), meas {t ∈ [0, τ χ fin ] , (ψ χ (t), ϕ χ (t)) / ∈ Ω 0 ∪ Ω π } ≤ T (f ) and hence, using (3.127) again, Λ(χ) ≤ -1 2 (τ χ fin -T (f )) + c f T (f ) (3.128)
with c f as in (3.9). This implies, according to (3.126), the second limit in (3.125). Case b: -π 2 < ϕ 0 < 0 and ϕ 1 = 0. If σ = 0, the solutions such that ϕ(0) = ϕ 0 , ϕ(τ fin ) = ϕ 1 must be in the region E, and the proof from case (a) applies, where ϕ 1 is replaced with zero. If σ > 0, the solutions on the stable manifolds {ψ = S 0 (ϕ)} and {ψ = S π (ϕ)} (see Figure 3.7) also qualify for this case, because ϕ reaches zero in finite time. Hence we have to examine the solutions such that (ψ(0), ϕ(0)) = (χ, ϕ 0 ) with χ ∈ [S 0 (ϕ 0 ), S π (ϕ 0 )] instead of the open interval; the solutions (ψ χ (.), ϕ χ (.)) to the Cauchy problem (3.54)-(3.122) are still well defined and depend continuously on χ; however, uniqueness of τ χ fin such that (3.123) holds for χ ∈ (S 0 (ϕ 0 ), S π (ϕ 0 )) but not for χ = S 0 (ϕ 0 ) or χ = S π (ϕ 0 ). If we call τ S 0 (ϕ 0 ) fin (resp. τ S π (ϕ 0 ) fin ) the first time t such that ϕ S 0 (ϕ 0 ) (t) = 0 (resp. ϕ S π (ϕ 0 ) (t) = 0), the solutions to be considered are these with initial condition (ψ(0), ϕ(0)) = (χ, ϕ 0 ), S 0 (ϕ 0 ) < χ < S π (ϕ 0 ) on the time interval [0, τ χ fin ] and these with initial condition (ψ(0), ϕ(0)) = (χ, ϕ 0 ), χ ∈ {S 0 (ϕ 0 ), S π (ϕ 0 )} on the time intervals [0, τ ], τ χ fin ≤ τ < +∞. With the first set of solutions, one reaches λ ∈ [Λ(S 0 (ϕ 0 )), Λ(S π (ϕ 0 ))]; the set of solutions with initial condition (S 0 (ϕ 0 ), ϕ 0 ) allows one to reach λ larger than Λ(S 0 (ϕ 0 )) as the value of τ 0 c(ψ S 0 (ϕ 0 ) (s), ϕ S 0 (ϕ 0 ) (s))ds varies from Λ(S 0 (ϕ 0 )) to +∞ as τ varies from τ S 0 (ϕ 0 ) (the smallest such value such that ϕS 0 (ϕ 0 )(τ ) = 0) to +∞; these solutions with initial condition (S π (ϕ 0 ), ϕ 0 ) allow one to reach λ smaller than Λ(S π (ϕ 0 )).

Case c: ϕ 0 = ϕ 1 = 0. It suffices to chose τ fin = | λ| and the solution to be the equilibrium

(0, 0) if λ ≥ 0 or the equilibrium (π, 0) if λ ≥ 0. Then (3.67) is satisfied because c(0, 0) = 1, c(π, 0) = -1. Case d: 0 < ϕ 0 < ϕ 1 < π 2 .
This is the other "generic" case, with case (a). According to Lemma 3.14, the solutions such that ϕ(0) = ϕ 0 , ϕ(τ fin ) = ϕ 1 must lie in one of the regions F , E, or F + or in the unstable manifolds U 0 or U π that separate them (see Figure 3.7). These solutions satisfy ϕ(0) = ϕ 0 , hence S 0 (ϕ 0 ) < ψ(0) < S π (ϕ 0 ). For any χ in the open interval (S 0 (ϕ 0 ), S π (ϕ 0 )), let (ψ χ (.), ϕ χ (.)) be the solution to the Cauchy problem (3.54)-(3.122). According to Lemma 3.14:

-If S 0 (ϕ 0 ) < χ < Z 0 b (ϕ 0 ) or Z 0 b (ϕ π ) < χ < S π (ϕ 0 ), t → ϕ χ (t) is first decreasing, then crosses the set of zeroes of b at some time t o : ψ χ (t o ) = Z 0 b (ϕ χ (t o )) or ψ χ (t o ) = Z π b (ϕ χ (t o ))
and is increasing for t larger than t o . Hence, since ϕ χ (t o ) < ϕ 0 < ϕ 1 , there is a unique τ χ fin (larger than

t o ) such that ϕ χ (τ χ fin ) = ϕ 1 . -If Z 0 b (ϕ 0 ) ≤ χ ≤ Z π b (ϕ 0 ), t → ϕ χ (t)
is monotonic increasing for positive times and cannot have a limit, hence it takes all the values between ϕ 0 and π 2 only once and there is a unique τ χ fin such that (3.123) holds. In both cases, Z 0 b (ϕ 1 ) < ψ χ (τ χ fin ) < Z π b (ϕ 1 ), hence b( ϕ χ (τ χ fin ) , ψ χ (τ χ fin ) ) > 0, and so there is a constant k > 0 such that |ϕ χ (t ′ ) -ϕ χ (t ′′ )| > k|t ′ -t ′′ | for t ′ , t ′′ in a neighborhood of τ χ fin . This implies that τ χ fin depends continuously on χ. This continuous dependence on χ allows us to define the continuous map Λ : (U 0 (ϕ 0 ), U π (ϕ 0 )) → R by (3.124); as in case (a), let us prove the following limits, sufficient to imply that Λ is onto:

lim χ→U 0 (ϕ 0 ) χ>U 0 (ϕ 0 ) Λ(χ) = +∞ , lim χ→U π (ϕ 0 ) χ<U π (ϕ 0 ) Λ(χ) = -∞ . (3.129)
This follows as in case (a): first we get

lim χ→U 0 (ϕ 0 ) χ>S 0 (ϕ 0 ) τ χ fin = lim χ→U π (ϕ 0 ) χ<S π (ϕ 0 ) τ χ fin = +∞
and then, with f = |ϕ 1 |, (3.127) holds for χ close to U 0 (ϕ 0 ) and (3.128) for χ close to U π (ϕ 0 ). Case e: 0 < ϕ 0 = ϕ 1 < π 2 . This is a similar to the previous case but degenerate in the sense that

τ χ fin = 0 if Z 0 b (ϕ 0 ) ≤ χ ≤ Z π b (ϕ 0 ).
The only nontrivial trajectories that display the same initial and final values of ϕ lie in the regions F or F + , and they join points on one side of the curve where b vanishes (ψ = Z 0 b (ϕ) or ψ = Z π b (ϕ)) to points on the other side. We have

Λ(χ)      > 0 if U 0 (ϕ 1 ) < χ < Z 0 b (ϕ 1 ) , = 0 if Z 0 b (ϕ 1 ) ≤ χ ≤ Z π b (ϕ 1 ) , < 0 if Z π b (ϕ 1 ) < χ < U π (ϕ 1 ) ;
(3.129) still holds and the same arguments prove that Λ is onto.

The Riemannian metric

The controlled Kepler equation, assuming the mass constant can be normalized to

d 2 q dt 2 = - q |q| 3 + u, (4.1)
where q = (q 1 , q 2 , q 3 ) is the position of the satellite and the thrust is bounded by |u| ≤ ǫ. The thrust can be decomposed in a moving frame u = u 1 F 1 + u 2 F 2 + u 3 F 3 e.g. the so-called radialorthoradial frame:

F 1 = q |q| , F 2 = F 3 ∧ F 1 and F 3 = q∧ q |q∧ q| .
The state of the system is described by an angle: the true longitude l and by five equinoctial elements x corresponding to first integrals of the uncontrolled motion. For instance, x = (P, e, h) where P is the semi-latus rectum of the osculating conic, e = (e x , e y ) is the eccentricity vector and h = (h x , h y ) is the inclination vector. We restrict the system to the elliptic domain, that is to the manifold X of elliptic trajectories of the Kepler equation X = {P > 0, |e| < 1}.

The system takes the form

dx dt = 3 i=1 u i F i (x, l) dl dt = w 0 (x, l) + g(x, l, u).
An important problem is to transfer the satellite between coplanar orbits, the corresponding subsystem is deduced by setting both the inclination h and the control u 3 to zero.

The energy minimization problem is studied in detail in [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF] and we present only the main results.

The control is rescaled using u = ǫv, |v| ≤ 1 to introduce the small parameter ǫ and we consider the energy minimization problem to transfer the system from (x 0 , l 0 ) to a terminal orbit x F . The terminal cumulated longitude is also fixed to l F . Parametrizing the trajectory by the cumulated longitude l, the system is written

dx dl = ǫ w 0 (x, l) + ǫg(x, l, u) 3 i=1 v i F i (x, l)
and the cost function to minimize is

ǫ 2 l F l 0 |v| 2 dl w 0 (x, l) + ǫg(x, l, u) .
In order to perform the analytic computation, we first relax the bound |v| ≤ 1. Indeed, for a fixed ǫ, the constraint will be fulfilled by a big enough final longitude l F .

Using the maximum principle [START_REF] Pontryagin | The mathematical theory of optimal processes Translated from the Russian by[END_REF], optimal trajectories are extremals, integral curves of the following Hamiltonian

H ǫ (x, l, p, v) = ǫ w 0 (x, l) + ǫg(x, l, u) (p 0 ǫ|v| 2 + 3 i=1 v i P i )
where p 0 ≤ 0 and P i = p, F i , i = 1, 2, 3. By controllability properties of the system we can restrict to the normal case p 0 < 0 and it can be normalized to -1 2ǫ . As a result, up to first order ǫ, we have the approximation

H ǫ (x, l, p, v) = ǫ w 0 (x, l) - 1 2 |v| 2 + 3 i=1 v i P i + o(ǫ).
In the computation of the averaged system, we can used the first order approximation

H(x, l, p, v) = 1 2 3 i=1 P i √ w 0 2 ,
since the trajectories are C 0 -closed [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF].

Definition 4.1. The averaged Hamiltonian is

H (x, p) = 1 2π 2π 0 H(x, l, p)dl.
Coplanar case We have H = 1 2ω 0 2 i=1 P 2 i and the averaged system is expressed in the coordinates (n, ρ, θ) where n = a -3/2 is the mean movement, a is the semi-major axis, ρ is the eccentricity and θ is the polar angle of the vector (e x , e y ) (hence ρ = e 2

x + e 2 y ),

P = 1 -ρ 2 n 2 3
, e x = ρ cos(θ), e y = ρ sin(θ), and we have Proposition 4.2. In coordinates (n, ρ, θ), the averaged Hamiltonian is

H 1 = 1 4n 5 3 18n 2 p 2 n + 5(1 -ρ 2 )p 2 ρ + (5 -4ρ 2 ) p 2 θ ρ 2
and H 1 is the Hamiltonian of the Riemannian metric

ds 2 = 1 9n 1 3 dn 2 + 2n 5 3 5(1 -ρ 2 ) dρ 2 + 2n 5 3 5 -4ρ 2 ρ 2 dθ 2 .
The coordinates (n, ρ, θ) are orthogonal coordinates.

Non Coplanar case

The complete Hamiltonian is H = 1 2 (P 2 1 + P 2 2 + P 2 3 ). As previously we use (n, ρ, θ) as coordinates and we make a polar representation of h,

h x = σ cos(Ω), h y = σ sin(Ω),
where the angle Ω is the longitude of the ascending node and σ = sin i 2 . Introducing θ = θ -Ω, the angle of the pericenter, and denoting

p θΩ = 2σ 2 σ 2 + 1 p θ + p Ω ,
we have Proposition 4.3. The averaged Hamiltonian of the non-coplanar transfer is

H = H 1 + H 2 with H 2 = (σ 2 + 1) 2 16n 5 3 × 1 + 4ρ 2 1 -ρ 2 (cos(θ)p σ + sin(θ) p θΩ σ ) 2 + (-sin(θ)p σ + cos(θ) p θΩ σ ) 2
and H is associated with a five-dimensional Riemannian metric.

Properties of the metric (coplanar case)

1. The metric associated to

H 1 g = 2 9n 1 3
dn 2 + 4n

5 3 5(1 -ρ 2 ) dρ 2 + 4n 5 3 5 -4ρ 2 ρ 2 dθ 2 ,
is isomorphic to g = dr 2 + r 2 (dφ 2 + G(φ)dθ 2 ) where r = 2 3/2 5 n 5 6 , φ = 1 c arcsin(ρ), G(φ) = 2. The metric g is Liouville integrable with a linear first integral and the geodesic flow can be integrated using elementary functions.

The perturbed case The system is written

dx dl = P (x, l, l ′ ) + 3 i=1 u i F i (x, l),
where P is the perturbation associated to the lunar perturbation, depending on an additional angular variable l ′ , e.g. the lunar longitude or the mean anomaly. The averaging procedure will produce an Hamiltonian which is the superposition of

• An averaged perturbation denoted H P .

• The averaged Hamiltonian H computed before and corresponding to the minimization problem.

This leads to the definition of a Zermelo navigation problem [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF][START_REF] Carathéodory | A Calculus of variations and partial differential equations of the first order[END_REF].

Definition 4.4. A Zermelo navigation problem on a n-dimensional Riemannian manifold (X , g) is a time minimal problem associated to the system

dx dl = F 0 (x) + n i=1 u i F i (x),
where F i form an orthonormal frame for the metric g and |u| ≤

The perturbations 4.3.1 Preliminairies

First of all, the perturbations lead to the definition of a vector field whose trajectories behavior can be roughly classified in the framework of properties of conservative systems in the large, introduced for the three-body problem by Poincaré [START_REF] Poincare | [END_REF] and see [START_REF] Nemytskii | Qualitative theory of differential equations[END_REF] for a modern presentation.

Definition 4.6. Let V be a smooth complete vector field on a manifold M and let x(t, x 0 ) be the solution starting at t = 0 from x 0 . The point x 0 is called Poisson-stable if for every neighbourhood U of x 0 and every T ≥ 0, there exists t 1 , t 2 ≥ T such that x(t 1 , x 0 ) and x(-t 2 , x 0 ) belong to U . The point x 0 is said to be departing if for each compact set K there exists T ≥ 0 such that if |t| ≥ T , x(t, x 0 ) / ∈ K.

Theorem 4.7. Let V be smooth complete conservative vector field on (M, θ), then almost every point is Poisson-stable or departing.

Hence Poisson stability corresponds to bounded motions. A more precise description was recently deduced from KAM theory [START_REF] Celletti | KAM stability for a three-body problem of the solar system[END_REF] which is briefly presented below. Coordinates The satellite position and velocity are represented by (q, q) and we denote (q, p) the standard symplectic coordinates. The motion of the satellite which is defined up to a proper normalization by Kepler Hamiltonian

H(q, p) = 1 2 |p| 2 - 1 |q| .
To analyze the effect of a perturbing force deriving from a potential R, one uses the Lagrange equations [START_REF] Pascoli | Eléments de mécanique céleste[END_REF][START_REF] Zarrouati | Trajectories spatiales[END_REF].

da dt = 1 n 2 a ∂R ∂τ dρ dt = 1 -ρ 2 n 2 a 2 ∂R ∂τ - 1 -ρ 2 na 2 ρ ∂R ∂θ di dt = cot(i) na 2 1 -ρ 2 ∂R ∂θ - 1 na 2 sin(i) 1 -ρ 2 ∂R ∂Ω dΩ dt = 1 na 2 sin(i) 1 -ρ 2 ∂R ∂i dθ dt = 1 -ρ 2 na 2 ρ ∂R ∂ρ - cot(i) na 2 1 -ρ 2 ∂R ∂i dτ dt = 2 n 2 a ∂R ∂a + 1 -ρ 2 n 2 a 2 ρ ∂R ∂ρ (4.2)
where i is the angle of inclination between the orbital plane of the satellite and the orbital plane of the Moon, Ω is the longitude of the ascending node, θ is the angle of the perigee and τ is the time of the perigee passage.

Remark It is useful to introduce the mean anomaly M to locate the satellite on its orbit. It is defined by the relation M = n(t -τ ). In this case, we set R(a, e, i, Ω, ω, τ ) = R(a, e, i, Ω, ω, M ) and the partial derivatives verify

∂R ∂τ = -n ∂ R ∂M , ∂R ∂a = ∂ R ∂a + na 2 dn dt (t -τ ).
The effect of the lunar perturbation on the satellite motions are well understood and we use the computations excerpted from [START_REF] Pascoli | Eléments de mécanique céleste[END_REF]. They are related to solar perturbation of the Moon. In this reference, they study Moon motion under the Sun perturbation, which we can adapt to the Earth-Moon-satellite case.

The lunar perturbation

Given a geocentric inertial frame of reference, the perturbing lunar potential can be expressed by

R(q, q ′ ) = μ′ 1 |q -q ′ | - q • q ′ r ′3 , (4.3)
where μ′ = µ ′ /m Earth is the standard gravitational parameter of the Moon divide by the mass of the Earth, q (resp. q ′ ) is the position vector of the satellite (resp. the Moon) and |q| (resp. |q ′ |) is denoted by r (resp. r ′ ). The potential (4.3) stands for the dynamics of the satellite of the two body problem Earth-satellite which is perturbed by the Moon.

Setting the reference plane as the orbital plane of the Moon, q ′ can be decomposed in terms of the osculating elements (x ′ = (n ′ , ρ ′ , θ ′ ), M ′ ) of the Moon where n ′ = μ′ /a ′3/2 is the mean movement, a ′ the semi-major axis, ρ ′ is the eccentricity, θ ′ the angle of the perigee and M ′ the mean anomaly.

In order to have a rough evaluation of the perturbation, we use a simplified model in [START_REF] Pascoli | Eléments de mécanique céleste[END_REF] based on the following assumptions : the eccentricity of the satellite ρ is small and the inclination i of the satellite with respect to the Moon orbital plane is small. We have,

1 |q -q ′ | = 1 r ′ 1 + r r ′ 2 -2 r r ′ cos(Ψ)
where Ψ is the angle between the two vectors q and q ′ . Assume the satellite on a low Earth orbit, then r ≪ r ′ and using Legendre polynomials P k

1 |q -q ′ | = 1 a ′ ∞ k=0 α k r a k a ′ r ′ k+1 P k (cos(Ψ))
where α = a a ′ ≪ 1. Using the approximation

1 |q-q ′ | = r r ′2 cos(Ψ) + 1 a ′ 1 + 1 2 α 2 r a 2 a ′ r ′ 3 (-1 + 3 cos 2 (Ψ)) + o(α 3 ), (4.4)
the perturbing potential expression becomes

R = μ′ 2a ′ α 2 r a 2 a ′ r ′ 3 (-1 + 3 cos 2 (Ψ)) + o(α 3 ) = n ′2 2n 4/3 r a 2 a ′ r ′ 3 (-1 + 3 cos 2 (Ψ)) + o(α 3 )
where the first term 1 a ′ in (4.4) has been removed since it does not depend on the satellite orbital elements (x = (n, e, i, θ, Ω), M ). Therefore,

Development of the terms

E = M + ∞ k=1 2 k J k (kρ) sin(kM ).
With a similar method (see for instance [START_REF] Vinti | Orbital and celestial mechanics[END_REF]), one obtains that

cos(E) = - e 2 + ∞ k=1 2 k 2 d dρ [J k (kρ)] cos(kM ).
From the relations,

a r = 1 n Ė r a = 1 -ρ cos(E),
we have the following expansions,

a ′ r ′ = 1 + 2 ∞ k=1 J k (kρ) cos(kM ), r a = 1 + ρ 2 2 - ∞ k=1 2ρ k 2 d dρ [J k (kρ)] cos(kM ).
In the sequel we use Maple software in the computations.

Development of the term cos 2 (Ψ)

The spherical trigonometry allows us to express the angle Ψ in terms of the orbital elements of the satellite and the Moon

cos(Ψ) = cos(θ + v) cos(θ ′ + v ′ -Ω) + cos(i) sin(θ + v) sin(θ ′ + v ′ -Ω),
where v and v ′ are respectively the true anomaly of the satellite and the Moon. We have

cos(v) = cos(E)-ρ 1-ρ cos(E) , sin(v) = √ 1-ρ 2 sin(E) 1-ρ cos(E) . Hence R(x, M, x ′ , M ′ ) = k,m∈Z C k,m (x, x ′ ) exp(I(kM + mM ′ )),
and more precisely

R(x, M, x ′ , M ′ ) = k∈Z 5
Ck (a, ρ, i, a ′ , ρ ′ ) exp(I k, ζ ).

where ζ = (Ω, θ, M, θ ′ , M ′ ).

Computations lead to

R(x, M, x ′ , M ′ ) = n ′2 n 4/3 1/4 1 + 3/2 ρ 2 + 3/2 ρ ′ 2 -6 sin 2 (i/2) -1/2 ρ cos (M ) + 3/4 ρ ′ cos M ′ -1/8 ρ 2 cos (2 M ) + 9/8 ρ ′ 2 cos 2 M ′ -3/4 ρρ ′ cos M -M ′ -3/4 ρρ ′ cos M + M ′ + 3/2 sin 2 (i/2) cos (2 θ + 2 M ) + 3/2 sin 2 (i/2) cos 2 θ ′ + 2 M ′ -2 Ω + 15 8 ρ 2 cos -2 Ω -2 θ + 2 θ ′ + 2 M ′ - 63 8 ρρ ′ cos 2 θ + 2 Ω + M -2 θ ′ -3 M ′ -9/4 ρ cos 2 θ + 2 Ω + M -2 θ ′ -2 M ′ + 9/8 ρρ ′ cos 2 θ + 2 Ω + M -2 θ ′ -M ′ + 3 4 - 3 2 sin 2 (i/2) - 15 8 ρ ′ 2 - 15 8 ρ 2 cos 2 θ + 2 Ω + 2 M -2 θ ′ -2 M ′ + 21 8 ρρ ′ cos 2 θ + 2 Ω + 3 M -2 θ ′ -3 M ′ + 3/4 ρ cos 2 θ + 2 Ω + 3 M -2 θ ′ -2 M ′ -3/8 ρρ ′ cos 2 θ + 2 Ω + 3 M -2 θ ′ -M ′ + 3/4 ρ 2 cos 2 θ + 2 Ω + 4 M -2 θ ′ -2 M ′ + 51 8 ρ ′ 2 cos -2 θ ′ -4 M ′ + 2 θ + 2 Ω + 2 M + 21 8 ρ ′ cos -2 θ ′ -3 M ′ + 2 θ + 2 Ω + 2 M -3/8 ρ ′ cos -2 θ ′ -M ′ + 2 θ + 2 Ω + 2 M + o(α) 3 + o(ρ) 3 + o(ρ ′ ) 3 . (4.5)
Definition 4.10. The double averaged potential is defined by

R M,M ′ (x, x ′ ) = 1 (2π) 2 2π 0 2π 0 R(x, M, x ′ , M ′ )dM dM ′ .
Then, the double averaged of the potential (4.5) is

R M,M ′ (x, x ′ ) = n ′2 4n 4/3 1 + 3/2 ρ 2 + 3/2 ρ ′ 2 -6 sin 2 (i/2)
(see [START_REF] Domingos | Third-body perturbation in the case of elliptic orbits for the disturbing body[END_REF] for more precise expansions). A simplified academic model is to set the eccentricity ρ ′ to zero, the inclination i to zero and to restrict the control to this plane, the lunar perturbing potential becomes

R M,M ′ (x, x ′ ) = n ′2 4n 4/3 (1 + 3 2 ρ 2 ). (4.6)

Computations

Shooting equation

Let H be an Hamiltonian vector field associated to the Zermelo navigation problem, z = (x, p), x ∈ X , denoting the state and adjoint vector and H being homogeneous of degree 1 in p. Fixing the initial and final state vectors (x 0 , x 1 ) and t F being the transfer time the shooting equation is defined by

S : p 0 → Π(exp(t F H(z 0 )) = x 1 ,
where z 0 = (x 0 , p 0 ), Π : (x, p) → x and p 0 can be normalized by homogeneity. is called the Jacobi equation. A Jacobi field is a non trivial solution δz = (δx, δp) of Jacobi equation and it is said to be vertical at time t if δx(t) = 0. Definition 4.12. We define the exponential mapping exp x 0 ,t (p 0 ) = Π(z(t, x 0 , p 0 ))

where p 0 can be restricted to the sphere |p 0 | = 1. If z = (x, p) is the reference extremal, a time t c > 0 is said to be conjugate to 0 if the mapping p 0 → exp x 0 ,t (p 0 ) is not of rank n -1 at t = t c (with n = dimX ) and the associated point x(t c ) is said to be conjugate to x 0 . We denote by t 1c the first conjugate time and C(x 0 ) is the conjugate locus formed by the set of first conjugate points.

Testing conjugacy An algorithm can be deduced which is implemented in the Hampath Code [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF] used in our numerical simulations. Let z(t) = (x(t), p(t)) be the reference extremal and consider the vector space of dimension n -1 generated by the Jacobi fields δz i = (δx i , δp i ), i = 1, ..., n -1 vertical at t = 0 and such that δp i (0) is orthogonal to p 0 . At a conjugate time t c , one has rank[δx 1 (t c ), ..., δx n-1 (t c )] < n -1 or equivalently, det[δx 1 (t c ), ..., δx n-1 (t c ), ẋ(t c )] = 0.

Hampath Code This code is used to . Integrate the Hamiltonian flow and compute the Jacobi fields along a given solution.

. Solve the shooting equation.

Computations

We are in the coplanar case and Lagrange equations (4.2) give using (4.6)

dn dt = dρ dt = 0 dθ dt = n 1/3 1 -ρ 2 ρ 3n ′2 ρ 4n 4/3 = 3n ′2 1 -ρ 2 4n .
This leads to the averaged Hamiltonian

H (n, ρ, θ) = p θ 3n ′2 √ 1-ρ 2 4n + λ 1 4n 5/3 18n 2 p 2 n + 5(1 -ρ 2 )p 2 ρ + (5 -4ρ 2 ) p 2 θ ρ 2 (4.7)
where λ is the scaling parameter of the control maximal magnitude.

Numerical results

This section is achieved by a series of numerical computations on the free system, on extremal trajectories and on conjugate loci. The simulations are computed thanks to the Hampath code [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF].

Free system The control is set to zero and the dynamical system is deduced from the Lagrange equations and the perturbative potential (4.6). The solutions are expressed in the (n(t), ψ(t), θ(t)) coordinates where φ is the angle such that ρ = sin(φ), and ψ(t) = π 2 -φ(t). 2 ) ≈ 2.8 • during one lunar revolution around the Earth.

Controlled system Due to the homogeneity of order 1 of the Hamiltonian (4.7), the time-minimum problem is considered and the adjoint vector is normalized. Extremal trajectories are computed for different values of λ. The shooting algorithm is performed to solve the boundary value problem by determining the initial adjoint vector p(0) and the optimal time t f . In the following figures, the perturbed case (λ = 10 -1 ) is represented in dash-dot line and is compared to the unperturbed case represented in solid line. The final points are indicated by cross markers. The first conjugate points, indicated by star markers, are computed thanks to the algorithm presented in the subsection 4.4.2 for which the time evolution of the determinant of the matrix (δx 1 (t) δx 2 (t) ẋ(t)) is presented in Figure 4.8. Figures 4.9 and 4.10 represent the projection of the extremal trajectories in (ψ, θ) coordinates starting from the initial point (ρ 0 , θ 0 ) = (0.60, π) in the unperturbed and perturbed case.

Conclusion

More general perturbations can be considered such as the J 2 -effect.

A brief description of the J 2 -effect [51]

The Earth is modelled by an homogeneous oblate ellipsoid of revolution whose axis of symmetry is identified to the axis of rotation passing through the pole denoted Oz and the position of the satellite can be represented in spherical coordinates (r, λ, φ), λ being the latitude and φ the longitude.

The perturbing potential in the normalized coordinates takes the form

R 2 = 1 2r R e r 2 J 2 (1 -3 sin 2 (φ))
where we have the relation sin(φ) = sin(i) sin(θ + v)

where v is the true anomaly. Hence the perturbing potential is given by

R 2 = 3 2 R 2 e J 2 a 3 a r 3 1 3 - 1 2 sin 2 (i) + 1 2 sin 2 (i) cos(2(θ + v)) , (4.8)
where R e is the mean Earth's equatorial radius and J 2 = 1.08263.10 -3 is a constant. The averaged perturbation computed with the formula

R 2 M = 1 2π 2π 0 RdM,
where M is the mean anomaly gives the following.

Proposition 4.13. The averaged perturbation associated to the J 2 -effect is described by the potential

R 2 M = 3 2 R 2 e J 2 a 3 (1 -ρ 2 ) 3 2 1 3 - 1 2 sin 2 (i) .

Additional perturbations [62]

In practise one may encounter other perturbations : solar perturbations and non conservative type of perturbations such as atmosphere drag and solar eclipses.

Extensions

Similar computations as in [START_REF] Domingos | Third-body perturbation in the case of elliptic orbits for the disturbing body[END_REF] will lead to more accurate model useful with longer time transfers. Finally the lunar perturbation and the J 2 -effect can be superposed and their effects numerically analysed for space mission where both effects have to be taken into account.

Chapter 5

Piecewise transfer from the Earth to L 1

In this chapter, we describe a method to compute a minimum-time low-thrust satellite transfer and implement it numerically. This method separates the transfer into an averaged and a non-averaged arc. Averaging simplifies numerical problems near the Earth, introduced by the low-thrust property of electrical engines (as discussed in sections 1.3.5, 2.3). The non-averaged arc is used to take into account the fact that averaging is no longer valid at a certain distance from the Earth (defined by solving for the radius in equation (1.18) for the chosen thrust magnitude u and thrust ratio ǫ ratio ). The approach in the numerical simulation is to construct the transfer by making use of overlapping optimization subproblems.

Transfer scheme

We wish to solve numerically a satellite transfer from an Earth-orbit to the L 1 Lagrange point (defined in 1.1.7), where the satellite is influenced by Earth, Moon and Sun gravity. In the planar inertial frame, q ∈ R 2 , the satellite dynamics take the form q = -µ E q -q E (q -q E ) 3 -µ M q -q m (q -q M ) 3 -µ S q -q S (q -q S ) 3 + u (5.1)

where q E , q M , q S are the positions of the Earth, Moon and Sun, respectively (when all these planets, and the satellite, are assumed to be co-planar), and µ E , µ M , µ S are the standard gravitational parameters of these planets, in the low-thrust case u ≤ ǫ.

Spheres of influence (SOI) are a concept of astrodynamics referring to the spheroid-shaped regions around a celestial body where the primary gravitational influence on an orbiting object is that body ( [START_REF] Carlson | An analytical solution to patched conic trajectories satisfying initial and final boundary conditions[END_REF], [START_REF] Battin | An introduction to the mathematics and methods of astrodynamics[END_REF], [START_REF] Szebehely | Theory of orbits in the restricted problem of three bodies[END_REF] ). The main benefit of the SOI concept is that it enables the partitioning of the transfer into a series of two-body orbits, which greatly simplifies mission analysis involving multiple planets: a transfer is broken down into portions where the arc takes the form (in the absence of control) of some conic, and one gravitational force is assumed to be acting upon the spacecraft at a time. We adapt this concept to the context of low-thrust problems, by defining the Earth SOI as the region of space where the low-thrust property (described using the 'thrust ratio' 99

ǫ ratio in equation (1.18)) is satisfied: u r 2 m µ E ≤ ǫ ratio . (5.2)
This means that the SOI is bounded by the ball of radius R f such that (for a chosen value of ǫ ratio and u , and the value µ E ) (5.4) From initial orbit, shoot to the chosen radius R f in minimum time with averaging in twobody motion perturbed by Moon and Sun gravity. The mechanical energy of the satellite at the final time t 1 f of this transfer (denoted by E f ) will be obtained numerically (we expect a low negative energy) 2. From L 1 , shoot 'backwards' in minimum time to the mechanical energy E f without averaging in (four-body) bicircular motion (by 'backwards', we mean that the final transfer time t f 2 of the transfer is going to be negative) .

R f = ǫ ratio µ E m u . ( 5
We choose to target the value E f in the second leg of the transfer rather than the radius R f or Cartesian position of the satellite because if the final position of the trajectory in subproblem 2 is reached with a positive energy, it may be impossible to find an initialization under which the adjoint vectors of subproblems 1 and 2 meet continuously for any time value. Targeting both the energy E f and radius R f in the second leg was more costly computationally than to target either one individually, and so we maintained the more significant parameter for the joining of the two trajectories, which is the energy. The dynamics (5.1) can be written in two different sets of coordinates:

Earth-centered Kepler orbital elements: we use (p, e x , e y , ℓ) (section 1.1.3), and, as in (1.16), we use the notation x = (p, e x , e y ) and the dynamics are given by

ẋ = u 1 G 1 (x, ℓ) + u 2 G 2 (x, ℓ) l = Q(x, ℓ) + u 1 g 1 (x, ℓ) + u 2 g 2 (x, ℓ). (5.5)
Synodic coordinates: these were denoted (see section 1.1.7) by (X, Y, Ẋ, Ẏ ) in (1.36)-(1.37) and this notation was mantained through 1.1.7, 1.1.8. Here, we write rather (X 1 , X 2 , X 3 , X 4 ):

X = (X 1 , X 2 , X 3 , X 4 ) stands for the notation (X, Y, Ẋ, Ẏ ) of sections 1.1.7, 1.1.8. (5.6)
The dynamics are given by (1.36), (1.37) under this change of notations, that reads

Ẋ = P 0 (X) + ǫ 2 i=1 P i (X)u i , (5.7)
with straightforward definitions of the vectors P 0 (X), P 1 (X), P 2 (X) from (1.36)-(1.37), detailed in (5.30).

Additional notation: we also need as an intermediatary the Cartesian coordinates in the Earthcentered frame introduced in section 1.1.7, see (1.25). We also change this notation:

ξ = (ξ 1 , ξ 2 , ξ 3 , ξ 4 ) stands for the notation ( X, Y , ˙ X, ˙ Y ) of sections 1.1.7. (5.8)
Of course, both sets of coordinates are valid everywhere, but we privilege one or the other depending on the region: when we are near the Earth, we prefer the Earth-centered equinoctical coordinates that clearly make the Moon and Sun appear as perturbations and will allow us to perform averaging with respect to the fast variable. When we are outside this zone, and in fact near L 1 , we prefer the second one that reflects the comparable contribution of the Earth and Moon, that almost cancel while the control thrust becomes preponderant in the vicinity of L 1 .

Remark 5.1. If we would consider (but it is not the case) trajectories where the spacecraft comes closer to the Moon, we would use Moon-centered orbital elements in that region.

To construct the time-minimal Earth-to-L 1 transfer, the method is to use the initialization (λ 0 , Λ 0 , t 1 f , t 2 f ) which arises from running the subproblems 1.) and 2.) separately (in forwards and backwards time, respectively), to obtain a trajectory of the transfer problem

Transfer           
( ẋ, l, Ẋ) : (5.5), (5.7) Initial condition : (x 0 , ℓ 0 ) = chosen initial orbit; (X 0 1 , X 0 2 , X 0 3 , X 0 4 ) = L 1 Final condition : (x(t 1 ), ℓ(t 1 )) = R f , T (x(t 1 ), l(t 1 ), λ(t 1 ), λℓ (t 1 )) = (X(t 2 ), Λ(t 2 )) t 1 -t 2 = t f minimize : t f (5.9) where the function T : (x, l, λ, λℓ ) → (X, Λ) is a diffeomorphism from 'averaged' (state and costate) Kepler coordinates (as we define in (5.48)-(5.49)) to non-averaged synodic coordinates in the state and costate.

In constructing this transfer numerically, we make use of a simplified version of the software T _3D developed by Thierry Dargent (which we adapt by including our own Lunar and Solar perturbations) to encode the dynamics of subproblem 1 numerically, while we constructed our own scripts for subproblem 2. The scripts are all prepared using matlab.

We now describe in detail the dynamics of subproblems 1 and 2, and the numerical methods used to solve them. We then show the details of our approach to solving (5.9) numerically and how the transformation function T is constructed.

Perturbed two-body problem (subproblem 1)

We stated in the previous section that the satellite motion in subproblem 1 is perturbed by the Sun and Moon perturbations, which gives rise to the dynamics (1.13)

ẋ = G(x, ℓ, t)(u + F ), l = Q(x, ℓ), (5.10)
where we express the third-body perturbation as a vector F in the radial-orthoradial frame, and apply this force vector to the Gauss equations (1.12) to describe the perturbation on each element. This approach was also used in [START_REF] Ronwell | Satellite lifetimes in nearly circular orbits for various Earth-Atmosphere models[END_REF]. Let us describe the computation of the perturbation due to the Moon and Sun. Since they are similar, we do only one, with a 'third body' that can be either the Moon or the Sun. The configuration of the Earth, satellite and the third body in the Earth-centered inertial frame is shown in figure 5.2 In this configuration, the third body's position is given by

--→ EB = r EB cos ωt sin ωt (5.11)
where ω is the angular speed of the third body and r EB is the length --→ EB . The satellite acceleration is given by q

= -µ q q 3 + µ b --→ BE --→ BE 3 - --→ EB -q --→ EB -q 3 , q ∈ S ⊆ R 2 (5.12)
where S is the SOI and µ b is the standard gravitational parameter of the third body (the term µ b --→ BE --→ BE 3 comes from the fact that the Earth experiences an acceleration due to the third body).

q = r 0 , --→ EB = r EB cos(-α) sin(-α) . (5.13)
This allows us to decompose the vector of acceleration of the satellite due to the third body from (5.12), as

F =   -µ r 2 + µ d 2 γ -1 + cos α W 1+γW + √ 1-γW µ d 2 γ -sin α W 1+γW + √ 1-γW   , (5.14)
where γ and W are defined by

γ = r r EB , (1 -2 cos αγ + γ 2 ) 3 = 1 -γW. (5.15)
Generally, ( [START_REF] Ronwell | Satellite lifetimes in nearly circular orbits for various Earth-Atmosphere models[END_REF], [START_REF] Pascoli | Eléments de mécanique céleste[END_REF]), the assumption is made that r << r EB , that is, that γ is a small value, and so the acceleration vector (5.14) is expressed up to first order in γ. Here, we do not make use of this simplification since the transfer in subproblem 1 approaches the boundary of the Earth SOI (as we have defined it) and so the ratio r/r EB does not necessarily remain small throughout the transfer. Note that the perturbing force in (5.14) is time-dependent. The vectors F Sun and F Moon are special cases of the perturbing force F for which α, µ, γ in (5. to denote often-repeated terms, and so the optimal Hamiltonian H 1 (x, ℓ, λ, λ ℓ , t) = H 1 (x, ℓ, λ, λ ℓ , u * , t) gives rise the optimal dynamics (the Hamilton's equations)

ṗ = p p + 2p p µ E u * q 1 + e x cos ℓ + e y sin ℓ ėx = e p x + p µ E (Au * q + u * r D sin ℓ) 1 + e x cos ℓ + e y sin ℓ ėy = e p y + p µ E (Bu * q -u * r D cos ℓ) 1 + e x cos ℓ + e y sin ℓ l = ℓ p + (1 + e x cos ℓ + e y sin ℓ) 2 p p µ E λp = p µ E -3λ p u * q + 3λ ℓ D 3 2K 2 p 2 + 0.5λe y (-Bu * q +u * r D cos ℓ) p + λe x (-Au * q -u * r D sin ℓ) 2p 1 + e x cos ℓ + e y sin ℓ + ∂s p ∂p λex = p µ E 2λppu * q cos ℓ D -2λ ℓ D 2 cos ℓ K 2 p - λe x (u * q D+cos ℓ(-Au * q +u * q D cos ℓ)) D + pe y cos ℓ(Bu * q -u * q D sin ℓ) D 1 + e x cos ℓ + e y sin ℓ + ∂s p ∂e x λey = p µ E 2λppu * q sin ℓ D -2λ ℓ D 2 sin ℓ K 2 p + λe x (Au * q -u * q D cos ℓ) sin ℓ D + λ ey Bu * q sin ℓ D -u * q 1 + sin ℓ 2 1 + e x cos ℓ + e y sin ℓ + ∂s p ∂e y λℓ = - 2Cλ p p p µ E u * q D 2 + 2Cλ ℓ p µ E D K 2 p + λ ex p µ E -ACu * q + Cu * q D cos ℓ -u * r D 2 cos ℓ + u * q D(1 + D) sin ℓ D 2 +λ ey   - BC p µ E u * q D 2 - p µ E u * q (1 + D) cos ℓ D + p µ E (Cu * q -u * r D) sin ℓ D   + ∂s p ∂ℓ ,
where the terms ∂s p ∂x i are fairly bulky and are described numerically in our simulation. We denote the optimal dynamics ( ẋ, l, λx , λℓ ) = J(x, ℓ, λ x , λ ℓ , t).

(5.24)

We have stated that equations of the form (5.10) where x ∈ S , are simplified using averaging. In the numerical simulation, the averaging is applied in the form of a Riemann sum, is given by i × 2π/n ave . Since the summand in (5.25) depends on time, there are two possible choices: to average with a 'constant time' in the Riemann sum (5.25), or to notice that the time t is a function of the angular variable, t(ℓ), and so introduce an 'averaged time' t into the variables x. While this use of the 'averaged time' is realistic as the Solar and Lunar perturbations are timevarying (and so variation of time with respect to ℓ is a consideration), maintaining time as an independent variable and ℓ as a state variable helps in simplifying the link between averaged and non-averaged dynamics, which is an important consideration for us, since we intend to 'shut off' averaging during the course of the trajectory's evolution in solving the transfer problem (5.9).

S = nave i=1 J(x, ℓ * i , λ x , λ ℓ , t) 1 Q(ℓ * i , x, t) 2π n ave , ℓ i-1 ≤ ℓ * i ≤ ℓ i , (5.25

Numerical solution of subproblem 1

To perform the transfer in subproblem 1 numerically, we choose the exact initial Earth-centered orbit and final values

       p 0 = 1.881e + 07, e 0 x = 0.1, e 0 y = 0, ℓ 0 = 0, p f 1 + e f x cos ℓ f + e f y sin ℓ f = R f ,
where the numerical value R f is given in (5.4) and (x f , ℓ f , λ f , λ f ℓ ) indicate the solutions x(t, x 0 , λ 0 , λ 0 ℓ ), ℓ(t, x 0 , λ 0 , λ 0 ℓ ), λ(t, x 0 , λ 0 , λ 0 ℓ ), λ ℓ (t, x 0 , λ 0 , λ 0 ℓ ) evaluated at the final time t 1 f , which are computed numerically using the function ode45. In this subproblem, we will in fact simplify the final condition by making a transfer to the semimajor axis a f = R f of a circular orbit e f

x = e f y = 0; since the chosen final orbit is circular, this transfer is equivalent to a transfer to the final radius R f , however we chose this more simple case since the condition on the final radius (of a general non-circular orbit) involves the 'averaged longitude' l in the averaged case, which increases the computation time unnecessarily. Applying the Pontryagin's principle to the time-minimal transfer        ( ẋ, l) : (5.5) (x(0), ℓ(0)) = (p 0 , e 0

x , e 0 y , ℓ 0 ) a(t f ) = R f , e f x = 0, e f y = 0 min t f (5.26) gives rise to the transversality condition λ f ℓ = 0. Also, we have the condition (1.51) on the associated maximum Hamiltonian H 1 , since the final time t 1 f is free; we replace this condition with the equivalent condition on the norm of the adjoint vector, λ 0 2 + (λ 0 ℓ ) 2 = 1. Thus we have constructed the shooting function for subproblem 1,

S 1 : (t 1 f , λ 0 , λ 0 ℓ ) →       a f -R f e f x e f y λ f ℓ λ 0 2 + (λ 0 ℓ ) 2 -1       . ( 5

.27)

To find the solution to subproblem 1, we find a suitable initial guess (for the final time and the initial adjoint variable) to solve S 1 (t 1 f , λ 0 , λ 0 ℓ ) = 0 as described in section 1.3. The orbital transfer in rotating coordinates Once we have determined the time-optimal transfer with averaging numerically, we determine the final numerical value of the specific orbital energy evaluated at the final time of this trajectory,

E f = - µ E 2a(t f 1 )
.

(5.28)

This numerical value is used in final conditions in subproblem 2.

Bicircular four-body problem (subproblem 2)

As stated in (5.6), we use the notation X = (X 1 , X 2 , X 3 , X 4 ) instead of (X, Y, Ẋ, Ẏ ) ; then the system (1.36), (1.37) reads

Ẋ =      X 3 X 4 2X 4 • ω b + X 1 • ω 2 b + µ S R 3 S X S -µ E (X 1 +µ M ) r 3 1 -µ M X 1 -µ E r 3 2 -µ S X 1 -X S r 3 3 -2X 3 • ω b + X 2 • ω 2 b -µ S R 3 S Y S -µ E X 2 r 3 1 -µ M X 2 r 3 2 -µ S X 2 -Y S r 3 3      +     0 0 ǫ 0     u 1 +     0 0 0 ǫ     u 2 (5.29)
giving the dynamics of the form (5.7), where

P 0 =      X 3 X 4 2X 4 • ω b + X 1 • ω 2 b + µ S R 3 S X S -µ E (X 1 +µ M ) r 3 1 -µ M X 1 -µ E r 3 2 -µ S X 1 -X S r 3 3 -2X 3 • ω b + X 2 • ω 2 b -µ S R 3 S Y S -µ E X 2 r 3 1 -µ M X 2 r 3 2 -µ S X 2 -Y S r 3 3      , P 1 =     0 0 ǫ 0     , P 2 =     ǫ     .
(5.30)

The time-minimal transfer with these dynamics has the form

       Ẋ : (5.7) |u| ≤ 1 X(0) = X 0 , E(X 1 (t f ), X 2 (t f ), X 3 (t f ), X 4 (t f )) = E f , min |t f | (5.31)
where X 0 is the L 1 Lagrange point (section 1.1.7, equations (1.33)-(1.35)), and we target the final mechanical energy, where the function E has the form

E(X 1 , X 2 , X 3 , X 4 ) = -X 2 3 -X 2 4 + 2X 2 X 3 ω b -2(c + X 1 )X 4 ω b -(c + X 1 ) 2 + X 2 2 ω 2 b 2 + µ E (c 1 + X 1 ) 2 + X 2 2 .
(5.32)

Applying the Pontryagin maximum principle to this subproblem leads to the Hamiltonian of the form H 2 (X, Λ, u) = Λ(P 0 (X) + ǫ The orbital transfer in subproblem 1, from Earth orbit to the target radius R f . For the solution (X(X 0 , Λ 0 , t 2 ), Λ(X 0 , Λ 0 , t 2 )) derived using ode45 on the dynamics (5.34)-(5.41), we use the function T -1 : X → x to map synodic to (averaged) equinoctical variables, and denote (p E , ex E , ey E , ℓ E ) = T -1 (X 1 (t 2 , X 0 , Λ 0 ), X 2 (t 2 , X 0 , Λ 0 ), X 3 (t 2 , X 0 , Λ 0 ), X 4 (t 2 , X 0 , Λ 0 )) (λ E p , λ E ex , λ E ex , λ E ℓ ) = T -1 (Λ X 1 (t 2 , X 0 , Λ 0 ), Λ X 2 (t 2 , X 0 , Λ 0 ), Λ X 3 (t 2 , X 0 , Λ 0 ), Λ X 4 (t 2 , X 0 , Λ 0 )).

We will set conditions at time time t 1 for the four state variables p(t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) = p E , e x (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) = e E x , e y (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) = e E y , ℓ(t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) ≡ ℓ E mod 2π, (5.45) as our first four final conditions (we require 10 final conditions because the unknown (t 1 , t 2 , λ 0 , Λ 0 ) is 10-dimensional) in order that the two trajectories 'join' in the (q, q) space. The next 4 conditions are on the adjoint vector:

λ p (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) = λ E p , λ ex (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) = λ E ex , λ ey (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) = λ E ey , λ ℓ (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) = λ E ℓ .

(5.46)

These conditions ensure that the adjoint vector is continuous on an extremal. We will use the condition λ 0 2 + (λ 0 ℓ ) 2 = 1 to replace the condition (1.51) for unknown final time. We have to add a condition that defines t 1 or t 2 ; we choose to keep t 2 to be the time where the mechanical energy reaches E f in the second arc: this has the advantage that it does not depend on the fast variable ℓ. Together, these conditions give rise to the shooting function associated to the time-optimal transfer (5.9):

S 3 : (λ 0 , Λ 0 , t 1 , t 2 ) →                 
p(t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) -p E e x (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) -e E x e y (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) -e E y mod(ℓ(t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ), 2π) -mod(ℓ E , 2π) λ p (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) -λ E p λ ex (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) -λ E ex λ ey (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) -λ E ey λ ℓ (t 1 , x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) -λ E ℓ µ E (-1+exE 2 +eyE 2 ) 2pE

-Ef λ 0 2 + (λ 0 ℓ ) 2 -1                  , (5.47)
where the solutions x(t, x 0 , ℓ 0 , λ 0 , λ 0 ℓ ), ℓ(t, x 0 , ℓ 0 , λ 0 , λ 0 ℓ ), λ(t, x 0 , ℓ 0 , λ 0 , λ 0 ℓ ), λ ℓ (t, x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) are computed using ode45 for the initial values x 0 , ℓ 0 given in (5.44) and evaluated at the (unknown) time t 1 and the notation mod(•, 2π) denotes the measure of the angle between 0 and 2π.

To find the optimal trajectory of transfer (5.9), we solve S 3 (λ 0 , Λ 0 , t 1 , t 2 ) = 0 for the initial satellite parameters given in tables A.1 -A.2 and the initial condition in (5.44). In figure 5.5 we show the 'initial' transfer: we plot the trajectories (x(t, x 0 , ℓ 0 , λ 0 , λ 0 ℓ ), ℓ(t, x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) for t ∈ [0, t 1 ] and X(t, x 0 , ℓ 0 , λ 0 , λ 0 ℓ ) for t ∈ [t 1 , t 1 -t 2 ] (in Cartesian coordinates) and their components where the initial values λ 0 ,λ 0 ℓ , Λ 0 and t 1 and t 2 are the 'initial guesses' obtained from solving the shooting equations S 1 = 0, S 2 = 0. This is effectively 'the transfer before solving S 3 = 0'.

Shown in figure 5.6 is the development of the orbital elements ((a, e x , e y ), ℓ) throughout the transfer, and then the plot of the complete satellite orbit in Cartesian coordinates ( X, Y ), after initializing using the solutions λ 0 ,λ 0 ℓ , Λ 0 , t 1 and t 2 of the shooting equation S 3 (λ 0 , Λ 0 , t 1 , t 2 ) = 0 (this is the solution of the transfer problem (5.9)).

Note that the solution of transfer problem (5.9) where averaging is used on one arc of the solution is not a 'true' time-optimal solution; however, if averaging is not used, then the use of the different dynamics (5.5) and (5.7) on two arcs of the solution is simply a canonical change of coordinates in the Hamiltonian of a time-optimal problem (performed at the time t = t 1 ) and the solution in this case is a time-optimal solution. To see how good the level of approximation brought about by averaging is, we plot in figure 5.7 the development of the coordinates a, e x , e y for the same transfer shown in figure 5.6, but performed both with and without averaging in the Kepler subproblem. It can be seen from this figure that the problem with averaging is very close to the time-optimal transfer. This closeness between the averaged and non-averaged solution could further be improved by using a more complex function than T to change from averaged to non-averaged variables which would also involve a change in the initial conditions between the averaged and non-averaged problem, as demonstrated in [START_REF] Dargent | Initial and final boundaries transformation when solving optimal control problem with averaging techniques and application to low thrust orbital transfer[END_REF]. Components (a, e x , e y , ℓ) of the solution of transfer (5.9) with averaging where Ψ : x → x is defined in equation (5.51) and T : x → X, where T = T 1 • T 2 for the function T 1 : x → ξ which maps the equinoctical coordinates to Cartesian (Earth-centered) coordinates, and the function T 2 : ξ → X which maps the Cartesian to synodic coordinates.

Specifically, the components of T 1 are given by The inverse function is

T -1 = G -1 • Ψ -1 • G • T -1 D T • DG -1 • DΨ • DG : X Λ → (x, l) ( λ, λℓ ) ,
where Ψ -1 is given by and the components of T 3 are (5.57)

T 3 ( X, Y , ˙ X, ˙ Y ) =          ( ˙ Y X- ˙ X Y ) 2 µ E ˙ Y 2 X- ˙ X ˙ Y Y µ E - X √ X 2 + Y 2 - ˙ X ˙ Y X+ ˙ X 2 Y µ E - Y √ X 2 + Y 2 tan -1 Y X          . ( 5
DT 3 =          2 ˙ Y ( ˙ Y X- ˙ XY ) µ E -2 ˙ X( ˙ Y X- ˙ XY ) µ E -2 Y ( ˙ Y X- ˙ XY ) µ E 2 X( ˙ Y X- ˙ XY ) µ E ˙ Y 2 µ E + X 2 ( X 2 + Y 2 ) ( 3/2) - 1 ( X 2 + Y 2 ) ( 1/2) - ˙ X ˙ Y µ E + X Y ( X 2 + Y 2 ) ( 3/2) - ˙ Y Y µ E (2 ˙ Y X- ˙ X Y ) µ E - ˙ X ˙ Y µ E + X Y ( X 2 + Y 2 ) ( 3/2) ˙ X 2 µ E + Y 2 ( X 2 + Y 2 ) ( 3/2) - 1 ( X 2 + Y 2 ) ( 1/2) - ˙ Y X+2 ˙ X Y µ E - ˙ X X µ E -Y ( X 2 + Y 2 ) X ( X 2 + Y 2 ) 0 0          .
(5.58)

To determine the Jacobian DΨ -1 , we make use of the definition (5.54) of Ψ -1 , giving

DΨ -1 =     1 0 0 0 0 1 0 0 0 0 1 0 0 ∂M ∂e ∂M ∂ω + 1 ∂M ∂ℓ     .
(5.59)

We will now compute ∂M ∂e , ∂M ∂ω and ∂M ∂ℓ . By differentiating the equations ( [START_REF] Battin | An introduction to the mathematics and methods of astrodynamics[END_REF]) In order to find the components DΨ, D T to construct the function T -1 , we simply need to determine the inverse of the matrices DΨ -1 , DT 4 , DT 3 in equations (5.58), (5.57), (5.71), which can be done analytically or numerically.

cos v = cos E -e 1

Chapter 6 Conclusion

The study set out to examine the effects of a low thrust in a number of optimal satellite transfers in different space contexts. The low-thrust property leads to complications both in an analytical and numerical sense in the solution of optimal transfers, which Sophie Geffroy's PhD [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] showed can be simplified through the process of averaging. We sought to exploit that simplification to determine optimal transfers in a variety of situations. An underlying idea is that influences in the spatial region in which the transfer takes place ultimately dictates whether a transfer truly satisfies the low-thrust property necessary for averaging. This separates the satellite motion into three distinct problems:

-Kepler two-body motion dominated by gravity of primary -Perturbed two-body motion dominated by primary's gravity but affected by a perturbing force (which plays a similar role in satellite acceleration as engine thrust)

-Multi-body motion dominated by more than one primary; here, the thrust also plays a dominant role and is no longer a perturbation Ultimately, the study had two goals 1. Contribute to existing studies of averaged transfers in the two-body problem 2. Answer the question 'how does averaging apply within the other two spheres in the space environment?'

We have already summarized our contributions for each type of transfer in chapter 2; together, they fulfilled the study's two research goals. In the first case, we considered the time-minimal twobody transfer to 'circular' orbits and discovered firstly that, while the Hamiltonian is not analytic as in the energy-minimal case (considered by Bonnard, Caillau and Dujol, [START_REF] Bonnard | Energy minimization of single input orbit transfer by averaging and continuation[END_REF]) averaging reduced the number of singular directions in the Hamiltonian to two, and also (again unlike the energyminimal case) that the elliptic domain is geodesically convex, which since the elliptic domain is the entire spatial region where averaging can be applied, is conclusive. In the Lunar-perturbed two-body problem (section 2.2), we made use of the 'double averaging' with respect to both the satellite and Moon angular variables. This double averaging has some precedent in the literature (e.g. Pascoli, [START_REF] Pascoli | Eléments de mécanique céleste[END_REF]), but we applied it within an optimal control problem 121 and showed that it lead to a Hamiltonian of the perturbed motion which relates the problem to a Zermelo problem, showing it has some useful implications for perturbed optimal transfers.

In the second case, in the multi-body problem we determined that it was necessary to identify the 'switching time' t 1 based on the criterion of specific orbital energy in order to develop two arcs, one where averaging was applied, and the second where it was not (section 2.3). We constructed such a transfer and demonstrated by a numerical experiment that the use of this switching time leads to an accuarate solution of a time-optimal transfer performed in a four-body context.

The implications to future work is that the method of averaging can be sucessfully adapted to a variety of low-thrust transfers rather than just those which take place in a two-body context. Since two-body motion is fairly restrictive when it comes to space missions, this means that benefits of averaging feasibly extended to numerical and analytical studies of more physical real-world cases.

There is also much scope for further studies: firstly, in our study of time-minimal transfers, we restricted to the case of 'transfer to circular orbits' (i.e. the Hamiltonian system is in 4 variables instead of 6). It would be of interest to lift this restriction to include transfers where we do specify a target value of ω, the periapsis argument. The problem could be extended even further to include out-of-plane transfers with a full (3-dimensional) engine thrust; however, this leads to a 10-dimensional Hamiltonian which is a significant obstacle for analytical and even numerical study.

For the lunar-perturbed transfer, a more extensive study of the conjugate points, building on the numerical evidence, would be interesting. The methods used in this study could be extended to include other perturbations such as the J 2 perturbation in a double-averaged, lunar-perturbed system, making a superposition of these perturbations. There is practical scope for this, since highly-elliptic (e.g. Molniya) orbits experience both perturbations simultaneously. This would involve mainly a refining of the simplifying assumptions and reference planes used within our study, since different perturbations are not generally measured with respect to the same frame of reference.

In studying the four-body transfer we restricted to the case that satellite mass is constant. This assumption reduced the number of variables but is not very realistic, especially for a mission of longer duration; to extend to a more physical case it would be useful to drop this assumption. The same technique of averaging up until a switching time could be performed for a transfer from a Lunar orbit to the L 1 Lagrange point using exactly the same techinques. Having carried out such a Lunar orbit-to -L 1 transfer in minimum time, it would then be quite feasible to use the shooting method to join the two trajectories, creating an Earth-to-Moon transfer with averaging around both of the two primary bodies. Finally, since the use of averaging in the time-minimal transfer in the four-body context had a conclusive result, it would also be of practical use to perform a numerical study comparing such transfers under different epochs (the reference date and time at which the configuration of the Sun and Moon is initially measured). It has been evidenced that the Sun plays a role in assisting fuel-minimum satellite transfers, and such a scheme of extensive simulations at different epochs would obtain data establishing whether the presence of the Sun also benefitted time-minimum transfers. The use of averaging in the transfer dynamics is especially helpful here, because it reduces the computation time, allowing for multiple tests to be carried out more quickly.

CONTENTS 3 . 2 . 1 1

 3211 Hamiltonian formalism, Pontryagin maximum principle . . . . . . . . . . 49 3.2.2 Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2.3 Averaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.2.3.1 Singularities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3 The analysis of the averaged systems for minimum energy . . . . . . . . . . . . . 54 3.3.1 The coplanar energy case . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.2 Transfer towards circular orbits . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 The analysis of the averaged systems for minimum time . . . . . . . . . . . . . . . 55 3.4.1 The Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.2 Singularities of the Hamiltonian in the single-input and two-input cases . . 56 3.4.3 The Hamiltonian flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.4 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.5 Comparison of time and energy cases . . . . . . . . . . . . . . . . . . . . . . . . 63 3.6 Conclusion and open problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.7 Proof of Proposition 3.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.8 Proof of Theorem 3.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Lunar perturbed orbital transfer 81 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.2 The Riemannian metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.3 The perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3.1 Preliminairies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3.2 The lunar perturbation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.4 Computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4.1 Shooting equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4.2 The geometric concept of conjugate point . . . . . . . . . . . . . . . . . . 90 4.4.3 Computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.4.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.5.1 A brief description of the J 2 -effect [51] . . . . . . . . . . . . . . . . . . . 96 4.5.2 Additional perturbations [62] . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.5.3 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Piecewise Earth to L 1 transfer 99 5.1 Transfer scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 Perturbed two-body problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.2.1 Numerical solution of subproblem 1 . . . . . . . . . . . . . . . . . . . . . 106 5.3 Bicircular four-body problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.3.1 Numerical solution of subproblem 2 . . . . . . . . . . . . . . . . . . . . . 109 5.4 Numerical solution of complete transfer . . . . . . . . . . . . . . . . . . . . . . . 112 5.4.1 Simulation for transfer (5.9) . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.4.2 Details of the construction of the function T . . . . . . . . . . . . . . . . 117 Controlled two-body (Newtonian) problem

Figure 1 . 1 :

 11 Figure 1.1: Arrangement of the Lagrange points L 1 to L 5 in the Earth-Moon system

Figure 1 . 2 :

 12 Figure 1.2:The arrangement of the four bodies in the bicircular problem (adapted from[START_REF] Gómez | II Fundamentals: the case of triangular libration points[END_REF]).

  Applying the maximum principle to the Gauss equations with the criteria minimum time or the quadratic u 2 dt, (sections 1.2.1, 1.3.1), we get the Hamiltonian o.d.e. (1.55) with one of two expressions of H, in (1.57) or (1.60).

Figure 2 . 1 :

 21 Figure 2.1: Behaviour of F Sun and F Moon with increasing r (after [56])

Theorem 3 . 6 .

 36 The geodesics of the averaged coplanar transfer towards circular orbits are straight lines in the domain X (see (3.10)) in suitable coordinates, namely

-e 2

 2 p e = ρ sin ψ , e = sin ϕ , -π 2 < ϕ < π 2 . (3.25) Equations (3.23) and (3.24) then yield

Proposition 3 . 8 .

 38 The maps L : C → R and M : C → R are real analytic away fromS = (ψ, ϕ), tan ψ = 1 + sin ϕ 2 cos ϕ ∪ (ψ, ϕ) , tan ψ = -1 + sin ϕ 2 cos ϕ . (3.34)

  [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems and bifurcations of vector fields[END_REF] where J(ψ, ϕ, E) vanishes for two distinct values of the angle E where it changes sign; these two values are given by cos E = R(ψ, ϕ) with R(ψ, ϕ) = cos ψ P (ψ, ϕ) , P (ψ, ϕ) = 2 cos ϕ sin ψ -sin ϕ cos ψ (3.39) (note that (3.37),(3.38) amount to R 1 = {|R(ψ, ϕ)| > 1}, R 2 = {|R(ψ, ϕ)| < 1}and S is the locus where R = ±1). Hence (3.32) yields

Proposition 3 . 10 .

 310 In the coordinates (n, ϕ, ψ, ρ) defined by (3.25), and after a time re-parametrization dt = n 1/3 dτ , (3.43) equation (3.41) (resp. equation (3.42)) becomes dψ dτ = a(ψ, ϕ) , dϕ dτ = b(ψ, ϕ) , dn dτ = -3n c(ψ, ϕ) , (3.44)

3 (resp. d dτ n 1 / 3

 313 43), one easily deduces (3.44) and the expressions (3.45) (resp. (3.46)) of a, b, c. Finally, (3.44) and (3.45) (resp. (3.44) and (3.46)) imply d dτ n 1/3 cos ψ L(ψ,ϕ)= n -1/cos ψ M (ψ,ϕ)= n 1/3 ), that implies (3.49) (resp. (3.50)) according to(3.43).

Figure 3 . 1 :

 31 Figure 3.1: Numerical plot (obtained using Matlab) of the stable and unstable manifolds (bold) and trajectories through a number of arbitrary initial values in C for the full control case. The other curve shown is ψ = Z b (ϕ).

Figure 3 . 2 :

 32 Figure 3.2: Numerical plot (obtained using Matlab) of the stable and unstable manifolds (bold) and trajectories through a number of arbitrary initial values in C for the tangential case. The other curve shown is ψ = Z b (ϕ).

Figure 3 . 3 :

 33 Figure 3.3: The phase portrait, for energy minimization, in the same coordinates as Figure 3.1 and 3.2.There are two lines of non isolated equilibria. The darker zone is all the points that can be reached in positive time from the line {ϕ = ϕ 0 }, with ϕ 0 rather close toπ 2 . The highest possible final value is ϕ is ϕ 1 max = ϕ 0 + 2/5 π.

  Continuity of Φ in (3.55), is, according to [40, chap. V, Theorem 2.1], guaranteed by uniqueness of solutions and continuity of a, b.

3 . 0 cos 2 E dE 1 -

 3021 Sign and zeroes of b. 3.1. Full control case (a, b given by (3.45)). On the one hand, one has b(0, ϕ) = 2 cos ϕ sin ϕ π π sin 2 ϕ cos 2 E , b(-π, ϕ) = -b(0, ϕ) . (3.71)

Figure 3 . 4 :

 34 Figure 3.4: Numerical plot of the function (3.73) on the interval [0, 1). Obtained with Maple 15.

. 4 ,+ e 2 √ 1 -e 2 , 3 . 2 .cos ϕ sin ϕ π π 0 cos 2 E 1 -

 42123221 the graph of the function e → γ(arcsin e) on [0, 1], i.e. e → b arctan -1 arcsin e . (3.73) On the one hand, this proves that b is smooth at points where it vanishes; on the other hand the derivative of b with respect to ψ is (see (3.72)) strictly positive at (Z b (ϕ), ϕ), ϕ > 0. This implies smoothness of Z b according to the inverse function theorem; and this extends to negative ϕ with Z b (0) = 0, hence point 3 of the proposition is satisfied with σ = 0; it is also easy to check that b(ψ, 0) only if ψ = 0 or ψ = π. Tangential thrust case (a, b given by (3.46)). In the region R 1 , one has b(ψ, ϕ) = (sign cos ψ) 2 sin 2 ϕ cos 2 E dE . (3.74)

Figure 3 .

 3 5 displays a numerical plot of the graph of the map ϕ → a(Z b (ϕ), ϕ) in the full control case; we also show ϕ → Z b (ϕ).

Figure 3 .

 3 6 displays a numerical plot of the graph of ϕ → a(Z b (ϕ), ϕ) in the tangential thrust case; we also show ϕ → Z b (ϕ) and ϕ → arctan -1+sin ϕ 2 cos ϕ to show that it is very close to Z b (ψ = arctan -1+sin ϕ 2 cos ϕ

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Plots (obtained with Matlab) of the maps ϕ → Z b (ϕ) (dashed line) and ϕ → a(Z b (ϕ), ϕ) (solid line) in the full control case.

7. 2

 2 tangential thrust case. Let us first compute some values of a and b on special lines 2 : according to (3.46) and (3.32)-(3.33),

τ fin 0 c

 0 (ψ ♯ (τ ), ϕ ♯ (τ ))dτ = τ fin 0 c(ψ + (τ ), ϕ + (τ ))dτ = τ fin 0 c(ψ(τ ), ϕ(τ ))dτ .(3.82)

  .86) and (3.58) may be replaced by b

  .87) From (3.84), (3.85) and (3.87), the relation (3.66) translates into

  .89) Also, (3.58) and (3.59) become: b(ψ, ϕ) = 0 ⇔ ϕ = 0 and ψ ∈ [-σ, σ] ∪ [π -σ, π + σ] or ϕ = 0 and ψ ∈ {Z 0 b (ϕ)} ∪ {Z π b (ϕ)} , (3.90) b(ψ, ϕ) > 0 if ϕ = 0 and σ < ψ < π + σ or ϕ = 0 and Z 0 b (ϕ) < ψ < Z π b (ϕ) .

Figure 3 . 7 :

 37 Figure 3.7: This picture reflects qualitatively the assumptions on a and b. It is provided as a help to follow the proof; a precise numerical drawing is provided in Figure 3.1 and Figure 3.2 for the specific expression of a, b in the case of full control or tangential thrust. The six invariant regions separated by the stable and unstable manifolds of (0, 0) and (π, 0) are shown. The other curves are ψ = Z 0 b (ϕ) and ψ = Z π b (ϕ), where b(ψ, ϕ) changes sign.

Proof. Points 1

 1 and 3 are obvious because b is negative in E and in the upper part of U 0 while it is positive in the upper part of S 0 except, if σ is nonzero, on the segment {ϕ = 0}, where it is zero. Let us prove point 2. In the region F , according to (3.59), b has the sign of ψ -Z 0 b (ϕ). Using differentiability of Z 0 b away from ϕ = 0 (see (3.63)), one may compute the derivative of ψ -Z 0 b (ϕ) with respect to time along a solution; it is is a(ψ, ϕ) -Z 0 b ′ (ϕ)b(ψ, ϕ), which, according to (3.60), is positive when ψ -Z 0 b (ϕ) = 0, i.e. when b(ψ, ϕ) = 0. Hence the region where b > 0 is positively invariant, this accounts for the behavior of solutions starting in {b > 0}, and no solution may stay on the locus where b = 0, this accounts for solutions that start in {b < 0}: either they stay in this part of F or they cross {b = 0} at one time and then remain in {b > 0}.

.98) Lemma 3 . 15 .

 315 Assume that a, b, c satisfy assumptions (3.53) to (3.62).

  3.111) and(3.99) and (3.100) hold: we only need to prove that (3.101) holds as well.First, "thicken" the curves ψ = Z 0 b (ϕ) and ψ = Z π b (ϕ) where b vanishes, using the flow Φ (see (3.55)):

  (3.119) Setting T (f ) = σ/a + 1 + (π-2ℓ)/m(f ), (3.115), (3.117), (3.118) and (3.119) imply (3.101).

25 2 sin 2

 2522 (cφ) 1+4 cos 2 (cφ) and c = 2 5 .

Definition 4 . 8 .Proposition 4 . 9 .

 4849 A solution of the systemİ = 0 φ = w(I),where (I, φ) are variables in (R d , T d ) is called quasi-periodic. There exists quasi-periodic solutions in the restricted circular planar three body problem.

r a and a ′ r ′

 r From the Kepler's equation the eccentric anomaly E satisfies E = M + ρ sin(E), the 2π-periodic function E → ρsin(E) can be expanded into Fourier series and E can be expressed asE = M + ∞ k=1 a k sin(kM )where a k = 2 π π 0 ρ sin(E) sin(kM )dM = 2 k J k (kρ), denoting J m (z) the Bessel functions of the first kind defined as J m (z) = 1 π π 0 cos(mθ -z sin(θ))dθ, (m ∈ Z, z ∈ C).

Figures 4. 1

 1

  and 4.2 yield two first integrals of the double averaged Hamiltonian system. The double average is taken with respect to the mean motion variables M and M ′ which correspond respectively to the satellite and the Moon. This integral over (M, M ′ ) ∈ [0, 2π] × [0, 2π] is computed with M ′ fixed with respect to M and the slow variables. The variation of θ(t) is equal to 3 4 n ′2 n (1 -ρ 2

Figure 4 . 2 :Figure 4 . 3 :

 4243 Figure 4.2: Evolution of n, ψ and θ of the double averaged (solid line) and the non averaged (dotted line) free system over ten lunar revolutions around the Earth.

Figure 4 . 4 :

 44 Figure 4.4: Evolution of adjoint vectors of extremal trajectories from the initial state point (n 0 , ρ 0 , θ 0 ) = (48.3, 0.60, π) to the final state point (n f , ρ f , θ f ) = (35.2, 0.10, π + 1). Final points are indicated.

Figure 4 . 5 :

 45 Figure 4.5: Evolution of state vectors of extremal trajectories from the initial state point (n 0 , ρ 0 , θ 0 ) = (48.3, 0.60, π) to (n f , ρ f , θ f ) = (35.2, 0.10, π + 1). Conjugate points are indicated (star markers).

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: Evolution of adjoint vectors of extremal trajectories from the initial state point (n 0 , ρ 0 , θ 0 ) = (48.3, 0.60, π) to the final state point (n f , ρ f , θ f ) = (35.2, 0.10, π + 1). Conjugate points are indicated.

detFigure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Rank condition for the determination of the first conjugate point for the perturbed case and the unpertubed one. The first zero of the determinant is the first conjugate time.

Figure 4 . 10 :

 410 Figure 4.10: Projection of extremal trajectories in the perturbed case (λ = 1) in (ψ, θ) coordinates starting from the same initial point (ρ 0 , θ 0 ) = (0.60, π). Conjugate points are indicated.

. 3 )

 3 The satellite motion described by (5.1) takes place both within and outside of the Earth SOI, since the transfer targets the L 1 -point, which lies in neither the SOI of the Earth or Moon (at this point the Moon and Earth gravity are precisely equal). Thus we will divide the transfer, as shown in the scheme depicted in figure5.1, into an arc within the Earth SOI and a second within the four-body problem as the satellite passes to the L 1 Lagrange point. Since for the controlled satellite, the radius r(t) is a function of time, in our computations this amounts to designating a particular time t * at which the satellite leaves the Earth SOI. For the chosen value of ǫ ratio = 4.116 × 10 -4 , and for the values of the (constant) mass m and thrust force F given in tables A.1 and A.2 used in these numerical simulations, this 'target radius' R f has the numerical value R f = 3.307546 × 10 7 .

Figure 5 . 1 :

 51 Figure 5.1: Possible transfer scheme from Earth orbit to L 1 showing the Earth and Moon SOIs Thus, we have established the transfer scheme 1. From initial orbit, shoot to the chosen radius R f in minimum time with averaging in twobody motion perturbed by Moon and Sun gravity. The mechanical energy of the satellite at the final time t 1 f of this transfer (denoted by E f ) will be obtained numerically (we expect a low negative energy)

Figure 5 . 2 :

 52 Figure 5.2: Configuration of Earth, satellite and third perturbing body in the Earth-centered perturbed 2body subproblem

  14) are replaced by α M , α S , respectively the Moon and Sun angular variables measured from epoch µ M , µ S , respectively the Moon and Sun standard gravitational parameters γ 1 , γ 2 , respectively the ratios r EM r and r ES r where r EM and r ES are the Earth-Moon and Earth-Sun distances.

  ) where ℓ ∈ I = [0, 2π] and the partition of I is Part = {[0, 2π/n ave ], [2(2π/n ave ), 3(2π/n ave )], . . . , [(n ave -1)(2π/n ave ), 2π]} , for n ave a numerical parameter chosen to fix the number of steps during averaging. The value ℓ * i

4 .

 4 In order to solve the problem, it is also necessary to obtain the initial angular displacement of the Moon and Sun, α 0 M and α 0 S by chosing an epoch. We choose the epoch 15 May 2014, 3:35 UTC. The spacecraft initial parameters as well as the parameters of the planets used in this numerical computation are given in tables A.1 and A.2. The optimal trajectory and the development of each of the orbital elements is shown in figure 5.4.

  Evolution of components in the Earth-orbit to R f transfer

Figure 5 . 3 :

 53 Figure 5.3: The orbital transfer in subproblem 1, from the Earth orbit to the target radius R f

2 i=1P

 2 i (X)u i ) -1, (5.33)Evolution of orbital elements in transfer from L 1 to the target mechanical energy E f .

Figure 5 . 4 :

 54 Figure 5.4: The orbital transfer of subproblem 2 from L 1 orbit to the target mechanical energy E f .

  Components of the 'initial' transfer

Figure 5 . 5 :

 55 Figure 5.5: Orbital transfer and (a, e x , e y , ℓ) components using the initial guess of time and adjoint vectors determined in subproblems 1 and 2.

Figure 5 . 6 :

 56 Figure 5.6: The orbital transfer (5.9) from Earth orbit to the Lagrange point L 1 with averaging on the first arc.

Figure 5 . 7 :

 57 Figure 5.7: Components a, e x , e y of Earth orbit to L 1 transfer, with and without averaging. Zoom on a shows oscillation of non-averaged a variable in Kepler arc.

T 1 (-c 1 +Finally, G - 1

 111 p, e x , e y , ℓ) X cos(ω b t) + Y sin(ω b t) Y cos(ω b t) -X sin(ω b t) ( ˙ X + ω b Y ) cos(ω b t) + ( ˙ Y -ω b X) sin(ω b t) ( ˙ Y -ω b X) cos(ω b t) -( ˙ X + ω b Y ) sin(ω b t)has the form G -1 (p, e, ω, ℓ) = (p, e cos ω, e sin ω, ℓ).

Ψ - 1 : 4 (X 1 , X 2 , X 3 , X 4 (c 1 +

 1412341 (p, e, ω, ℓ) T -1 is the composition T -1 = T 3 • T 4 of the functions T 4 : X → ξ which maps synodic orbital elements to Earth-centered Cartesian elements and T 3 : ξ → x which maps Cartesian elements to equinoctical variables.The components of T 4 are given byT X 1 ) cos(ω b t) -X 2 sin(ω b t) X 2 cos(ω b t) + (c 1 + X 1 ) sin(ω b t) (-ω b X 2 + X 3 ) cos(ω b t) -(ω b (c 1 + X 1 ) + X 4 ) sin(ω b t) (ω b (c 1 + X 1 ) + X 4 ) cos(ω b t) + (-ω b X 2 + X 3 ) sin(ω b t)

. 56 )

 56 The gradient DΨ -1 • D T -1 used to construct the function T is computed by determining the Jacobians D T -1 = D(T 3 • T 4 ) = DT 3 • DT 4 , where sin(ω b t) -ω b cos(ω b t) cos(ω b t) -sin(ω b t) ω b cos(ω b t) -ω b sin(ω b t) sin(ω b t) cos(ω b t)

-e cos E , sin v = √ 1 -e 2 sin E 1 -√ 1 -e 2 ( 1 -e 2 ) ( 1 + 1 -√ 1 -e 2 ( 1 -e 2 ) ( 1 + e cos v) 2 dv de - ( 2 +e cos v) sin v √ 1 -e 2 ( 1 + 1 -e 2 ( 1 -e 2 ) ( 1 + e cos v) 2 ( 1 - 2 √ 1 -e 2 ( 1 -e 2 )(1+e cos v) 2 + 1 √ 1 -e 2 ( 1

 111212111212121211212121212121121 e cos E ,(5.60) which relate the true and eccentric anomalies, we obtaindE = -sin v √ 1 -e 2 (1 + e cos v) de + √ 1 -e 2 1 + e cos v dv,(5.61)where, from the expression M = e -e sin E, thendM = (1 -e cos E)dE -sin Ede (5.62)and substituting in dE from (5.61) givesdM = e cos v) 2 dv -(2 + e cos v) sin v √ 1 -e 2 (1 + e cos v) e 2 (1 -e 2 ) (1 + e cos v) 2 dv dv -(2 + e cos v) sin v √ 1 -e 2 (1 + e cos v) 2 de dv (5.64) = √ 1 -e 2 (1 -e 2 )(1 + e cos v)2 (5.65)and∂M ∂e = e cos v) sin v √ 1 -e 2 (1 + e cos v) e cos v) 2 .(5.67)Using the relation v = ℓ -ω, then ∂M ∂ℓ = e 2 (1 -e 2 ) (1 + e cos v) 2 , constructed the function T = T • G -1 • Ψ • G DG -1 • DΨ -1 • DG • D T -1 : (x, l) ( λ, λℓ ) → X Λ .

  [START_REF] Arnold | Mathematical methods of classical mechanics[END_REF]

Table 1 .

 1 1.1.

	Orbital altitude	Perturbation
	LEO	asymmetry of the Earth
	160km < alt < 2000 km	atmospheric drag
		tidal deformations of oblateness
		Earth's magnetic field
	MEO	atmospheric drag
	2000 km < alt < 35786 km third body perturbations
		Earth's magnetic field
	HEO	third to n-body perturbations
	35786km < alt	light pressure
		cosmic particles

1: Perturbations associated to orbital altitudes

  The considerations in the present section are valid both in the full control case and in the "tangential thrust" case. of notation we use H e , H t to represent both the Hamiltonians and the averaged Hamiltonians, although the inputs into these functions are different). These may be re-computed in terms of H 1 , H 2 . Unlike when H 1 , H 2 are used in the computation, using the Hamiltonian lifts H 1 , H 2 allows for an explicit expression of the averaged Hamiltonians H e (x, p), H t (x, p). Making the change of variables x 0 = Ξ(e, ω, l) -with Ξ deduced from x 0 = (E -e sin E)/n, (3.13) and (3.8)-in the integral, and using the facts that ∂Ξ/∂l = 1/Ω(x, l) and H

	The respective averaged Hamiltonians are
		H e (x, p) =	n 2π	0	2π/n	H e (x, p, x 0 )dx 0	(3.16)
		H t (x, p) =		n 2π	0	2π/n	H t (x, p, x 0 )dx 0 .	(3.17)
	(for ease then, using (3.12),
	H e (x, p) =	(1 -e 2 ) 3/2 2π	0	2π	i=1,2

Remark 3.2 (Tangential thrust). If the normal component u n is forced to be zero, there is a single term in the sums in

(3.15)

, and these equations become H e = H

2

1 , H t = H 1 . e (x, p, Ξ(e, ω, l)) = H e (x, p, l) , H t (x, p, Ξ(e, ω, l)) = H t (x, p, l) ,

  [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF] contains two terms, and we obtain

	H t (n, e, p n , p e ) =	(1 -e 2 ) 3/2 2π	0	2π	2 i=1	h i (n, e, p n , p e , v) 2	dv (1 + e cos v) 2 ,	(3.23)
	In the tangential thrust case it only contains h 1 -see remark 3.2-and we get (the superscript 1
	in H 1 t denotes single input):
	H 1

t (n, e, p n , p e ) = (1 -e 2 ) 3/2 2π 2π 0 h 1 (n, e, p n , p e , v) dv (1 + e cos v) 2 . (3.24)

  1. Observe that F 0 represents the current of magnitude |F 0 | g . If |F 0 | g < 1, this defines a Finsler metric. Definition 4.5. If we apply the Maximum principle to the previous optimal problem this defines an Hamiltonian which is homogeneous in p. Conversely, one can associate to the Hamiltonian H = H P + λ H a Zermelo navigation problem, where λ is a scaling parameter associated to the maximal control magnitude.

  4.4.2 The geometric concept of conjugate point Definition 4.11. Let z = (x, p) be a reference extremal solution of H on [0, t F ].

			The variational
	equation	δz(t) = d	-→ H (z(t))δz(t)

  Evolution of n, ψ and θ of the double averaged (solid line) and the non averaged (dotted line) free system over one lunar revolution around the Earth.

						n
	48.06
	48.04
	48.02
	48
	47.98
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
					t (nbs revolutions Moon)
						ψ
	1.562
	1.561
	1.56
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
					t (nbs revolutions Moon)
						θ
	1.6
	1.5
	1.4
	1.3
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
					t (nbs revolutions Moon)
	47.98 48 48.02 48.04 48.06 Figure 4.1: 0	1	2	3	4	5 n	6	7	8	9	10
					t (nbs revolutions Moon)
						ψ
	1.562
	1.561
	1.56
	0	1	2	3	4	5	6	7	8	9	10
					t (nbs revolutions Moon)
						θ
	2.5
	2
	1.5
	0	1	2	3	4	5	6	7	8	9	10
					t (nbs revolutions Moon)

The journal version unfortunately contains some misprints in equations (79)-(80); they are corrected here.

ACKNOWLEDGEMENTS

-To Bernard Bonnard, but for whom this thesis would never have happened in the first place -To Jean-Baptiste Pomet, for perserverence, patience and perfectionism APPENDIX A. TABLES Bibliography

Chapter 4

Lunar perturbation of the metric associated to the averaged orbital transfer This chapter is a reproduction of the paper 'Lunar and J2 perturbations of the metric associated to the averaged orbital transfer.' which is authored by B. Bonnard, H. Henninger and J. Rouot, [START_REF] Bonnard | Lunar and J 2 perturbations of the metric associated to the averaged orbital transfer[END_REF] .

Introduction

Recent space missions like lunar Smart-1 mission, Boeing orbital transfer, using electric propulsion are innovative design feature to reduce launch costs and lead to the analyse of the low thrust controlled Kepler equation using averaging techniques in optimal control. Pioneering work in this direction associated to the energy minimization problem are due to Edelbaum [START_REF] Edelbaum | Optimum low-thrust rendezvous and station keeping[END_REF][START_REF] Edelbaum | Optimum power-limited orbit transfer in strong gravity fields[END_REF], Epenoy-Geffroy [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] and more recently to Bonnard-Caillau [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF][START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF]. Under some simplifying assumption they lead to the definition of a Riemannian distance between Keplerian orbits, and this is a preliminary step in computing the time minimal or find mass maximizing solutions using numerical continuation techniques ( [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF], [START_REF] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation : resolution par homotopie différentielle[END_REF]).

The objective of this article is to analyse the deformation of this metric taking into account the lunar perturbation which affect a wide range of missions. Again in the framework of the continuation techniques, we shall make simplifying assumptions. The main point is to deduce from the averaged system the qualitative policy to make the transfer and to initialize the shooting algorithm.

Making such assumptions leads to the analysis of a Zermelo navigation problem defined by an Hamiltonian which is a deformation of the Hamiltonian associated to the Riemannian metric and complete analysis of the transfer is made using continuation about trajectories computations and conjugate and cut analysis.

The organization of this article is the following. In section 4.2, we recall the computations and properties of the Riemannian metric based on [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF]. In section 4.3, we present lunar perturbation and we describe the averaged system. In section 4.4, we give numerical simulations computations on extremal trajectories and on conjugate loci. In section 4.5, we conclude by considering more general perturbations. where x p = G(x, ℓ) • (F Sun + F Moon ), and the components are

We denote the sum s p = λ • x p = λ p p p + λ ex e p x + λ ey e p y . Applying Pontryagin's principle (theorem 1.1) to subproblem 1 (for which the dynamics have the form (1.14) where G is the matrix associated to (1.12) in planar form) gives rise to a Hamiltonian of the form

for which the optimal controls are

where (for brevity) we use the coefficients for which the optimal control is

and the optimal Hamiltonian of the form H 2 (X, Λ) = H 2 (X, Λ, u * ) gives rise to the Hamilton's equations

(5.35)

(5.36)

(5.37)

3

where X S , Y S and r 1 , r 2 , r 3 are as defined in (1.36), (1.37), and (for simplicity) we denote

(5.42)

Numerical solution of subproblem 2

To perform the transfer in subproblem 2 numerically, we choose the exact initial values

where the numerical value of E f is the final specific orbital energy described by equation (5.28) for the optimal trajectory of subproblem 1. From the Pontryagin maximum principle, as described in section 1.3.2, our final condition (5.43) gives rise to the three transversality conditions

where (X f , Λ f ) denote the solutions X(t, X 0 , Λ 0 ), Λ(t, X 0 , Λ 0 ) evaluated at the final time t = t 2 f , which are determined by solving the Hamilton's equations numerically using ode45 (for the numerical values of X 0 given in the equation (5.43)). Finally, the condition (1.51) for the case when the final time t f is free on the associated maximum Hamiltonian is replaced with the equivalent condition on the norm of the adjoint vector, Λ(t 2 0 ) = 1. Thus for choices of initial and final orbit, the shooting function for subproblem 2 is

To find the solution to subproblem 2, we find a suitable initial guess for the final time and initial adjoint vector to solve the shooting equation S 2 (t f 2 , Λ 0 ) = 0. The spacecraft initial parameters are given in tables A.1 and A.2. The development of each of the orbital elements, and the optimal trajectory in figure 5.4 Note that the 'spiraling' of the trajectory before converging to specific orbital energy E f signifies that the orbit has become elliptic, which is what we would expect for it to target the low energy E f .

Details of the construction of the function T

In section 1.3.5 we noted that the average x is the first asymptotic approximation to x; that is, x = x + ǫ for some 'small' ǫ. In fact, we will further assume that (p, ēx , ēy ) = (p, e x , e y ) (5.48) for the averaged optimal trajectory of subproblem 1.

Using the fact that lim u→0 l = L, we take L = l where L is the mean longitude. That is,

where M is the mean anomaly, which can be expressed in terms of e and v by

(5.50)

The inverse of M in (5.50) for a particular e cannot be computed explicitly, but the implicit function M -1 such that v = M -1 (M, e) can be computed numerically using Newton's method. Thus in what follows, in converting from l to ℓ, we will make use of the expression ℓ = v + ω, in the form

a numerical expression of ℓ. Since ℓ is best described using the variables (e, ω) rather than (e x , e y ). Thus we will make use of the change of variables G : (p, e x , e y , ℓ) → (p, e, ω, ℓ), where G(p, e x , e y , ℓ) = (p, e 2

x + e 2 y , tan -1 ey ex , ℓ). We construct the function

where, by the definitions in (5.48), (5.49)

i.e. Ψ truly sends averaged coordinates to 'non-averaged' coordinates, where averaging is defined by (5.48), (5.49).

The function T (x, l, λ, λℓ ) = (X, Λ)

is best described as a composition of functions, and has the form