{'fraction_non_alphanumeric': 0.060561299852289516, 'fraction_numerical': 0.01965300553288436, 'mean_word_length': 3.7637447823494337, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time Versus Energy in the Averaged Optimal Coplanar Kepler Transfer towards Circular Orbits

The aim of this note is to compare the averaged optimal coplanar transfer towards circular orbits when the costs are the transfer time transfer and the energy consumption. While the energy case leads to analyze a 2D Riemannian metric using the standard tools of Riemannian geometry (curvature computations, geodesic convexity), the time minimal case is associated to a Finsler metric which is not smooth. Nevertheless a qualitative analysis of the geodesic flow is given in this article to describe the optimal transfers. In particular we prove geodesic convexity of the elliptic domain.


Introduction
We consider the controlled Kepler equation describing orbital transfers with low thrust engines, that we normalize asq = − q q 3 + u; (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: 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 andq 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 distribution spanned by ∂/∂q 1 , ∂/∂q 2 in the original cartesian coordinates). In these coordinates, the free motion isẋ = 0,l = Ω(x, l); there may be a control term inl 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 [8,9], [12,11], [4]. 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 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 [4]. Moreover this result is still true if the thrust is oriented only in the tangential direction [5].
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 [3], these norms are not associated with inner products on these tangent spaces -this defines a Finsler metric [2], 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 2 recalls the equations and the computation of the average Hamiltonians; section 3 recalls the results from [4,5] on the minimum energy problem; section 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 5 explains why that proof fails in the minimum energy problem, which is consistent with the non-convexity mentioned in [4].

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 thaṫ x ∈ X, (u 1 , u 2 , ..., u ) ∈ B ⊂ R 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 (2) is the function 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 [15] (see also [6,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, 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 (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: 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: 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 + 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 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 ω .
The equations become: The angle v is the true anomaly v = l − ω.
In these coordinates, the elliptic domain is Remark 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: The fact that negative values of e are allowed comes from identifying (−e, ω) with (e, ω + π), or, equivalently, considering that (e x , e y ) (see (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}.
where u 1 , u 2 stand for u n , u t , x = (n, e, ω), the vectors F 1 , F 2 are readily obtained from (7a)-(7c), and 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 x 0 = (E − e sin E)/n , 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 becomeṡ 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):

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 which is approximated for small ε by For this system, we consider the following minimization problems: Applying the Pontryagin maximum principle leads to the following respective Hamiltonians (normal case in the energy minimization problem), where the lifts H i , defined by (14), are periodic with respect to x 0 with period 2π/n.
Remark 2 (Tangential thrust) If the normal component u n is forced to be zero, there is a single term in the sums in (15), and these equations become H e = H 2 1 , H t = H 1 . The considerations in the present section are valid both in the full control case and in the "tangential thrust" case.
The respective averaged Hamiltonians are  (13) and (8)-in the integral, and using the facts that ∂Ξ/∂l = 1/Ω(x, l) and H e (x, p, Ξ(e, ω, l)) = H e (x, p, l) , H t (x, p, Ξ(e, ω, l)) = H t (x, p, l) , then, using (12), Remark 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 forṁ 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 are the Hamiltonian lifts of F i (x). More complicated extremals are related to the switching surface Σ : 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 Σ) [7].

The analysis of the averaged systems for minimum energy
First of all we recall the results from the energy case [4] . The energy minimization problem is expressed as 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 1 In the coordinates (n, e, ω) the averaged Hamiltonian (up to a positive scalar) is given by where the singularity e = 0 corresponds to circular orbits. In particular (n, e, ω) are orthogonal coordinates for the Riemannian metric associated to H, namely Further normalizations are necessary to capture the main properties of the averaged orbital transfer.

Transfer towards circular orbits
As noticed in Remark 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 (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: Theorem 1 The geodesics of the averaged coplanar transfer towards circular orbits are straight lines in the domain X (see (10)) in suitable coordinates, namely with c = 2/5. Since c < 1, the domain is not (geodesically) convex and the metric is not complete.
Remark 4 (Tangential thrust) The properties of theorem 1 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 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 ϕ 4 The analysis of the averaged systems for minimum time

The Hamiltonian
We compute H t according to (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) .
The right-hand side does not depend on the cyclic variable ω, see Remark 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 (7), we get 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 (19). In the full control case (both tangential and normal control), the sum in (19) contains two terms, and we obtain In the tangential thrust case it only contains h 1 -see remark 2-and we get (the superscript 1 in H 1 t denotes single input): 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 Equations (23) and (24) then yield with L and M some functions C → R, where C is the cylinder The expressions of L and M are, taking the eccentric anomaly E as the variable of integration instead of v (see (13); in particular, dv/(1 + e cos v) = dE/ q − e 2 ) and restricting the revert of integration from [0, 2π] to [0, π] because the integrand depends on cos E only: and In the sequel we take advantage of this double homogeneity with respect to ρ and n.

Singularities of the Hamiltonian in the single-input and two-input cases
According to (26) and (27), the Hamiltonians H t and H 1 t , have the same degree of smoothness as, respectively the maps L and M .
Proposition 3 The maps L : C → R and M : C → R are real analytic away from They are both continuously differentiable on C, but their differentials are not locally Lipschitz-continuous on the set S; we have the following moduli of continuity of the differentials: in a neighborhood of a point ξ = (ψ, ϕ) ∈ S, Proof The set S is the set of points (ψ, ϕ) such that I(ψ, ϕ, E) vanishes for some value of E; hence the integrand in (29) is real analytic on C \ S and so is L. The degree of regularity (35) for L at points in S is given in [3].
Let us now treat M . It turns out that S is also the border between the region where the sign of J(ψ, ϕ, E) does not depend on E and the region 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 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 [3], 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 (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: with H t given by (23) and H 1 t by (24). Specifically, we establish geodesic convexity of the elliptic domain X (see (10)), i.e. any two points in X can be joined by a extremal curve. This is contained in the following result: Theorem 2 (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 (41) (resp. of (42)) such that (n(0), e(0)) = (n 0 , e 0 ) and (n(T ), e(T )) = (n 1 , e 1 ).
In order to ease the proof, let us write (41) and (42) in other coordinates.
The first two equations in (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(τ ): 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 (44).
Let us now gather some properties of the maps a, b, c, i.e. of the differential equation (44), that are valid both for a, b, c given by (45) and for a, b, c given by (46); they contain all the information to prove Theorem 2.
Proposition 5 The maps a, b, c given by (45) satisfy the following properties with σ = 0. The maps a, b, c given by (46) satisfy the same properties with σ = arctan 1 2 .
2. Uniqueness of solutions. The following differential equation on C: This defines a flow Φ from an open subset of C × R to C such that the above unique solution is 3. Sign and zeroes of b. There exists a continuous map continuously differentiable on the open interval (0, π 2 ), such that 4. Sign and zeroes of a. One has 5. Hyperbolic saddle point at (0, 0). The maps a and b are smooth in a neighborhood of (0, 0) and 6. Values of c at equilibria.
7. Stable and unstable manifolds of (0, 0). There exists continuous maps continuously differentiable on the open interval (0, π 2 ), and a number σ with such that the stable and unstable manifolds of (0, 0) are described by Furthermore, the zeroes of b are positioned with respect to the stable and unstable manifolds so that the maps S, U, Z b satisfy: Proof See Appendix A.
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 5, the differential equation (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 5. Theorem 2 will be easily deduced from Theorem 3.

Simulations
A numerical simulation of the phase portrait of the differential equation (54) (or the first two equations in (44)) is displayed in Figure 1 in the "full control case" where a and b are given by (45) and in Figure 2 in the "tangential thrust case" where a and b are given by (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 5; 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 1) but not in the tangential thrust case (Figure 1) where it comprises a segment of the ψ-axis.
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 Appendix B) where all trajectories go "up" (ϕ is monotone increasing), one (called F + in Appendix B) where all trajectories go "down", and four other regions (called E, E , E + and E + in Appendix B) 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. The generic figure 7 is a drawing used to support that proof, that figures in an illustrative manner the features contained in the assumptions of Theorem 3, and that can also be observed in the numerical simulations of the two cases that we are really interested in (Theorem 2).

Comparison between the minimum-energy and minimum-time cases from the convexity point of view
In section 3 we recalled some results from [4] (and previous work by the same authors); in particular, Theorem 1 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 minimum-energy case. Using the coordinates from Theorem 1 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 (20) and can be written as follows in the coordinates (n, ϕ, p n , p ϕ ) that result from the symplectic change of coordinates e = sin ϕ, p ϕ = √ 1 − e 2 p e . The Hamiltonian equations can be writteṅ n = 12 n −2/3 (3n p n ) , d dt (3n p n ) = 5n −5/3 2(3n p n ) 2 + 5p 2 ϕ , ϕ = 10 n −5/3 p ϕ ,ṗ ϕ = 0 .
It is easy to describe the solutions of these equations on the cylinder C (see (28)). There are two lines of equilibria at ψ = 0 and ψ = π and ϕ + 2 5 arctan is a first integral (it is smooth at ψ = π 2 ). These solutions are drawn on Figure 3. It is clear that, on a solution, the maximum possible variation of the variable ϕ is 2/5 π; this implies that, if |ϕ 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 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 [4,5], 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.

APPENDIX A Proof of Proposition 5
Let us prove that the seven points in Proposition 5 are satisfied by a, b, c given by (45) (full control case) and also by a, b, c given by (46) (tangential thrust case). (29)  This yields identities (53) with a, b, c given either by (45) or by (46).

Uniqueness of solutions.
This follows from the classical Cauchy-Lipschitz theorem away from S (see Proposition 3). On S, -in the full control case (a, b given by (45)), as seen in [3], the regularity properties (35) of the right hand side of (54) guarantee the existence and uniqueness of solutions to the Cauchy problem (Kamke uniqueness Theorem [13, chap. III, Th. 6.1]), -in the tangential thrust case (a, b given by (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. [10]).
Continuity of Φ in (55), is, according to [13, chap. V, Theorem 2.1], guaranteed by uniqueness of solutions and continuity of a, b.
Since (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 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 .
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.

Hyperbolic saddle. Smoothness around the origin follows from Proposition 3. 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 (74)).
In the full control case, The eigenvalues of the Jacobian 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 (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 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 A short computation shows that a(0, ϕ) = cos ϕ sin ϕ π π 0 cos 2 E dE 1 − sin 2 ϕ cos 2 E , a(π, ϕ) = −a(0, ϕ) , 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 becausē ϕ(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 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 (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 (66). 7.2 tangential thrust case. Let us first compute some values of a and b on special lines. Making the change of variables y = cos E one gets The unstable manifold comprises the equilibrium point (0, 0) and two solutions that tend towards it as time τ tends to −∞, tangent, according to (77), to the line {ψ = 1/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 (79) and (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.

B Proof of Theorem 3
Theorem 3 and this section are independent of the rest of the paper: here we only refer to the seven conditions ranging from equation (53) to equation (66).
We also use the " " and "+" notation to denote the transformations in C: Let us make further constructions and remarks on the conditions (53)-(66) before proceeding with the proof per se.
Stable and unstable manifolds. Equation (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. [14, section 8.3]). Thus it has a stable manifold S 0 and an unstable manifold U 0 ; these are curves passing through (0, 0) tangent to the corresponding eigenvectors. Their existence is a consequence of (61) but point 7 assumes a more specific description.
(91) Fig. 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 1 and Figure 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.
We now state and prove two preliminary lemmas and give the proper proof of Theorem 3.
Lemma 2 Assume that a, b satisfy assumptions (53) to (60), and consider a solution t → (ψ(t), ϕ(t)) of (54) defined on [0, τ fin ], τ fin > 0. For any number f , 0 < f < π 2 , let Lemma 3 Assume that a, b, c satisfy assumptions (53) to (62). For any number f , 0 < f < π 2 , there are two neighborhoods Ω 0 and Ω π of (0, 0) and (π, 0) respectively and, a number T (f ) ≥ 0 such that no solution t → (ψ(t), ϕ(t)) of (54) may cross both Ω 0 and Ω π , and any solution t → (ψ(t), ϕ(t)) defined on the time interval [0, 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. 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: and, if σ = 0,  (53) and (103), we get that, for ε small enough, a does not vanish on K ε .
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 1, 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 -This is obvious if the solution is an equilibrium or is on the unstable manifold U 0 , that do not cross K (see (110)). -In E: among the rectangles in (103) 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 (106) and a(0, σ) < 0 according to (60)) andφ is positive in E, a solution may only exit through the top or left-hand edge; if it exits through the top edge {ϕ =kε}, it will not enter again because ϕ will remain larger thankε. 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 (116). 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 2, 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 (116) for solutions that remain in F . -A solution in the upper part of the stable manifold S 0 also satisfies (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 (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 (111), ψ(t) varies monotonically in [t 1 , t 2 ], and so its variation is at most σ: a(t 2 − t 1 ) < σ. This yields 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 (112). Hence, for any solution, Using Lemma 2, 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 (113) then implies Setting T (f ) = σ/a + 1 + (π−2 )/m(f ), (115), (117), (118) and (119) imply (101).
Proof of Theorem 3. From Lemma 1, the "+" symmetry allows one to interchange ϕ 1 and ϕ 0 while the " " symmetry changes their sign: we may assume in the proof without loss of generality. We distinguish four cases.